logo标签列表

使用numpy进行向量操作

pythonnumpymatrixvectorization
3
我有三个numpy数组：
X：一个3073 x 49000的矩阵
W：一个10 x 3073的矩阵
y：一个49000 x 1的向量
y中包含0到9之间的值，每个值代表W中的一行。
我想将X的第一列添加到由y的第一个元素给出的W中的行中。即，如果y的第一个元素为3，则将X的第一列添加到W的第四行。然后将X的第二列添加到由y的第二个元素给出的W中的行中，以此类推，直到X的所有列都添加到由y指定的W中的行中，这意味着总共添加了49000行。
"W[y] += X.T"对我来说不起作用，因为这不会将多个向量添加到W中的一行。
请注意：我只在寻找向量化的解决方案。即无需使用for循环。
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]])

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

因此，所需的结果为：

W = [[17, 22, 16, 16],
     [ 7, 11, 14, 17],
     [ 8,  2, 10, 11]]
看起来非常类似于将加法向量化为由另一个数组索引的数组。 - Divakar
1
“no loop”是一个速度问题还是一个编程挑战问题？ - hpaulj
这是一个编程挑战问题，其动机在于速度。就我而言，这是一个编程挑战，但我不能使用循环的原因是为了练习编写更高效的代码。 - Skeppet
4个回答
3

向量化方法

方法 #1

基于这个答案，下面是使用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

方法二

您可以充分利用布尔索引和np.einsum，以简洁的向量化方式完成工作 -

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

循环方法

方法 #3

由于您需要从每个唯一的y中选择和累加来自X的大量列，因此在性能方面，运行一个循环比较好，其复杂度等于这些唯一的y's的数量，这似乎最多等于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()])
感谢您的解决方案。但是，它仍然比for循环慢...特别是在原始问题中提到的矩阵大小上尤为明显。 - Skeppet
@Skeppet 对此并不感到意外。请查看修改后的for循环版本for_loop_v2，在性能方面略有改进，请点击here。 - Divakar
@Skeppet 在循环中添加了低复杂度的代码。 - Divakar
http://stackoverflow.com/a/28205888/901925 提到，如果 y 中存在空缺，bincount 的解决方案会引发错误。 - hpaulj
@hpaulj 感谢您的编辑和指出错误案例！此外，我的编辑应该处理好了越界/错误情况。 - Divakar
3

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

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是我几年前对类似问题（针对1D数组）给出的答案。

但当处理大型数组时：

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维上迭代也可能具有速度优势。或者最好是在10维上，如Divakar所示。
对于小的测试案例a,b,c=3,4,5，add.at解决方案最快，其次是Divakar的bincount和einseum。对于较大的a,b,c=10,1000,20000，add.at变得非常慢，bincount是最快的。
相关SO答案 https://stackoverflow.com/a/28205888/901925（指出bincount要求完全覆盖y）。 https://dev59.com/5F0a5IYBdhLWcg3wsKe2#30041823（其中Divakar再次表明bincount最优！）
https://github.com/numpy/numpy/issues/5922 - ufunc.at 的性能太慢，比 10x。开发人员已经讨论了 .at 和 bincount 的相对速度。这部分归因于一个的通用性与另一个的专业用途。 - hpaulj
哇，那是一篇非常详尽的帖子。非常感谢！关于矩阵大小，它们不是我的选择。我只是在跟随大学课程做练习。 - Skeppet
好的，一个类。这很有道理。我看过至少一个早期的SO问题，其中包含3073和49000。 - hpaulj
1
这将实现您想要的效果：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]]
谢谢您清晰明了的回答！很抱歉我的解释不够清楚。我已经编辑了问题并附上了一个例子。然而，使用您的解决方案仍然存在从“X”添加多列的问题。 - Skeppet
不好意思，就像我之前说的一样，你提出的解决方案仍存在问题。如果同一个数字在“y”中出现多次（对于49000行来说肯定会出现），你提出的解决方案只会将由该数字的最后一次出现所指示的“X”列添加。我会再次更新问题让你看到区别。 - Skeppet
正如您所看到的，它适用于W中的第二行，因为该行仅在y中出现一次，但对于第0行和第1行则不适用。 - Skeppet
0
虽然我不能说这是非常符合Python风格的，但它是一个解决方案（我认为）：
for column in range(x.shape[1]):
    w[y[column]] = x[:,column].T
抱歉，不能使用for循环。需要一个向量化的解决方案。还是非常感谢你的努力！ - Skeppet
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接