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.
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:
tmp) is created to hold the sum of b[1] and ones(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.