Numerical stability

The last lesson left a queasy feeling: arithmetic can be flawless and the conclusion still wrong. There the culprit was a hidden variable. Here it is a hidden limitation, and it sits one level deeper — inside the machine itself.

Every derivation in this course assumed real numbers with infinite precision. The math your code actually runs does not — it uses floating point, a finite approximation with a finite range and about seven honest digits. Most of the time the gap is invisible. But ML lives exactly where it bites: exp of a large logit is infinity, log of a probability that rounded to zero is minus-infinity, and the product of a thousand small probabilities underflows to zero. A naive softmax or cross-entropy then returns NaN and poisons the whole training run. The reassuring news: the fixes are a small, standard toolkit — and once you have seen them you will spot the subtract-the-max, the log-space, and the +1e-9 everywhere in real code.

Two ways floating point fails

A float has finite range and finite precision:

• Out of range. Above roughly 3.4 × 10^38 (for 32-bit float) a value overflows to inf; below about 1.2 × 10^-38 it underflows to 0. Both are unrecoverable — once you have inf or 0, the information is gone.
• Lost precision. A 32-bit float keeps only ~7 significant digits (16-bit far fewer). Subtracting two nearly-equal numbers can wipe out every meaningful digit — catastrophic cancellation.

ML’s favorite operations — exp, log, products of probabilities — live right at these edges.

Trap 1: overflow in softmax — subtract the max

Softmax exponentiates its inputs. If the logits are large, exp overflows and the whole thing collapses to NaN. The fix uses a free identity: softmax is shift-invariant, softmax(x) = softmax(x − c) for any constant c. Pick c = max(x) so the largest exponent is exp(0) = 1 and nothing overflows:

import numpy as np
x = np.array([1000.0, 1001.0, 1002.0])

naive = np.exp(x) / np.exp(x).sum()           # exp(1000) overflows to inf
z = x - x.max()                                # shift so the largest is 0
stable = np.exp(z) / np.exp(z).sum()

print("naive softmax :", naive)
print("stable softmax:", np.round(stable, 4))

naive softmax : [nan nan nan]
stable softmax: [0.09   0.2447 0.6652]

Same math, but the naive version returns NaN (it computed inf / inf) while the shifted version is exact. Every framework’s softmax does the subtraction for you.

Trap 2: underflow — why ML lives in log space

Now the opposite end. A sequence’s probability is a product of many per-token probabilities, each less than 1. Multiply a thousand of them and the result underflows to 0 — and log(0) = -inf, so your loss explodes. The fix is to never form the product: work in log space, where a product becomes a sum.

log( p_1 · p_2 · ... · p_n )  =  log p_1 + log p_2 + ... + log p_n

A sum of a thousand log-probabilities (each a modest negative number) is perfectly stable, while their product is not. This is the deep reason ML uses log-likelihood, cross-entropy, and log-probabilities everywhere — not mathematical taste, numerical survival.

Trap 3: the log-sum-exp trick

Sometimes you need log( Σ exp(x_i) ) — the normalizer of a softmax, the denominator of a log-likelihood. Computing the exp directly overflows. The log-sum-exp identity factors the max back out, exactly like the softmax shift:

log Σ exp(x_i)  =  m + log Σ exp(x_i − m)        where m = max(x_i)

import numpy as np
x = np.array([1000.0, 1001.0, 1002.0])
m = x.max()
print("naive log-sum-exp :", np.log(np.exp(x).sum()))     # log(inf) = inf
print("stable log-sum-exp:", round(float(m + np.log(np.exp(x - m).sum())), 4))

naive log-sum-exp : inf
stable log-sum-exp: 1002.4076

The naive form overflows to inf; the stable form gives the right answer (1002.41). scipy.special.logsumexp and torch.logsumexp are exactly this.

In one breath

• Floating point has finite range (overflow → inf, underflow → 0) and finite precision (~7 digits in fp32); ML’s exp/log/probability-products sit right at these edges.
• Softmax overflow: large logits make exp blow up; subtract the max first — softmax is shift-invariant, so the result is unchanged and nothing overflows.
• Underflow: products of many small probabilities round to 0, so ML works in log space (a product becomes a sum) — the real reason log-likelihood/cross-entropy are everywhere.
• The log-sum-exp trick, log Σ exp(x) = m + log Σ exp(x − m), computes sums-of-exps safely; watch for catastrophic cancellation when subtracting nearly-equal numbers.
• Use the framework’s stable ops (log_softmax, cross_entropy, logsumexp); never hand-roll the naive form — and it’s worse in fp16/bf16.

Practice

A question to carry forward

Notice the move that rescued every trap in this lesson. Underflow vanished the moment we stopped multiplying probabilities and started adding their logs; overflow vanished the moment we factored the max out and back in. In each case we never touched the answer — we changed the representation, did the work where it was easy, then changed back. log → add → exp is that round trip in miniature.

Hold that shape, because the final stretch of this course is built on it at a far grander scale. There is a change of representation so powerful that it turns one of the most expensive operations in computing — convolution — into a plain multiplication, and it is the front door to all of audio ML. What is the Fourier transform, why is any signal secretly a sum of sine waves, and how does moving from the time domain to the frequency domain, doing the easy thing, and moving back become a tool hiding inside spectrograms, rotary positional embeddings, and fast long-convolution models?