sumjuliabroadcastin-placearrayofarrays

Julia "strange" behaviour using fill() and .+=

I observe an unexpected behaviour for ".+=" in my code (it's probably just me, I'm rather new to Julia). Consider the following example:

julia> b = fill(zeros(2,2),1,3)
       1×3 Array{Array{Float64,2},2}:
       [0.0 0.0; 0.0 0.0]  [0.0 0.0; 0.0 0.0]  [0.0 0.0; 0.0 0.0]

julia> b[1] += ones(2,2)
       2×2 Array{Float64,2}:
       1.0  1.0
       1.0  1.0

julia> b
       1×3 Array{Array{Float64,2},2}:
       [1.0 1.0; 1.0 1.0]  [0.0 0.0; 0.0 0.0]  [0.0 0.0; 0.0 0.0]

julia> b[2] .+= ones(2,2)
       2×2 Array{Float64,2}:
       1.0  1.0
       1.0  1.0

julia> b
       1×3 Array{Array{Float64,2},2}:
       [1.0 1.0; 1.0 1.0]  [1.0 1.0; 1.0 1.0]  [1.0 1.0; 1.0 1.0]

As it can be seen, the last command changed not only the value of b[2] but also of b[3], while b[1] remains the same as before (*), as we can confirm running:

julia> b[2] .+= ones(2,2)
       2×2 Array{Float64,2}:
       2.0  2.0
       2.0  2.0

julia> b
       1×3 Array{Array{Float64,2},2}:
       [1.0 1.0; 1.0 1.0]  [2.0 2.0; 2.0 2.0]  [2.0 2.0; 2.0 2.0]

Now, using simply "+=" instead I can obtain the behaviour I would have expected also for ".+=", that is:

julia> b = fill(zeros(2,2),1,3); b[2]+=ones(2,2); b
       1×3 Array{Array{Float64,2},2}:
       [0.0 0.0; 0.0 0.0]  [1.0 1.0; 1.0 1.0]  [0.0 0.0; 0.0 0.0]

Can anyone explain me why does it happen? I can use of course just +=, or maybe something different from an Array of Arrays, but since I'm striving for speed (I have a code that needs to perform these operations millions of times, and on much larger matrices) and .+= is considerably faster I would like to understad if I can still exploit this feature. Thank you all in advance!

EDIT: (*) apparently only because b[1] was not zero. If I run:

julia> b = fill(zeros(2,2),1,3); b[2]+=ones(2,2);
julia> b[1] .+= 10 .*ones(2,2); b
       [10.0 10.0; 10.0 10.0]  [1.0 1.0; 1.0 1.0]  [10.0 10.0; 10.0 10.0]

you can see that only the zero-values are changed. This beats me.

Solution

  • This happens because of the combination of several factors. Let's try and make things clearer.

    First, b = fill(zeros(2,2),1,3) does not create a new zeros(2,2) for each element of b; instead it creates one 2x2 array of zeros, and sets all elements of b to that unique array. In short, this line behaves equivalently to

    z = zeros(2,2)
b = Array{Array{Float64,2},2}(undef, 1, 3)
for i in eachindex(b)
    b[i] = z
end

    therefore, modifying z[1,1] or any of the b[i,j][1,1] would modify the other values as well. To illustrate this:

    julia> b = fill(zeros(2,2),1,3)
1×3 Array{Array{Float64,2},2}:
 [0.0 0.0; 0.0 0.0]  [0.0 0.0; 0.0 0.0]  [0.0 0.0; 0.0 0.0]

# All three elements are THE SAME underlying array
julia> b[1] === b[2] === b[3]
true

# Mutating one of them mutates the others as well
julia> b[1,1][1,1] = 42
42

julia> b
1×3 Array{Array{Float64,2},2}:
 [42.0 0.0; 0.0 0.0]  [42.0 0.0; 0.0 0.0]  [42.0 0.0; 0.0 0.0]

    Second, b[1] += ones(2,2) is equivalent to b[1] = b[1] + ones(2,2). This implies a succession of operations:

    1. a new array (let's call it tmp) is created to hold the sum of b[1] and ones(2,2)
    2. b[1] is rebound to that new array, thereby losing its connection to z (or all other elements of b.

    This is a variation on the classical theme that although both involve = signs in their notations, mutation and assignment are not the same thing. Again, to illustrate:

    julia> b = fill(zeros(2,2),1,3)
1×3 Array{Array{Float64,2},2}:
 [0.0 0.0; 0.0 0.0]  [0.0 0.0; 0.0 0.0]  [0.0 0.0; 0.0 0.0]

# All elements are THE SAME underlying array
julia> b[1] === b[2] === b[3]
true

# But that connection is lost when `b[1]` is re-bound (not mutated) to a new array
julia> b[1] = ones(2,2)
2×2 Array{Float64,2}:
 1.0  1.0
 1.0  1.0

# Now b[1] is no more the same underlying array as b[2]
julia> b[1] === b[2]
false

# But b[2] and b[3] still share the same array (they haven't be re-bound to anything else)
julia> b[2] === b[3]
true

    Third, b[2] .+= ones(2,2) is a different beast altogether. It does not imply re-binding anything to a newly created array; instead, it mutates the array b[2] in place. It effectively behaves like:

    for i in eachindex(b[2])
    b[2][i] += 1  # or b[2][i] = b[2][i] + 1
end

    Neither b itself nor even b[2] is re-bound to anything, only elements of it are modified in place. And in your example this affects b[3] as well, since both b[2] and b[3] are bound to the same underlying array.