Project

General

Profile

Bug #15233

Speeding Up Matrix Powers

Added by octonion (Christopher Long) 8 months ago. Updated 8 months ago.

Status:
Assigned
Priority:
Normal
Target version:
-
ruby -v:
ruby 2.5.1p57 (2018-03-29 revision 63029) [x86_64-linux-gnu]
[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

History

Updated by duerst (Martin Dürst) 8 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) 8 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) 8 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) 8 months ago

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

Also available in: Atom PDF