使用 numpy 进行向量运算

Question

我有三个 numpy 数组：

X：3073 x 49000 矩阵
W：10 x 3073 矩阵
y：49000 x 1 向量

y 包含 0 到 9 之间的值，每个值代表W 中的一行。

我想将X 的第一列添加到W 中由y 中的第一个元素给出的行中。 IE。如果y 的第一个元素是3，则将X 的第一列添加到W 的第四行。然后将X的第二列添加到y中第二个元素给定的W中的行中，以此类推，直到X中的所有列都添加到@指定的W中的行中987654338@，表示总共增加了49000行。

W[y] += X.T 对我不起作用，因为这不会在W 的一行中添加多个向量。

请注意：我只是在寻找矢量化解决方案。 IE。没有 for 循环。

编辑：为了澄清，我将添加一个小矩阵大小的示例，该示例改编自萨尔瓦多·达利（Salvador Dali）的以下示例。

In [1]: import numpy as np

In [2]: a, b, c = 3, 4, 5

In [3]: np.random.seed(0)

In [4]: X = np.random.randint(10, size=(b,c))

In [5]: W = np.random.randint(10, size=(a,b))

In [6]: y = np.random.randint(a, size=(c,1))

In [7]: X
Out[7]: 
array([[5, 0, 3, 3, 7],
       [9, 3, 5, 2, 4],
       [7, 6, 8, 8, 1],
       [6, 7, 7, 8, 1]])

In [8]: W
Out[8]: 
array([[5, 9, 8, 9],
       [4, 3, 0, 3],
       [5, 0, 2, 3]])

In [9]: y
Out[9]: 
array([[0],
       [1],
       [1],
       [2],
       [0]])

In [10]: W[y.ravel()] + X.T
Out[10]: 
array([[10, 18, 15, 15],
       [ 4,  6,  6, 10],
       [ 7,  8,  8, 10],
       [ 8,  2, 10, 11],
       [12, 13,  9, 10]])

In [11]: W[y.ravel()] = W[y.ravel()] + X.T

In [12]: W
Out[12]: 
array([[12, 13,  9, 10],
       [ 7,  8,  8, 10],
       [ 8,  2, 10, 11]])

问题是BOTH将 X 中的第 0 列和第 4 列添加到 W 中的第 0 行，以及将 X 中的第 1 列和第 2 列都添加到 W 中的第 1 行。

因此，期望的结果是：

W = [[17, 22, 16, 16],
     [ 7, 11, 14, 17],
     [ 8,  2, 10, 11]]

Answer 1

矢量化方法

方法#1

基于this answer，这里是使用np.bincount 的矢量化解决方案-

N = y.max()+1
id = y.ravel() + np.arange(X.shape[0])[:,None]*N
W[:N] += np.bincount(id.ravel(), weights=X.ravel()).reshape(-1,N).T

方法 #2

您可以充分利用boolean indexing 和np.einsum 以简洁的矢量化方式完成工作-

N = y.max()+1
W[:N] += np.einsum('ijk,lk->il',(np.arange(N)[:,None,None] == y.ravel()),X)

---

循环方法

方法#3

由于您要从X 中选择并添加大量列，每个唯一的y，在性能方面，使用complexity 等于此类唯一@987654337 的数量运行循环可能会更好@，似乎在max 等于W 中的行数，而在您的情况下，这只是10。因此，循环只有 10 次迭代，还不错！这是实现这些愿望的实现 -

for k in range(W.shape[0]):
    W[k] += X[:,(y==k).ravel()].sum(1)

方法#4

您可以引入np.einsum 来进行按列求和并得到最终输出，如下所示 -

for k in range(W.shape[0]):
    W[k] += np.einsum('ij->i',X[:,(y==k).ravel()])

Answer 2

首先直接循环解决方案作为参考：

In [65]: for i,j in enumerate(y):
    W[j]+=X[:,i]
   ....:     

In [66]: W
Out[66]: 
array([[17, 22, 16, 16],
       [ 7, 11, 14, 17],
       [ 8,  2, 10, 11]])

add.at 解决方案：

In [67]: W=W1.copy()
In [68]: np.add.at(W,(y.ravel()),X.T)
In [69]: W
Out[69]: 
array([[17, 22, 16, 16],
       [ 7, 11, 14, 17],
       [ 8,  2, 10, 11]])

add.at 进行无缓冲计算，绕过阻止W[y.ravel()] += X.T 工作的缓冲。它仍然是迭代的，但循环已移至编译代码。这不是真正的矢量化，因为应用程序的顺序很重要。一行X.T 的加法取决于前几行的结果。

https://stackoverflow.com/a/20811014/901925 是我几年前对类似问题（对于一维数组）给出的答案。

但是在处理大型数组时：

X: a 3073 x 49000 matrix
W: a 10 x 3073 matrix
y: a 49000 x 1 vector 

这可能会遇到速度问题。请注意W[y.ravel()] 与X.T 的大小相同（您为什么选择这些需要转置的大小？）。它是一个副本，而不是一个视图。所以已经有时间惩罚了。

bincount已经在之前的问题中提出过，我认为它更快。 Making for loop with index arrays faster（bincount 和 add.at 解决方案）

在小的 3073 维度上进行迭代也可能具有速度优势。或者像Divakar 演示的那样在尺寸为 10 的维度上更好。

---

对于小型测试用例a,b,c=3,4,5，add.at 解决方案最快，其次是Divakar's bincount 和einseum。对于更大的a,b,c=10,1000,20000，add.at 变得非常慢，bincount 是最快的。

相关的 SO 答案

https://stackoverflow.com/a/28205888/901925（注意bincount 需要完全覆盖y）。

https://stackoverflow.com/a/30041823/901925（其中Divakar 再次显示bincount 规则！）

Answer 3

虽然我不能说这是非常 Pythonic，但它是一个解决方案（我认为）：

for column in range(x.shape[1]):
    w[y[column]] = x[:,column].T

Answer 4

这将实现你想要的：X + W[y.ravel()].T

---

要查看这确实有效，这里有一个可重现的示例：

import numpy as np
np.random.seed(0)
a, b, c = 3, 5, 4  # you can use your 3073, 49000, 10 later

X = np.random.rand(a, b)
W = np.random.rand(c, a)
y = np.random.randint(c, size=(b, 1))

现在你的矩阵是：

[[ 0.0871293   0.0202184   0.83261985]
 [ 0.77815675  0.87001215  0.97861834]
 [ 0.79915856  0.46147936  0.78052918]
 [ 0.11827443  0.63992102  0.14335329]]

[[3]
 [0]
 [3]
 [2]
 [0]]

[[ 0.5488135   0.71518937  0.60276338  0.54488318  0.4236548 ]
 [ 0.64589411  0.43758721  0.891773    0.96366276  0.38344152]
 [ 0.79172504  0.52889492  0.56804456  0.92559664  0.07103606]]

而W[y.ravel()] 为您提供“由 y 中的第一个元素给出的 W”。通过转置它，您将得到一个准备添加到 X 的矩阵：

[[ 0.11827443  0.0871293   0.11827443  0.79915856  0.0871293 ]
 [ 0.63992102  0.0202184   0.63992102  0.46147936  0.0202184 ]
 [ 0.14335329  0.83261985  0.14335329  0.78052918  0.83261985]]