Kahan Summation Algorithm

Innumerical analysis, the Kahan summation algorithm , also known as compensated summation ,significantly reduces thenumerical error in the total obtained by adding asequence of finite-precision floating-point numbers , compared to the obvious approach. This is done by keeping a separate running compensation (a variable to accumulate small errors).

In particular, simply summing n numbers in sequence has a worst-case error that grows proportional to n , and aroot mean square error that grows as for random inputs (the roundoff errors form arandom walk).With compensated summation, the worst-case error bound is independent of n , so a large number of values can be summed with an error that only depends on the floating-pointprecision.

Thealgorithm is attributed toWilliam Kahan.Similar, earlier techniques are, for example, Bresenham's line algorithm , keeping track of the accumulated error in integer operations (although first documented around the same time) and the delta-sigma modulation (integrating, not just summing the error).

The algorithm [ edit ]

Inpseudocode, the algorithm is:

<b>function</b> KahanSum(input)
    <b>var</b> sum = 0.0                    // Prepare the accumulator.
    <b>var</b> c = 0.0                      // A running compensation for lost low-order bits.
    <b>for</b> i = 1 <b>to</b> input.length <b>do</b>     // The array <i>input</i> has elements indexed input[1] to input[input.length].
        <b>var</b> y = input[i] - c         // <i>c</i> is zero the first time around.
        <b>var</b> t = sum + y              // Alas, <i>sum</i> is big, <i>y</i> small, so low-order digits of <i>y</i> are lost.
        c = (t - sum) - y            // <i>(t - sum)</i> cancels the high-order part of <i>y</i>; subtracting <i>y</i> recovers negative (low part of <i>y</i>)
        sum = t                      // Algebraically, <i>c</i> should always be zero. Beware overly-aggressive optimizing compilers!
    <b>next</b> i                           // Next time around, the lost low part will be added to <i>y</i> in a fresh attempt.
    <b>return</b> sum

Worked example [ edit ]

This example will be given in decimal. Computers typically use binary arithmetic, but the principle being illustrated is the same. Suppose we are using six-digit decimal floating-point arithmetic, sum has attained the value 10000.0, and the next two values of input[i] are 3.14159 and 2.71828. The exact result is 10005.85987, which rounds to 10005.9. With a plain summation, each incoming value would be aligned with sum , and many low-order digits would be lost (by truncation or rounding). The first result, after rounding, would be 10003.1. The second result would be 10005.81828 before rounding and 10005.8 after rounding. This is not correct.

However, with compensated summation, we get the correct rounded result of 10005.9.

Assume that c has the initial value zero.

y = 3.14159 - 0.00000             <i>y = input[i] - c</i>
  t = 10000.0 + 3.14159
    = 10003.14159                   But only six digits are retained.
    = 10003.1                       Many digits have been lost!
  c = (10003.1 - 10000.0) - 3.14159 This <b>must</b> be evaluated as written! 
    = 3.10000 - 3.14159             The assimilated part of <i>y</i> recovered, vs. the original full <i>y</i>.
    = -0.0415900                    Trailing zeros shown because this is six-digit arithmetic.
sum = 10003.1                       Thus, few digits from <i>input(i</i>) met those of <i>sum</i>.

The sum is so large that only the high-order digits of the input numbers are being accumulated. But on the next step, c gives the error.

y = 2.71828 - (-0.0415900)        The shortfall from the previous stage gets included.
    = 2.75987                       It is of a size similar to <i>y</i>: most digits meet.
  t = 10003.1 + 2.75987             But few meet the digits of <i>sum</i>.
    = 10005.85987                   And the result is rounded
    = 10005.9                       To six digits.
  c = (10005.9 - 10003.1) - 2.75987 This extracts whatever went in.
    = 2.80000 - 2.75987             In this case, too much.
    = 0.040130                      But no matter, the excess would be subtracted off next time.
sum = 10005.9                       Exact result is 10005.85987, this is correctly rounded to 6 digits.

So the summation is performed with two accumulators: sum holds the sum, and c accumulates the parts not assimilated into sum , to nudge the low-order part of sum the next time around. Thus the summation proceeds with "guard digits" in c , which is better than not having any, but is not as good as performing the calculations with double the precision of the input. However, simply increasing the precision of the calculations is not practical in general; if input is already in double precision, few systems supplyquadruple precision, and if they did, input could then be in quadruple precision.

