如何在numpy中将多个矩阵乘法向量化？

Question

为了更好地理解我的意思，我有两个数据点：

x_0 = np.array([0.6, 1.4])[:, None]
x_1 = np.array([2.6, 3.4])[:, None]

以下是一个2x2矩阵:

y = np.array([[2, 2], [2, 2]])

如果我执行x_0.T @ y @ x_0，我会得到array([[ 8.]])。类似地，x_1.T @ y @ x_1返回array([[ 72.]])。

但是有没有一种方法可以在不使用循环的情况下一次性执行这两个计算？显然，在这里加速可以忽略不计，但是我处理的数据点比这里呈现的要多得多。

- Keir Simmons

---

x_是如何存储的？ - Divakar

它们是列堆叠的。每一列包含一个像素的RGB值，共有n列（对应n个像素）。存储在一个np.array中。 - Keir Simmons

Answer 1

将 x_0、x_1等列进行堆叠，得到 x，我们可以使用np.einsum函数 -

np.einsum('ji,jk,ki->i',x,y,x)

通过使用np.einsum和矩阵乘法的组合 -

np.einsum('ij,ji->i',x.T.dot(y),x)

正如先前所述，x 假定是按列堆叠的，就像这样：

x = np.column_stack((x_0, x_1))

运行时测试 -

In [236]: x = np.random.randint(0,255,(3,100000))

In [237]: y = np.random.randint(0,255,(3,3))

# Proposed in @titipata's post/comments under this post
In [238]: %timeit (x.T.dot(y)*x.T).sum(1)
100 loops, best of 3: 3.45 ms per loop

# Proposed earlier in this post
In [239]: %timeit np.einsum('ji,jk,ki->i',x,y,x)
1000 loops, best of 3: 832 µs per loop

# Proposed earlier in this post
In [240]: %timeit np.einsum('ij,ji->i',x.T.dot(y),x)
100 loops, best of 3: 2.6 ms per loop

- Divakar

---

谢谢，这似乎有效。我已经查看了“einsum”的文档，但还不理解它的工作原理，尽管它似乎非常强大。我会继续研究它。同样，是否可能在没有“einsum”的情况下完成这个操作？ - Keir Simmons

@KeirSimmons，能否解答我在问题下提出的疑问？ - Divakar

@KeirSimmons 如果需要额外的内存要求，一个更简单的方法是 - (x.T.dot(y)*x.T).sum(1)。 - Divakar

Answer 2

基本上，您想对所有拥有的x执行操作(x.T).dot(A).dot(x)。

x_0 = np.array([0.6, 1.4])[:, None]
x_1 = np.array([2.6, 3.4])[:, None]
x = np.hstack((x_0, x_1)) # [[ 0.6  2.6], [ 1.4  3.4]]

简单来说，对于所有的x_i，你可以将它们与y相乘，如下：

[x_i.dot(y).dot(x_i) for x_i in x.T]
>> [8.0, 72.0]

当然，这种方法并不是太高效。但是，你可以使用技巧，先将x与y进行点积，然后再乘回去并对列求和，即手动执行点积。这将使计算速度更快：

x = x.T
(x.dot(y) * x).sum(axis=1)
>> array([  8.,  72.])

请注意，我首先转置矩阵，因为我们想将y的列乘以每行的x。

- titipata