Floating point, cancellation, mixed precision, conditioning and backward error

Before this lesson: lesson 08 (the SVD), lesson 11 (norms) and lesson 12 (QR and least squares). Time: about 85 minutes, plus the notebook.

---

Why this matters #

Every number in a training run is rounded, and bf16 keeps only 8 significant bits (2–3 decimal digits). Rounding explains failures that look like bugs: a loss that becomes NaN in fp16, a running sum that stops growing, a variance that comes out negative, regression coefficients that change completely when one data point is added. Two ideas organize all of it. Conditioning says how sensitive the problem is: you lose about $\log_{10}\kappa$ digits, whatever algorithm you use. Stability says whether the algorithm adds more error than that.

1. Predict first #

2. The picture #

A 3-digit decimal computer. Suppose a calculator keeps 3 significant digits: $1.23 \times 10^5$, $4.56 \times 10^{-2}$. The gap between neighbours grows with size (0.01 near 1, 1000 near $10^5$), but the relative gap stays about the same. Rounding to the nearest number has a relative error of at most $u = \tfrac12 \times 10^{-2} = 0.005$.

Binary formats work the same way, with powers of 2 instead of powers of 10. Between $2^e$ and $2^{e+1}$ the numbers are evenly spaced; at each power of 2 the gap doubles.

One fp16 number, bit by bit. fp16 uses 16 bits: 1 sign bit $s$, 5 exponent bits $E$ and 10 fraction bits $f$. A normal number is

The 15 is the bias. It lets the stored $E$ be an unsigned integer from 0 to 31 while the true exponent $e = E - 15$ can be negative. Decode the bits 1 10001 1010000000:

1. Sign: $s = 1$, so the number is negative.
2. Exponent: $E = 10001_2 = 16 + 1 = 17$, so $e = 17 - 15 = 2$.
3. Fraction: $f = 0.1010000000_2 = 2^{-1} + 2^{-3} = 0.625$, so the significand is $(1.f)_2 = 1.625$. The leading 1 is not stored (the hidden bit), so fp16 has $p = 10 + 1 = 11$ bits of precision.
4. Value: $x = -1.625 \times 2^2 = -6.5$.

Two exponent values are special. $E = 0$ gives zero and the subnormal numbers (no hidden bit). $E = 31$ gives $\pm\infty$ (if $f = 0$) and NaN (if $f \ne 0$). So the largest fp16 number has $E = 30$ and all fraction bits 1: $(2 - 2^{-10}) \times 2^{15} = 65504$.

Interactive figure: loading…

Interactive figure: loading…

Cancellation on the 3-digit computer. The true values $1.2349$ and $1.2151$ are stored as $1.23$ and $1.22$, each with a relative error below 0.4%. The subtraction $1.23 - 1.22 = 0.01$ is exact, but the true difference is $0.0198$: the result is wrong by 49%. The leading digits cancel, and the old rounding errors become the main part of what is left.

A matrix as an amplifier. $A$ maps the unit circle to an ellipse with semi-axes $\sigma_1 \ge \sigma_2$. Solving $Ax = b$ runs the map backwards. If $b$ lies along the long axis, $x$ is small ($\|x\| = \|b\|/\sigma_1$). A change $\delta b$ along the short axis gives a large change $\|\delta x\| = \|\delta b\|/\sigma_2$. So the relative error can grow by $\sigma_1/\sigma_2 = \kappa$. A flat ellipse (large $\kappa$) means a nearly singular matrix.

Interactive figure: loading…

3. Definitions and main results #

Floating-point numbers #

A binary floating-point number is $x = \pm(1.f_1f_2\dots f_{p-1})_2 \times 2^e$: a significand with $p$ bits of precision (the leading 1 is not stored) and an exponent $e$ in a fixed range.

format	bits s/e/f	$\varepsilon = 2^{1-p}$
fp64 (double)	1/11/52	$2^{-52} \approx 2.2\times10^{-16}$
fp32	1/8/23	$2^{-23} \approx 1.2\times10^{-7}$
fp16	1/5/10	$2^{-10} \approx 9.8\times10^{-4}$
bf16	1/8/7	$2^{-7} \approx 7.8\times10^{-3}$

format	largest	smallest normal
fp64 (double)	$1.8\times10^{308}$	$2.2\times10^{-308}$
fp32	$3.4\times10^{38}$	$1.2\times10^{-38}$
fp16	$65504$	$6.1\times10^{-5}$
bf16	$3.4\times10^{38}$	$1.2\times10^{-38}$

