Multiplying N matrices - symbolic calculation

Question

What is the most efficient (quickest) way to multiply 20 identical 6x6 matrices (M)?

N = 20
w = sy.Symbol("w");v = sy.Symbol("v");p = sy.Symbol("p");q = sy.Symbol("q");c = 1;n = 1;nc = 1
M = np.array([[w*p*q,w*q,0,0,0,0], 
              [0,0,v,0,0,0], 
              [0,0,0,nc,0,c], 
              [0,0,0,0,v,0], 
              [w,w,v,nc,0,c],
              [0,0,0,n,0,1]])
Mi = np.array([[w*p*q,w*q,0,0,0,0],
               [0,0,v,0,0,0],
               [0,0,0,nc,0,c],
               [0,0,0,0,v,0], 
               [w,w,v,nc,0,c],
               [0,0,0,n,0,1]])
for l in range(N-1):
    M = np.dot(M, Mi)
difZ = sy.diff(Z2,w)
expr = w*(np.divide(difZ,Z2))
Z_lamda = sy.lambdify([w,v,p,q], expr, "numpy")

https://stackoverflow.com/questions/11838352/multiply-several-matrices-in-numpy — Employee, Jan 20 '19 at 12:38
The last 3 lines of your code sample are useless - you haven't defined `Z2`. Your array, `M` is object dtype, containing `sympy` objects. Math using object dtype arrays is much slower than math using `numpy` numeric arrays, and somewhat hit-or-miss, depending on whether it can delegate the actions to methods of the elements. — hpaulj, Jan 20 '19 at 17:34

tel · Accepted Answer · 2019-01-20T17:43:51.020

For your special use case, I'd recommend using numpy.linalg.matrix_power (which wasn't mentioned in the linked question).

Timings

Here's the setup code I used:

import numpy as np
import sympy as sy
sy.init_printing(pretty_print=False)

N = 20
w = sy.Symbol("w");v = sy.Symbol("v");p = sy.Symbol("p");q = sy.Symbol("q");c = 1;n = 1;nc = 1
M = np.array([[w*p*q,w*q,0,0,0,0], 
              [0,0,v,0,0,0], 
              [0,0,0,nc,0,c], 
              [0,0,0,0,v,0], 
              [w,w,v,nc,0,c],
              [0,0,0,n,0,1]])
Mi = M.copy()

and here's some timings comparing your original iterative dot approach to matrix_power:

%%timeit
M = Mi.copy()
for _ in range(N-1):
    M = np.dot(M, Mi)
# 527 ms ± 14.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%%timeit
np.linalg.matrix_power(Mi, N)
# 6.63 ms ± 96.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

So matrix_power is about 80X faster.

Added bonus: `matrix_power` works better with arrays of Sympy expressions

For whatever reason, matrix_power seems to work better with Sympy than the iterative dot method. The expressions in the resulting array will be more simplified with fewer terms. Here's how you can count the terms in the result array:

import numpy as np
import sympy as sy

def countterms(arr):
    return np.sum([len(e.args) for e in arr.flat])

N = 20
w = sy.Symbol("w");v = sy.Symbol("v");p = sy.Symbol("p");q = sy.Symbol("q");c = 1;n = 1;nc = 1
M = np.array([[w*p*q,w*q,0,0,0,0], 
              [0,0,v,0,0,0], 
              [0,0,0,nc,0,c], 
              [0,0,0,0,v,0], 
              [w,w,v,nc,0,c],
              [0,0,0,n,0,1]])
Mi = M.copy()

for _ in range(N-1):
    M = np.dot(M, Mi)

Mpow = np.linalg.matrix_power(Mi, N)

print("%d terms total in looped dot result\n" % countterms(M))
print("%d terms total in matrix_power result\n" % countterms(Mpow))

Output:

650 terms total in looped dot result

216 terms total in matrix_power result

In particular, print(Mpow) runs much, much faster than print(M).

`M` is an object dtype array. `np.dot` works because the elements themselves implement basic math, * and +. It is not using the fast numeric BLAS routines. `matrix_power` gains speed by doing repeated matrix squares, reducing the total number of calls to `np.dot`. — hpaulj, Jan 20 '19 at 17:24
`matrix_power` is written in Python, using `np.dot` (or `matmul` in non-object cases). No special `C` code. — hpaulj, Jan 20 '19 at 17:37
@hpaulj Thanks for the info. I updated my answer to match, and included some timings. `matrix_power` is a heck of a lot faster, even without any C deep magic. "repeated matrix squares, reducing the total number of calls to `np.dot`": is that something like the [peasant multiplication algorithm](https://en.wikipedia.org/wiki/Ancient_Egyptian_multiplication)? — tel, Jan 20 '19 at 17:46
`matrix_power(M,4)` produces the same thing as `temp=np.dot(M,M); M4=np.dot(temp,temp)`. `matrix_power(M,5)` would be the `M4` dotted once more with `M`. — hpaulj, Jan 20 '19 at 17:54

Multiplying N matrices - symbolic calculation

1 Answers1

Timings

Added bonus: matrix_power works better with arrays of Sympy expressions

Added bonus: `matrix_power` works better with arrays of Sympy expressions