# Kahan summation in `sum`?

I was under the impression that `sum(itr)` uses Kahan summation, but after inspecting the Julia repo it appears that it calls `mapreduce(identity, add_sum, itr)`. And `add_sum` is basically just `+` with some promotions for small integer types. Does `sum(itr)` not use Kahan summation, or is there a code path that I’m missing?

Looking at the docs for `numpy.sum`, it appears that they use pairwise summation (not Kahan summation) to improve the output precision.

Long ago some code for Kahan summation used to live in `Base`, but now it’s in a separate Git repo on Github:

So `Base` doesn’t use compensated summation.

2 Likes

You may be interested by

3 Likes

Here’s the history, which includes a good high-level comparison between the naive, pairwise and Kahan summation strategies: RFC: use pairwise summation for sum, cumsum, and cumprod by stevengj · Pull Request #4039 · JuliaLang/julia · GitHub

Pairwise summation recursively divides the array in half, sums the halves recursively, and then adds the two sums. As long as the base case is large enough (here, n=128 seemed to suffice), the overhead of the recursion is negligible compared to naive summation (a simple loop). The advantage of this algorithm is that it achieves O(sqrt(log n)) mean error growth, versus O(sqrt(n)) for naive summation, which is almost indistinguishable from the O(1) error growth of Kahan compensated summation.

4 Likes

Compensated algorithms have a catch that’s not mentioned in that repo’s README, though: they usually assume that the exponent range of the relevant FP format is large enough, which may not be the case. This means that the compensated algorithms may produce `NaN` when a naive algorithm would, more correctly, produce an infinity. For example, Ogita, Rump & Oishi say (in the compensated dot product paper):

We assume that no overflow occurs, but allow underflow.

So I’m pretty sure this catch also applies to simple compensated summation.

1 Like

this doesn’t apply to compensated sumation (unless you disable subnormals). Addition of floating point numbers never underflows.

1 Like

Imprecision of sum(::Generator) · Issue #30421 · JuliaLang/julia · GitHub has a kahan summation implementation.

Yeah, you can only do the pairwise recursion if you can index into the object at arbitrary indices (which generators can’t do in general). But if you have an array, then pairwise summation is the obvious answer that balances performance and accuracy — it’s why others do the same thing.

2 Likes

Just to clarify, it sounds like `sum` performs pairwise summation for objects that support linear indexing?

Yes, but not just linear indexing and not just `sum` — sufficiently large AbstractArrays (and lazy broadcasts) use a recursive divide and conquer strategy for most reductions. The size cutoff varies by operator.

This is why we have both `reduce` (whose order of traversal is unspecified) and `foldl`/`foldr`.

4 Likes