Sum of float64 vector gives slightly incorrect answer

General Usage Performance
41

This is normal, use eg isapprox with a low tolerance (eg sqrt(eps()), 50*eps(), depending on algorithm, YMMV) for testing.

42

This is apparently an old idea, but it turns out not to help in general. It was analyzed in Higham (1993) for a variety of orderings, who wrote:

… this “+/–” method does not reduce the amount of cancellation—it simply concentrates all the cancellation in one step. … In summary, the +/– method appears to have no advantages over the other methods considered here, and in cases where there is heavy cancellation in the sum it can be expected to be the least accurate method.

3 Likes
43

Yes, thanks for pointing this out.

However, it can help when combined with compensation techniques. For example, Priest’s doubly compensated summation is guaranteed to achieve a 2 ulp accuracy by sorting its operands (in contrast to most other doubly compensated algorithms, which only triple the available precision).

2 Likes
44

You may like using KahanSummation.jl


julia> using KahanSummation

julia>  a = [0.1, 0.1, 100.0];

julia> sum(a), sum(reverse(a))
(100.2, 100.19999999999999)

julia> sum_kbn(a), sum_kbn(reverse(a))
(100.2, 100.2)
1 Like
45

I know this has been a long and wide-ranging discussion, but I feel like it’s incomplete without linking to Stefan’s post about the ability to sum a list of 2046 elements to any floating point value you so choose just based upon its ordering. Because sometimes it’s helpful to see the absolute worst case while you’re chasing the best.

6 Likes