2

将 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")
4

1 回答 1

0

对于您的特殊用例，我建议使用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_powerSympy 似乎比迭代dot方法更有效。结果数组中的表达式将使用更少的项更加简化。以下是计算结果数组中项的方法：

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

于 2019-01-20T15:00:43.557 回答