对于非矢量化版本,我有一系列(2, 2)
形矩阵(即形状为(n, 2, 2)
的ndarray),并且需要按顺序乘以它们(矩阵乘法),这意味着n
顺序矩阵乘法。 这是一个最小示例的样子
def get_matrix_product_eig_val(J):
# J holds the sequence of matrices to multiply and has shape=(n, 2, 2)
M = np.identity(2, dtype=np.double)
for i in range(n_gens):
M = np.matmul(M, J[i])
eig_val, eig_vec = np.linalg.eig(M) # eig_val is what I'm interested in
return eig_val
现在我有一个k
这样的矩阵序列数组(形状(k, n, 2, 2)
的ndarray),并且需要为每个k
条目做相同的顺序矩阵积。 天真的方法是做
# Now J has shape=(k, n, 2, 2)
for i in range(k):
eig_vals[i] = get_matrix_product_eig_val(J[i, :, :, :])
有没有办法摆脱循环并以矢量化的方式做到这一点?
笔记:
1)n
预计在~100的数量级内。k
可以是 ~100 到 ~10,000 之间的任何位置
2)有人建议用np.linalg.multi_dot
替换内环。这实际上大大减慢了速度
3)我可以看到3x3
矩阵的遥远未来应用,但2x2
的特定解决方案很好。无论哪种方式,所有矩阵都是正方形的,并且具有相同的维度
for i in range(n - 1):
J[:, i + 1, :, 0] = np.broadcast_to((J[:, i+1, 0, 0])[..., None], (k, 2)) * J[:, i, :, 0] +
np.broadcast_to((J[:, i+1, 1, 0])[..., None], (k, 2)) * J[:, i, :, 1]
J[:, i + 1, :, 1] = np.broadcast_to((J[:, i+1, 0, 1])[..., None], (k, 2)) * J[:, i, :, 0] +
np.broadcast_to((J[:, i+1, 1, 1])[..., None], (k, 2)) * J[:, i, :, 1]
eig_vals, _ = np.linalg.eig(J[:, -1, :, :])