numpy：广播矩阵乘以数组答案

【问题标题】：numpy: broadcast matrix multiply accross arraynumpy：广播矩阵乘以数组
【发布时间】：2014-04-10 05:00:21
【问题描述】：

我有一个 3xN 数组，概念上是一个 N 3 向量数组，我想构造矩阵将给定的3x3 矩阵与大批。有没有一种好的方法可以以矢量化的方式做到这一点？

目前，我的问题是3xN，但我以后可能需要考虑3xNxM（或更多）。

循环方法

U=numpy.rand( [3,24] )

R=numpy.eye(3) # placeholder

for i in xrange( U.shape[1]):
    U[:,i]=numpy.dot( R, U[:,i] )

【问题讨论】：

在这里查看我的答案：stackoverflow.com/a/22081723/553404
@MrE 这个答案很好，对于我的体型 (3xN)，我避免了你在那里提到的转置。

标签： python arrays numpy matrix

【解决方案1】：

使用 np.einsum 函数，您甚至可以解决多维问题：

U = np.random.rand(3,24,5) 
R = np.eye(3,3)
result = np.einsum( "ijk,il", U,R )

符号有点棘手：您首先给出的字符串说明了数组维度的索引；因此对于 U 而言，索引 ijk 是每个维度的运行索引。它遵循爱因斯坦求和约定，因此字符串中具有相同字母的索引将被求和。详情请阅读ScipyDocs。我敢肯定，在您的情况下，点会更快，因为开销更少，并且可能会使用一些 blas 例程，但是正如您所说，您想扩展到更多维度，这可能是要走的路。

【讨论】：

在实际的基准测试中，einsum 经常领先。特别是在这些情况下，涉及小尺寸的收缩。我想 BLAS 是通过 gufunc 机制调用的；也就是说，对每个微小的子矩阵进行单独的 BLAS 调用。在这种情况下，BLAS 会因为 O(N) 开销而受到影响，这很快就会超过 einsum 的 O(1) 开销。

【解决方案2】：

在这种情况下，您只需调用 np.dot(R, U) 即可：

import numpy as np

np.random.seed(0)

U = np.random.rand(3,24) 
R = np.random.rand(3,3)

result = np.empty_like(U)

for i in range( U.shape[1]): 
    result[:,i] = np.dot(R, U[:,i])

print result

print np.allclose(result, np.dot(R, U))

对于(3xNxM) 的情况，您可以重新整形为 (3x(N.M)), dot 并将结果重新整形，类似于我的回答 here

【讨论】：

我认为尽可能使用 np.dot 实际上是可行的方法，因为它经常使用经过高度优化的非常快速的库，例如 Blas 或 Lapack。