Feature #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) about 1 year 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 1 year 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) about 1 year 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 1 year ago
Cool. I'll definitely have a look, in a few days probably as I'm travelling right now
Updated by jeremyevans0 (Jeremy Evans) 3 months 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