乘N矩阵符号计算

Question

将20个相同的6x6矩阵（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")

Answer 1

对于您的特殊用例，我建议使用^{}（链接问题中没有提到）。你知道吗

时间安排

以下是我使用的设置代码：

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()

下面是一些时间安排，将您最初的迭代dot方法与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)

所以matrix_power大约快80倍。你知道吗

额外的好处：matrix_power在Sympy表达式数组中工作得更好

不管出于什么原因，matrix_power似乎比迭代dot方法更适合于symphy。结果数组中的表达式将更加简化，只需较少的项。下面是如何计算结果数组中的项：

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))

输出：

650 terms total in looped dot result

216 terms total in matrix_power result

特别是，print(Mpow)比print(M)运行得快得多。你知道吗