# 🧊 Floating Point, Catastrophic Cancellation and Mixed Precision

By Better Wise

Modules 1 and 2 gave you the objectives (KL, cross-entropy, mutual information) and the solvers (QR, Cholesky, SVD, Tikhonov). Conditioning told you how much a problem amplifies input error. This lesson is about the other half of the error budget: the arithmetic itself rounds, and some expressions throw away almost every digit they were given. A backward-stable solver run on a cancellation-prone residual still returns garbage, and a bf16 softmax can silently tie logits that were never tied.

By the end you should be able to:

- Read a float format as exponent plus mantissa and convert that into a unit roundoff.
- Recognize cancellation-prone expressions and rewrite them in stable form.
- Choose a summation strategy whose error grows like $\sqrt{n}u$ rather than $nu$.
- Decide which operations in a training loop can drop to bf16 or fp16 and which must stay fp32.

## What a float stores, and the error it promises

Every format splits its bits into sign, exponent (range) and mantissa (relative precision). With $p$ stored mantissa bits, one correctly rounded operation on exact inputs satisfies

$$
u = 2^{-p-1}, \qquad \text{fp32: } u \approx 6\times10^{-8}
$$

| Format | Exponent bits | Mantissa bits | Unit roundoff $u$ | Approx. range |
| --- | --- | --- | --- | --- |
| float64 | 11 | 52 | $1.1\times10^{-16}$ | $10^{\pm308}$ |
| float32 | 8 | 23 | $6.0\times10^{-8}$ | $10^{\pm38}$ |
| float16 | 5 | 10 | $4.9\times10^{-4}$ | $6\times10^{-5}$ to $65504$ |
| bfloat16 | 8 | 7 | $3.9\times10^{-3}$ | $10^{\pm38}$ |

bfloat16 keeps 8 exponent bits but only 7 mantissa bits, so its unit roundoff is about _____, roughly three decimal digits. Note the shape of the promise: the bound is *relative* to the result. When the result is tiny compared with the inputs, the same relative bound on the inputs becomes an enormous relative error on the output.

## Catastrophic cancellation: the subtraction is innocent

Subtraction of two floats is exact when they are close. The damage was done earlier. If $\hat{x} = x + \delta_x$ with $|\delta_x| \lesssim u|x|$, then

$$
\frac{|\hat{x}-\hat{y} - (x-y)|}{|x-y|} \lesssim u\,\frac{|x|+|y|}{|x-y|}
$$

Subtracting two nearly equal numbers is dangerous because the _____ error in the inputs stays fixed while the result shrinks, so the ratio explodes.

**Worked example: the quadratic formula.** For $x^2 - 10^8 x + 1 = 0$, the small root is $10^{-8}$. In fp32, $b^2 - 4ac = 10^{16} - 4$ and the relative perturbation $4\times10^{-16}$ sits far below $u$, so the discriminant rounds to exactly $10^{16}$, the square root to $10^{8}$, and the small root to $0$. Fix it with the identity $x_1 x_2 = c/a$:

```python
disc = b*b - 4*a*c
x_large = (-b - np.sign(b)*np.sqrt(disc)) / (2*a)   # no cancellation
x_small = c / (a * x_large)                          # product of roots
```

**Worked example: variance.** For $x_i = 10^4 + \varepsilon_i$ with $\varepsilon \sim \mathcal{N}(0,1)$, the variance is $1$ but $\mathbb{E}[x^2] \approx 10^8$. In fp32 the absolute error on $\mathbb{E}[x^2]$ is about $10^8 \times 6\times10^{-8} = 6$, so `E[x**2] - E[x]**2` returns a variance wrong by hundreds of percent. Use a two-pass computation or Welford's online update.

**Worked example: log and exp.** `log(1 + x)` for $x = 10^{-9}$ loses the answer entirely in fp32, so use `np.log1p(x)`. The same idea stabilizes the KL term $p\log(p/q)$ when $p \approx q$: compute it as `p * log1p((p - q) / q)`, and prefer fp64 for the small-difference regime.

