The function "sum"

[ami...@gmail.com]

Dear all,

Some (serious?) issue was raised on a mailing list concerning the "sum" functions involved in computing the sum of elements of an array.

It seems that we obtain different results for the three following procedures:

- using sum over the whole matrix
- using sum over the columns/rows of the matrix and summing these results within a loop
- computing the sum with a double loop over all the elements of the matrix

Attached is a file with the different procedures and their results for a particular matrix.

What is wrong with the implementation of the sum function?
This could have important consequences with iterative procedures.

Many thanks,

[Stefan Karpinski]

This is expected. Floating-point addition is non-associative and these codes add values in different orders. Julia's built-in sum function uses recursive pairwise summation, which is the most accurate of these three approaches but nearly as fast as linear scanning. A slower but more accurate algorithm is compensated summation. Even compensated summation is not completely accurate; completely accurate floating-point summation is tricky.

[ami...@gmail.com]

Dear Stefan,

Thank you for your quick answer. That's what I had in mind. I just thought, very naively I must admit, that sum(a) would have computed the total sum in the same order as one of the other 4 approaches.
I'll have a look at compensated summation.

Thanks again,

[Ivar Nesje]

Seems like this thread is missing a reference to the standard library function sum_kbn.

[Stefan Karpinski]

Right, wasn't sure if we kept that around or not.

[Tamas Papp]

If one needs completely accurate floating point summation, I guess they
can just convert to Rational{BigInt} and sum them -- it will be
expensive, but accurate.

But I am wondering why anyone would expect "complete accuracy" from
floating point anyway.

Best,

Tamas

On Sat, Nov 15 2014, Stefan Karpinski <ste...@karpinski.org> wrote:

> This is expected. Floating-point addition is non-associative and these
> codes add values in different orders. Julia's built-in sum function uses
> recursive pairwise summation, which is the most accurate of these three
> approaches but nearly as fast as linear scanning. A slower but more
> accurate algorithm is compensated summation

> <http://en.wikipedia.org/wiki/Kahan_summation_algorithm>. Even compensated

[Dahua Lin]

The current implementation of the sum goes in this way,

For arrays longer than a certain threshold (1024, 2048?), it recursively divide the array into two halves, and sum different halves respectively.  As the part to be summed gets shorter, it resorts to a four-way sequential summing algorithm, i.e. x[i], x[i+1], x[i+2], x[i+3] are summed to four accumulators.

I guess better implementations would be developed, and we may switch to an SIMD-based summing algorithm in future. 

So, never count on the `sum` function to sum the values in any specific order. 

Dahua