How Does Kahan Summation Minimize Rounding Errors?

The Kahan compensated summation algorithm reduces accumulated rounding errors when adding long sequences of floating-point numbers by tracking and reintroducing the lost low-order bits at each step. Standard floating-point addition discards small digits when adding a small number to a much larger running total, causing errors to grow with the length of the sequence. By maintaining a separate compensation variable that captures this lost precision, Kahan summation keeps the total error nearly independent of the number of terms summed.

The Problem with Naive Floating-Point Summation

In classical computing, floating-point numbers are represented with a finite number of significand bits, such as 24 bits in IEEE 754 single precision or 53 bits in double precision. When two numbers of vastly different magnitudes are added, the smaller number must be shifted right to align its exponent with the larger number. During this alignment, the least significant bits of the smaller number are shifted past the end of the significand and rounded away.

When summing a sequence of \(n\) numbers \(x_1, x_2, \dots, x_n\), naive summation computes:

\[S_i = S_{i-1} + x_i\]

If \(S_{i-1}\) grows much larger than \(x_i\), each addition loses a small amount of precision. Over \(n\) additions, the worst-case rounding error grows proportionally to \(O(n \epsilon)\), where \(\epsilon\) is the machine epsilon. For long sequences, this accumulated drift can severely degrade numerical accuracy.

How Kahan Compensated Summation Works

Introduced by William Kahan in 1965, compensated summation uses a running compensation variable \(c\) to store the low-order bits lost during each addition. At each step, the algorithm subtracts the previous compensation from the next input value, adds the result to the running sum, and then recovers the rounding error.

For each element \(x_i\) in the sequence, the algorithm performs four operations starting with \(S = 0\) and \(c = 0\):

\[y = x_i - c\] \[t = S + y\] \[c = (t - S) - y\] \[S = t\]

Here, \(y\) is the current input adjusted by the previous error. When \(t = S + y\) is computed, the low-order bits of \(y\) are lost if \(S\) is much larger than \(y\). The expression \((t - S)\) recovers the high-order part of \(y\) that actually survived the addition into \(S\). Subtracting \(y\) from \((t - S)\) isolates the negated lost low-order bits, storing them in \(c\) for the next iteration.

Error Bounds and Practical Considerations

With Kahan summation, the worst-case error bound improves from \(O(n \epsilon)\) to:

\[\left| S_{\text{computed}} - \sum_{i=1}^{n} x_i \right| \le (2\epsilon + O(n \epsilon^2)) \sum_{i=1}^{n} |x_i|\]

Because \(\epsilon^2\) is extremely small (around \(10^{-32}\) for double precision), the error remains effectively constant and bounded by a small multiple of \(\epsilon\) even for millions of terms.

One important caveat in classical computing is compiler optimization. Aggressive floating-point optimizations, such as algebraic reassociation (for example, -ffast-math in C/C++ compilers), can simplify \((S + y) - S - y\) to \(0\), completely eliminating the compensation. To preserve accuracy, Kahan summation requires strict IEEE 754 evaluation order without associative reordering.