沿具有不同 bin 边缘的轴有效计算 3D numpy 数组的直方图

Question

问题描述

我有一个 3D numpy 数组，表示为 data，形状为 N x R x C，即 N 个样本，R 行和 C 列。我想为样本和行的每个组合获取沿列的直方图。然而，固定长度 S 的 bin 边缘（参见 numpy.histogram 中的参数 bins）在不同行将不同，但在样本之间共享。考虑这个例子来说明，对于第一个样本 (data[0])，其第一行的 bin 边缘序列与其第二行的 bin 边缘序列不同，但与第二个样本的第一行的 bin 边缘序列相同 (data[1] ）。因此，所有 bin 边缘序列都存储在形状为 R x S 的 2D numpy 数组中，表示为bin_edges。

我的问题是如何有效地计算直方图？

一个有效但缓慢的解决方案

使用numpy.histogram，我能够想出一个可行但相当慢的解决方案，如下面的代码 sn-p 所示

```
Get dummy data

    N: number of samples
    R: number of rows (or kernels)
    C: number of columns (or pixels)
    S: number of bins
```
import numpy as np

N, R, C, S = 100, 50, 1000, 10
data = np.random.randn(N, R, C)

# for each row/kernel, pool pixels of all samples
poolsamples = np.swapaxes(data, 0, 1).reshape(R, -1)
# use quantiles as bin edges
percentiles = np.linspace(0, 100, num=(S + 1))
bin_edges = np.transpose(np.percentile(poolsamples, percentiles, axis=1))


```
A working but slow solution of getting histograms along column
```
hist = np.empty((N, R, S))
for idx in np.arange(R):
    bin_edges_i = bin_edges[idx, :]
    counts = np.apply_along_axis(
        lambda a: np.histogram(a, bins=bin_edges_i)[0],
        1, data[:, idx, :])
    hist[:, idx, :] = counts

可能的方向

• 花式 numpy 重塑以避免使用 for 循环
• 此问题源于为通过训练的神经网络转发的每个图像提取低端特征。因此，如果直方图的提取可以嵌入到 TensorFlow 图中，最终在 GPU 上进行，那就再理想不过了！
• 我注意到一个 python 包fast-histogram 声称比numpy.histogram 快7-15 倍。然而，一维直方图函数只能采用 bin 数量而不是实际 bin 位置
• numexpr?

我很想听听任何意见！提前致谢！

Answer 1

使用2D 版本的np.searchsorted : searchsorted2d -

def vectorized_app(data, bin_edges):
    N, R, C = data.shape
    a = np.sort(data.reshape(-1,C),1)
    b = np.repeat(bin_edges[None],N,axis=0).reshape(-1,bin_edges.shape[-1])

    idx = searchsorted2d(a,b)
    idx[:,0] = 0
    idx[:,-1] = a.shape[1]
    out = (idx[:,1:] - idx[:,:-1]).reshape(N,R,-1)
    return out

运行时测试-

In [591]: N, R, C, S = 100, 50, 1000, 10
     ...: data = np.random.randn(N, R, C)
     ...: 
     ...: # for each row/kernel, pool pixels of all samples
     ...: poolsamples = np.swapaxes(data, 0, 1).reshape(R, -1)
     ...: # use quantiles as bin edges
     ...: percentiles = np.linspace(0, 100, num=(S + 1))
     ...: bin_edges = np.transpose(np.percentile(poolsamples, percentiles, axis=1))
     ...: 

In [592]: %timeit org_app(data, bin_edges)
1 loop, best of 3: 481 ms per loop

In [593]: %timeit vectorized_app(data, bin_edges)
1 loop, best of 3: 224 ms per loop

In [595]: np.allclose(org_app(data, bin_edges), vectorized_app(data, bin_edges))
Out[595]: True

超过2x 加速。

仔细观察发现，所提出的矢量化方法的瓶颈在于排序本身 -

In [594]: %timeit np.sort(data.reshape(-1,C),1)
1 loop, best of 3: 194 ms per loop

我们需要这种排序才能使用searchsorted。