np.bincount 1 行，矢量化多维平均

Question

我正在尝试使用 numpy 对操作进行矢量化，我在我分析过的 python 脚本中使用它，发现此操作是瓶颈，因此需要进行优化，因为我将运行很多次。

操作是在两部分的数据集上进行的。首先，一大组 (n) 不同长度的一维向量（最大长度 Lmax），其元素是从 1 到 maxvalue 的整数。向量集排列在大小为 (num_samples,Lmax) 的二维数组 data 中，每行中的尾随元素都归零。第二部分是一组标量浮点数，一个与每个向量相关联，我有一个计算并且取决于它的长度和每个位置的整数值。标量集被制成一维数组 Y，大小为 num_samples。

所需的操作是 形成 n 个样本的 Y 的平均值，作为 (value,position along length,length) 的函数。

整个操作可以在 matlab 中使用 accumarray 函数向量化：通过使用 3 个与 data 大小相同的二维数组，其元素是相应的值、位置和所需最终数组的长度索引：

sz_Y = num_samples;
sz_len = Lmax 
sz_pos = Lmax 
sz_val = maxvalue
ind_len = repmat( 1:sz_len      ,1         ,sz_samples);
ind_pos = repmat( 1:sz_pos      ,sz_samples,1         );
ind_val = data
ind_Y   = repmat((1:sz_Y)',1         ,Lmax      );
copiedY=Y(ind_Y);
mask = data>0; 
finalarr=accumarray({ind_val(mask),ind_pos(mask),ind_len(mask)},copiedY(mask), [sz_val sz_pos sz_len])/sz_val;

我希望用 np.bincounts 来模拟这个实现。但是，np.bincounts 在两个相关方面不同于 accumarray：

两个参数必须具有相同的一维大小，并且

没有选择输出数组形状的选项。

在 accumarray 的上述用法中，索引列表 {ind_val(mask),ind_pos(mask),ind_len(mask)} 是用作索引元组的 1x3 数组的一维元胞数组，而在 np.bincounts 中它必须是一维据我了解，标量。我希望 np.ravel 可能有用，但我不确定如何在这里使用它来做我想做的事。我从 matlab 来到 python，有些东西不能直接翻译，例如冒号运算符，以与 ravel 相反的顺序 ravels。所以我的问题是 我如何使用 np.bincount 或任何其他 numpy 方法来实现此操作的高效 python 实现 .

编辑： 为了避免浪费时间：对于这些具有复杂索引操作的多维索引问题，推荐的方法是只使用 cython 来显式实现循环吗？

EDIT2: 替代 Python 实现我刚刚想出。
这是一个重型 ram 解决方案：

首先预计算：
使用长度的索引单位（即长度 1 =0）制作一个 4D 布尔数组，大小为 (num_samples,Lmax+1,Lmax+1,maxvalue) ，保存 Y.

中每个值满足条件的位置

ALLcond=np.zeros((num_samples,Lmax+1,Lmax+1,maxvalue+1),dtype='bool')
for l in range(Lmax+1):
    for i in range(Lmax+1):
        for v in range(maxvalue+!):
            ALLcond[:,l,i,v]=(data[:,i]==v) & (Lvec==l)`

其中 Lvec=[len(row) for row in data]。然后使用 np.where 获取这些索引并初始化一个 4D 浮点数组，您将在其中分配 Y:

的值

[indY,ind_len,ind_pos,ind_val]=np.where(ALLcond)
Yval=np.zeros(np.shape(ALLcond),dtype='float')

现在在我必须执行操作的循环中，我用两行计算它：

Yval[ind_Y,ind_len,ind_pos,ind_val]=Y[ind_Y]
Y_avg=sum(Yval)/num_samples

这使直接循环实现的速度提高了 4 倍左右。我期待更多。也许，这是 Python 头脑消化的更具体的实现。欢迎任何更快的建议:)

Answer 1

一种方法是将 3 "indices" 转换为线性索引，然后应用 bincount。 Numpy 的 ravel_multi_index 与 MATLAB 的 sub2ind 本质上是一样的。所以移植的代码可能是这样的：

shape = (Lmax+1, Lmax+1, maxvalue+1)
posvec = np.arange(1, Lmax+1)

ind_len  = np.tile(Lvec[:,None], [1, Lmax])
ind_pos  = np.tile(posvec,       [n,    1])
ind_val  = data
Y_copied = np.tile(Y[:,None],    [1, Lmax])

mask = posvec <= Lvec[:,None]  # fill-value independent
lin_idx = np.ravel_multi_index((ind_len[mask], ind_pos[mask], ind_val[mask]), shape)
Y_avg = np.bincount(lin_idx, weights=Y_copied[mask], minlength=np.prod(shape)) / n
Y_avg.shape = shape

这是假设 data 的形状为 (n, Lmax)，Lvec 是 Numpy 数组等。您可能需要稍微调整代码以摆脱逐一错误。

有人可能会争辩说 tile 操作不是很有效，也不是很 "numpythonic"。带 broadcast_arrays 的东西可能不错，但我想我更喜欢这种方式：

shape = (Lmax+1, Lmax+1, maxvalue+1)
posvec = np.arange(1, Lmax+1)

len_idx  = np.repeat(Lvec, Lvec)
pos_idx  = np.broadcast_to(posvec, data.shape)[mask]
val_idx  = data[mask]
Y_copied = np.repeat(Y, Lvec)

mask = posvec <= Lvec[:,None]  # fill-value independent
lin_idx = np.ravel_multi_index((len_idx, pos_idx, val_idx), shape)
Y_avg = np.bincount(lin_idx, weights=Y_copied, minlength=np.prod(shape)) / n
Y_avg.shape = shape

注意 broadcast_to 是在 Numpy 1.10.0 中添加的。

np.bincount 1 行，矢量化多维平均

np.bincount for 1 line, vectorized multidimensional averaging

python

arrays

matlab

numpy

cython