mean! may return incorrect results

Question

I recently had the following question regarding the function mean! of the Statistics.jl package.

The following bug is reported regarding the behavior of mean!. As indicated, the mean! function does not properly consider that its arguments may alias each other. In some such cases the result from mean! is not correct:

julia> let a = [1 2 3]
       mean!(a, a)
   end
1×3 Array{Int64,2}:
0  0  0

julia> let a = [1 2 3]
       mean!(copy(a), a)
   end
1×3 Array{Int64,2}:
1  2  3

julia> versioninfo()
Julia Version 1.5.3
Commit 788b2c77c1 (2020-11-09 13:37 UTC)
Platform Info:
 OS: macOS (x86_64-apple-darwin18.7.0)
 CPU: Intel(R) Core(TM) i9-9980HK CPU @ 2.40GHz
 WORD_SIZE: 64
 LIBM: libopenlibm
 LLVM: libLLVM-9.0.1 (ORCJIT, skylake)

However, I think that this behavior is normal since viewing the definition of mean!, the result of the operation mean!(r, v) is written in r. Therefore it seems logical to me that if you use the same object as variable r and as variable v, the result is unpredictable.

I have seen that this also happens with the sum! function.

Someone can tell me if I am right or indeed, there is something that I am not understanding.

Answer 1

mean! behaves the way you observe because it internally calls sum!.

Now sum! behaves this way for the following reason. It was designed to perform summation without making any allocations. Therefore the first thing sum! does is initializing the target vector to 0 (a neutral element of summation). After this is done your a vector contains only 0s, and thus later you get all 0s also.

However, indeed it would make sense that the sum! (and similar) functions docstring should mention that the target should not alias with the source. Here is another example of the same you have observed:

julia> x = [1 2 3
            4 5 6
            7 8 9]
3×3 Matrix{Int64}:
 1  2  3
 4  5  6
 7  8  9

julia> y = view(x, :, 1)
3-element view(::Matrix{Int64}, :, 1) with eltype Int64:
 1
 4
 7

julia> sum!(y, x)
3-element view(::Matrix{Int64}, :, 1) with eltype Int64:
  5
 11
 17

julia> x
3×3 Matrix{Int64}:
  5  2  3
 11  5  6
 17  8  9