高效的 2d cumsum答案

【问题标题】：Efficient 2d cumsum高效的 2d cumsum
【发布时间】：2015-10-19 01:03:13
【问题描述】：

假设我有一个这样的数组

>>> a = np.arange(1,8).reshape((1,-1))
>>> a
array([[1, 2, 3, 4, 5, 6, 7]])

我想为a 中的每个项目创建一个“接下来的 4 个项目的总和”。也就是说，我的预期输出是

1,       2,      3, 4, 5, 6, 7, 8
1+2,     2+3,     ...
1+2+3    2+3+4    ...
1+2+3+4  2+3+4+5  ...

即一个矩阵，包含

1, 2, 3, 4, 5, 0, 0, 0
3, 5, 7, 9, 11,0, 0, 0
6, 9, 12,15,18,0, 0, 0
10,14,18,21,26,0, 0, 0

由于无法对最后 3 项正确执行 cumsum 操作，我希望那里有一个 0。我知道怎么做一个cumsum。其实数组是

a[:4].cumsum().reshape((-1,1)); a[1:5].cumsum().reshape((-1,1))...

水平堆叠。但是，我不知道如何以有效的方式做到这一点。这样做的好的矢量化 numpy 方式是什么？我也对scipy 软件包持开放态度，只要它们在效率或可读性方面优于numpy。

【问题讨论】：

标签： python arrays numpy scipy cumsum

【解决方案1】：

一种可能的方法是将滚动窗口方法与cumsum()结合使用。

例如：

from numpy.lib.stride_tricks import as_strided

a = np.arange(1, 9) # the starting array
slice_length = 4

那么你可以写：

arr = as_strided(a, (slice_length, len(a)), (a.strides[0], a.strides[0])).cumsum(axis=0)

这可以帮助您完成大部分工作，但要填写剩余的 0 值，您可以使用 slice 和 assign 来获得所需的输出：

arr[:, (1-slice_length):] = 0

那么你就有了数组：

>>> arr
array([[ 1,  2,  3,  4,  5,  0,  0,  0],
       [ 3,  5,  7,  9, 11,  0,  0,  0],
       [ 6,  9, 12, 15, 18,  0,  0,  0],
       [10, 14, 18, 22, 26,  0,  0,  0]])

我不知道是否有任何方法可以使用 NumPy 中的一种矢量化方法（即没有切片）准确地产生您想要的输出。（accumulateat，有点像reduceat，添加到 NumPy 的 ufunc 中可能是一件有趣的事情……）

【讨论】：

现在，我比另一个更喜欢这个答案，因为我可以简单地做其他计算，例如将cumsum 更改为cumprod。但是，对于不同的起始数组，它似乎会分解，例如a = np.ones((1,8))。不知道发生了什么，需要先了解一下strides。
strides 只是 NumPy 必须在内存中跳转以达到行或列中的下一个值的字节数。我的代码 sn-p 适用于 1D 输入 - 如果 a 是 2D 的一行，您需要将其更改为 as_strided(a, (slice_length, a.shape[1]), (a.strides[1], a.strides[1]))。

【解决方案2】：

您可以使用一种称为summed area table 的技术的更简单变体有效地进行计算，该技术在图像处理应用程序中也称为积分图像。首先，您计算并存储总面积表，即第一行的完整总和，并在前面添加 0：

a = np.arange(1, 8)
cs = np.concatenate(([0], np.cumsum(a)))

您现在可以将“下一个n 项目的累积和”创建为cs[:n] - cs[:-n]：

>>> for n in range(1, 5):
...     print n, '-->', (cs[n:] - cs[:-n])[:4]
...
1 --> [1 2 3 4]
2 --> [3 5 7 9]
3 --> [ 6  9 12 15]
4 --> [10 14 18 22]

您需要将它们正确地排列成您想要的形状，但是一旦完成原始计算，您就可以用一个减法计算输出的每个项目，这几乎是可以达到的效率。

【讨论】：

【解决方案3】：

你可以像这样使用broadcasting -

In [53]: a
Out[53]: array([ 4, 13,  4, 18,  1,  2, 11, 15])

In [54]: WSZ = 4 # Window size

In [55]: idx = np.arange(WSZ)[:,None] + np.arange(a.size-WSZ+1) # Broadcasted indices

In [56]: a[idx].cumsum(axis=0) # Index into "a" & perform cumsum along axis-0
Out[56]: 
array([[ 4, 13,  4, 18,  1],
       [17, 17, 22, 19,  3],
       [21, 35, 23, 21, 14],
       [39, 36, 25, 32, 29]], dtype=int32)

如果需要，用零填充 -

In [57]: np.lib.pad(a[idx].cumsum(0),((0,0),(0,WSZ-1)),'constant',constant_values=0)
Out[57]: 
array([[ 4, 13,  4, 18,  1,  0,  0,  0],
       [17, 17, 22, 19,  3,  0,  0,  0],
       [21, 35, 23, 21, 14,  0,  0,  0],
       [39, 36, 25, 32, 29,  0,  0,  0]], dtype=int32)

【讨论】：