Accuracy [ edit ]

A careful analysis of the errors in compensated summation is needed to appreciate its accuracy characteristics. While it is more accurate than naive summation, it can still give large relative errors for ill-conditioned sums.

Suppose that one is summing n values x _i , for i = 1, ... , n . The exact sum is

(computed with infinite precision).

With compensated summation, one instead obtains , where the error is bounded by

where ε is themachine precision of the arithmetic being employed (e.g. ε ≈ 10 ⁻¹⁶ for IEEE standarddouble-precision floating point). Usually, the quantity of interest is therelative error , which is therefore bounded above by

In the expression for the relative error bound, the fraction Σ| x _i |/|Σ x _i | is thecondition number of the summation problem. Essentially, the condition number represents the intrinsic sensitivity of the summation problem to errors, regardless of how it is computed.The relative error bound of every (backwards stable) summation method by a fixed algorithm in fixed precision (i.e. not those that usearbitrary-precision arithmetic, nor algorithms whose memory and time requirements change based on the data), is proportional to this condition number.An ill-conditioned summation problem is one in which this ratio is large, and in this case even compensated summation can have a large relative error. For example, if the summands x _i are uncorrelated random numbers with zero mean, the sum is arandom walk, and the condition number will grow proportional to . On the other hand, for random inputs with nonzero mean the condition number asymptotes to a finite constant as . If the inputs are allnon-negative, then the condition number is 1.

Given a condition number, the relative error of compensated summation is effectively independent of n . In principle, there is the O( nε ² ) that grows linearly with n , but in practice this term is effectively zero: since the final result is rounded to a precision ε , the nε ² term rounds to zero, unless n is roughly 1/ ε or larger.In double precision, this corresponds to an n of roughly 10 ¹⁶ , much larger than most sums. So, for a fixed condition number, the errors of compensated summation are effectively O ( ε ), independent of  n .

In comparison, the relative error bound for naive summation (simply adding the numbers in sequence, rounding at each step) grows as multiplied by the condition number.This worst-case error is rarely observed in practice, however, because it only occurs if the rounding errors are all in the same direction. In practice, it is much more likely that the rounding errors have a random sign, with zero mean, so that they form a random walk; in this case, naive summation has aroot mean square relative error that grows as multiplied by the condition number.This is still much worse than compensated summation, however. However, if the sum can be performed in twice the precision, then ε is replaced by ε ² , and naive summation has a worst-case error comparable to the O( nε ² ) term in compensated summation at the original precision.

By the same token, the Σ| x _i | that appears in above is a worst-case bound that occurs only if all the rounding errors have the same sign (and are of maximal possible magnitude).In practice, it is more likely that the errors have random sign, in which case terms in Σ| x _i | are replaced by a random walk, in which case, even for random inputs with zero mean, the error grows only as (ignoring the nε ² term), the same rate the sum grows, canceling the factors when the relative error is computed. So, even for asymptotically ill-conditioned sums, the relative error for compensated summation can often be much smaller than a worst-case analysis might suggest.

Further enhancements [ edit ]

Neumaierintroduced an improved version of Kahan algorithm, which he calls an "improved Kahan–Babuška algorithm", which also covers the case when the next term to be added is larger in absolute value than the running sum, effectively swapping the role of what is large and what is small. Inpseudocode, the algorithm is:

<b>function</b> NeumaierSum(input)
    <b>var</b> sum = 0.0
    <b>var</b> c = 0.0                       // A running compensation for lost low-order bits.
    <b>for</b> i = 1 <b>to</b> input.length <b>do</b>
        <b>var</b> t = sum + input[i]
        <b>if</b> |sum| >= |input[i]| <b>then</b>
            c += (sum - t) + input[i] // If <i>sum</i> is bigger, low-order digits of <i>input[i]</i> are lost.
        <b>else</b>
            c += (input[i] - t) + sum // Else low-order digits of <i>sum</i> are lost.
        <b>endif</b>
        sum = t
    <b>next</b> i
    <b>return</b> sum + c                    // Correction only applied once in the very end.

