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

【问题标题】：np.bincount for 1 line, vectorized multidimensional averaging1 行的 np.bincount，矢量化多维平均
【发布时间】：2016-03-31 23:09:58
【问题描述】：

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

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

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

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

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，有些东西不能直接翻译，例如以相反顺序散开的冒号运算符。所以我的问题是如何使用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 负责人来说，这是一个更切实的实现方式来消化。欢迎任何更快的建议:)

【问题讨论】：

看看meshgrid，它有助于生成这样的二维数组，这里有一个简单的例子：stackoverflow.com/questions/36300023/…numpy对这些二维数组的操作自动是vectorized（没有循环）。
正是我的想法。我可以使用np.tile（或您建议的np.meshgrid）轻松地制作索引数组。我没有看到的部分是如何在numpy 设置中执行，这就是accumarray 在上面的示例中所做的。您不能直接将np.meshgrid 的输出输入到np.bincount，因为np.bincount 只需要一维数组。 np.ravel 将破坏进入np.bincount 的索引信息，所以没有用。是否可以在 Cython 中显式地执行循环……如果数组操作被证明过于复杂，人们会建议这样做吗？无论如何感谢您的建议。

标签： python arrays matlab numpy cython

【解决方案1】：

一种方法是将 3 个“索引”转换为线性索引，然后应用 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 中添加。

【讨论】：