Numerical Stability and Conditioning
Dive into Deep Learning · §25.5
Why the math is right but the loss is NaN
floating point · stable softmax · cancellation · conditioning.
The math is right; the loss is NaN
Motivation
Every proof in this chapter was over \mathbb{R}. Your GPU computes over a finite, gappy imitation. When an answer is wrong, two questions split the blame:
- Did the algorithm solve a nearby problem? (backward error)
- Do nearby problems have wildly different answers? (conditioning)
The fixes are reformulations: max-subtraction, log-space, Welford, ridge.
01
Floating point
a number system with gaps
A number system with gaps
Floating point
A float is base-2 scientific notation with a fixed digit budget:
x = (-1)^s\,(1.m_1\ldots m_p)_2\; 2^{e}.
The exponent buys range; the mantissa fixes relative precision. Between powers of two the spacing is constant, so the absolute gap doubles at every power, and each format ends at an overflow threshold.
Machine epsilon \varepsilon_{\text{mach}} = 2^{-p} is the gap from 1 to its successor: \mathrm{fl}(x) = x(1+\delta), |\delta| \le \tfrac12 \varepsilon_{\text{mach}}.
fp8: E4M3 and E5M2
Floating point
Since 2022 hardware also offers fp8, in two flavors with a division of labor:
E4M3 keeps digits: \varepsilon = 0.125 (about one decimal), max = 448: for weights and activations. E5M2 trades a mantissa bit for fp16’s full range: max = 57344, smallest normal 6.1\times10^{-5}, at \varepsilon = 0.25: for gradients, which need range.
At \varepsilon_{\text{mach}} = 0.125 there is no slack left: fp8 training must pair every tensor with a scale factor. (The ml_dtypes package provides both formats.)
Where the thresholds are
Floating point
Because e^x turns additive scale into multiplicative scale, a modest logit overflows: fp32 dies at x \approx 88.7, fp16 at x \approx 11.1.
1 + eps != 1 : True
1 + eps/2 == 1 : True
gap between adjacent float32 values near 1 : 1.1920929e-07
gap between adjacent float32 values near 2^20 : 0.125
float16: exp(x) overflows for x > 11.09
float32: exp(x) overflows for x > 88.72
fp16 gradients below 6\times10^{-5} vanish; updates of relative size below \varepsilon_{\text{mach}}/2 round to no update at all. Both limits bite mixed-precision training (next slide).
Both fp16 failure modes, both rescues
Floating point
The two failure modes, and their two escapes, in mixed-precision training:
- A true gradient of 10^{-8} underflows an fp16 backward pass to an exact 0.0; multiplying the loss by 2^{14} before differentiating (and unscaling after) shifts the whole gradient chain into representable territory and recovers 1.000\times10^{-8}.
- A healthy update of relative size 10^{-4} is lost to round-to-nearest in fp16 (w - \eta g = w exactly); an fp32 master copy of the weights accepts it without complaint.
Loss scaling is underflow management; master weights are rounding management: this is what the library’s mixed-precision utilities automate, and what bfloat16’s fp32-sized exponent makes unnecessary.
02
Softmax & cross-entropy
the one-line bug, and the shift that fixes it
Softmax overflows; subtract the max
Stable softmax
The most common stability bug is one line: \mathrm{softmax} exponentiates logits, so any logit past 88.7 makes the numerator inf and the ratio NaN. But softmax is shift-invariant:
\mathrm{softmax}(\mathbf{z} - c\mathbf{1}) = \mathrm{softmax}(\mathbf{z}),
so, for finite logits, shift by c = \max_i z_i: every exponent is \le 0, the denominator sits in [1,n], and exponential overflow is avoided.
naive, logits z : [0.09003058 0.24472846 0.66524094]
naive, logits z + 100: [nan nan nan]
stable, logits z + 100: [0.09003057 0.24472846 0.66524094]
naive and stable agree where both work: False
Log-sum-exp: an exact, safe identity
Stable softmax
The softmax normalizer earns its own operator. The same shift rewrites it exactly:
\mathrm{lse}(\mathbf{z}) = \log\textstyle\sum_j e^{z_j}
= c + \log\textstyle\sum_j e^{z_j - c},
\qquad \max_j z_j \le \mathrm{lse}(\mathbf{z}) \le \max_j z_j + \log n.
Logits near 1000 overflow even float64; in log space they are effortless:
naive log(sum(exp(z))) : inf
stable log_sum_exp(z) : 1002.4076
log softmax : [-2.4075928 -1.4075928 -0.40759277]
probabilities sum to 1 : 1.000013
A soft maximum, within \log n of the true max, and the reason naive Bayes sums logs instead of multiplying probabilities.
Pass logits, not probabilities
Stable softmax
Cross-entropy is computable straight from logits with one stable lse:
-\log\mathrm{softmax}(\mathbf{z})_y = \mathrm{lse}(\mathbf{z}) - z_y.
The via-probabilities route forces the loss through the representable range of probabilities and fails in one of three ways, depending on the library: subnormal noise before inf, inf outright, or a silent clip:
gap CE from logits CE via probabilities
20 20.0000 20.0000
60 60.0000 60.0000
103 103.0000 inf
104 104.0000 inf
How the from-probabilities route fails
Stable softmax
The label is the unlikely class, so the true loss is the logit gap. From logits it is exact at every gap; via probabilities it fails:
MXNet. The softmax does not linger in the subnormal range: e^{-t} underflows to exactly 0 already at gap 103, and the loss reads inf at gaps 103 and 104.
The lesson generalizes: losses and likelihoods should live in log space from birth; convert to probabilities last.
03
Catastrophic cancellation
subtracting near-equal numbers
Subtraction annihilates digits
Cancellation
Subtracting nearly equal numbers is exact, yet it strips the leading digits they agreed on, exposing the trailing noise. Relative error is amplified by \tfrac{|a| + |b|}{|a - b|}\,u, which blows up precisely when a \approx b. In float32, 1 + 10^{-8} rounds to 1, so \log(1+x) returns 0, but log1p is exact:
float32 rounds 1 + x to : 1.0
log(1 + x) = 0.0 log1p(x) = 1e-08
a - b in float32 : 2.3841858e-07 (true value 3.0e-07)
amplification (|a|+|b|)/|a-b| ~ 8.4e+06
Standard victims: \log(1{+}x), e^x{-}1 near 0 (log1p, expm1); 1{-}\cos x; the quadratic formula near a double root. Reformulate; don’t add bits.
Welford beats E[x²] − E[x]²
Cancellation
The one-pass variance formula \mathbb{E}[x^2] - \mathbb{E}[x]^2 subtracts two numbers near \mu^2 to get \sigma^2, amplification \mu^2/\sigma^2. Welford keeps a running mean and centered sum of squares, so nothing large is ever subtracted:
m_k = m_{k-1} + \frac{x_k - m_{k-1}}{k},
\qquad M_k = M_{k-1} + (x_k - m_{k-1})(x_k - m_k).
Mean 10^9, true variance 1, 10^5 samples, all in float64:
naive E[x^2] - E[x]^2 : -256.000000
Welford, one pass : 1.000257
two-pass reference : 1.000257
The naive answer is pure amplified noise; here it even comes out negative (-256), a variance below zero, its sign hostage to the summation order. This is what BatchNorm avoids with running moments.
Summation order is an algorithm
Cancellation
The naive variance’s noise changed sign between NumPy builds. That is summation error at work. Adding n floats left to right commits one rounding per addition, worst case \approx n\,u; the repairs reorganize the additions:
left-to-right O(n\,u) · pairwise (sum halves recursively) O(u \log n) (what NumPy’s sum does, blocking and all) · Kahan (carry each rounding in a second accumulator) O(u), independent of n
Welford composes with either: the pairwise merge rule is exactly how running moments are combined across devices.
04
Conditioning
backward error, forward error, and κ
Backward and forward error
Conditioning
Forward error is what you want: \|\hat{\mathbf{x}} - \mathbf{x}\|. Backward error judges the algorithm: the smallest input perturbation for which \hat{\mathbf{x}} is exactly right. The condition number converts one into the other:
\frac{\|\hat{\mathbf{x}} - \mathbf{x}\|}{\|\hat{\mathbf{x}}\|}
\le \kappa(\mathbf{A})\,\varepsilon.
correct digits \approx format digits -\,\log_{10}\kappa(\mathbf{A}). A backward-stable float64 solve carries \approx 16; \kappa = 10^k costs you k.
Hilbert matrices: the rule of thumb, verified
Conditioning
\kappa of the Hilbert matrix H_{ij} = 1/(i{+}j{-}1) grows exponentially. Solving \mathbf{H}\mathbf{x} = \mathbf{b} with \mathbf{x} = \mathbf{1}, the surviving digits track the rule of thumb row by row:
n kappa log10 kappa forward error correct digits backward error
4 1.6e+04 4.2 4.1e-14 13.4 0.0e+00
6 1.5e+07 7.2 1.4e-10 9.8 1.3e-16
8 1.5e+10 10.2 6.1e-08 7.2 5.7e-17
10 1.6e+13 13.2 8.7e-05 4.1 1.0e-16
12 1.8e+16 16.2 3.2e-01 0.5 9.8e-17
The backward error never leaves the 10^{-16} floor; the solver is blameless at every row. The matrix, not the algorithm, amplifies the error.
Normal equations square the pain
Conditioning
Solving least squares via \mathbf{A}^\top\mathbf{A}\,\mathbf{w} =
\mathbf{A}^\top\mathbf{b} replaces \kappa(\mathbf{A}) with its square:
\kappa(\mathbf{A}^\top\mathbf{A}) = \kappa(\mathbf{A})^2.
With \kappa(\mathbf{A}) = 10^5, that is five extra digits lost versus an SVD/QR solve on \mathbf{A} directly:
kappa(A) = 1.0e+05 kappa(A^T A) = 1.0e+10
normal equations: relative error 7.7e-08 (7.1 correct digits)
SVD (lstsq) : relative error 5.5e-14 (13.3 correct digits)
This is why lstsq exists and numerical libraries solve least squares by QR or SVD, never by forming \mathbf{A}^\top\mathbf{A}.
Ridge regularization as preconditioning
Conditioning
Adding \lambda\|\mathbf{w}\|^2 lifts every eigenvalue of \mathbf{A}^\top\mathbf{A} by \lambda, so \kappa = \tfrac{\sigma_1^2 + \lambda}{\sigma_n^2 + \lambda}\downarrow 1: the valley rounds into a bowl, and the same \lambda pays twice, accurate solves and fast gradient descent. Unlike a true preconditioner, ridge changes the minimizer: it shrinks \mathbf{w}_\lambda toward \mathbf{0}.
iterations / kappa: [5.89 5.9 5.94 5.99 6.06 6.15 6.19 6.25 6. ]
Recap
Wrap-up
- Floating point: relative precision \varepsilon_{\text{mach}}, gaps that double, overflow thresholds (e^x dies at x\approx 88.7 in fp32).
- Stable softmax: subtract the max; log-sum-exp is exact; compute cross-entropy from logits as \mathrm{lse}(\mathbf{z}) - z_y.
- Cancellation: reformulate, don’t add bits (
log1p, Welford).
- Conditioning: forward \le \kappa \times backward error; normal equations square \kappa, ridge lowers it.
\kappa sets both solve accuracy and gradient-descent speed, and ridge is the one knob that helps both. The fixes throughout are reformulations rather than extra bits.