"s/e/f" counts the sign, exponent and fraction bits; the precision is $p$ = fraction bits + 1. Below the smallest normal number come subnormal numbers with fewer significant bits; the smallest positive fp16 number is $2^{-24} \approx 6\times10^{-8}$.

• Machine epsilon $\varepsilon = 2^{1-p}$ is the gap between 1 and the next larger float: $2^{-52}$ in double. It is not the smallest relative gap: the relative gap varies between $\varepsilon/2$ and $\varepsilon$.
• Spacing. In $[2^e, 2^{e+1})$ the floats are evenly spaced with gap $\varepsilon\,2^e = 2^{e-p+1}$. So the gap near $x$ lies between $\tfrac12\varepsilon|x|$ and $\varepsilon|x|$: the absolute gap grows with $|x|$, the relative gap stays near $\varepsilon$. Example: the gap at 1024 in double is $2^{-42} \approx 2.3\times10^{-13}$.
• Unit roundoff $u = \varepsilon/2$ ($2^{-53} \approx 1.1\times10^{-16}$ in double) bounds the relative error of rounding to nearest: $\operatorname{fl}(x) = x(1 + \delta)$ with $|\delta| \le u$ (no overflow or underflow).
• Standard model. IEEE 754 rounds the exact result of $+, -, \times, /$ and $\sqrt{\ }$: $\operatorname{fl}(x \circ y) = (x \circ y)(1 + \delta)$, $|\delta| \le u$. One operation is always accurate; trouble comes from how many operations combine.

Why $0.1 + 0.2 \ne 0.3$. $0.1 = (0.0\overline{0011})_2$ has no finite binary expansion, so it is rounded when stored ($0.1000000000000000055\ldots$), and so are $0.2$ and $0.3$. The sum of the stored $0.1$ and $0.2$ rounds to $0.3000000000000000444\ldots$, while the stored $0.3$ is $0.2999999999999999888\ldots$: two neighbouring doubles. The cause is that these decimals are not representable in binary, not the size of $\varepsilon$.

Cancellation #

Subtracting nearly equal floats is exact (Sterbenz: if $y/2 \le x \le 2y$, then $x - y$ is a float). The damage comes from the errors already in the inputs. If $\hat x = x(1 + \delta_1)$ and $\hat y = y(1 + \delta_2)$ with $|\delta_i| \le u$, then $(\hat x - \hat y) - (x - y) = x\delta_1 - y\delta_2$, so

The factor $(|x| + |y|)/|x - y|$ is the condition number of subtraction. When $x \approx y$ it is huge: this is catastrophic cancellation.

The fix: rewrite so that no nearly equal numbers are subtracted.

• Algebra: $\sqrt{x + 1} - \sqrt{x} = 1/(\sqrt{x + 1} + \sqrt{x})$; $1 - \cos x = 2\sin^2(x/2)$.
• Library functions for small $x$: log1p(x) $= \log(1 + x)$, expm1(x) $= e^x - 1$.
• Quadratic $ax^2 + bx + c$: compute the root of larger magnitude with the sign that adds ($-b$ and $-\operatorname{sign}(b)\sqrt{b^2 - 4ac}$ have the same sign), then the other root from Vieta's formula $x_1x_2 = c/a$.
• Variance: $\frac1N\sum x_i^2 - \bar x^2$ subtracts two large, nearly equal numbers when the mean is large compared with the spread. Use two passes, $\frac1N\sum(x_i - \bar x)^2$, or Welford's one-pass update of a running mean and a running sum of squared deviations.

Interactive figure: loading…

Conditioning against stability, in one example. The small root of $x^2 - 10^8x + 1 = 0$ is about $10^{-8}$. It is a well-conditioned function of the coefficients: a relative change of $b$ or $c$ changes the root by the same relative amount (condition number 1). Still, the textbook formula $(10^8 - \sqrt{10^{16} - 4})/2$ gives $7.45 \times 10^{-9}$ in double, 25% wrong. The formula is an unstable algorithm for a well-conditioned quantity: the formula itself creates the cancellation. The rewrite $2/(10^8 + \sqrt{10^{16} - 4})$ changes the algorithm, not the problem.

Mixed precision #

