CODEWITHSUNDEEP
numpy
Ziqi Liu
Ziqi Liu

Reputation: 3171

numpy large matrix multiplication optimize

I have to do a iterative calculation with large matrix: R(t) = M @ R(t-1), where M is n x n, and R is n x 1

if I write this:

for _ in range(iter_num):
    R = M @ R

I suppose it will be very slow, because it has to copy and create new array each time. Is that any way to optimize this? (maybe do it inplace?)

Upvotes: 1

Views: 1718

Answers (4)

Daniel F
Daniel F

Reputation: 14399

Using explicit spectral decomposition my be useful if iter_num is large compared to n (assuming np.lialg.matrix_power doesn't do this already) and M is invertible:

def mat_pow(a, p):
    vals, vecs = np.linalg.eig(a)
    return vecs @ np.diag(vals**p) @ vecs.T

mat_pow(M, iter_num) @ R

If M is symmetric, you could use the even faster np.linalg.eigh

Upvotes: 1

Tai
Tai

Reputation: 7984

Using numpy.linalg.multi_dot

np.linalg.multi_dot([M] * iter_num + [R])

([M] * iter_num creates a list of references to M.)

Some thoughts mentioned in the documentation,

(multi_dot) Compute the dot product of two or more arrays in a single function call, while automatically selecting the fastest evaluation order.

and

Think of multi_dot as:

def multi_dot(arrays): return functools.reduce(np.dot, arrays)

Note OP's method is actually quite fast. See Paul Panzer's answer for more timing results.

Thanks for Paul Panzer's suggestion for using reference rather than view.

Upvotes: 2

Paul Panzer
Paul Panzer

Reputation: 53029

A few timings to show that OP's approach is actually quite competitive:

>>> import functools as ft
>>> kwds = dict(globals=globals(), number=1000)
>>> R = np.random.random((200,))
>>> M = np.random.random((200, 200))
>>> def f_op(M, R):
...     for i in range(k):
...         R = M@R
...     return R
... 
>>> def f_pp(M, R):
...     return ft.reduce(np.matmul, (R,) + k * (M.T,))
... 
>>> def f_ag(M, R):
...     return np.linalg.matrix_power(M, k)@R
... 
>>> def f_tai(M, R):
...     return np.linalg.multi_dot([M]*k+[R])
... 
>>> k = 20
>>> repeat('f_op(M, R)', **kwds)
[0.14156094897771254, 0.1264056910004001, 0.12611976702464744]
>>> repeat('f_pp(M, R)', **kwds)
[0.12594187198556028, 0.1227772050187923, 0.12045996301458217]
>>> repeat('f_ag(M, R)', **kwds)
[2.065609384997515, 2.041590739012463, 2.038702343008481]
>>> repeat('f_tai(M, R)', **kwds)
[3.426795684004901, 3.4321794749994297, 3.4208814119920135]
>>>
>>> k = 500
>>> repeat('f_op(M, R)', **kwds)
[3.066054102004273, 3.0294102499901783, 3.020273027010262]
>>> repeat('f_pp(M, R)', **kwds)
[2.891954762977548, 2.8680382019956596, 2.8558325179910753]
>>> repeat('f_ag(M, R)', **kwds)
[5.216210452985251, 5.1636185249954, 5.157578871003352]

Upvotes: 3

Aguy
Aguy

Reputation: 8059

Would this work for you?

R_final = np.linalg.matrix_power(M, iter_num) @ R

It seems like you are doing M @ M @ M @ ... @ M @ R, which can be cast to M ** iter_num @ R

Upvotes: 1

Related Questions