沿轴的 2D 数组的 Numpy 总和 = 1,浮动范围

Numpy sum of 2D array along axis=1, floating range

我想在第二个轴上对二维数组求和,但范围是可变的。没有向量化它是:`

import numpy as np

nx = 3
ny = 5
a = np.ones((nx, ny))
left_bnd  = np.array([0, 1, 0])
right_bnd = np.array([2, 2, 4])

b = np.zeros(nx)
for jx in range(nx):
    b[jx] = np.sum(a[jx, left_bnd[jx]: right_bnd[jx]])

print(b)

输出 b 是 [2. 1. 4.] 我喜欢矢量化循环,有点

b = np.sum(a[:, left_bnd[:]: right_bnd[:], axis=1)

为了加快计算速度,因为我的“n”一般是1e6几个。不幸的是我找不到合适的工作语法。

在 for 循环中手动求和的 jitted numba 实现大约快 100 倍。在 numba 函数内使用带有切片的 np.sum 只有一半的速度。此解决方案假定所有切片都在有效范围内。

为基准测试生成足够大的样本数据

import numpy as np
import numba as nb
np.random.seed(42)    # just for reproducibility

n, m = 5000, 100
a = np.random.rand(n,m)
bnd_l, bnd_r = np.sort(np.random.randint(m+1, size=(n,2))).T

numba 合影。请确保通过 运行 函数对编译的热代码进行基准测试至少两次。

@nb.njit
def slice_sum(a, bnd_l, bnd_r):
    b = np.zeros(a.shape[0])
    for j in range(a.shape[0]):
        for i in range(bnd_l[j], bnd_r[j]):
            b[j] += a[j,i]
    return b
slice_sum(a, bnd_l, bnd_r)

输出

# %timeit 1000 loops, best of 5: 297 µs per loop
array([ 4.31060848, 35.90684722, 38.03820523, ..., 37.9578962 ,
        3.61011028,  6.53631388])

在 python 循环中使用 numpy(这是一个很好、简单的实现)

b = np.zeros(n)
for j in range(n):
    b[j] = np.sum(a[ j, bnd_l[j] : bnd_r[j] ])
b

输出

# %timeit 10 loops, best of 5: 29.2 ms per loop
array([ 4.31060848, 35.90684722, 38.03820523, ..., 37.9578962 ,
        3.61011028,  6.53631388])

验证结果是否相等

np.testing.assert_allclose(slice_sum(a, bnd_l, bnd_r), b)

这是一个纯 numpy 解决方案,其速度接近已发布的 numba 解决方案。它利用 reduceat 但设置非常复杂。

def slice_sum_np(a, left_bnd, right_bnd):
    nx, ny = a.shape
    linear_indices = np.c_[left_bnd, right_bnd] + ny * np.arange(nx)[:,None]
    sums = np.add.reduceat(a.ravel(), linear_indices.ravel())[::2]
    # account for reduceat special case
    sums[left_bnd >= right_bnd] = 0
    return sums