使用 Numpy（无循环）沿轴相乘矩阵（点积）答案

【问题标题】：Multiplying matrices (dot product) along an axis with Numpy (without loops)使用 Numpy（无循环）沿轴相乘矩阵（点积）
【发布时间】：2017-12-13 11:47:06
【问题描述】：

我正在与 Numpy 合作解决图像处理问题，我正在尝试避免循环并执行以下操作：

我有一个 Dims NxNxKxK 的矩阵 M（这是一个矩阵 KxK 的矩阵 NxN），并且对于每一行，我希望将行中的所有 N 个矩阵（KxK）相乘（点积）。因此，如果我在整个 M（所有行）上执行此操作，我会得到矩阵 (KxK) 的向量 V (Nx1)，其中 V[i] 包含 M[i,0]xM[i,1 ]x...xM[i,N-1]。

我使用循环实现了这个问题的解决方案，但我想不出没有循环的方法。

实现（带循环）：

    a = np.array([[1,1,1], [1,1,1], [1,1,1]])
    mat = np.array([[a,a,a,a], [a*2,a*2,a*2,a*2], [a*3,a*3,a*3,a*3],
                    [a*4,a*4,a*4,a*4]])  # the original matrix
    N, N, k, k = mat.shape
    result = np.ones((N, k, k))  # resulting matrix
    for i in range(N):
        k = functools.reduce(np.dot, mat[i,:])
        result[i,:] = k
    print(result)

【问题讨论】：

在您的示例中，N 为 4。在您的实际应用中，N 的大小是否有上限？
你的行是否总是由相同的矩阵组成，即[a,a,a,a]、[a*2,a*2,a*2,a*2]等？
@WarrenWeckesser - 在我们的示例中，您是对的 N=4，但一般来说，否 - 没有限制。
我建议保持原样。
原则上，您可以执行np.einsum('ab,bc,cd,de,ef,fg,gh,hi', *factors[:8]) 之类的操作，但事实证明它比reduce 方法慢得多。

标签： numpy matrix product multiplication

【解决方案1】：

以下使用reduce，但不是 N 上的循环：

mat = mat.swapaxes(0, 1)
result = functools.reduce(lambda a, b: np.einsum('ijk,ikl->ijl', a, b), mat[:])

einsum 符号'jk,kl->jl' 表示矩阵乘法，索引i 表示它应该对第一个索引的每个值进行。 mat[0] 或mat[1] 的第一个索引实际上是mat 的第二个索引（列索引），所以写起来，乘法发生在mat 的每一列中。您希望它在每一行中完成，因此使用swapaxes。

这比for循环版本快还是慢取决于N和k的相对大小。 np.dot 方法经过高度优化，但如果 N 上的循环很长，einsum 可能会获胜。一些%timeit 的结果：

N=100，k=2，for-loop 版本需要 7.5 ms，einsum 版本需要 4.31 ms。
N=100，k=20，for-loop 版本耗时 27.3 ms，einsum 版本耗时 153 ms。

因此，在特定情况下会有适度的收益，而在大多数其他情况下会出现重大损失。但是您没有要求“有效”的解决方案，您要求的是“无循环”，所以它是（“无循环”！=“更快”）。正如 Divakar 在评论中建议的那样，您最好保持代码不变。

【讨论】：