Training in bf16 or fp16 halves memory traffic and runs on fast tensor cores, but it needs three safeguards. The first one comes from a simple fact: a 16-bit running sum stops growing.

Why a 16-bit sum stops growing, step by step. fp16 has $p = 11$.

1. In $[2^e, 2^{e+1})$ the fp16 numbers are $2^{e-p+1} = 2^{e-10}$ apart.
2. Below 2048 (in $[1024, 2048)$) the gap is $2^{0} = 1$. From $2048 = 2^{11}$ on, the gap is $2^{1} = 2$: the numbers are $\dots, 2046, 2047, 2048, 2050, 2052, \dots$
3. Up to 2048, each $s + 1$ is a float, so the sum is exact.
4. $2048 + 1 = 2049$ lies exactly halfway between 2048 and 2050. A tie goes to the neighbour whose significand ends in 0 (round to even): $2048 = (1.0000000000)_2 \times 2^{11}$ ends in 0, and $2050 = (1.0000000001)_2 \times 2^{11}$ ends in 1. So $2048 + 1$ rounds to $2048$.
5. From now on, every addition gives 2048 again. The sum is stuck.

bf16 has $p = 8$. The gap at $256 = 2^8$ is $2^{8-7} = 2$, so $256 + 1$ rounds back to 256: bf16 sums of ones stop at 256. In general, adding $t$ to $s$ is lost when $|t|$ is less than half the gap at $s$ (this loss is called swamping). fp32 has $p = 24$, so the same sum of ones only stops at $2^{24} \approx 1.7 \times 10^7$.

Interactive figure: loading…

The three safeguards:

1. fp32 accumulation. Matmuls take 16-bit inputs but accumulate in fp32. Large reductions (softmax sums, normalization statistics, the loss) also run in fp32. The error bound for a sum of $n$ terms is about $(n-1)u\sum|x_i|$; with bf16's $u = 2^{-8}$ this bound exceeds the sum itself once $n$ passes a few hundred.
2. fp32 master weights. An update $\eta g$ vanishes when $|\eta g|$ is below half the gap at $w$. In bf16 the gap above $w = 1$ is $2^{-7} \approx 0.0078$, so $1 + 10^{-3}$ rounds back to 1. The optimizer keeps an fp32 copy of the weights, applies updates there, and casts down for the forward and backward pass.
3. Loss scaling (fp16). fp16 cannot represent magnitudes below about $6\times10^{-8}$, and it loses precision below $6.1\times10^{-5}$. Small gradients underflow: $10^{-8}$ becomes $0$. Multiply the loss by $S$ (a power of 2, so the scaling itself is exact). By the chain rule every gradient is multiplied by $S$ and moves into range ($10^{-8}\cdot 2^{16} \approx 6.6\times10^{-4}$). Unscale in fp32 before the update. Dynamic loss scaling (PyTorch GradScaler starts at $S = 2^{16}$) halves $S$ and skips the step when an inf or NaN appears (overflow above 65504), and doubles $S$ after 2000 steps without one.

Conditioning #

The relative condition number of a problem $x \mapsto f(x)$ is the worst amplification of a small relative input change:

For a scalar differentiable $f$: $\kappa_f(x) = |x f'(x)/f(x)|$.

For an invertible matrix $A$, the condition number is $\kappa(A) = \|A\|\,\|A^{-1}\|$ (an induced norm). In the 2-norm, $\|A\|_2 = \sigma_1$ and $\|A^{-1}\|_2 = 1/\sigma_n$ (because $A^{-1} = V\Sigma^{-1}U^\top$), so

$\kappa_2(Q) = 1$ for orthogonal $Q$, $\kappa(cA) = \kappa(A)$ for $c \ne 0$, and $\kappa = \infty$ for a singular $A$. For a rectangular $A$ with full column rank, $\kappa_2(A) = \sigma_1/\sigma_n$ as well. The determinant does not measure conditioning: $0.1\,I_{100}$ has $\det = 10^{-100}$ and $\kappa = 1$.

Theorem (error bound for $Ax = b$) #

Let $A$ be invertible, $Ax = b \ne 0$, and $A(x + \delta x) = b + \delta b$. Then

and in the 2-norm the bound is attained.

Proof. Each line says which fact it uses.

