使用Numpy(无循环)将轴乘以矩阵(点积)

Dan*_*lah 7 product numpy matrix multiplication

我正在与Numpy一起处理图像处理问题,我正在尝试避免循环并执行以下操作:

我有一个Dims矩阵M NxNxKxK(它是矩阵KxK的矩阵NxN),对于每一行,我希望乘以(点积)行中的所有N个矩阵(KxK).因此,如果我在完整的M(所有行)上执行此操作,我得到矩阵的向量V(Nx1)(KxK),其中V [i]保持M [i,0] xM [i,1]的点积]×... 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)
Run Code Online (Sandbox Code Playgroud)

小智 2

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

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

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

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

  • 当 N=100、k=2 时,for 循环版本需要 7.5 ms,einsum 版本需要 4.31 ms。
  • 当 N=100、k=20 时,for 循环版本需要 27.3 ms,einsum 版本需要 153 ms。

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