Floating point, conditioning and numerical stability

Problem

Add 0.1 to itself 100,000 times in float32, one term at a time, and the answer is 9998.5566 rather than 10,000. Compute the variance of four readings near 10⁹ with the textbook shortcut E[x²] − E[x]² in float64 and it comes out negative. Evaluate a softmax over the logits 1000, 999, 998 and every probability is NaN. None of these is a bug in the code that ran: each is ordinary rounding, amplified by the way the formula was arranged.

Machine learning runs into all of them because it works at scale and at low precision. Sums run over millions of terms, logits and log-probabilities span huge ranges, statistics are accumulated over streams, and training increasingly happens in float16, bfloat16 and smaller, where each rounding is a thousand times coarser than in float32. Knowing which formulas lose accuracy, and what to use instead, is the difference between a model that trains and one that silently produces garbage.

Intuition

A floating-point number is a significand with a fixed number of bits times a power of two. Within one power of two (one binade) the representable numbers are evenly spaced, and the spacing doubles at every power of two. So the absolute gap between neighbours grows with the size of the number, while the relative gap stays roughly constant. That relative gap is the whole story: every operation is exact up to a relative error of at most half of it.

The numbers come from toy_floats in the page's helpers. float32 is the same picture with 23 fraction bits instead of 2, so its relative gap is 2⁻²³ ≈ 1.2 × 10⁻⁷ instead of ¼.

With that picture, two distinct questions decide whether a computed answer can be trusted. Conditioning belongs to the problem: if the inputs move by a tiny relative amount, how far can the exact answer move? Stability belongs to the algorithm: does it add more error than the unavoidable rounding of its inputs would cause anyway? A stable algorithm on a well-conditioned problem gives an accurate answer. On an ill-conditioned problem even a perfect algorithm cannot, and an unstable algorithm can ruin a well-conditioned one.

Mechanics

The rounding model

IEEE 754 arithmetic rounds the exact result of each +, −, ×, ÷ and √ to the nearest representable number. With p significand bits (counting the implicit leading 1), that gives the standard model

fl(x ∘ y) = (x ∘ y)(1 + δ),   |δ| ≤ u = 2^−p,

where u is the unit roundoff: 2⁻²⁴ ≈ 6.0 × 10⁻⁸ for float32 and 2⁻⁵³ ≈ 1.1 × 10⁻¹⁶ for float64. The machine epsilon that torch.finfo(dtype).eps reports is the gap above 1.0, which is 2u. Two consequences follow at once. Rounding is relative, so a quantity is only known to about log₁₀(1/u) significant digits: 7 in float32, 16 in float64, 3 in float16 and 2 to 3 in bfloat16. And floating-point addition is not associative: (0.1 + 0.2) + 0.3 = 0.6000000000000001 but 0.1 + (0.2 + 0.3) = 0.6 in float64. Any reduction whose order changes, across threads, GPU blocks or devices, changes the last bits of the result, which is why bitwise reproducibility needs a fixed order.

Summation

Adding n numbers one at a time rounds after every addition, and each rounding is relative to the running total, not to the term. The error bound for this recursive summation is

|ŝ − s| ≤ (n − 1) u Σ|x_i| + O(u²),

linear in n. Once the total is large, small terms lose most of their bits on the way in. Three standard fixes:

• Pairwise (cascade) summation adds neighbours, then sums of neighbours, in a tree. Each term passes through only ⌈log₂ n⌉ additions, so the bound becomes ⌈log₂ n⌉ u Σ|x_i|. PyTorch's CPU sum does this: it accumulates blocks of at least 16 elements into four levels of partial sums so that only numbers of similar size are added together.
• Compensated (Kahan) summation keeps a second variable holding the low-order bits each addition dropped and feeds them back into the next term. The bound is 2u Σ|x_i| + O(nu²)Σ|x_i|, essentially independent of n, for four operations per term instead of one.
• A wider accumulator: sum float16 or float32 inputs in float32 or float64. scikit-learn's _safe_accumulator_op runs sums and dot products in float64, and GPU kernels accumulate half-precision products in float32.

Both traces from f32_naive_sum and f32_kahan_sum, which round to float32 after every operation with struct. Between powers of two the total's gap is fixed, so each addition of 0.1 rounds the same way and the error grows linearly; at each power of two (8,192, 16,384, 32,768, 65,536) the gap doubles and the direction or size of the drift changes.

