## Bug #15233

### Speeding Up Matrix Powers

**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

### History

#### Updated by duerst (Martin DÃ¼rst) 4 months 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) 4 months 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
index 62852bdad0..fd9115cfad 100644
--- 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[0] == 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) 4 months 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) 4 months ago

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