1. $A\,\delta x = \delta b$. Subtract $Ax = b$ from $A(x + \delta x) = b + \delta b$.
2. $\delta x = A^{-1}\delta b$. $A$ is invertible.
3. $\|\delta x\| \le \|A^{-1}\|\,\|\delta b\|$. An induced norm satisfies $\|Mv\| \le \|M\|\,\|v\|$ (lesson 11); use it with $M = A^{-1}$.
4. $\|b\| \le \|A\|\,\|x\|$. The same inequality with $M = A$, since $b = Ax$.
5. $\dfrac{1}{\|x\|} \le \dfrac{\|A\|}{\|b\|}$. Divide line 4 by $\|x\|\,\|b\|$; both are positive because $b \ne 0$.
6. $\dfrac{\|\delta x\|}{\|x\|} \le \|A\|\,\|A^{-1}\|\,\dfrac{\|\delta b\|}{\|b\|} = \kappa(A)\,\dfrac{\|\delta b\|}{\|b\|}$. Multiply lines 3 and 5 (both sides are nonnegative).
7. Sharpness in the 2-norm. Lines 3 and 4 are equalities for the right vectors. Take $b = u_1$: then $x = v_1/\sigma_1$ and $\|b\| = \|A\|\,\|x\|$. Take $\delta b = \eta\,u_n$: then $\delta x = \eta\,v_n/\sigma_n$ and $\|\delta x\| = \|A^{-1}\|\,\|\delta b\|$. The ratio is exactly $\sigma_1/\sigma_n$. $\blacksquare$

A $2 \times 2$ check. $A = \operatorname{diag}(1, 0.01)$ has $\sigma_1 = 1$, $\sigma_2 = 0.01$, $\kappa_2 = 100$.

• $b = (1, 0)$, so $x = (1, 0)$. Take $\delta b = (0, 0.01)$: then $\delta x = (0, 1)$. The relative change of $b$ is $0.01$; the relative change of $x$ is $1 = 100 \times 0.01$. The bound is attained.
• Same $b$, but $\delta b = (0.01, 0)$: then $\delta x = (0.01, 0)$, a relative change of $0.01$. The same size of $\delta b$ has a 100 times smaller effect: the direction matters.
• $b = (0, 1)$, so $x = (0, 100)$, and $\delta b = (0, 0.01)$: $\delta x = (0, 1)$, a relative change of $0.01$ again. The worst case needs $b$ along $u_1$ (small $x$) and $\delta b$ along $u_n$.

Perturbing $A$ instead of $b$

If $(A + \delta A)(x + \delta x) = b$, then $A\,\delta x + \delta A\,(x + \delta x) = 0$, so $\delta x = -A^{-1}\delta A\,(x + \delta x)$ and

$$\frac{\|\delta x\|}{\|x + \delta x\|} \le \kappa(A)\,\frac{\|\delta A\|}{\|A\|}.$$

To first order, $\|\delta x\|/\|x\| \lesssim \kappa(A)\,\|\delta A\|/\|A\|$. So $\kappa(A)$ also measures the sensitivity to changes in the matrix.

Theorem: $\kappa_2(A^\top A) = \kappa_2(A)^2$ #

Let $A$ have full column rank and SVD $A = U\Sigma V^\top$. Then $A^\top A = V\Sigma^\top\Sigma V^\top = V\operatorname{diag}(\sigma_i^2)V^\top$. This is an SVD (and an eigendecomposition) of the symmetric positive definite $A^\top A$, so its singular values are $\sigma_i^2$ and $\kappa_2(A^\top A) = \sigma_1^2/\sigma_n^2 = \kappa_2(A)^2$. $\blacksquare$

Meaning for the normal equations. Any method that works with $A^\top A$ meets condition number $\kappa^2$. With $\kappa(A) = 10^6$, QR keeps about 10 digits and the normal equations about 4. For $\kappa(A) \gtrsim 1/\sqrt{u} \approx 10^8$, the computed $A^\top A$ can be exactly singular (the Läuchli example of lesson 12). The digits are lost when $A^\top A$ is formed, before any solver runs.

Backward error and stability #

The picture first. A solver returns $\hat x$ instead of the exact $x$. Ask a different question: for which data is $\hat x$ the exact answer? Put $\hat b = A\hat x$. Then $\hat x$ solves $A\hat x = \hat b$ exactly. If $\hat b$ is very close to $b$, the solver did its job well: it solved a nearby problem exactly. The relative size of that change of the data is the backward error.