For many sequences of numbers, both algorithms agree, but a simple example due to Petersshows how they can differ. For summing in double precision, Kahan's algorithm yields 0.0, whereas Neumaier's algorithm yields the correct value 2.0.

Higher-order modifications of better accuracy are also possible. For example a variant suggested by Klein, which he called a second-order "iterative Kahan–Babuška algorithm". Inpseudocode, the algorithm is:

<b>function</b> KleinSum(input)
    <b>var</b> s = 0.0
    <b>var</b> cs = 0.0
    <b>var</b> ccs = 0.0
    <b>for</b> i = 1 <b>to</b> input.length <b>do</b>
        <b>var</b> t = s + input[i]
        <b>if</b> |s| >= |input[i]| <b>then</b>
            c = (s - t) + input[i]
        <b>else</b>
            c = (input[i] - t) + s
        <b>endif</b>
        s = t
        t = cs + c
        <b>if</b> |cs| >= |c| <b>then</b>
            cc = (cs - t) + c
        <b>else</b>
            cc = (c - t) + cs
        <b>endif</b>
        cs = t
        ccs = ccs + cc
    <b>next</b> i
    <b>return</b> s + cs + ccs

Alternatives [ edit ]

Although Kahan's algorithm achieves error growth for summing n numbers, only slightly worse growth can be achieved bypairwise summation: onerecursively divides the set of numbers into two halves, sums each half, and then adds the two sums.This has the advantage of requiring the same number of arithmetic operations as the naive summation (unlike Kahan's algorithm, which requires four times the arithmetic and has a latency of four times a simple summation) and can be calculated in parallel. The base case of the recursion could in principle be the sum of only one (or zero) numbers, but toamortize the overhead of recursion, one would normally use a larger base case. The equivalent of pairwise summation is used in many fast Fourier transform (FFT) algorithms and is responsible for the logarithmic growth of roundoff errors in those FFTs.In practice, with roundoff errors of random signs, the root mean square errors of pairwise summation actually grow as .

Another alternative is to use arbitrary-precision arithmetic , which in principle need no rounding at all with a cost of much greater computational effort. A way of performing exactly rounded sums using arbitrary precision is to extend adaptively using multiple floating-point components. This will minimize computational cost in common cases where high precision is not needed.Another method that uses only integer arithmetic, but a large accumulator, was described by Kirchner and Kulisch;a hardware implementation was described by Müller, Rüb and Rülling.

Possible invalidation by compiler optimization [ edit ]

In principle, a sufficiently aggressiveoptimizing compiler could destroy the effectiveness of Kahan summation: for example, if the compiler simplified expressions according to theassociativity rules of real arithmetic, it might "simplify" the second step in the sequence

t = sum + y;

c = (t - sum) - y;

to

c = ((sum + y) - sum) - y;

and then to

c = 0;

thus eliminating the error compensation.In practice, many compilers do not use associativity rules (which are only approximate in floating-point arithmetic) in simplifications, unless explicitly directed to do so by compiler options enabling "unsafe" optimizations,although theIntel C++ Compiler is one example that allows associativity-based transformations by default.The originalK&R C version of the C programming language allowed the compiler to re-order floating-point expressions according to real-arithmetic associativity rules, but the subsequentANSI C standard prohibited re-ordering in order to make C better suited for numerical applications (and more similar toFortran, which also prohibits re-ordering),although in practice compiler options can re-enable re-ordering, as mentioned above.

Support by libraries [ edit ]

In general, built-in "sum" functions in computer languages typically provide no guarantees that a particular summation algorithm will be employed, much less Kahan summation.

[ citation needed ]

TheBLAS standard forlinear algebra subroutines explicitly avoids mandating any particular computational order of operations for performance reasons,and BLAS implementations typically do not use Kahan summation.

The standard library of thePython computer language specifies an fsum function for exactly rounded summation, using theShewchuk algorithmto track multiple partial sums.

In theJulia language, the default implementation of the sum function doespairwise summation for high accuracy with good performance,but an external library provides an implementation of Neumaier's variant named sum_kbn for the cases when higher accuracy is needed.

In theC♯ language, HPCsharp nuget package implements the Neumaier variant andpairwise summation: both as scalar, data-parallel usingSIMD processor instructions, and parallel multi-core.

See also [ edit ]

• Algorithms for calculating variance , which includes stable summation

References [ edit ]