Array cumulation

General Usage
1

Hello,

I am trying to, quickly, accumulate an array without over-riding previous allocations. In particular consider the following example:

B = zeros(2,2);
A = [0.1 0.3; 0.2 0.5];
idx = [1 1; 1 1];

I would like to add to B, according to idx all the elements in A. Since idx all maps to the first element, I would hope to get:

B = [1.1 0; 0 0];

but with the following command I get:

@views B[idx] .+= A
B = [0.5 0; 0 0];

I.e. only the last allocation is recorded, overwriting everything else. Is there a way out of this?

Thank you!

2

You may want to opt for a for-loop:

for (i, Aᵢ) in zip(idx, A)
    B[i] += Aᵢ
end

julia> B
2×2 Array{Float64,2}:
 1.1  0.0
 0.0  0.0
1 Like
3

I was trying to avoid a loop because I have to do this many many times. I was under the impression there was a quick way of doing this, but maybe not?

4

For-loops are very fast in Julia. If you need to do it many times, you can just write a small function like so:

function cumulate!(B, A, idx)
    for (i, Aᵢ) in zip(idx, A)
        B[i] += Aᵢ
    end
end

julia> cumulate!(B, A, idx)

julia> B
2×2 Array{Float64,2}:
 1.1  0.0
 0.0  0.0
4 Likes
5

I think Julia is specifically preventing you from doing this. It wants B .+= A to mean B .= B .+ A with the right hand side evaluated first, and since this contains views of the same data as on the left, it makes copies to be safe. I’m not sure I can always predict what it’s going to do in such cases, but I do think that anything which relies on broadcasting being done sequentially (or worse, in a particular order) is fishy. In this case, one way to check is by timing it – the 15MB allocation tells you it’s in fact making two copies, although I’m not sure what the second one is:

julia> N = 10^3; A = rand(N,N); B = zeros(N,N); idx = ones(Int, N,N);

julia> @time B[idx] .= A;
  0.003550 seconds (6 allocations: 208 bytes)

julia> @time B[idx] .+= A;  # has B[idx] on the right
  0.007897 seconds (7 allocations: 7.630 MiB)

julia> @time @views B[idx] .+= A;  # has aliased view on the right
  0.007404 seconds (12 allocations: 15.259 MiB)

julia> @time B[idx];
  0.004902 seconds (2 allocations: 7.629 MiB)

julia> @time @view B[idx];
  0.000582 seconds (3 allocations: 128 bytes)

As @stillyslalom says, if you do want it done sequentially, then you want a loop. There is no reason this need be slower – broadcasting is itself ultimately just loops.

2 Likes