AI Engineering Foundations 13: Numerical Stability - Floating Point Is a Leaky Abstraction

Forum topic · 小凯 · 2026-10-09

Summary

This is a detailed Chinese-language review of Lesson 13 (Numerical Stability) from the open-source course 'AI Engineering from Scratch' (58k GitHub stars, MIT license). The post explains how silent training failures happen: logits becoming inf, gradients becoming NaN, and precision loss from float16 training without loss scaling. It covers IEEE 754 representation (float32's 23-bit mantissa vs float16's 10 bits), the three failure modes (overflow, underflow, catastrophic cancellation - with a 19.2% relative error example), and the max-subtraction trick that makes log-sum-exp and softmax numerically safe. It also audits the course's own demos: the Python demo code runs in float64 so the promised NaN/inf failures never occur unless inputs are cast to float32, and several exercises have flawed verification targets (machine epsilon, float16 zeroing ratios). Highlights include Exercise 1, where the naive variance formula E[x^2]-E[x]^2 produces -65536.0 in float32 while Welford's algorithm yields 0.6667 accurately, and practical guidance on bfloat16 for training, loss scaling, gradient clipping, and normalization layers as built-in numerical stabilizers.

This article is part of a Chinese series reviewing the open-source course AI Engineering from Scratch (58k GitHub stars, MIT license). Lesson 13 covers numerical stability.

Key points

How training silently dies

Three realistic scenarios: (1) logits become inf at step 9,001 and all gradients become NaN by step 9,002 - no error, just a NaN loss; (2) accuracy is 2% below the paper because you trained in float16 without gradient scaling while the paper used float32; (3) a hand-written cross-entropy returns inf once logits exceed 100, because exp(100) overflows.
Floating point is a "leaky abstraction": you think you're computing with real numbers, but you're hopping on a sparse grid.

Arithmetic on a grid

IEEE 754 floats have sign, exponent, and mantissa. float32: ~7 decimal digits of precision; float16: ~3 digits, max value 65,504 - dangerously small for ML.
0.1 + 0.2 != 0.3 because 0.1 is non-terminating in binary. Never compare floats with ==; use abs(a-b) < epsilon.
Three ways a grid computation dies: overflow (becomes inf, exp() is the usual culprit), underflow (small numbers silently become zero, killing gradients), and catastrophic cancellation. Example: 1.0000001 - 1.0000000 in float32 gives 1.1920929e-07 instead of 1e-7 - 19.2% relative error in one subtraction.

Two lines of code that save training

The core identity: log(sum(exp(x))) = max(x) + log(sum(exp(x) - max(x))). After subtracting the max, the largest exponent term is exp(0) = 1 (no overflow) and the sum is at least 1 (no log(0)). Mathematically identical; numerically, it turns a minefield into a highway.
bfloat16 keeps float32's 8-bit exponent (range up to 3.4e38) while sacrificing precision: training needs range, inference needs precision, so train in bf16, often infer in float16.
Pure float16 training hits gradient underflow; the fix is loss scaling (multiply loss by 1024, divide gradients back before the weight update). Gradient clipping by norm (transformer convention: 1.0) rescales without changing direction. Batch norm and layer norm also act as numerical stabilizers by keeping activations in a safe range.

The reviewer's audit of the course's own demos

The Build It demos run in Python's math module (float64), so the promised failures never occur: e^100 = 2.7e43 and e^502 = 1.0e218 fit comfortably in float64's 1.8e308 limit. Naive and stable implementations match to six decimals. To see the documented nan/inf, you must cast inputs to np.float32 yourself - the chosen "dangerous" inputs (100, 500) are safe for float64 and only fatal for float32.
The softmax section contains leftover author reasoning in the prose, including an inequality written backwards (2.69e43 < 3.4e38); the final conclusion is still correct.
Exercise 2 is ambiguous: "the smallest float32 x with 1.0 + x == 1.0" literally yields 1.4e-45 (the grid floor); interpreted as "smallest x where 1+x != 1", you get 5.9604652e-08 - still only half of eps (1.1920929e-07) because 1 + 2^-24 ties between grid points and round-to-even swallows it. Three readings, three numbers, none equal to eps.
Exercise 5 (fraction of gradients that become zero in float16) depends on an unstated sampling distribution: uniform over [1e-9, 1e-3] gives 0.004% zeros; log-uniform (closer to real gradients) gives 24.5%, rescued by 1024x loss scaling.
Minor errata: exp(88.7) measured as 3.33e38 (not 3.40e38); the 0.1+0.2 example's sum wasn't itself rounded to float32 (0.300000011920929 vs the quoted exact sum); an opening anecdote misattributes the accuracy loss to "32-bit accumulated rounding error" when it was the float16-vs-float32 gap; bfloat16 simulation uses truncation while modern frameworks use round-to-nearest.

What the course gets right

Exercise 1 is among the series' best: computing the variance of [1000000, 1000001, 1000002], the naive formula E[x^2] - E[x]^2 returns -65536.0 in float32 (negative! true value 0.6667), while Welford's algorithm yields 0.6666666865 (3e-8 relative error).
The log-sum-exp proof is airtight; the catastrophic cancellation figure is accurate to the third decimal; hard boundaries (65504, 5.96e-8, format tables) all check out; and the bf16-vs-float16 tradeoff is presented more even-handedly than most materials.

Takeaway

Three numbers worth running yourself: 1.4e-45 (grid floor), 5.96e-08 (tie rounding), 1.19e-07 (finfo's eps). And keep the Exercise 1 screenshot of -65536.0 - once you've seen variance come out negative, you never forget that floating point promises real numbers but delivers points on a grid.