Numerical Stability: Where Neural Nets Quietly Break

Floating point fails in three ways: numbers get too big, numbers get too small, and subtracting two nearly equal numbers destroys all the precision in the answer. The third is the dangerous one, because it produces a plausible wrong number rather than an infinity, and this page works one that returns exactly zero for a dataset with obvious variance.

Three ways floating point loses

The first two are visible in a log as inf or nan. The third produces a number that looks fine and is not, which is why it gets its own section here and no error message anywhere.

A variance that comes out as zero

Variance has two textbook formulas that are algebraically identical. One of them is a trap. Take three values: x = [10000.0, 10001.0, 10002.0], in fp32.

True population variance:
  mean = 10001
  deviations = [-1, 0, +1]
  variance = (1 + 0 + 1) / 3 = 0.6667

Now the one-pass formula, E[x^2] - E[x]^2:

  10000^2 = 100,000,000
  10001^2 = 100,020,001
  10002^2 = 100,040,004
  sum     = 300,060,005

In fp32 the ULP near 3e8 is 32, because 3e8 sits between
2^28 and 2^29 and fp32 has 24 significand bits:

  2^29 * 2^-24 = 2^5 = 32

  300,060,005 rounds to 300,060,000
  / 3 = 100,020,000        <- E[x^2], already wrong by 1.67

  E[x]^2 = 10001^2 = 100,020,001
  ULP near 1e8 is 8, so this rounds to 100,020,000

  variance = 100,020,000 - 100,020,000 = 0.0

Zero. Not approximately zero, exactly zero, for a dataset whose variance is obviously two-thirds. Every digit of the answer lived in the part that rounding removed, and the subtraction of two hundred-million-scale numbers cancelled everything that was left.

The two-pass version, which is what every library actually uses:

mean = 10001.0                 exact in fp32
d    = [-1.0, 0.0, 1.0]        exact
var  = (1.0 + 0.0 + 1.0) / 3 = 0.66666667

Correct, because nothing large was ever subtracted
from anything else large.

The general rule this illustrates: subtract early, while the numbers are still small relative to their difference, and never late. It is the same principle as subtracting the maximum before exponentiating in softmax — do the dangerous operation where the representation has room.

This matters directly in layer normalisation, which computes a variance over the residual stream at every layer, and where the values are large and close together exactly when the model is behaving oddly.

A sum of ones that stops at 2,048

fp16 represents integers exactly only up to 2^11 = 2048. Above that, the gap between representable values is 2, and adding 1 to 2048 lands exactly halfway between 2048 and 2050. Round-half-to-even picks 2048.

total = fp16(0.0)
for i in range(1_000_000):
    total = total + fp16(1.0)

Result: 2048.0

At total = 2047: 2047 + 1 = 2048, representable, fine.
At total = 2048: ULP is 2. 2049 is not representable.
                 It is midway between 2048 and 2050.
                 Round to even -> 2048.
                 Every subsequent addition does the same.

The loop runs 997,952 more times and changes nothing.

This is why every framework accumulates reductions in fp32 even when the inputs are fp16 or bf16. A mean over a 4,096-element hidden state, a loss summed over a batch, an attention row summed over 32,000 keys — all of these are reductions long enough to stall in a 16-bit accumulator, and the fix is one dtype argument that most people never see because the library set it.

The same effect is worse in fp8: E4M3 represents integers exactly only up to 16. Nobody accumulates in fp8; the format exists for the operands of a matmul whose accumulator is fp32.

The guards inside every framework

None of these change what is being computed. Each one changes where in the number line the computation happens, so that the representation can hold it.

Non-determinism is not instability

Floating-point addition is not associative. Change the order of a sum and the answer can change, and on a GPU the order depends on how the work happened to be scheduled.

In fp64, which has 53 significand bits:

  a =  1e16
  b = -1e16
  c =  1.0

  (a + b) + c = 0 + 1        = 1.0
  a + (b + c) = 1e16 + -1e16 = 0.0

Why: 2^53 = 9.007e15, so the gap between representable
values near 1e16 is 2. Therefore -1e16 + 1 rounds
straight back to -1e16, and the 1 is gone before the
outer addition ever sees it.

Same three numbers. Same operations. Different answer.

A reduction over 4,096 elements on a GPU is split across threads, summed in a tree whose shape depends on the block size, and the block size depends on the batch size, the kernel the library selected, and sometimes on what else was resident on the device. None of that is a bug, and all of it changes the last bits of the result.

Those last bits reach the output. Two logits that differ by 1e-7can swap places, the argmax picks a different token, and from that point the two generations diverge completely — which is the mechanism behind temperature 0 not being reproducible, and why running the same request in a batch of 1 and a batch of 32 can give different text.

The distinction worth holding: this is non-determinism, not instability. The computation is correct to within its precision every time; there is simply more than one correct answer at that precision. Instability is when the error grows without bound, which is what the variance example above shows. The two get confused constantly, and they need opposite responses — non-determinism is addressed by fixing the execution configuration if you need bit-exactness, and instability is addressed by changing the algorithm.

Diagnosing a NaN

1. Find the first NaN, not the tenth. NaN propagates, so by the time the loss is NaN the origin is many operations back. PyTorch’s anomaly detection mode makes the backward pass raise at the operation that produced it, at a substantial speed cost. Turn it on for one run.
2. Check whether it is overflow or a bad input. Log the maximum absolute value of activations per layer for a few steps before the failure. A steady climb toward 65,504 is fp16 overflow. A sudden jump from nowhere is a bad batch or a corrupt sample.
3. Check the denominators. Any division, any log, any sqrt of a computed quantity. A normalisation by a count that can be zero — an empty mask, a filtered batch — is the most common source in data-handling code as opposed to model code.
4. Check the masks. Attention masks are usually implemented by adding a large negative number to masked positions. If that constant is -inf and an entire row is masked, the softmax over that row is exp(-inf - -inf), which is NaN. Use a large finite negative value such as -1e9 in fp32, and note that -1e9 is itself out of range in fp16.
5. Reproduce in fp32. If the same run is stable in fp32, the problem is precision and the fix is a guard. If it is unstable in fp32 too, the problem is the maths and precision was hiding it.

Rules that prevent most of it

• Work in log space for anything multiplicative. Products of probabilities underflow within a few dozen terms. Sums of log probabilities do not. This is why every sequence scorer reports log probabilities and not probabilities.
• Accumulate in a wider type than you multiply in. The hardware already does this for matmuls; do it yourself for any reduction you write by hand.
• Never subtract two large nearly-equal numbers. If an algorithm requires it, there is almost always an algebraically equivalent form that does not — two-pass variance, log1p and expm1 for small arguments, the quadratic formula rearranged for the near-cancelling root.
• Prefer bf16 to fp16 for training, fp16 to bf16 for inference. Training needs range because gradients are small; inference needs precision because activations are well-scaled and there is no accumulation across steps.
• Add epsilon inside the square root, not outside. sqrt(var + eps) is stable; sqrt(var) + eps has an infinite derivative at zero and will produce a NaN gradient even though the forward pass looks fine.