**Worked example: log-sum-exp.** In fp16, `exp(20)` overflows the $65504$ maximum. Subtract the max logit first; the shift cancels exactly in the ratio, so the result is unchanged mathematically and enormous numerically. Before exponentiating, subtract the _____ logit so that every exponent is at most zero.

## Reductions: error grows with the number of terms

Naive left-to-right summation has worst case $nu\sum|x_i|$ and random-walk behavior around $\sqrt{n}\,u\sum|x_i|$. For $n = 10^6$ in fp32 that is already a few digits. Pairwise (tree) summation brings this to $O(u\log n)$; compensated summation brings it near $2u$ independent of $n$:

```python
def kahan(xs):
    s = 0.0; c = 0.0
    for x in xs:
        y = x - c
        t = s + y
        c = (t - s) - y
        s = t
    return s
```

The practical consequence is that reduction order changes results at the last digits, so bitwise comparison across hardware or batch shapes is not a meaningful correctness test. It also means norms, means and batch statistics are exactly where low-precision accumulation quietly costs you accuracy.

## Mixed precision: storage, compute, accumulation

The useful mental split is three separate decisions: what format stores the tensor, what format the multiply runs in, and what format the accumulation runs in. Tensor cores do bf16 or fp16 multiplies with fp32 accumulation, which is why matmul is the safest place to lower precision.

```mermaid
flowchart TD
  A["Operation in the loop"] --> B{"Reduction or near-equal subtraction?"}
  B -->|"yes"| C["🛡️<br/>fp32 or fp64"]
  B -->|"no"| D{"Range risk: exp, norms?"}
  D -->|"yes"| E["bf16 compute<br/>fp32 accumulate"]
  D -->|"no"| F["fp16 compute<br/>plus loss scaling"]
```

| Stage | Safe choice | Why |
| --- | --- | --- |
| Matmul, convolution | bf16 or fp16 inputs, fp32 accumulate | wide accumulator absorbs per-product rounding |
| Master weights | fp32 | updates are smaller than the low-precision ulp near 1 |
| Softmax, log, exp, logsumexp | fp32, max-shifted | cancellation and range |
| Loss and its reduction | fp32 | small terms vanish against large ones |
| Norms, variance, means | fp32, pairwise or Welford | $nu$ growth in long reductions |
| Optimizer state | fp32 | second moments span many orders of magnitude |
| Epsilon constants | format-aware | `1e-8` flushes to zero in fp16 |

**Worked audit: a contrastive loss in bf16.** You compute InfoNCE logits in bf16 and take log-softmax. Near magnitude 12, the bf16 ulp is $12 \times 2^{-7} \approx 0.06$, so two logits $0.01$ apart become the same number, the softmax returns uniform probabilities, and the gradient with respect to that pair is zero. The loss curve looks fine; the representation stops learning.

```python
import torch
torch.tensor([12.0, 12.01], dtype=torch.bfloat16)
# tensor([12., 12.], dtype=torch.bfloat16): the gap is gone
```

**Worked audit: weight updates.** A weight at $1.0$ with a gradient step of $10^{-4}$ gets rounded to $1.0$ in both bf16 (ulp $0.0039$) and fp16 (ulp $0.001$). Keeping an fp32 master copy fixes it; the low-precision tensor is only a compute view. fp16 needs _____ because gradients below $6\times10^{-5}$ underflow to zero; bf16 does not, since its exponent range matches fp32. Optimizer state and master weights stay in _____ because updates are tiny compared to the weights.

## Guardrails: spending your precision budget

- Do not test stability on "nice" random data. Generate inputs with an offset of $10^6$ and a signal of $1$; cancellation bugs only appear with a large mean and a small spread.
- When a solver returns a bad answer, check the residual computation before blaming the factorization. A residual computed in the working precision cannot certify accuracy below $u\|A\|\|x\|$.
- Add `eps` at a scale that is representable in the target format, and add it to the squared norm, not the norm.
- Compensated summation can be optimized away by aggressive fast-math flags; verify, do not assume.
- Lower precision in the forward pass, but keep the loss, the reductions and the update in fp32. That single rule prevents most silent mixed-precision failures.
- If two runs differ only in the last three digits, suspect reduction order, not a bug.

