Adding up ten million numbers

A floating-point sum looks like the most innocent computation in a program. It is computed one addition at a time, $\hat s_k = \operatorname{fl}(\hat s_{k-1} + x_k)$, and every addition rounds its result to the nearest representable number. Each rounding is tiny, but a long sum performs millions of them, and their errors do not have to cancel.

For recursive summation, the classic bound is¹

where $u$ is the unit roundoff: $2^{-53} \approx 1.1 \times 10^{-16}$ in double precision and $2^{-24} \approx 6.0 \times 10^{-8}$ in single precision. The bound is pessimistic because it assumes every rounding error points the same way. It still grows with $n$, and for positive terms the partial sum grows too, so later terms are added to an ever larger total and lose more of their low-order bits.

Single precision is where this stops being academic. Sensor data, graphics, and machine-learning code store float32 arrays as a matter of course, and $10^7$ terms is an ordinary array.

I drew $n$ values uniformly from $[0, 1)$ as float32 and summed them three ways, comparing each result with math.fsum, which returns the correctly rounded sum of the same values:

• recursive: one addition after another, as a plain loop does;
• pairwise: np.sum, which splits the array into blocks and adds partial sums in a tree;
• compensated: Kahan’s algorithm, below.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

import numpy as np

def kahan_sum(values):
    total = np.float32(0.0)
    compensation = np.float32(0.0)
    for value in values:
        y = value - compensation
        t = total + y
        compensation = (t - total) - y
        total = t
    return total

Relative errors of the single-precision results:

At ten million terms the recursive sum is about 4,999,406 and off by 127, which is 255 units in the last place. The compensated sum was correctly rounded at every size: no single-precision number lies closer to the exact sum. Pairwise summation was correctly rounded at every size but one, where it was a single unit in the last place away.

In double precision the same experiment is quieter: the recursive sum of $10^7$ terms is off by a relative $1.8 \times 10^{-14}$, and both other methods match math.fsum exactly. The problem has not gone away; it has moved to longer sums and to data with cancellation.

Kahan’s algorithm keeps a second variable for the part of each addition that rounding throws away.² In exact arithmetic, line 9 computes $(t - \text{total}) - y = 0$. In floating-point arithmetic, $t - \text{total}$ is the part of $y$ that actually made it into the sum, so subtracting $y$ leaves the negative of what was lost. Line 7 adds it back to the next term before that term is added to the total.

The bound for compensated summation no longer grows with $n$ at first order:

Pairwise summation reaches a similar effect differently. Adding numbers in a balanced tree means that each term passes through about $\log_2 n$ additions instead of up to $n - 1$, so the factor $(n - 1)$ in the first bound becomes $\lceil \log_2 n \rceil$.

Compensated summation performs four floating-point operations per term instead of one, and the loop carries a dependency through compensation that limits how much a processor can overlap. In interpreted Python the difference disappears in the interpreter’s overhead; in compiled code it is real but often hidden by memory traffic, because a long sum is limited by how fast the data arrive.

The practical advice is short. Prefer library reductions, which usually sum pairwise. Write a compensated loop when you must accumulate in a loop, in single precision, or over many terms. Never assume that a sum is accurate because each addition is.