N个矩阵相乘——符号计算
Multiplying N matrices - symbolic calculation
将 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")
对于您的特殊用例,我建议使用 numpy.linalg.matrix_power
(链接问题中未提及)。
时间
这是我使用的设置代码:
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
方法更适合 Sympy。结果数组中的表达式将用更少的项更简化。以下是计算结果数组中的项的方法:
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)
运行得快得多。
将 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")
对于您的特殊用例,我建议使用 numpy.linalg.matrix_power
(链接问题中未提及)。
时间
这是我使用的设置代码:
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
方法更适合 Sympy。结果数组中的表达式将用更少的项更简化。以下是计算结果数组中的项的方法:
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)
运行得快得多。