Interactive figure: loading…

Now the definitions. For a computed $\hat x$:

• the forward error is $\|\hat x - x\|/\|x\|$: how wrong the answer is. It is usually unknown, because $x$ is unknown.
• the backward error is the smallest relative change of the data that makes $\hat x$ exact. It is computable from the residual $r = b - A\hat x$.
  • If only $b$ may change: $\hat b - b = -r$, so the backward error is $\|r\|/\|b\|$.
  • If only $A$ may change: $\min\{\|\delta A\|_2/\|A\|_2 : (A + \delta A)\hat x = b\} = \|r\|_2 / (\|A\|_2\,\|\hat x\|_2)$. This is the usual definition for linear solvers.
• An algorithm is backward stable if its backward error is $O(u)$ for every input. Householder QR and SVD solvers are backward stable; Gaussian elimination with partial pivoting is backward stable in practice. The normal equations are not backward stable for least squares.

Why the smallest change of $A$ has size $\|r\|/\|\hat x\|$

1. Any $\delta A$ with $(A + \delta A)\hat x = b$ satisfies $\delta A\,\hat x = b - A\hat x = r$.
2. So $\|r\| = \|\delta A\,\hat x\| \le \|\delta A\|\,\|\hat x\|$, that is, $\|\delta A\| \ge \|r\|/\|\hat x\|$.
3. The rank-1 matrix $\delta A = r\hat x^\top/\|\hat x\|^2$ reaches this value: $\delta A\,\hat x = r$, and the 2-norm of $uv^\top$ is $\|u\|\,\|v\|$, so $\|\delta A\|_2 = \|r\|\,\|\hat x\|/\|\hat x\|^2 = \|r\|/\|\hat x\|$.
4. Divide by $\|A\|_2$. $\blacksquare$

Forward error from backward error. The computed $\hat x$ is the exact solution for the data $\hat b = b - r$. Apply the error bound theorem with $\delta b = -r$:

The same holds to first order for changes in $A$. This gives the rule of thumb:

Conditioning belongs to the problem, stability to the algorithm, and the error is roughly their product. In double precision you lose about $\log_{10}\kappa$ of the 16 digits. A small residual does not mean an accurate $x$ when $\kappa$ is large.

Answers to section 1. (1) $s = 2048$. From 2048 on, the gap between fp16 numbers is 2, so $2048 + 1$ is a tie and rounds back to 2048 (ties go to the even significand). (2) About $16 - 10 = 6$ digits (forward error $\approx \kappa u \approx 10^{-6}$). Yes, the residual is tiny, about $u\,\|A\|\,\|x\|$, because the solver is backward stable.

4. Worked example #

Two nearly parallel columns.

Change $b_2$ in its fifth digit, $b' = (2, 2.0002)$. Then $x' = (0, 2)$.

• $\det A = 10^{-4}$ and $A^{-1} = \begin{pmatrix} 10001 & -10000 \\ -10000 & 10000\end{pmatrix}$.
• $A$ is symmetric positive definite, so its singular values are its eigenvalues, the roots of $\lambda^2 - 2.0001\lambda + 10^{-4} = 0$: $\sigma_1 \approx 2.00005$ and $\sigma_2 \approx 5.0\times10^{-5}$. So $\kappa_2(A) \approx 4.0\times10^4$.
• Bound: $\|\delta b\|/\|b\| = 10^{-4}/2.83 = 3.5\times10^{-5}$; times $\kappa$ gives $1.41$. Actual: $\|\delta x\|/\|x\| = \|(-1, 1)\|/\|(1, 1)\| = 1$. The bound holds and is nearly sharp.
• Normal equations would face $\kappa_2(A^\top A) = \kappa_2(A)^2 \approx 1.6\times10^9$.
• In double, np.linalg.solve(A, b) has relative error $2.2\times10^{-12}$, inside $\kappa u \approx 4\times10^4 \times 1.1\times10^{-16} = 4.4\times10^{-12}$.

In a regression the two columns are two almost identical features (say, temperature at two neighbouring grid points). The data cannot separate their effects, so the two coefficients follow the noise.

The same example, read backwards. Suppose a solver had returned $\hat x = (0, 2)$ for the original $b$. Step by step:

1. Residual: $r = b - A\hat x = (2 - 2,\; 2.0001 - 2.0002) = (0, -10^{-4})$.
2. Backward error in $b$: $\|r\|/\|b\| = 10^{-4}/2.83 = 3.5\times10^{-5}$. With changes in $A$: $\|r\|/(\|A\|_2\|\hat x\|) = 10^{-4}/(2.00005 \times 2) = 2.5\times10^{-5}$.
3. Forward error: $\|\hat x - x\|/\|x\| = 1$, that is, 100% wrong.
4. Check the bound: $\kappa \times 3.5\times10^{-5} \approx 1.41 \ge 1$. ✓

So an answer that is completely wrong is the exact answer of a problem whose data differs by 0.0035%. If $b$ comes from a measurement with 4 digits, $\hat x$ is as good an answer as $x$.

Cancellation and its fix. $f(x) = (1 - \cos x)/x^2$ tends to $1/2$ as $x \to 0$. It is very well conditioned there: $f(x) \approx \tfrac12 - x^2/24$, so $\kappa_f(x) = |xf'(x)/f(x)| \approx x^2/6$.

$x$	naive $(1 - \cos x)/x^2$	rewrite $2\sin^2(x/2)/x^2$
$10^{-4}$	$0.4999999969612645$ (rel. error $5\times10^{-9}$)	$0.4999999995833334$
$10^{-6}$	$0.5000444502911705$ (rel. error $9\times10^{-5}$)	$0.4999999999999583$
$10^{-8}$	$0$ (rel. error 1)	$0.5$

At $x = 10^{-8}$, $\cos x = 1 - 5\times10^{-17}$ rounds to exactly 1 (the gap below 1 is $1.1\times10^{-16}$), so all information is gone before the subtraction. The rewrite has no subtraction and is correct to about $10^{-16}$ for every $x$. The problem was well conditioned; only the formula was unstable.

5. In code #

NumPy shows the facts of this lesson. NumPy has no bfloat16 type, so a three-line helper rounds float32 values to bf16 by cutting the last 16 bits with round-to-nearest-even.

import numpy as np
from decimal import Decimal

def bf16(x):
    """Round to the nearest bfloat16 (ties to even); the result is a float32 that bf16 can store."""
    b = np.asarray(x, dtype=np.float32).view(np.uint32).astype(np.uint64)
    b = (b + 0x7FFF + ((b >> 16) & 1)) & 0xFFFF0000
    return b.astype(np.uint32).view(np.float32)

for t in (np.float64, np.float32, np.float16):
    fi = np.finfo(t)
    print(t.__name__, fi.eps, fi.max, fi.smallest_normal)   # ε, largest, smallest normal
print("bfloat16 eps", bf16(1 + 2**-7) - 1)            # 0.0078125 = 2^-7
print(np.spacing(1.0), np.spacing(1024.0))            # 2.2e-16, 2.3e-13: the gap grows with |x|
print(Decimal(0.1 + 0.2))                             # 0.3000000000000000444...
print(Decimal(0.3))                                   # 0.2999999999999999888...: a different double

bits = np.array([-6.5], dtype=np.float16).view(np.uint16)[0]
print(f"{bits:016b}")                                 # 1100011010000000 = 1 | 10001 | 1010000000

s = np.float16(0)
for _ in range(10_000):
    s = np.float16(s + np.float16(1))
print(s)                                              # 2.048e+03 (= 2048): swamping
print(bf16(bf16(1.0) + bf16(1e-3)))                   # 1.0: the update is lost
print(np.float16(1e-8), np.float16(1e-8 * 2**16))     # 0.0 (underflow), 0.000655 (after scaling)
print(bf16(1e-8))                                     # 1.0011718e-08: bf16 has the range of fp32

The cancellation table:

import numpy as np

np.set_printoptions(precision=16)
x = np.array([1e-4, 1e-6, 1e-8])
print((1 - np.cos(x)) / x**2)          # [0.4999999969612645 0.5000444502911705 0.]: wrong digits
print(2 * np.sin(x / 2)**2 / x**2)     # [0.4999999995833334 0.4999999999999583 0.5]: all correct

Conditioning, and a tiny backward error with a large forward error:

import numpy as np
from scipy.linalg import hilbert

A = np.array([[1., 1.], [1., 1.0001]])
b = np.array([2., 2.0001])
print(np.linalg.cond(A), np.linalg.cond(A.T @ A))     # 4.0e4 and 1.6e9 = κ²
print(np.linalg.solve(A, b))                          # [1. 1.]
print(np.linalg.solve(A, [2., 2.0002]))               # [0. 2.] (up to ~1e-16)

