Math

Numerical Stability

Floating point is not arithmetic — it is arithmetic with a budget of about sixteen digits, and the failures come from spending that budget carelessly rather than from any bug in your logic.

Shrinking a step size reduces the error of the approximation and increases the error from rounding, so total error is U-shaped and there is a best step size beyond which smaller is worse

Why Does This Exist?

Your loss becomes NaN at step 3,000 and the model is unrecoverable. Nothing in your code changed. There is no exception, no stack trace pointing at a line, and the run before this one was fine.

That is the failure this page prevents, and its causes are a short list: a value exceeded its format's ceiling, a subtraction destroyed every meaningful digit, or a logarithm met a zero. Each has a standard fix that is well known and easy to omit.

Three things you already use exist purely because of this:

log_softmax rather than log(softmax(x)) — the log-sum-exp trick, and the reason frameworks fuse softmax into cross-entropy.
The epsilon in every normalisation layer — (x - μ) / sqrt(σ² + ε), guarding a division.
Loss scaling in mixed precision — multiplying the loss so small gradients survive float16.

None of these are arbitrary constants. They are each answers to a specific arithmetic failure.

Think of It Like This

A calculator with a sixteen-digit window

Imagine every number you compute must fit in a window showing sixteen significant digits, plus an exponent. Large and small numbers are fine — the exponent handles those. What you cannot do is hold sixteen digits of a big number and details far below them.

Now subtract two nearly equal numbers: 1.2345678901234567 minus 1.2345678901234000. Both filled the window. Their difference is about $5.67 \times 10^{-14}$ , and only the last few digits contributed — the leading fifteen cancelled and carried no information into the answer. You started with two sixteen-digit numbers and ended with two or three meaningful digits.

Nothing overflowed and nothing warned you. You simply spent your precision on digits that then cancelled, and the survivors are mostly rounding noise. That is catastrophic cancellation, and it is the failure that produces wrong answers rather than obvious errors.

The other failure is simpler: the exponent has limits too. Push past the largest representable number and you get infinity, and infinity minus infinity is not a number.

How It Actually Works

Three failure modes, one fix each.

Overflow, and the log-sum-exp trick

float32 tops out near $3.4 \times 10^{38}$ , reached by $e^{x}$ at only $x \approx 88$ . float16 tops out at 65,504, reached at $x \approx 11$ . Logits routinely exceed both.

Once a value is inf, the next operation usually produces nan: inf/inf, inf - inf, 0 × inf. And nan propagates through everything it touches, so one bad value poisons every parameter on the backward pass.

The fix rests on an exact identity. For any constant $m$ :

\log \sum_i e^{x_i} = m + \log \sum_i e^{x_i - m}

Choose $m = \max_i x_i$ . Now the largest exponent is $e^0 = 1$ , so nothing can overflow — and because the identity is exact, this is not an approximation. Softmax gets the same treatment: subtracting the maximum from every logit leaves the output unchanged, because the constant cancels between numerator and denominator.

This is why log_softmax is a single fused operation rather than a composition. Computing log(softmax(x)) forms probabilities first, and a small probability underflows to exactly 0, after which log(0) is $-\infty$ . Folding the logarithm in algebraically never forms the probability at all.

Catastrophic cancellation

Subtracting nearly equal numbers discards their agreeing leading digits and promotes rounding noise into the result. Three places this bites in ML:

Variance as $E[X^2] - (E[X])^2$ . With large means both terms agree in their leading digits, so the difference is noise — and can come out negative, which is impossible for an average of squares.
Numerical gradient checks. $f(x+h) - f(x-h)$ is a subtraction of nearly equal values, which is why shrinking $h$ makes the check worse.
Any "compute then subtract the mean" pattern on uncentred data.

The fixes are structural rather than parametric. Welford's algorithm computes variance without ever forming large squares. np.linalg.norm factors out the largest component before squaring. log1p(x) computes $\log(1+x)$ accurately for tiny $x$ , and expm1(x) does the same for $e^x - 1$ — both exist precisely to avoid a cancelling addition.

The U-shaped error curve

The diagram shows the pattern that governs every step-size choice. Truncation error — how wrong your approximation is — shrinks as the step shrinks. Round-off error — how much precision the arithmetic destroys — grows as the step shrinks, because you are subtracting ever-closer values and dividing by an ever-smaller number.

Total error is therefore U-shaped, with a minimum. For a central-difference derivative in float64 that minimum sits near $h \approx 10^{-5}$ . Smaller is not better; it is worse, and this surprises people every time.

Where epsilon goes

(x - μ) / sqrt(σ² + ε) puts the epsilon inside the square root, added to the variance. Writing sqrt(σ²) + ε instead is subtly different and worse: it fails to protect the square root's derivative, which is $1/(2\sqrt{\sigma^2})$ and blows up as the variance approaches zero. The gradient becomes enormous even though the forward value looked fine.

Typical values are $10^{-5}$ for float32 and $10^{-3}$ for float16, and they are dtype-dependent for a reason — an epsilon smaller than the format's own resolution does nothing at all.

Worked example

Softmax of the logits $[1000, 1001, 1002]$ .

Naive. $e^{1000}$ overflows even float64, whose ceiling is about $1.8 \times 10^{308}$ . So the numerator is inf, the denominator is inf, and the result is nan.

Stable. Subtract the maximum, 1002, giving shifted logits $[-2, -1, 0]$ :

