# Kahan summation in \`sum\`?
**URL:** https://discourse.julialang.org/t/kahan-summation-in-sum/102723
**Category:** General Usage
**Tags:** float, precision
**Created:** [August 11, 2023, 4:18pm UTC](https://discourse.julialang.org/t/kahan-summation-in-sum/102723 "2023-08-11T16:18:41Z")
**Posts on this page:** 14
**Page:** 1
### Author:  [@CameronBieganek](https://discourse.julialang.org/u/CameronBieganek)
#### Post date: [August 11, 2023, 4:18pm UTC](https://discourse.julialang.org/t/kahan-summation-in-sum/102723/1 "2023-08-11T16:18:41Z")
I was under the impression that `sum(itr)` uses [Kahan summation](https://en.wikipedia.org/wiki/Kahan_summation_algorithm), but after inspecting the Julia repo it appears that [it calls](https://github.com/JuliaLang/julia/blob/eb4416b16b8a865376da5c76451a4d60516e2c4a/base/reduce.jl#L535) `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](https://numpy.org/doc/stable/reference/generated/numpy.sum.html) for `numpy.sum`, it appears that they use pairwise summation (not Kahan summation) to improve the output precision.
---
### Author:  [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [August 11, 2023, 4:21pm UTC](https://discourse.julialang.org/t/kahan-summation-in-sum/102723/2 "2023-08-11T16:21:09Z")
Long ago some code for Kahan summation used to live in `Base`, but now it’s in a separate Git repo on Github:
> **[GitHub - JuliaMath/KahanSummation.jl: Sum and cumulative sum using the...](https://github.com/JuliaMath/KahanSummation.jl)**
>
> Sum and cumulative sum using the Kahan-Babuska-Neumaier algorithm - GitHub - JuliaMath/KahanSummation.jl: Sum and cumulative sum using the Kahan-Babuska-Neumaier algorithm
So `Base` doesn’t use compensated summation.
---
### Author:  [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)
#### Post date: [August 11, 2023, 4:22pm UTC](https://discourse.julialang.org/t/kahan-summation-in-sum/102723/3 "2023-08-11T16:22:52Z")
You may be interested by
> **[GitHub - JuliaMath/AccurateArithmetic.jl: Calculate with error-free,...](https://github.com/JuliaMath/AccurateArithmetic.jl)**
>
> Calculate with error-free, faithful, and compensated transforms and extended significands. - GitHub - JuliaMath/AccurateArithmetic.jl: Calculate with error-free, faithful, and compensated transform...
---
### Author:  [@mbauman](https://discourse.julialang.org/u/mbauman)
#### Post date: [August 11, 2023, 4:26pm UTC](https://discourse.julialang.org/t/kahan-summation-in-sum/102723/4 "2023-08-11T16:26:24Z")
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](https://github.com/JuliaLang/julia/pull/4039)
> 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.
---
### Author:  [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [August 11, 2023, 4:30pm UTC](https://discourse.julialang.org/t/kahan-summation-in-sum/102723/5 "2023-08-11T16:30:33Z")
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.
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [August 11, 2023, 4:33pm UTC](https://discourse.julialang.org/t/kahan-summation-in-sum/102723/6 "2023-08-11T16:33:48Z")
this doesn’t apply to compensated sumation (unless you disable subnormals). Addition of floating point numbers never underflows.
---
### Author:  [@cstjean](https://discourse.julialang.org/u/cstjean)
#### Post date: [August 11, 2023, 4:34pm UTC](https://discourse.julialang.org/t/kahan-summation-in-sum/102723/7 "2023-08-11T16:34:14Z")
>
>
> It seems that \`sum\` uses the naive sequential sum algorithm for generators. With… large vectors, it eventually saturates, and yields an incorrect answer:
>
> \`\`\`julia
> julia\> N = 100000000; aa = rand(Float32, N);
>
> julia\> mean((x for x in aa))
> 0.16777216f0
>
> julia\> mean(aa)
> 0.500059f0
> \`\`\`
>
> I have a real-world case where it causes an alarmingly large difference:
>
> \`\`\`julia
> julia\> mean(skipmissing(Umat))
> 1.0638367f0 V
>
> julia\> mean(collect(skipmissing(Umat)))
> 3.1320891f0 V
> \`\`\`
>
> As @simonbyrne pointed out \[on discourse\](https://discourse.julialang.org/t/imprecision-of-mean-over-iterators-of-large-vectors/18759), \`sum(::Array)\` already uses a smarter algorithm. It could presumably be used for generators, too.
[Imprecision of sum(::Generator) · Issue #30421 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/issues/30421#issuecomment-449129195) has a kahan summation implementation.
---
### Author:  [@mbauman](https://discourse.julialang.org/u/mbauman)
#### Post date: [August 11, 2023, 4:40pm UTC](https://discourse.julialang.org/t/kahan-summation-in-sum/102723/8 "2023-08-11T16:40:31Z")
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.
---
### Author:  [@CameronBieganek](https://discourse.julialang.org/u/CameronBieganek)
#### Post date: [August 11, 2023, 5:14pm UTC](https://discourse.julialang.org/t/kahan-summation-in-sum/102723/9 "2023-08-11T17:14:43Z")
Just to clarify, it sounds like `sum` performs pairwise summation for objects that support linear indexing?
---
### Author:  [@mbauman](https://discourse.julialang.org/u/mbauman)
#### Post date: [August 11, 2023, 5:30pm UTC](https://discourse.julialang.org/t/kahan-summation-in-sum/102723/10 "2023-08-11T17:30:21Z")
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`.
---
### Author:  [@e3c6](https://discourse.julialang.org/u/e3c6)
#### Post date: [December 28, 2024, 2:39pm UTC](https://discourse.julialang.org/t/kahan-summation-in-sum/102723/11 "2024-12-28T14:39:10Z")
Should this be mentioned in the docs for `sum`? That is, that the implementation can often give better accuracy than naively summing the elements one by one?
---
### Author:  [@e3c6](https://discourse.julialang.org/u/e3c6)
#### Post date: [December 28, 2024, 2:47pm UTC](https://discourse.julialang.org/t/kahan-summation-in-sum/102723/12 "2024-12-28T14:47:45Z")
Also a related question: is `LinearAlgebra.dot` using similar compensated algorithms for improved accuracy (like `sum`)?
---
### Author:  [@giordano](https://discourse.julialang.org/u/giordano)
#### Post date: [December 28, 2024, 6:15pm UTC](https://discourse.julialang.org/t/kahan-summation-in-sum/102723/13 "2024-12-28T18:15:54Z")
As I also said in your cross-post on Slack, [`LinearAlgebra.dot(::AbstractArray{T}, ::AbstractArray{T}) where {T<:Union{Float32,Float64}}`](https://github.com/JuliaLang/LinearAlgebra.jl/blob/a622302e59a99606936d33dab9da0b51e085e95f/src/blas.jl#L389-L399) calls the corresponding routines in the currently used BLAS library, so that’s a question for the BLAS library you’re using, there’s no promise from the Julia side to do that.
---
### Author:  [@mbauman](https://discourse.julialang.org/u/mbauman)
#### Post date: [December 29, 2024, 1:04am UTC](https://discourse.julialang.org/t/kahan-summation-in-sum/102723/14 "2024-12-29T01:04:03Z")
> [@mbauman](#):
>
> 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).
This isn’t true! My mind couldn’t imagine the possibility, but a recursive approach with iteration alone is totally doable… and there’s a WIP for it, too!
>
>
> Currently, \`mapreduce\` uses a pairwise reduction order for arrays, but switches …to \`mapfoldl\` order for other iterators. This has the unfortunate effect that \[pairwise summation\](https://en.wikipedia.org/wiki/Pairwise\_summation) is only used for arrays, and other iterators get much less accurate floating-point sums (and related reductions). For example, passing an array through a generator or an \`Iterators.filter\` would suddenly make sums less accurate.
>
> I had long been under the mistaken impression that pairwise reduction required random access to an iterator, but @mikmoore pointed out in #52365 that this is not the case.
>
> This WIP PR changes \`mapreduce\` to use a pairwise order by default for arbitrary iterators. I've only done light testing so far, but it should make summation about equally accurate for arrays and other iterators, and the performance seems about the same as the old \`mapfoldl\` fallback.
>
> More testing and benchmarking required, but I wanted to post some code to get the ball rolling. (Should not be a breaking change, in theory, since we explicitly document that the associativity of \`mapreduce\` is implementation-dependent.)
>
> Note that if you want to try out this code on an older version of Julia, just define the following \`foldl\`-like functions
> \`\`\`jl
> import Base: \_InitialValue, mapreduce\_empty\_iter, reduce\_empty\_iter, \_xfadjoint, pairwise\_blocksize, mapreduce\_first, Generator, MappingRF, mapfoldl\_impl
> mapreduce\_empty\_iter(f::F, op::OP, itr) where {F,OP} = reduce\_empty\_iter(\_xfadjoint(op, Generator(f, itr))...)
> mapfoldp(f, op, itr; init=\_InitialValue()) = mapreduce\_impl(f, op, init, itr)
> foldp(op, itr; init=\_InitialValue()) = mapfoldp(identity, op, itr; init=init)
> \`\`\`
> and copy the new \`mapreduce\_impl\` and \`\_mapreduce\_impl\` methods from the PR (the chunk of the diff starting at \`macro \_repeat\`).
>
> Closes #30421.