## Feature #15233  ### Speeding Up Matrix Powers

Status:
Closed
Priority:
Normal
Target version:
-
[ruby-core:89452]

Description

I've been looking at SymPy's slow matrix power speed, and noticed that replacing Ruby's matrix power code with the equivalent code from SymPy makes it significantly faster. As this is a recursively defined function, presumably it can be made even faster with a proper iterative version.

Base:

``````require 'matrix'

Q = Matrix[[0,0,1,0,2,0],
[0,0,0,1,0,0],
[1,0,0,1,0,2],
[0,2,1,0,0,0],
[1,0,0,0,0,0],
[0,0,1,0,0,0]]

r = Vector[1,0,0,0,0,0]

p = (Q**100000000*r).sum

puts p.to_s.size
``````

Output:

```time ruby main.rb
35950264

real    1m42.250s
user    1m41.050s
sys     0m1.184s
```

Modified:

``````require 'matrix'

def pow(a,num)
if (num==1)
return a
end
if (num%2)==1
return a*pow(a,num-1)
end
ret = pow(a,(num/2))
ret*ret
end

Q = Matrix[[0,0,1,0,2,0],
[0,0,0,1,0,0],
[1,0,0,1,0,2],
[0,2,1,0,0,0],
[1,0,0,0,0,0],
[0,0,1,0,0,0]]

r = Vector[1,0,0,0,0,0]

p = (pow(Q,100000000)*r).sum

puts p.to_s.size
``````

Output:

```time ruby main.rb
35950264

real    1m9.489s
user    1m8.661s
sys     0m0.828s
```

#### Updated by duerst (Martin Dürst)about 2 years ago

The code doesn't look very Ruby-ish. What about rewriting it as follows:

``````def pow(a, num)
if num==1
a
elsif num.odd?
a*pow(a, num-1)
else
ret = pow(a, num/2)
ret*ret
end
end
``````

#### Updated by mame (Yusuke Endoh)about 2 years ago

• Assignee set to marcandre (Marc-Andre Lafortune)
• Status changed from Open to Assigned

Interesting. The current implementation accumulates the value in LSB-to-MSB order, and SymPy's implementation does in MSB-to-LSB because. The former is slower than the latter because the latter multiplies the original matrix (which is smaller than an accumulated matrix for each digit).

marcandre (Marc-Andre Lafortune), could you review this?

``````diff --git a/lib/matrix.rb b/lib/matrix.rb
--- a/lib/matrix.rb
+++ b/lib/matrix.rb
@@ -1086,12 +1086,12 @@ def **(other)
return self.class.identity(self.column_count) if other == 0
other = -other
end
-      z = nil
-      loop do
-        z = z ? z * x : x if other == 1
-        return z if (other >>= 1).zero?
-        x *= x
+      z = self
+      (other.bit_length - 2).downto(0) do |i|
+        z *= z
+        z *= self if other[i] == 1
end
+      return z
when Numeric
v, d, v_inv = eigensystem
v * self.class.diagonal(*d.each(:diagonal).map{|e| e ** other}) * v_inv
``````

#### Updated by CaryInVictoria (Cary Swoveland)about 2 years ago

The O(ln n) method could be written as follows.

```def pow1(m, n)
return m if n == 1
p = pow1(m, n/2)
n.even? ? p*p : p*p*m
end

t = Time.now
p = (pow(Q,100000000)*r).sum % 1_000_000
#=> 109376
Time.now-t
#=> 53.2

t = Time.now
p1 = (pow1(Q,100000000)*r).sum % 1_000_000
#=> 109376
Time.now-t
#=> 54.4
```

#### Updated by marcandre (Marc-Andre Lafortune)about 2 years ago

Cool. I'll definitely have a look, in a few days probably as I'm travelling right now

#### Updated by jeremyevans0 (Jeremy Evans)over 1 year ago

• Backport deleted (2.3: UNKNOWN, 2.4: UNKNOWN, 2.5: UNKNOWN)
• ruby -v deleted (ruby 2.5.1p57 (2018-03-29 revision 63029) [x86_64-linux-gnu])
• Tracker changed from Bug to Feature

#### Updated by Anonymous about 4 hours ago

• Status changed from Assigned to Closed

Applied in changeset git|a83a51932dbc31b549e11b9da8967f2f52a8b07c.

[ruby/matrix] Optimize **

Avoiding recursive call would imply iterating bits starting from
most significant, which is not easy to do efficiently.
Any saving would be dwarfed by the multiplications anyways.
[Feature #15233]

Also available in: Atom PDF