几个 2x2 矩阵序列的向量化乘法答案

【问题标题】：Vectorized multiplication of several sequences of 2x2 matrices几个 2x2 矩阵序列的向量化乘法
【发布时间】：2019-06-12 18:50:19
【问题描述】：

对于未矢量化的版本，我有一系列(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 的特定解决方案很好。无论哪种方式，所有矩阵都是正方形的，并且具有相同的维度

【问题讨论】：

您在寻找np.linalg.multi_dot吗？
multi_dot 按顺序进行乘法运算。它所添加的只是一个测试，看看以不同的顺序进行乘法是否有任何优势。如果所有数组都具有相同的形状，则使用它没有任何优势。
你能把这个问题看作是对 (k,l,2,2) 数组进行n 顺序乘法吗？ matmul 可以将一个 (k,l,2,2) 数组与另一个数组相乘。 A1 @ A2 @ A3 ...
如果您的目标是性能，您可以这样做 stackoverflow.com/a/51062850/4045774 如果小矩阵的维度始终为 2x2，则可能会有进一步的改进。只需用一个真实的完整工作示例更新您的问题...

标签： python numpy vectorization

【解决方案1】：

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, :, :])

【讨论】：