matrix - Julia 中 cumsum 的倒数

Question

矩阵Y定义为

Y = cumsum(cumsum(X,dims=1), dims=2)

例如，

julia> X = [1 4 2 3; 2 4 5 2; 4 3 4 1; 2 5 4 2];
julia> Y = cumsum(cumsum(X,dims=1), dims=2)
4x4 Matrix{Int64}:
1  5  7 10
3 11 18 23
7 18 29 35
9 25 40 48

我想从 Y 复制矩阵 X。函数 diff 似乎很有帮助。但是，如下所示，我们无法重现 X 的第一行和第一列。

julia> diff(diff(y, dims=1), dims=2)
3x3 Matrix{Int64}:
4 5 2
3 4 1
5 4 2

所以，我连接零。然后，它起作用了。

julia> y00 = vcat(zeros(5)',hcat(zeros(4), y))
5x5 Matrix{Int64}:
0  0  0  0  0
0  1  5  7 10
0  3 11 18 23
0  7 18 29 35
0  9 25 40 48

julia> diff(diff(y00, dims=1), dims=2)
4x4 Matrix{Int64}:
1  5  7 10
3 11 18 23
7 18 29 35
9 25 40 48

但我认为连接需要时间和内存。
从 Y 重现 X 有什么更好的主意吗？

上下文

我想将上面的矩阵 X 和 Y 展开为任意维数组。例如，我想根据给定的三维数组重建一个三维数组 X

Y = cumsum( cumsum( cumsum(X, dims=1), dims=2), dims=3)

Answer 1

当同时需要速度和简洁时，很难击败像 Tullio.jl 这样强大的 Julia 包.这是一个比@DanGetz 最快的解决方案快大约 4 倍的单行代码。

using Tullio
cumdiff(Y) = @tullio X[i,j] = Y[i,j] - Y[i,j-1] - Y[i-1,j] + Y[i-1,j-1]

使用 100×100 矩阵进行基准测试可得出:

X = rand(0:100,100,100)
Y = cumsum(cumsum(X,dims=1), dims=2)
@btime cumdiff($Y)
@btime decumsum3($Y)
  4.957 μs (17 allocations: 464 bytes)
  21.300 μs (2 allocations: 78.17 KiB)

修复:上面的代码使用了预定义的 X 而不是创建一个新的。这是在下面修复的，加速更像是 3.5 倍而不是 4 倍。

function cumdiff(Y)
    X = similar(Y)
    X[1] = Y[1]
    for i = 2:size(Y,1) X[i,1] = Y[i,1] - Y[i-1,1] end
    for j = 2:size(Y,2) X[1,j] = Y[1,j] - Y[1,j-1] end
    @tullio X[i,j] = Y[i,j] - Y[i,j-1] - Y[i-1,j] + Y[i-1,j-1]
end

@btime cumdiff($Y)
@btime decumsum3($Y)
  6.000 μs (4 allocations: 78.23 KiB)
  21.300 μs (2 allocations: 78.17 KiB)