Bug #15233: Speeding Up Matrix Powers

• Status: Assigned
• Priority: Normal
• Assignee: marcandre (Marc-Andre Lafortune)
• Target version: -
• ruby -v: ruby 2.5.1p57 (2018-03-29 revision 63029) [x86_64-linux-gnu]
• Backport: 2.3: UNKNOWN, 2.4: UNKNOWN, 2.5: UNKNOWN

[octonion (Christopher Long)]

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

[duerst (Martin Dürst)]

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

[mame (Yusuke Endoh)]

• 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

[CaryInVictoria (Cary Swoveland)]

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

[marcandre (Marc-Andre Lafortune)]

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