Accuracy issue with summation of `@view`

[sjkelly]

This was originally noticed in:
#28848 (comment)

I don't know if this is somewhat related to this issue, but I also found an accuracy issue which I think is more troublesome because it also affects DataFrames.jl

here is the MVP

using Distributions
using Random

# It doesn't happen always, here is a seed where this happens.
rng = MersenneTwister(630);
v = rand(rng, Normal(zero(Float32), one(Float32)), 1000)
sa = @view v[collect(1:end)]

# View (as SubArray) vs Vector
sum(sa) ≈ sum(v) # false
# They are different! and worse part is that the view version is less accurate! (according to Kahan compensated summation)

IndexStyle(v) isa IndexLinear # true
IndexStyle(sa) isa IndexCartesian #true

# They are dispatched to different implementations in base/reduce.jl > _mapreduce 

Another MWE is:

julia> using Statistics, Random

julia> Random.seed!(1);

julia> x = rand(Float32, 10^6);

julia> mean(x)
0.49974534f0

julia> mean(@view x[:])
0.49974534f0

julia> mean(@view x[collect(1:end)])
0.49974778f0

@mbauman points out that the latter in this MWE hits the naive iterator: foldl(+, x)/length(x) which leads to the poor floating point accuracy.

cc @el-oso

[mikmoore]

Would probably be resolved by #52397 (or easily resolved after it).

[stevengj]

Yes, this line should resolve it: https://github.com/JuliaLang/julia/pull/52397/files#diff-97ad8b900c30b39398c32d78754c404f4a6df9a153d45284b2c11ab54837deb9R524 (changing the fallback for CartesianIndex arrays)

[el-oso]

Thanks for the fast response!

[el-oso]

I just tested master and got this

 julia> sum(sa) ≈ sum(v)
true

Great!

[stevengj]

Note (a) that the new pairwise sum is not merged yet and (b) the ≈ operator only checks whether half of the significant digits match.

[mbauman]

Interestingly, the particular MWEs in the OP aren't fully resolved by either pairwise algorithms or with a unification of the IndexLinear/IndexCartesian implementations. These views are using arbitrary Arrays of collected indices, which means that — regardless of the implementation — if we allow sum(v) to SIMD, then sum(@view v[collect(begin:end)]) will give a different answer because it cannot SIMD.

Pairwise re-associations will indeed help with larger arrays (the default pairwise split is greater than the 1000 elements in the first example), but the additional SIMD reassociations will still be likely to produce different values. In addition to views like these, this could also theoretically manifest in sum(A; dims=1) vs. sum(A'; dims=2)' for a Matrix{Float64} A; it's possible for the former to SIMD in a way that effectively adds more pairwise splits, while the latter cannot.