## Deck for spaced review

Front: What does unit roundoff u bound?
Back: The relative error of one correctly rounded operation on exact inputs, with u = 2^(-p-1) for p stored mantissa bits.

Front: Unit roundoff of float32?
Back: About 6e-8, roughly 7 decimal digits of relative precision.

Front: bfloat16 versus float16: where does each trade precision?
Back: bfloat16 keeps fp32's 8-bit exponent and drops to 7 mantissa bits; float16 keeps 10 mantissa bits but overflows at 65504.

Front: Root cause of catastrophic cancellation?
Back: The inputs already carry absolute error of order u|x|, and subtracting near-equal values divides that fixed error by a tiny difference.

Front: Stable way to get the small root of a quadratic?
Back: Compute the large root with a sign-aware formula, then use the product of roots: x_small = c / (a * x_large).

Front: Why is E[x^2] - E[x]^2 unstable for large means?
Back: Both terms are about 1e8 while the variance is O(1), so the subtraction exposes the rounding error of the large terms.

Front: Two stable ways to compute a variance?
Back: Welford's online update, or a two-pass algorithm that subtracts the mean first.

Front: Why subtract the max in softmax and logsumexp?
Back: It makes every exponent at most zero, preventing overflow while leaving the normalized result unchanged.

Front: What are log1p and expm1 for?
Back: Accurate computation of log(1+x) and exp(x)-1 near zero, where the naive forms lose precision.

Front: Where is float16 typically safe inside a training loop?
Back: In matrix multiplies and convolutions on tensor cores, with fp32 master weights and fp32 accumulation.

Front: Why keep a float32 master copy of weights in mixed precision?
Back: Small updates are lost if added directly to bfloat16 or float16 weights; the fp32 copy preserves them.

Front: What does loss scaling protect against in float16 training?
Back: Underflow of tiny gradients to zero; scale the loss up before backward and unscale before the optimizer step.

Front: Why is summing many small terms in float32 risky?
Back: Rounding error grows with the number of additions; pairwise or Kahan summation keeps the error near u instead of n*u.

Front: What is the first question to ask about a suspicious numerical result?
Back: What is the condition number of the underlying problem, and how much precision does the algorithm lose relative to it?

Unit roundoff for p stored mantissa bits is _____.
Catastrophic cancellation turns a relative error of order u into an absolute error of order _____, which then gets divided by the tiny difference.
In mixed precision, matrix multiplies run in low precision while _____ stays in float32.

## Precision Habits That Keep Pipelines Honest

- Treat unit roundoff as a budget, not a footnote. If your algorithm needs to resolve a difference of size $\delta$ from inputs of size $M$, you need roughly $M/\delta$ dynamic range, and fp32 gives you about $10^7$ relative units before answers become noise.
- Compute the condition number before you blame the code. A stable algorithm on an ill-conditioned problem still returns an ill-conditioned answer. Precision work and conditioning work are separate conversations.
- For variances, means of squares, and softmax denominators, prefer the shifted or online form by default. The naive formula is not "slightly worse"; it can return negative variances and infinities.
- In mixed precision, decide per operation: low precision for throughput-heavy matrix multiplies, fp32 for accumulation, master weights, and any reduction whose result feeds a sensitive branch.
- Common mistake: assuming bfloat16 is always safer because its exponent range matches fp32. The smaller mantissa makes small weight updates disappear sooner, so master weights matter even more.

## Closing the Loop on Finite Precision

Floating point error is relative to magnitude. Cancellation converts that relative error into absolute error and then divides by a small difference, which is why subtraction of near-equal quantities is the usual suspect when a model silently breaks. Stable rewrites are local and cheap: sign-aware quadratic roots, centered or Welford variance, max-shifted softmax and logsumexp, and `log1p`/`expm1` near zero.