x_hat = np.array([0., 2.])                            # a wrong answer for the original b
r = b - A @ x_hat
print(f"{np.linalg.norm(r) / np.linalg.norm(b):.1e}")  # 3.5e-05: backward error in b; forward error is 1

H = hilbert(10)                                       # κ ≈ 1.6e13
x_true = np.ones(10)
b = H @ x_true
x = np.linalg.solve(H, b)                             # LU with partial pivoting
backward = np.linalg.norm(b - H @ x) / (np.linalg.norm(H, 2) * np.linalg.norm(x))
forward = np.linalg.norm(x - x_true) / np.linalg.norm(x_true)
print(f"cond {np.linalg.cond(H):.1e}  backward {backward:.0e}  forward {forward:.0e}")
# cond 1.6e+13  backward 1e-16  forward 4e-04: tiny residual, only about 3-4 correct digits

6. Where it shows up #

• Mixed-precision training. torch.autocast runs matmuls in bf16 or fp16 and keeps numerically sensitive ops such as softmax, layer norm and the losses in fp32; GradScaler is for fp16. fp8 formats (E4M3, E5M2) have even less range and precision, so fp8 training uses per-tensor or per-block scale factors: the loss-scaling idea, applied to every tensor.
• Stable formulas in ML code. Log-sum-exp subtracts the maximum before exp (an overflow fix); log1p, expm1 and softplus avoid cancellation in losses and log-probabilities; normalization layers and streaming statistics use two-pass or Welford variance.
• ML for science. Physical fields have large means and small anomalies: temperature near 288 K with anomalies of about 1 K, surface pressure near $10^5$ Pa. A variance computed as $\overline{x^2} - \bar x^2$ in fp32 loses about 5 of its 7 digits; work with anomalies. Kernel (Gram) matrices of smooth kernels are extremely ill-conditioned, which is why Gaussian processes add a small jitter $\delta I$. Discretized differential operators have $\kappa$ growing like $h^{-2}$ with the grid spacing $h$, which slows iterative solvers (lesson 14).
• Regression diagnostics. A large $\kappa(X)$ means some combination of features is nearly zero. Rescale columns, remove collinear features, or use ridge: $X^\top X + \lambda I$ has condition number $(\sigma_1^2 + \lambda)/(\sigma_n^2 + \lambda)$.

7. Common mistakes #

• Taking $\varepsilon$ for the smallest positive float. $\varepsilon$ is the gap at 1; the smallest positive double is about $5\times10^{-324}$ (a subnormal). Also, books use "machine precision" for both $\varepsilon$ and $u = \varepsilon/2$: say which one you mean.
• Blaming the problem for an unstable formula, or the algorithm for an ill-conditioned problem. Ask both questions: how large is $\kappa$, and does the algorithm add error beyond $\kappa u$?
• Reading a small residual as an accurate solution. The residual measures backward error; the forward error can be $\kappa$ times larger.
• Judging near-singularity by $\det A$. Use $\kappa$, or $\sigma_{\min}/\sigma_{\max}$.
• Treating loss scaling as the bf16 fix while sums run in 16 bits. bf16's problem is precision, not range: accumulate in fp32.

8. Check yourself #

1. What is the gap between doubles near $1000$?
   Answer

   $1000 \in [2^9, 2^{10})$, so the gap is $\varepsilon \cdot 2^9 = 2^{-43} \approx 1.14\times10^{-13}$.
2. Rewrite $\sqrt{x + 1} - \sqrt{x}$ for large $x$. How bad is the naive formula at $x = 10^{12}$?
   Answer

   $1/(\sqrt{x + 1} + \sqrt{x})$. The true value is $\approx 5.0\times10^{-7}$; the naive formula has relative error $\approx 7.6\times10^{-6}$ (the two square roots agree in their first 13 digits), the rewrite $\approx 10^{-16}$.
3. What are $\kappa_2(\operatorname{diag}(10, 0.1))$ and $\kappa_2(0.1\,I_{100})$?
   Answer

   $10/0.1 = 100$, and $1$: all singular values of $0.1\,I$ are equal, although $\det = 10^{-100}$.