e^{-2} = 0.135335,\quad e^{-1} = 0.367879,\quad e^{0} = 1

Their sum is $1.503214$ , so:

\text{softmax} = \frac{[0.135335,\; 0.367879,\; 1]}{1.503214} = [0.09003,\; 0.24473,\; 0.66524]

These sum to 1 ✓, and every intermediate value stayed between 0 and 1.5.

The log-sum-exp of the original logits follows from the identity:

\log\sum_i e^{x_i} = 1002 + \log(1.503214) = 1002 + 0.407606 = 1002.407606

A perfectly ordinary number, computed without any intermediate exceeding 1.5. And note the shift changed nothing: the softmax of $[1000,1001,1002]$ is identical to the softmax of $[0,1,2]$ , because only the differences matter.

Show Me the Code

1import numpy as np2
3
4def softmax_stable(x: np.ndarray) -> np.ndarray:5    shifted = x - x.max()                  # largest becomes 0, so exp() is safe6    e = np.exp(shifted)7    return e / e.sum()8
9
10logits = np.array([1000.0, 1001.0, 1002.0])11print(np.round(softmax_stable(logits), 5))            # -> [0.09003 0.24473 0.66524]12print(round(float(logits.max() + np.log(np.exp(logits - logits.max()).sum())), 6))13# -> 1002.407606   the log-sum-exp14
15with np.errstate(over="ignore", invalid="ignore"):16    print(np.exp(logits) / np.exp(logits).sum())       # -> [nan nan nan]17
18# Cancellation: the one-pass variance formula on a large mean.19big = np.array([1e9 + 4, 1e9 + 7, 1e9 + 13])20print((big**2).mean() - big.mean() ** 2, round(float(big.var()), 4))  # -> wrong, 13.555621
22# log1p keeps precision where log(1+x) cannot.23x = 1e-1624print(np.log(1 + x), np.log1p(x))                     # -> 0.0 1e-16

The stable softmax and the 1002.4076 match the hand calculation, and the naive path returns nan from the same input. The last two blocks are cancellation and its two standard remedies.

Watch Out For

Hunting a NaN at the line where it surfaced

nan propagates, so by the time your loss prints nan the cause is usually many operations upstream and possibly several steps earlier. Debugging at the point of discovery is the slow path.

The efficient approach is to make the failure loud at its origin. torch.autograd.set_detect_anomaly(True) runs the backward pass with checks and raises at the operation that first produced a non-finite gradient, including the forward traceback. It is slow, so enable it only while debugging — but it converts a mystery into a line number.

Failing that, bisect with assertions: assert torch.isfinite(x).all() after each block narrows the origin quickly.

The usual culprits, roughly in order: exponentiating an unbounded logit; log of something that reached zero — a probability, a variance, a count; dividing by a norm or standard deviation that collapsed; sqrt of a value that went slightly negative through rounding; and a float16 gradient that overflowed in a long chain-rule product.

One more, easy to miss: nan fails every comparison, including nan != nan. So if loss > threshold is False for nan and a guard written that way lets it straight through. Use torch.isnan or math.isfinite explicitly.

Adding epsilon in the wrong place, or at the wrong scale

An epsilon is a specific guard against a specific division or logarithm, and moving it changes what it protects.

sqrt(var + eps) protects both the division and the square root's derivative. sqrt(var) + eps protects only the division — the gradient of sqrt at a near-zero variance is $1/(2\sqrt{\text{var}})$ , which explodes, so you get a finite forward pass and an enormous gradient. That is a real and confusing bug, and the two expressions look interchangeable.

Scale matters too, and it is dtype-dependent. float16's machine epsilon is about $10^{-3}$ , so an epsilon of $10^{-8}$ rounds away to nothing and protects nothing at all. Mixed-precision training needs larger epsilons than float32 — which is why frameworks change the default rather than exposing one constant.

Over-large epsilons cause the opposite problem: an epsilon comparable to your actual variance biases the normalisation, so activations are systematically under-scaled and the layer partly stops doing its job.

Two habits. Clamp rather than add where the quantity is genuinely bounded — probs.clamp(min=1e-7) before a log states the intent more clearly than adding. And prefer the fused operation whenever one exists: log_softmax, logsumexp, log1p, cross_entropy on logits. Each was written by someone who worked out where the epsilon goes, and using it means you do not have to.

The Quick Version

Floating point carries about 16 significant digits in float64 and 7 in float32. Both the range and the precision run out.
exp overflows float32 at a logit of ~88 and float16 at ~11. Then inf becomes nan, and nan spreads everywhere.
Log-sum-exp: subtract the maximum before exponentiating. Exact, not an approximation.
Softmax is unchanged by subtracting a constant from every logit, which is what makes the trick free.
Use log_softmax and cross_entropy on raw logits. log(softmax(x)) underflows to log(0).
Catastrophic cancellation: subtracting nearly equal numbers destroys the leading digits. Source of negative variances and failed gradient checks.
Use Welford for variance, np.linalg.norm for norms, log1p and expm1 for near-zero arguments.
Total error is U-shaped in step size. For float64 central differences the best $h$ is about $10^{-5}$ ; smaller is worse.
Epsilon goes inside the square root: sqrt(var + eps). Outside leaves the derivative unprotected.
Epsilon scale is dtype-dependent — $10^{-8}$ does nothing in float16.
Debug nan with set_detect_anomaly(True), and remember nan fails every comparison silently.