python中的无环3D矩阵乘法答案

【问题标题】：A loopless 3D matrix multiplication in pythonpython中的无环3D矩阵乘法
【发布时间】：2011-07-17 17:18:27
【问题描述】：

我希望在 python (numpy) 中执行以下操作。

Matrix A is M x N x R
Matrix B is N x 1 x R

矩阵乘法 AB = C，其中 C 是 M x 1 x R 矩阵。本质上，A（其中的 R）的每个 M x N 层都是矩阵独立地乘以 B 中的每个 N x 1 向量。我确信这是一个单线。我一直在尝试使用 tensordot()，但我似乎给了我意想不到的答案。

我已经在 Igor Pro 中编程近 10 年了，现在我正在尝试将它的页面转换为 python。

【问题讨论】：

标签： python matrix numpy operation

【解决方案1】：

对不起，死灵术，但这个答案可以大大改善，使用无价的 np.einsum。

import numpy as np

D,M,N,R = 1,2,3,4
A = np.random.rand(M,N,R)
B = np.random.rand(N,D,R)

print np.einsum('mnr,ndr->mdr', A, B).shape

请注意，它有几个优点：首先，速度快。 np.einsum 通常优化得很好，但此外，np.einsum 足够聪明，可以避免创建 MxNxR 临时数组，而是直接对 N 执行收缩。

但也许更重要的是，它的可读性很强。毫无疑问，这段代码是正确的；你可以让它变得更复杂而没有任何麻烦。

请注意，如果您愿意，可以简单地从 B 和 einsum 语句中删除虚拟“D”轴。

【讨论】：

我看到 np.dot() 也可以做一些多维操作，但是在一些奇怪的规则下工作。你有这方面的知识吗？
我的知识可以最好地概括为“改用 einsum”。可能有点冗长，但在我看来，“显式胜于隐式”从未如此适用。

【解决方案2】：

numpy.tensordot() 是正确的做法：

a = numpy.arange(24).reshape(2, 3, 4)
b = numpy.arange(12).reshape(3, 1, 4)
c = numpy.tensordot(a, b, axes=[1, 0]).diagonal(axis1=1, axis2=3)

编辑：第一个版本有问题，这个版本计算了更多它应该计算的，并丢弃了大部分。也许在最后一个轴上使用 Python 循环是更好的方法。

另一个编辑：我得出的结论是numpy.tensordot() 不是最好的解决方案。

c = (a[:,:,None] * b).sum(axis=1)

会更有效率（虽然更难掌握）。

【讨论】：

感谢您这么快回复我。至少，它会让我开始。我将尝试编写的大部分代码都是面向矩阵运算的......所以我真的应该尝试了解这里发生了什么。话虽如此，代码中有两个部分让我感到困惑。首先是 tensordot 内部的调用“轴”。我对这实际上做了什么感到困惑，因为我希望 a x b 应该只给你 c （就像我上面描述的那样）而不声明任何特别的东西。也许一旦我理解了这一点，我就会明白为什么有必要使用 .diagonal。
这很聪明......我不知道我需要多长时间才能做这样的事情（看起来你创建了一个新轴来进行乘法和，然后基本上重新组合）非常感谢您的宝贵时间，非常感谢！
如果第一个维度是矩阵列表 (R)，第二个维度是公共矩阵维度 (N) 即 a 的形状为 (4, 3, 2)，则更容易想到这一点) 和 b 是 (4, 3, 1)。然后乘法运算变为 (a * b).sum(axis=1)。 (a * b) 是矩阵中每个元素的行和列的乘积，需要对每一行求和才能得到最终的矩阵。就像您手动处理多个矩阵一样。

【解决方案3】：

另一种方法（对于像我这样不熟悉爱因斯坦符号的人来说更容易）是np.matmul()。重要的是在最后两个索引中具有匹配的维度 ((M, N) x (N, 1))。为此使用np.transpose() 示例：

M, N, R = 4, 3, 10
A = np.ones((M, N, R))
B = np.ones((N, 1, R))

# have the matching dimensions at the very end
C = np.matmul(np.transpose(A, (2, 0, 1)), np.transpose(B, (2, 0, 1))) 
C = np.transpose(C, (1, 2, 0))

print(A.shape)
# out: #(4, 3, 10)
print(B.shape)
# out: #(3, 1, 10)
print(C.shape)
# out: #(4, 1, 10)

【讨论】：