使用numpy在Python中计算多个矩阵的乘积

Question

假设您有n个方阵A1，...，An。有没有一种简洁的方法来乘这些矩阵？据我所知，Numpy中的点积只接受两个参数。一种显而易见的方法是定义一个调用自身并获取结果的函数。有没有更好的方法来完成它？

- NNsr

Answer 1

这可能是一个相对较新的特性，但我喜欢它：

A.dot(B).dot(C)

或者，如果你拥有一个长链，你可以这样做：

reduce(numpy.dot, [A1, A2, ..., An])

更新:

这里有更多关于reduce的信息，点击这里。下面是一个可能会有帮助的例子。

>>> A = [np.random.random((5, 5)) for i in xrange(4)]
>>> product1 = A[0].dot(A[1]).dot(A[2]).dot(A[3])
>>> product2 = reduce(numpy.dot, A)
>>> numpy.all(product1 == product2)
True

2016年更新： 从Python 3.5开始，有一个新的矩阵乘法符号，@：

R = A @ B @ C

- Bi Rico

---

谢谢您的回复。第一个选项可以正常工作；但是第二个选项不行；或者至少我没能让它工作。您能否再详细解释一下或者给一个例子？非常感谢。 - NNsr

我经常遇到这个问题，最终写了一个辅助函数。希望这可以成为NumPy的一部分：def xdot(*args): return reduce(np.dot, args)。（翻译者注：该段英文已被完整翻译为中文，无需再翻译。） - rd11

只是添加一条注释，当A、B和C的类型为numpy.ndarray时，这个函数可以正常工作。这可能也适用于其他类型，但我没有进行过检查。 - OfLettersAndNumbers

NameError: 名称 'reduce' 未定义 - rsc05

在Python 3中，reduce被移动到了functools模块中，可以使用from functools import reduce来导入。 - Bi Rico

你能否编辑一下，将 R = A @ B @ C 放在答案顶部（我无法编辑，因为建议队列已满）。然后当搜索这个问题时，duckduckgo 将会预览它。现在，它显示的是 A.dot(B).dot(C)。 - Archisman Panigrahi

Answer 2

更新一个老问题:

截至2014年11月13日，现在有一个np.linalg.multi_dot函数，它正好做你想要的事情。 它还具有优化调用顺序的好处，尽管在您的情况下这并不是必要的。

请注意，此功能从numpy版本1.10开始提供。

- perimosocordiae

---

这应该是未来任何人到达此处的首选答案。 - SajidSalim

“@”和“np.linalg.multi_dot”之间是否存在速度差异？ - kanso37

我使用 A_list = [np.random.random(100, 100) for i in range(3)] 创建了一个数组列表，通过运行 %timeit np.linalg.multi_dot(A_list) 和 %timeit A_list[0] @ A_list[1] @ A_list[2] 来进行简单测试。看起来第二种方法（在我的机器上为100微秒对比85微秒）比第一种方法表现更好，但我不能确定这在一般情况下是否正确。我也想知道如何使用递归地将列表推广到第二种方法。 - Daniel Casas-Orozco

Answer 3

如果你预先计算所有矩阵，那么你应该使用优化方案进行矩阵链乘法。请查看此维基百科文章。

- Florian Brucker

---

谢谢您的评论，但我认为对于方阵来说这并不重要。对吧？ - NNsr

@Nikandish：没错，我在你原始答案中漏掉了那部分。 - Florian Brucker

Answer 4

另一种实现此目的的方式是使用einsum，它为NumPy实现了爱因斯坦求和约定。

简要解释此约定与该问题的关系：当您将多个矩阵乘积写成一个大的乘积总和时，您会得到类似以下的结果：

P_im = sum_j sum_k sum_l A1_ij A2_jk A3_kl A4_lm

P是您的产品结果，A1、A2、A3和A4是输入矩阵。请注意，您需要对出现两次的索引（即j、k和l）进行求和。由于该特性在物理学、向量微积分和其他领域中经常出现，因此存在一个名为einsum的NumPy工具可以实现该功能。

在上面的示例中，您可以使用如下方式计算您的矩阵乘积：

P = np.einsum( "ij,jk,kl,lm", A1, A2, A3, A4 )

在这里，第一个参数告诉函数要应用于参数矩阵的哪些索引，然后所有重复出现的索引都会求和，得到所需的结果。

请注意，计算效率取决于几个因素（因此最好只是进行测试）：

• 为什么numpy的einsum比numpy的内置函数更慢？
• 为什么numpy的einsum比numpy的内置函数更快？

- Wrzlprmft

Answer 5

A_list = [np.random.randn(100, 100) for i in xrange(10)]
B = np.eye(A_list[0].shape[0])
for A in A_list:
    B = np.dot(B, A)

C = reduce(np.dot, A_list)

assert(B == C)

- Nicolas Barbey