向量化 Numpy 中的操作答案

【问题标题】：Vectorize an operation in Numpy向量化 Numpy 中的操作
【发布时间】：2016-05-15 09:35:58
【问题描述】：

我正在尝试在不使用循环的情况下在 Numpy 上执行以下操作：

我有一个维度为 N*d 的矩阵 X 和一个维度为 N 的向量 y。 y 包含从 1 到 K 的整数。
我正在尝试获取大小为 K*d 的矩阵 M，其中 M[i,:]=np.mean(X[y==i,:],0)

我可以在不使用循环的情况下实现这一点吗？

使用循环，它会变成这样。

import numpy as np

N=3
d=3 
K=2 

X=np.eye(N)
y=np.random.randint(1,K+1,N)
M=np.zeros((K,d))
for i in np.arange(0,K):
    line=X[y==i+1,:]
    if line.size==0:
        M[i,:]=np.zeros(d)
    else:
        M[i,:]=mp.mean(line,0)

提前谢谢你。

【问题讨论】：

K == N 吗？ y 的值是否唯一？
如果你显示一些代码会很酷。
如果您将输入写成我们可以运行的短代码 sn-p 会有所帮助，例如y = np.random.random_integers(low=1, high=K, size=N) 和 X = ? 并定义 K、N 等等。添加一些使用循环计算输出的代码会更有帮助，这样我们就可以清楚地了解您要实现的目标。
X 在您的实际情况中也是单位矩阵还是它有随机数？
与@Divakar 类似的问题：由于您将X 的元素与索引进行比较，X 通常可以包含什么，它的dtype 是什么？感谢您更新问题，现在更清楚了

标签： python numpy vectorization

【解决方案1】：

代码基本上是从 X 中收集特定行并添加它们，我们在 np.add.reduceat 中有一个 NumPy 内置函数。因此，重点关注这一点，以矢量化方式解决它的步骤可能如下所列 -

# Get sort indices of y
sidx = y.argsort()

# Collect rows off X based on their IDs so that they come in consecutive order
Xr = X[np.arange(N)[sidx]]

# Get unique row IDs, start positions of each unique ID
# and their counts to be used for average calculations
unq,startidx,counts = np.unique((y-1)[sidx],return_index=True,return_counts=True)

# Add rows off Xr based on the slices signified by the start positions
vals = np.true_divide(np.add.reduceat(Xr,startidx,axis=0),counts[:,None])

# Setup output array and set row summed values into it at unique IDs row positions
out = np.zeros((K,d))
out[unq] = vals

【讨论】：

【解决方案2】：

这解决了问题，但创建了一个中间 K×N 布尔矩阵，并且不使用内置的均值函数。在某些情况下，这可能会导致更差的性能或更差的数值稳定性。我让类标签的范围从0 到K-1 而不是1 到K。

# Define constants
K,N,d = 10,1000,3

# Sample data
Y = randint(0,K-1,N) #K-1 to omit one class to test no-examples case
X = randn(N,d)

# Calculate means for each class, vectorized 

# Map samples to labels by taking a logical "outer product"
mark = Y[None,:]==arange(0,K)[:,None] 

# Count number of examples in each class    
count = sum(mark,1)

# Avoid divide by zero if no examples
count += count==0

# Sum within each class and normalize
M = (dot(mark,X).T/count).T

print(M, shape(M), shape(mark))

【讨论】：