numpy中的多维累积和答案

【问题标题】：Multidimensional cumulative sum in numpynumpy中的多维累积和
【发布时间】：2018-05-22 08:42:50
【问题描述】：

我希望能够计算一个大型 n 维 numpy 数组的累积和。最终数组中每个元素的值应该是索引大于或等于当前元素的所有元素的总和。

二维：xᶦʲ = ∑xᵐⁿ ∀ m ≥ i 和 n ≥ j

3D：xᶦʲᵏ = ∑xᵐⁿᵒ ∀ m ≥ i 且 n ≥ j 且 o ≥ k

2D 示例：

1 1 0       2  1  0
1 1 1  ->   5  3  1
1 1 1       8  5  2

1 2 3       6  5  3
4 5 6  ->  21 16  9
7 8 9      45 33 18

3D 示例：

1 1 1       3   2   1
1 1 1       6   4   2
1 1 1       9   6   3

1 1 1       6   4   2
1 1 1  ->  12   8   4
1 1 1      18  12   6

1 1 1       9   6   3
1 1 1      18  12   6
1 1 1      27  18   9

【问题讨论】：

标签： python arrays numpy

【解决方案1】：

沿最后一个轴翻转，沿同一轴进行 cumsum，将其翻转回来，最后沿倒数第二个轴向上 cumsum 直到第一个轴 -

def multidim_cumsum(a):
    out = a[...,::-1].cumsum(-1)[...,::-1]
    for i in range(2,a.ndim+1):
        np.cumsum(out, axis=-i, out=out)
    return out

示例 2D 案例运行 -

In [107]: a
Out[107]: 
array([[1, 1, 0],
       [1, 1, 1],
       [1, 1, 1]])

In [108]: multidim_cumsum(a)
Out[108]: 
array([[2, 1, 0],
       [5, 3, 1],
       [8, 5, 2]])

示例 3D 案例运行 -

In [110]: a
Out[110]: 
array([[[1, 1, 1],
        [1, 1, 1],
        [1, 1, 1]],

       [[1, 1, 1],
        [1, 1, 1],
        [1, 1, 1]],

       [[1, 1, 1],
        [1, 1, 1],
        [1, 1, 1]]])

In [111]: multidim_cumsum(a)
Out[111]: 
array([[[ 3,  2,  1],
        [ 6,  4,  2],
        [ 9,  6,  3]],

       [[ 6,  4,  2],
        [12,  8,  4],
        [18, 12,  6]],

       [[ 9,  6,  3],
        [18, 12,  6],
        [27, 18,  9]]])

【讨论】：

当然，谢谢！为了完整起见，我的 4D 用例的解决方案是：a[...,::-1].cumsum(-1)[...,::-1].cumsum(-2).cumsum(-3).cumsum(-4)
@user2685230 使用基于相同想法的通用解决方案进行编辑。
这太神奇了，真的应该被称为魔方：D

【解决方案2】：

对于那些想要左上角最小的“numpy-like”cumsum 的人：

def multidim_cumsum(a):
    out = a.cumsum(-1)
    for i in range(2,a.ndim+1):
        np.cumsum(out, axis=-i, out=out)
    return out

改编自@Divakar（感谢他！）

【讨论】：

完美，无价！

【解决方案3】：

这是一个通用的解决方案。我是按照描述来的，而不是示例，即垂直显示的顺序是自上而下而不是自下而上：

import itertools as it
import functools as ft

ft.reduce(np.cumsum, it.chain((a[a.ndim*(np.s_[::-1],)],), range(a.ndim)))[a.ndim*(np.s_[::-1],)]

或就地：

for i in range(a.ndim):
    b = a.swapaxes(0, i)[::-1]
    b.cumsum(axis=0, out=b)

【讨论】：

您的就地版本似乎有错误？ b 在后续迭代中被覆盖。
@shaunc "[...] bug? bis overwritten [...]" ???这不就是in-place的意思吗？
我的意思是——假设 a.ndim > 1。在每次迭代中，b 不依赖于先前计算的 b 值，但它会从 a 获取一个新副本，而忽略之前的工作。