Catastrophic cancellation

Subtracting two nearly equal numbers is itself exact, but it removes their shared leading digits and leaves only the trailing ones, which are exactly the digits earlier roundings corrupted. If a and b each carry a relative error u, then a − b carries a relative error up to u(|a| + |b|)/|a − b|, which is enormous when a ≈ b. The textbook variance Σx²/n − (Σx/n)² is the classic case: both terms are about the square of the mean, their difference is the much smaller variance, and every digit of the variance smaller than u times the mean squared is lost. Subtract the mean first (two passes), or update the mean and the sum of squared deviations together as each point arrives (Welford's algorithm), and the large common part never enters a subtraction.

Conditioning and backward stability

The condition number of a problem is the worst-case ratio of the relative change in the answer to the relative change in the input. For solving Ax = b it is κ(A) = ‖A‖‖A⁻¹‖ = σ_max/σ_min, met on the linear algebra toolkit page. An algorithm is backward stable if its computed answer is the exact answer for a slightly perturbed input, with a relative perturbation of order u. Then the two ideas combine into the rule of numerical analysis:

forward error ≲ κ × backward error ≈ κ u.

LU with partial pivoting, Cholesky and Householder QR are backward stable in practice, so their error is as small as the problem allows: you lose about log₁₀ κ digits, and no algorithm can do better in that precision. Forming XᵀX for least squares squares κ before a solver ever sees it, which is why the QR route is preferred. The residual b − Ax̂ measures the backward error, so a tiny residual proves the algorithm did its job, not that x̂ is accurate.

Stable rewrites in machine learning

Most of the numerically delicate code in a deep learning framework is one of a few rewrites, each keeping a large intermediate from overflowing or a small one from vanishing:

• logsumexp and softmax: with m = max_i x_i, log Σ e^x_i = m + log Σ e^{x_i − m}. Every shifted exponent is at most 0, so nothing overflows, and one term equals 1, so the sum is between 1 and n and the log never sees 0. Softmax is e^{x − m} / Σ e^{x − m}, and log-softmax is x − logsumexp(x), computed directly instead of as log(softmax(x)), whose small probabilities would underflow to log 0.
• log1p and expm1: for small z, 1 + z rounds away z's low digits before the log sees them; log1p(z) computes log(1 + z) without forming 1 + z.
• softplus, log(1 + eˣ): overflows for large x and loses everything for very negative x. The stable form switches regimes: eˣ when x is very negative, log1p(eˣ), then x + e⁻ˣ, then just x. PyTorch's Softplus returns x once βx exceeds a threshold (20 by default).
• losses on logits: binary cross-entropy through sigmoid and log loses accuracy once σ(z) rounds to 0 or 1. On the logit z directly, −log σ(z) = max(−z, 0) + log1p(e^−|z|), which is what BCEWithLogitsLoss and cross_entropy compute. See entropy and KL and loss functions.
• epsilons in denominators: normalization layers divide by √(var + ε) and Adam by √v̂ + ε, so a zero variance or second moment does not divide by zero. The ε must be representable in the dtype the division runs in: Adam's default 10⁻⁸ is below float16's smallest positive number, 6.0e-08, and rounds to 0.

The pieces of PyTorch and scikit-learn that implement these rules:

Worked example

Computed by hilbert_solve with the page's lu_factor and lu_solve; κ from its singular values. The gap between forward and backward error is the condition number at work.

Implementation

import math

def logsumexp(xs):
    m = max(xs)
    if math.isinf(m):
        return m                       # all -inf, or an inf present
    return m + math.log(sum(math.exp(x - m) for x in xs))

def softmax(xs):
    m = max(xs)
    exps = [math.exp(x - m) for x in xs]
    total = sum(exps)
    return [v / total for v in exps]

def softplus(x):                       # log(1 + e^x), cutoffs as in scikit-learn
    if x <= -37: return math.exp(x)
    if x <= -2:  return math.log1p(math.exp(x))
    if x <= 18:  return math.log(1 + math.exp(x))
    if x <= 33.3: return x + math.exp(-x)
    return x

def bce_with_logits(z, y):             # -[y log s(z) + (1 - y) log(1 - s(z))]
    return (1 - y) * z - min(z, 0.0) + math.log1p(math.exp(-abs(z)))

def kahan_sum(xs):
    total = carry = 0.0
    for x in xs:
        y = x - carry
        t = total + y
        carry = (t - total) - y        # the low bits that did not fit in t
        total = t
    return total

def welford(xs):
    n = mean = m2 = 0
    for x in xs:
        n += 1
        delta = x - mean
        mean += delta / n
        m2 += delta * (x - mean)
    return mean, m2 / (n - 1)

welford([1e9 + 4, 1e9 + 7, 1e9 + 13, 1e9 + 16])   # (1000000010.0, 30.0)

The same habits in PyTorch and scikit-learn:

import torch
import torch.nn.functional as F

logits = torch.tensor([1000., 999., 998.])
torch.logsumexp(logits, dim=0)               # 1000.4076, not inf
F.log_softmax(logits, dim=0)                 # not softmax(...).log()
F.binary_cross_entropy_with_logits(z, y)     # not F.binary_cross_entropy(torch.sigmoid(z), y)
F.cross_entropy(logits_batch, labels)        # takes logits, applies log_softmax inside

x = torch.randn(10_000_000, dtype=torch.float16)
x.sum(dtype=torch.float32)                   # accumulate in a wider type
torch.finfo(torch.float16).eps, torch.finfo(torch.bfloat16).eps   # 0.000977, 0.0078

torch.testing.assert_close(actual, expected)                  # per-dtype rtol/atol defaults
torch.testing.assert_close(a, b, rtol=1e-4, atol=1e-6)         # or explicit ones

torch.linalg.cond(A)                         # κ before trusting a solve
torch.linalg.lstsq(X, y)                     # QR/SVD based, not solve(X.T @ X, X.T @ y)

from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
for batch in batches:                        # streaming mean and variance, float64 accumulators
    scaler.partial_fit(batch)

Tradeoffs

• Precision against speed and memory. float64 makes most of these problems disappear and costs 2× the memory and, on GPUs, often 30× or more the time. Deep learning instead runs low precision where errors average out (matrix products with wide accumulators) and keeps float32 for the parts that accumulate: reductions, normalization statistics, optimizer state and the loss.
• Compensated sums cost operations. Kahan summation does four floating-point operations per term and a serial dependency; pairwise summation is nearly free and parallel, which is why libraries default to it. Compilers allowed to reassociate (-ffast-math) can optimize the compensation away entirely.
• Stable rewrites change the graph. A fused log_softmax or a loss on logits is both faster and more accurate than composing the pieces, but only if the model exposes logits; putting a sigmoid or softmax inside the model and a log in the loss throws the stability away.
• Determinism against throughput. Fixing the reduction order makes results bitwise reproducible, but rules out atomics and some fast kernels. Differences in the last few bits between runs are expected; differences larger than the tolerances below are bugs.
• Tolerances must match the dtype. Compare with a relative tolerance a few times the dtype's epsilon plus a small absolute floor for values near zero. assert_close's defaults:

Check yourself

Connections

References

• Goldberg. What Every Computer Scientist Should Know About Floating-Point Arithmetic. ACM Computing Surveys 23(1), 5–48, 1991.
• Higham. Accuracy and Stability of Numerical Algorithms, 2nd ed. SIAM, 2002: chapters 2 (the rounding model), 4 (summation) and 7–9 (conditioning and LU).
• Trefethen and Bau. Numerical Linear Algebra. SIAM, 1997: lectures 12–15 on conditioning and backward stability.
• Kahan. Further Remarks on Reducing Truncation Errors. Communications of the ACM 8(1), 40, 1965.
• Welford. Note on a Method for Calculating Corrected Sums of Squares and Products. Technometrics 4(3), 419–420, 1962.
• Chan, Golub and LeVeque. Algorithms for Computing the Sample Variance: Analysis and Recommendations. The American Statistician 37(3), 242–247, 1983.
• Mächler. Accurately Computing log(1 − exp(−|a|)), Assessed by the Rmpfr Package. 2012, the note scikit-learn's log1pexp follows.
• Blanchard, Higham and Higham. Accurately Computing the Log-Sum-Exp and Softmax Functions. IMA Journal of Numerical Analysis 41(4), 2311–2330, 2021.
• IEEE. 754-2019, Standard for Floating-Point Arithmetic. 2019.
• PyTorch v2.15.0: logsumexp, _log_softmax, BCEWithLogitsLoss, cascade_sum, _DTYPE_PRECISIONS.
• scikit-learn 1.9.1: log1pexp, _incremental_mean_and_var, _safe_accumulator_op.