Mixed precision is a placement decision, not a global switch. Keep fp32 where accumulation, master weights, and sensitive reductions live; spend low precision where tensor cores buy throughput and the algorithm is already averaging over many terms. Next, Spectral Methods: Eigenbases, Filtering and Smoothing carries the same conditioning lens into the eigenbasis, where filtering and smoothing become diagonal operations.

---

## Contents

Part of [🎲 Mathematics for Machine Learning 5: Information Theory, Matrix Methods and Numerical Stability](https://everwise.study/courses/mathematics-for-machine-learning-5-information-theory-matrix-methods-and.md)

- **🧮 Information-Theoretic Training Objectives**
  - [📉 Entropy, Cross-Entropy and KL: One Family of Losses](https://everwise.study/courses/mathematics-for-machine-learning-5-information-theory-matrix-methods-and/entropy-cross-entropy-kl-one-family.md)
  - [🔗 Mutual Information as an Objective: InfoNCE and Contrastive Bounds](https://everwise.study/courses/mathematics-for-machine-learning-5-information-theory-matrix-methods-and/mutual-information-objective-infonce-contrastive-bounds.md) (sign-in required)
  - [🪞 Forward vs Reverse KL: Mode Covering, Mode Seeking and Failure Modes](https://everwise.study/courses/mathematics-for-machine-learning-5-information-theory-matrix-methods-and/forward-vs-reverse-kl-mode-covering-seeking.md)
  - [✅ Checkpoint: Information-Theoretic Objectives](https://everwise.study/courses/mathematics-for-machine-learning-5-information-theory-matrix-methods-and/checkpoint-information-theoretic-objectives.md)
- **🔢 Solving and Regularizing Linear Systems**
  - [🧱 QR and Cholesky: Factorizations That Solve](https://everwise.study/courses/mathematics-for-machine-learning-5-information-theory-matrix-methods-and/qr-cholesky-factorizations-that-solve.md)
  - [🔍 Full SVD: Pseudoinverse, Truncation and Rank Control](https://everwise.study/courses/mathematics-for-machine-learning-5-information-theory-matrix-methods-and/full-svd-pseudoinverse-truncation-rank-control.md)
  - [⚖️ Conditioning, Ill-Posed Problems and Tikhonov Regularization](https://everwise.study/courses/mathematics-for-machine-learning-5-information-theory-matrix-methods-and/conditioning-ill-posed-tikhonov-regularization.md)
  - [✅ Checkpoint: Linear Systems in Practice](https://everwise.study/courses/mathematics-for-machine-learning-5-information-theory-matrix-methods-and/checkpoint-linear-systems-in-practice.md)
- **🌊 Floating Point, Spectral Methods and Variational Inference**
  - [🧊 Floating Point, Catastrophic Cancellation and Mixed Precision](https://everwise.study/courses/mathematics-for-machine-learning-5-information-theory-matrix-methods-and/floating-point-catastrophic-cancellation-mixed-precision.md) — this lesson
  - [🌈 Spectral Methods: Eigenbases, Filtering and Smoothing](https://everwise.study/courses/mathematics-for-machine-learning-5-information-theory-matrix-methods-and/spectral-methods-eigenbasis-filtering-smoothing.md)
  - [🎯 First Look at Variational Inference: The ELBO as a KL Objective](https://everwise.study/courses/mathematics-for-machine-learning-5-information-theory-matrix-methods-and/first-look-variational-inference-elbo-kl.md) (premium content)
  - [✅ Checkpoint: Precision, Spectra and Variational Bounds](https://everwise.study/courses/mathematics-for-machine-learning-5-information-theory-matrix-methods-and/checkpoint-precision-spectra-variational-bounds.md) (premium content)
- [Course cheatsheet](https://everwise.study/courses/mathematics-for-machine-learning-5-information-theory-matrix-methods-and/cheatsheet.md) (premium content)