4. Why does bf16 training usually need no loss scaling, but still needs fp32 accumulation?
   Answer

   bf16 has 8 exponent bits like fp32, so small gradients do not underflow. It has only 8 bits of precision, so long sums swamp ($256 + 1 = 256$) and updates below half the gap at $w$ vanish: accumulate and update in fp32.
5. A solver returns $\hat x$ with $\|b - A\hat x\|/(\|A\|\|\hat x\|) = 10^{-16}$, and $\kappa(A) = 10^{13}$. Is $\hat x$ accurate?
   Answer

   The backward error is tiny (the solver is backward stable), but the forward error can be up to $\kappa \times 10^{-16} = 10^{-3}$: only about 3 correct digits. The Hilbert matrix of size 10 shows exactly this ($4\times10^{-4}$).

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. What are machine epsilon $\varepsilon$ and the unit roundoff $u$ in double precision?

   Show answer

2. Decode the fp16 bits 0 10000 1000000000 (sign, exponent, fraction; bias 15).

   Show answer

3. What is catastrophic cancellation?

   Show answer

4. What does bf16 trade compared with fp16?

   Show answer

5. Why do 16-bit matmuls and reductions accumulate in fp32?

   Show answer

6. Write $\kappa_2(A)$ with singular values.

   Show answer

7. State the error bound for $Ax = b$ when $b$ changes by $\delta b$.

   Show answer

8. What is $\kappa_2(A^\top A)$, and why does it matter?

   Show answer

9. What is the backward error of a computed solution $\hat x$ of $Ax = b$?

   Show answer

10. Conditioning versus stability: which belongs to the problem, which to the algorithm, and how do they combine?

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "Why does fp16 training need loss scaling and bf16 does not?" — fp16 has 5 exponent bits (smallest subnormal $6\times10^{-8}$, max 65504), so small gradients underflow; bf16 has fp32's 8 exponent bits. Both need fp32 accumulation and master weights.
• "What is the condition number, and why should an ML engineer care?" — $\sigma_{\max}/\sigma_{\min}$; you lose about $\log_{10}\kappa$ digits, $\kappa(X^\top X) = \kappa(X)^2$, and gradient descent on a quadratic converges at a rate set by the Hessian's $\kappa$.
• "Your loss became NaN after switching to fp16. What do you check?" — overflow above 65504 in activations or logits, missing or badly tuned loss scaling, softmax or loss computed in 16 bits, and constants such as an $\epsilon = 10^{-8}$ that become 0 in fp16.

10. Practice #

The notebook for this lesson runs in Colab. It has five exercises, each with a check cell and a folded solution:

1. Find machine epsilon of float16, float32 and float64 by halving until $1 + e = 1$, and compare with np.finfo.
2. Rewrite $(1 - \cos x)/x^2$ and $\sqrt{x + 1} - \sqrt{x}$ stably, and test them against exact values from Python's decimal module.
3. Write a stable quadratic solver and compare it with the textbook formula on $x^2 - 10^8x + 1$.
4. Write Welford's one-pass variance and compare it with $\overline{x^2} - \bar x^2$ on data with a large mean.
5. Write condition_number and backward_error, check $\kappa(A^\top A) = \kappa(A)^2$ on random matrices, and compare backward and forward error on the Hilbert matrix.

11. Watch and read #

Floating Point Numbers - ComputerphileComputerphile · 9 min

Harvard AM205 video 0.4 - Finite-precision arithmeticChris Rycroft · Harvard AM205 (numerical methods) · 31 min

1 4 1 The condition number of a matrixAdvanced LAFF · UT Austin, Advanced Linear Algebra (week 1) · 8 min

Lecture 10: Survey of Difficulties with Ax = bGilbert Strang · MIT 18.065, Matrix Methods in Data Analysis (2018) · 50 min

• Driscoll and Braun, Fundamentals of Numerical Computation, free online: §1.1 floating-point numbers, §1.2 conditioning, §1.3 algorithms, §1.4 stability and §2.8 matrix condition numbers.
• David Goldberg, "What every computer scientist should know about floating-point arithmetic" (1991), section "Rounding Error" (relative error and ulps, cancellation).
• Nick Higham, "What is backward error?" (2020); for the full theory, Trefethen and Bau, Numerical Linear Algebra, Lectures 12–15 and 18.
• Micikevicius et al., "Mixed Precision Training" (ICLR 2018), §3.1–3.3 (fp32 master weights, loss scaling, fp32 accumulation).