Project

General

Profile

Actions

Feature #15233

closed

Speeding Up Matrix Powers

Added by octonion (Christopher Long) over 5 years ago. Updated over 3 years ago.

Status:
Closed
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) over 5 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) over 5 years ago

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

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) over 5 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) over 5 years ago

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

Actions #5

Updated by jeremyevans0 (Jeremy Evans) over 4 years ago

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

Updated by Anonymous over 3 years 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]

Actions

Also available in: Atom PDF

Like0
Like0Like0Like0Like0Like0Like0