Numerical Methods & Stability

Floating point, conditioning, and the numerical-stability tricks that keep ML code from exploding — taught from zero, starting with why 0.1 + 0.2 isn't 0.3, and built up to the log-sum-exp trick, gradient checking, and mixed-precision training.

intermediate#numerics#stability#floating-point

Start here — what this is really about

It is 3 a.m. Your model trained beautifully for six hours and then the loss printed nan. Nothing in your code changed. The math on the whiteboard is correct. The model is dead anyway.

That is what this page is about. Every idea below exists because of one uncomfortable fact:

A computer does not do arithmetic with real numbers. It does arithmetic with approximations of real numbers, and it rounds after every single operation.

Almost always, that rounding is far too small to notice. Occasionally, it is the difference between a working model and nan. Numerical methods is the craft of telling those two situations apart — and of writing the second kind of code so it behaves like the first.

How to read this page

It starts from first principles: no symbols until you already understand the idea they stand for. Flip the Depth switch at the top for formal notation, error bounds, and edge cases — and the panels open automatically once you've finished the prerequisite (Linear Algebra). Every number on this page was produced by actually running the code, not by hand-waving.

Why this topic pays for itself

This is the least glamorous page in the math track and quite possibly the most useful. Almost nobody debugs a nan by rereading their model architecture — they debug it by knowing the four or five things that produce nan, all of which are on this page. It is also a favourite interview topic precisely because it separates people who have shipped training runs from people who have only read about them.

Part 1 — How a computer actually stores a number

The fixed budget of digits

A computer stores each number in a fixed number of bits — usually 32 or 64. A fixed number of bits means a fixed budget of digits. And a fixed budget of digits means most numbers simply cannot be stored exactly; the machine keeps the closest one it can store and throws away the rest.

is the scheme every computer uses to do this.

Think of it like scientific notation on a very small index card:

You are allowed to write a number only as d.dddddd × 10^k — exactly seven digits, then a power of ten. You can write the mass of the sun and the mass of an electron, because the exponent k moves the scale wherever you need. What you cannot do is write both of them and their difference: with only seven digits of room, adding a gram to the sun changes nothing on your card. The card is not broken. It is simply full.

Precision and range are two different budgets

A floating-point number spends its bits on two separate jobs. The exponent bits buy range — how enormous or how minuscule a number you can name at all. The significand bits buy precision — how many digits of detail you keep once you get there. These trade off against each other, and nearly every numerical disaster in machine learning is one of the two budgets running out. Overflow and underflow are range failures. Cancellation and stalled sums are precision failures.

The famous one: 0.1 + 0.2

Why 0.1 + 0.2 is not 0.3

Computers store numbers in binary, and in binary the fraction 0.1 does not terminate — the same way 1/3 does not terminate in decimal (0.3333…). So the machine stores the nearest number it can represent. Here is what a 64-bit float really holds when you type 0.1:

you type:  0.1
it stores: 0.1000000000000000055511151231257827021181583404541015625

Same for 0.2:

you type:  0.2
it stores: 0.200000000000000011102230246251565404236316680908203125

Now add the two stored values. Their exact sum is a hair above 0.3 — and the nearest storable number to that sum is not the same as the nearest storable number to 0.3 itself. So:

0.1 + 0.2 gives 0.30000000000000004
0.3 on its own gives 0.29999999999999998889776975374843...
Therefore 0.1 + 0.2 == 0.3 is False

Nothing went wrong. Two roundings simply landed on different neighbours.

Run it yourself — the exact stored values are printed, so you can see there is no mystery, only truncation:

Python · runs in your browser

What this does: Prints the exact decimal value your computer actually stores when you write 0.1, 0.2 and 0.3, then shows why their sum misses 0.3 by one tiny step — and why adding 0.1 ten times does not give exactly 1.0 either.

The rule this buys you

Never test floating-point numbers with ==. Test whether they are close: abs(a - b) < 1e-9, or np.allclose(a, b), or torch.allclose(a, b, atol=1e-6). Every equality assertion you have ever seen fail mysteriously in a test suite is this. The one safe exception is comparing against a value you stored rather than computed.

Try to recall

Your test asserts loss == 0.3 and it fails, printing 0.30000000000000004. Is this a bug in your loss function?

Hint: Think about what the machine can and cannot store exactly.

The gaps between numbers are not evenly spaced

Here is the single most important picture on this page. Because a float keeps a fixed number of significant digits, the gap between one storable number and the next grows with the size of the number. Near 1 the gaps are microscopic. Near a billion they are not.

How far apart are the numbers a computer can store?— interactive, drag & zoom

Loading chart…

Measured with numpy.spacing. Near 1.0 a float32 resolves about one ten-millionth; near 100 million the gap has grown to 8, so float32 literally cannot tell 100000000 from 100000004. A float64 near 10 quadrillion has a gap of 2 — it cannot represent odd numbers up there at all. This is why adding a small number to a large one so often does nothing.

That gap has a name: one . And the gap measured right next to 1.0 has its own famous name: .

Read machine epsilon as -digits you get-

Machine epsilon is roughly $10^{-(\text{number of reliable decimal digits})}$ . float32's $\epsilon \approx 1.19\times10^{-7}$ means about 7 decimal digits of trust. float64's $\epsilon \approx 2.22\times10^{-16}$ means about 16. Every claim on the rest of this page is really a claim about how fast you burn through those 7 or 16 digits.

Python · runs in your browser

What this does: Asks numpy to report the real limits of each floating-point format — its precision near 1.0, its largest and smallest values — and then shows the gaps widening as the numbers get bigger, which is the plot above in numbers.

In float32, why does 1e8 + 1 give back exactly 1e8?

Part 2 — The three ways floating point bites

Rounding on its own is harmless. It becomes visible in exactly three ways, and every numerical bug you will ever meet is one of them.

Bite 1 — Rounding that piles up

One rounding costs you the sixteenth digit. A billion roundings can cost you the third.

Think of it like a cashier who rounds every transaction to the nearest cent:

Round one transaction and nobody notices. Round ten million transactions in the same direction and the drawer is meaningfully short. Long-running sums — a loss accumulated over an epoch, a running mean over a dataset — are exactly this.

The worst version is the stalled accumulator: once your running total is large enough that the next increment falls inside one ULP, the sum stops growing entirely, no matter how many more terms you add.

Python · runs in your browser

What this does: Adds 1.0 over and over in float32 and watches the running total freeze at 16777216, because past that point adding 1 no longer reaches the next representable float — a stalled accumulator, the classic way a long summation silently stops working.

Bite 2 — Overflow and underflow (the range budget runs out)

and are the range failures, and the exponential function is where they nearly always start.

exp() runs out of room absurdly early

Your intuition says float32 holds numbers up to $3.4\times10^{38}$ , which sounds like an enormous amount of headroom. But $e^x$ reaches that ceiling at $x \approx 88.7$ . A logit of 90 — a perfectly ordinary number for a confident model — is enough to produce infinity. In float16 the ceiling is 65 504, and $e^x$ hits it at $x \approx 11.1$ .

Where exp(x) runs out of room— interactive, drag & zoom

Loading chart…

The two dashed lines are the largest finite value each format can hold. exp(x) crosses the float16 line at x = 11.09 and the float32 line at x = 88.72 — both verified with numpy. Anything above the line is stored as infinity, and infinity minus infinity, or infinity divided by infinity, is NaN. This single chart explains most NaN losses in classification models.

Underflow is quieter and often nastier: a probability of $10^{-40}$ becomes exactly 0.0, and the very next line takes its logarithm and hands you $-\infty$ . Nothing warned you, because zero is a perfectly legal number.

Try to recall

A model outputs a probability that underflows to 0.0, and the loss is computed as -log(p). What do you see, and what is the standard fix?

Hint: What is the logarithm of zero?

Bite 3 — Catastrophic cancellation (the precision budget vanishes at once)

This is the subtle one, and the one that separates people who know the topic from people who have merely heard of it.

does not create error. It reveals error that was already there.

Think of it like two people guessing the weight of a truck:

Two surveyors estimate a loaded truck at 30 000 kg and 30 001 kg. Each is accurate to about 10 kg — excellent, a relative error of 0.03%. Now subtract to get the weight of the cargo: 1 kg, with an uncertainty of 10 kg. The answer is not merely imprecise, it is meaningless, and it could easily be negative. Both inputs were superb. The subtraction destroyed them.

The quadratic formula gets a root 25% wrong

Solve $x^2 + 10^{8}x + 1 = 0$ in float64, the format with sixteen good digits. The two roots are approximately $-10^{8}$ and $-10^{-8}$ .

Use the formula everyone memorised, $x = \frac{-b + \sqrt{b^2 - 4ac}}{2a}$ , for the small root:

$b^2 - 4ac = 10^{16} - 4$ . But the gap between float64 neighbours near $10^{16}$ is 2, so the $-4$ barely registers.
$\sqrt{b^2 - 4ac} \approx 10^{8}$ , agreeing with $b$ in essentially all sixteen digits.
$-b + \sqrt{\cdot}$ subtracts two numbers that agree in sixteen digits. Every meaningful digit cancels.
The computer returns $-7.45\times10^{-9}$ . The true root is $-1.00\times10^{-8}$ .

That is a 25% error, in the format with sixteen digits of precision. And here is the fix — multiply the formula by $\frac{-b-\sqrt{b^2-4ac}}{-b-\sqrt{b^2-4ac}}$ to get an algebraically identical expression with no dangerous subtraction:

x = \frac{2c}{-b - \sqrt{b^{2} - 4ac}}

$2c$ on top, and $-b - \sqrt{\cdot}$ on the bottom — an addition of two same-signed numbers, since $b$ is positive here. Same-signed addition never cancels.
The two expressions are equal in exact arithmetic — this is algebra, not an approximation.
Evaluated in floating point, this one returns exactly $-1\times10^{-8}$ : every digit correct.
Why it matters: this is the whole discipline in miniature. The same mathematics, written two ways, gives 2 correct digits or 16. Which formula you type is a numerical decision, not a cosmetic one.

The machine-learning version of this trap is the "computational" variance formula, $\mathrm{Var}(x) = E[x^2] - (E[x])^2$ , which is still printed in textbooks:

Python · runs in your browser

What this does: Computes the variance of four numbers two ways — the textbook E[x²]-E[x]² shortcut and the honest two-pass formula — on data that sits far from zero, and shows the shortcut returning a negative variance, which is mathematically impossible.

The pattern to recognise

Cancellation happens whenever you subtract two quantities that are close to each other and large compared to their difference. Once you know the shape, you see it everywhere: $\log(1+x)$ for small $x$ , $e^x - 1$ for small $x$ , $1 - \cos\theta$ for small $\theta$ , numerical derivatives, sample variances, and computing a Euclidean distance as $\sqrt{\lVert a\rVert^2 + \lVert b\rVert^2 - 2a\cdot b}$ — which can go negative under the square root and hand you nan on your own embedding-similarity code.

Every one of those has a purpose-built, cancellation-free replacement. Two you should simply memorise:

Python · runs in your browser

What this does: Shows log1p and expm1 — the built-in cancellation-free versions of log(1+x) and exp(x)-1 — recovering digits that the naive formulas throw away for tiny inputs, including a case where the naive answer is exactly zero and the correct answer is not.

Which of these subtractions is at risk of catastrophic cancellation?

Part 3 — Conditioning vs stability: the problem or the algorithm?

When a computation gives a bad answer there are exactly two suspects, and telling them apart is the central skill of this whole subject.

A bad answer is either the problem-s fault or the method-s fault

Conditioning is a property of the problem: how much the true answer moves when the input wobbles a little. Stability is a property of the algorithm: how much extra error your particular method adds on top. An ill-conditioned problem cannot be rescued by better code — the answer genuinely is that sensitive. An unstable algorithm on a well-conditioned problem is your bug, and it is fixable.

Think of it like weighing a letter on a bathroom scale:

The problem — how heavy is this letter? — is perfectly well-conditioned; the true answer is stable and definite. The algorithm — step on the scale, note it, step on holding the letter, subtract — is catastrophically unstable, because it computes a few grams as the difference of two 80-kilogram readings. Same question, different method, and the method is what fails. Now imagine instead balancing a pencil on its tip and asking where it will fall: no measuring technique on earth saves you, because the problem is ill-conditioned.

A 2x2 system where a 0.0035% nudge flips the answer completely

Solve $A\mathbf{x} = \mathbf{b}$ with two nearly parallel lines:

A = \begin{bmatrix} 1 & 1 \\ 1 & 1.0001 \end{bmatrix}, \qquad \mathbf{b} = \begin{bmatrix} 2 \\ 2.0001 \end{bmatrix}

Solving gives $\mathbf{x} = [1, 1]$ — clean and exact.
Now nudge the last entry of $\mathbf{b}$ from $2.0001$ to $2.0002$ . That is a relative change of about $0.0035\%$ — smaller than a measurement error, far bigger than any rounding error.
Solve again: $\mathbf{x} = [0, 2]$ .

The input moved by 0.0035% and the answer moved by 100%. No algorithm was harmed in the making of this disaster — solving a 2×2 system is about as stable as computation gets. The problem is the villain: two nearly parallel lines have an intersection point that slides enormously when either line tilts by a hair.

The number that measures this sensitivity is the . A problem with a large condition number is called .

The bound is not theoretical hand-waving — you can watch it happen. Below, the notoriously ill-conditioned Hilbert matrices (entry $i,j$ is $1/(i+j+1)$ ), solved for a known answer of all ones, at sizes 2 through 12:

Error grows exactly as fast as the condition number says it will— interactive, drag & zoom

Loading chart…

Hilbert matrices of size 2, 4, 6, 8, 10 and 12, solved in float64 with numpy for a right-hand side whose exact answer is all ones. The measured error (solid) tracks the predicted bound (dashed) across fifteen orders of magnitude. At the far right the condition number is 1.6e16 — larger than 1 over machine epsilon — and the answer is 13% wrong despite a perfectly good solver. The problem, not the code, is the limit.

Try to recall

Your linear solve returns garbage. You rewrite it three different ways and get three different garbage answers. What should you compute before writing a fourth version?

Hint: Is it the algorithm or the problem?

Part 4 — Numerical differentiation, and the U-curve everyone meets

Here is where the two enemies — truncation and rounding — meet head-on in one formula you will actually use: gradient checking, the standard way to verify that a hand-written backward pass agrees with the math.

The idea is the definition of the derivative with the limit removed:

f'(x) \;\approx\; \frac{f(x+h) - f(x)}{h}

$f'(x)$ — the true slope at $x$ , the quantity you are trying to recover.
$h$ — a small step you choose. The whole art is in this one number.
$f(x+h) - f(x)$ — how much the function moved over that step.
dividing by $h$ — rise over run, the ordinary slope you met in Calculus.
$\approx$ instead of $=$ — because we stopped short of the limit $h \to 0$ . That shortfall is truncation error, and it shrinks as $h$ shrinks.
Why it matters: this is how you check an autograd implementation. If your analytic gradient and this estimate disagree by more than about $10^{-6}$ relatively, your backward pass has a bug.

Naively you would make $h$ as small as possible. Try it and something surprising happens:

Gradient checking: the U-curve of choosing h— interactive, drag & zoom

Loading chart…

Estimating the derivative of sin at x = 1 in float64, h shrinking from left to right. Both curves fall as truncation error shrinks, bottom out, then climb again as cancellation in f(x+h) - f(x) takes over. Forward differences are best near h = 1e-8 (about sqrt of machine epsilon); central differences reach roughly a thousand times lower error at h = 1e-5. Making h smaller than the sweet spot actively makes your gradient check worse.

Two errors pulling in opposite directions

Shrinking $h$ reduces truncation error — you are closer to a true limit. But it increases rounding error, because $f(x+h)$ and $f(x)$ become nearly identical numbers and subtracting them is catastrophic cancellation, then you divide the wreckage by a tiny number, magnifying it. Somewhere in between is the best you can do, and it is nowhere near as small as your instincts suggest.

Python · runs in your browser

What this does: Runs a real gradient check on a small function — comparing the analytic derivative against forward and central differences at several step sizes — so you can see the central difference win by a thousandfold and see the error start climbing again when h gets too small.

You are gradient-checking a backward pass. Which setup gives the most trustworthy comparison?

Part 5 — The stability toolkit ML code actually uses

Everything so far has been diagnosis. Here is the treatment: the handful of transformations that keep real training runs alive. Each one rewrites a formula into an algebraically identical form that is numerically far better behaved — exactly the quadratic-formula move, applied to the operations deep learning does constantly.

The log-sum-exp trick (the single most important one)

Softmax turns a vector of scores into probabilities, and it does it with exp — the function you just watched overflow at 88.7.

Softmax with large logits, before and after

Take logits [1000, 1001, 1002]. Nothing about them is pathological — they are just large.

The naive way:

$e^{1000}$ , $e^{1001}$ , $e^{1002}$ — all three overflow to inf (the ceiling is $e^{88.7}$ ).
The denominator is inf + inf + inf = inf.
Each output is inf / inf, which is nan. All three probabilities are nan.

The stable way — subtract the largest logit first:

Largest is $1002$ . Subtract it from all: [-2, -1, 0].
$e^{-2}, e^{-1}, e^{0}$ = [0.135, 0.368, 1.0]. No overflow — the biggest exponent is now exactly zero, so the biggest term is exactly 1.
Divide by their sum: [0.0900, 0.2447, 0.6652].

And here is the crucial part: that is the exactly correct answer. It is bit-identical to the softmax of [0, 1, 2], because softmax only ever depended on the differences between logits. Subtracting a constant from every logit changes nothing mathematically, and everything numerically.

Python · runs in your browser

What this does: Runs softmax on large logits both ways — naive exp-then-normalize, and the max-subtraction version — showing the naive one return all NaN while the stable one returns the exact right probabilities.

Stay in log-space: never compute log(softmax(x))

Even a stable softmax can underflow a small probability to zero, and cross-entropy immediately takes its log.

Python · runs in your browser

What this does: Takes a confident prediction and computes its log-probability two ways — log of the softmax, and log_softmax computed directly — showing the first collapse to negative infinity while the second returns the correct large finite number.

The practical rule, in framework terms

Hand raw logits to the loss function; never hand it probabilities.
• Use torch.nn.CrossEntropyLoss (or F.cross_entropy), not F.softmax followed by torch.log and NLLLoss.
• Use BCEWithLogitsLoss, not Sigmoid followed by BCELoss.
• Use F.log_softmax when you need log-probabilities, not torch.log(F.softmax(x)).
These fused versions do the max-subtraction internally and never materialise a probability that can underflow. This is the single highest-value habit on this page.

The epsilon in every denominator

You have seen a mysterious eps in the signature of every normalization layer and optimizer. Now you know exactly what it is for.

\hat{x} = \frac{x - \mu}{\sqrt{\sigma^{2} + \epsilon}}

$x - \mu$ — subtract the mean: centre the values on zero.
$\sigma^2$ — the variance of the batch or the layer. It can be exactly zero, when every value in the group happens to be identical — a dead channel, a constant feature, a batch of one.
$\epsilon$ — a tiny constant, typically $10^{-5}$ , added inside the square root.
$\sqrt{\sigma^2 + \epsilon}$ in the denominator — with $\epsilon$ present, the divisor can never be zero, so the layer returns a large finite number instead of inf or nan.
Why it matters: without it, one degenerate batch anywhere in training kills the run permanently — a nan propagates into the weights on the very next optimizer step and never leaves. The $\epsilon$ costs nothing when $\sigma^2$ is healthy (it is swamped) and saves the run when it is not. Adam's $\epsilon$ in $\hat{m}/(\sqrt{\hat{v}}+\epsilon)$ does the identical job for a parameter whose gradient has been zero for a while.

Try to recall

Someone sets LayerNorm eps to 0 to be -more mathematically pure- and training runs fine for hours before producing nan. What happened?

Hint: What must have been true of one group of activations at that moment?

Part 6 — Deliberately using less precision

Everything above treats lost precision as a hazard. Modern training treats it as a budget to spend: half-precision numbers are half the memory and several times the throughput on current accelerators. The engineering question is which bits you can afford to lose.

Where each format spends its bits— interactive, drag & zoom

Loading chart…

float16 and bfloat16 are both 16 bits and split them completely differently. bfloat16 keeps all 8 exponent bits of float32 — so it has the same enormous range and a float32 value converts to it without overflowing — but keeps only 7 mantissa bits, roughly 2 to 3 decimal digits. float16 keeps 3 more digits of precision and pays for it with a range that stops at 65504. In training, range turns out to matter far more than precision, which is why bfloat16 won.

Why bfloat16 beat float16 for training

Gradients span an enormous range of magnitudes and cluster near zero. float16 underflows to zero below about $6\times10^{-8}$ and overflows above 65 504, so a raw float16 training run loses small gradients and blows up on large activations. bfloat16 has float32's full range and simply carries fewer digits — and it turns out neural networks tolerate coarse gradients far better than they tolerate missing ones. Noise averages out over a batch; zeros do not.

Loss scaling, the fix that made float16 training work

Suppose a gradient value is $10^{-8}$ — small but real, and cumulatively meaningful over thousands of steps.

Store it in float16: the smallest representable magnitude is about $6\times10^{-8}$ , so it becomes exactly 0. The signal is gone, not approximated.
Now multiply the loss by $S = 1024$ before calling backward. By the chain rule every gradient is scaled by the same $S$ , so this one becomes $1.02\times10^{-5}$ — comfortably representable in float16.
Before the optimizer step, divide the gradients by $S$ in float32. You recover $1.001\times10^{-8}$ — the original value back, to three good digits.

The scale factor cancels exactly in the math and moves the whole gradient distribution up out of the underflow zone in the arithmetic. Frameworks pick $S$ automatically, raising it while things are fine and halving it whenever an inf appears (torch.cuda.amp.GradScaler).

g \;=\; \frac{1}{S}\,\nabla_\theta\!\left(S \cdot L\right)

$L$ — the loss; $\nabla_\theta$ — its gradient with respect to the parameters.
$S$ — the loss-scale factor, a power of two (so multiplying and dividing by it are exact, changing only the exponent and never touching the mantissa).
$S \cdot L$ inside — scale up before backpropagating, so every intermediate gradient in the backward pass is $S$ times larger and clears the underflow floor.
$\frac{1}{S}$ outside — scale back down after, in float32, before the weights are updated.
$=$ — the two sides are exactly equal in real arithmetic, because differentiation is linear. The entire trick is a no-op mathematically and a rescue numerically.
Why it matters: this is the difference between float16 training that diverges and float16 training that matches float32 accuracy. It is also why you must call scaler.step(optimizer) rather than optimizer.step() — the scaler is what unscales, checks for inf, and skips the step if the batch overflowed.

Python · runs in your browser

What this does: Demonstrates loss scaling on a single tiny gradient value — showing it vanish to exactly zero in float16, then survive when multiplied by 1024 first and divided back out afterwards.

Why does mixed-precision training keep a float32 master copy of the weights?

The NaN hunting checklist

When a run dies, work down this list. It is ordered by how often each cause is the culprit.

Your loss just became NaN

Learning rate too high. By far the most common cause and nothing to do with numerics — the parameters diverged, activations exploded, exp overflowed. Drop it 10× and see if the problem disappears before investigating anything else.
A raw exp, log, sqrt, or division you wrote yourself. Search the diff for them. Is there a max-subtraction before the exp? An eps in the denominator? A clamp before the log? Is anything under a sqrt provably non-negative?
A probability handed to a loss instead of logits. log(softmax(x)) or BCELoss(sigmoid(x)) — replace with the fused version.
A zero-variance normalization group, or a division by a count that can be zero (an empty mask, a class absent from the batch, a padded sequence).
float16 overflow. Switch to bfloat16 or check that the grad scaler is enabled and actually stepping.
Bad input data. A nan or inf already present in a feature, propagating in from the dataloader. Assert on it at the boundary: assert torch.isfinite(x).all().

To localise it, enable torch.autograd.set_detect_anomaly(True) — slow, but it points at the exact backward op that first produced the nan — and log grad_norm every step. The step where the norm spikes by orders of magnitude is the step before the death.

Where this shows up

Every item here is one of the ideas above, wearing work clothes:

CrossEntropyLoss takes logits, not probabilities — the log-sum-exp trick, fused.
The eps in LayerNorm, BatchNorm, RMSNorm, Adam, and cosine similarity — a guarantee that a denominator is never zero.
Attention masks use $-10^{9}$ rather than $-\infty$ — avoiding inf - inf and 0 × inf.
GradScaler and bfloat16 defaults — range budgeting under mixed precision.
Ridge's $\lambda I$ and feature standardization — conditioning fixes as much as statistical ones.
Gradient clipping — containment for the numerics you have not fixed.
torch.allclose in every test suite — because == on floats is a bug.
Non-reproducible results across GPUs — non-associative addition under a different reduction order.

Explain it yourself

Explain to a friend, without formulas, why subtracting two nearly equal numbers is dangerous — and then explain the difference between a problem being ill-conditioned and an algorithm being unstable, using the bathroom-scale picture. If you cannot say which of the two a bad answer is, reread Part 3.

Recap — the key ideas

A computer stores a fixed budget of significant digits, so the gaps between storable numbers grow with magnitude; machine epsilon ( $\approx 10^{-7}$ for float32, $\approx 10^{-16}$ for float64) is that gap near 1.
Floating point breaks three familiar rules: 0.1 + 0.2 != 0.3, == is unsafe, and addition is not associative — which is why reduction order changes results.
The three bites are accumulated rounding (stalled sums), overflow/underflow (a range failure, usually via exp), and catastrophic cancellation (subtracting near-equals wipes out every significant digit).
Conditioning is the problem's sensitivity ( $\kappa$ , and you lose about $\log_{10}\kappa$ digits); stability is the algorithm's own added error. Ill-conditioned problems need reformulating, not better code.
Numerical derivatives have a U-curve: too-small $h$ is worse, not better. Use central differences with $h \approx 10^{-5}$ and expect agreement to about $10^{-11}$ .
The toolkit: subtract the max before exp, stay in log-space, put an eps in every denominator, use log1p/expm1, two-pass variance, clip gradients, accumulate in float32.
Mixed precision spends precision on purpose: bfloat16 for range, float32 master weights and accumulation, loss scaling to lift gradients out of the underflow zone.

Practice — and how to make it stick

Learn it the way that actually works

• Retrieval practice: before scrolling up, try to name the three bites and one fix for each. Pulling them from memory beats rereading them.
• Spaced repetition: mark this topic complete and it joins your Review queue, resurfacing right before you would have forgotten it.
• Interleaving: mix these problems with Linear Algebra (condition numbers are singular-value ratios) and Optimization (an ill-conditioned Hessian is why momentum exists) rather than grinding this page alone.

By hand: without running anything, predict what (1e16 + 1) - 1e16 and 1e16 + (1 - 1e16) each return, and explain the difference in one sentence.
Find the cancellation: the expression $\sqrt{x+1} - \sqrt{x}$ loses all its digits for large $x$ . Rewrite it with no subtraction of near-equals (multiply by the conjugate), then check both versions at $x = 10^{12}$ .
Break a solver on purpose: build a Hilbert matrix at $n=12$ , solve for a known answer of all ones, and compare the error to np.linalg.cond times machine epsilon.
Gradient-check something real: implement a small function and its analytic gradient, then sweep $h$ from $10^{-1}$ to $10^{-14}$ with central differences and plot the U-curve yourself.

Start with the lab below — edit and run it right here, and if you get stuck or hit an error, ask Ada on the right: she can see your code and your terminal output.

Practice labPython · runs in your browser · Ada can see your terminal

Your task: Make the unstable code stable. Run it first to see the two failures: a softmax that returns NaN, and a variance formula that returns a negative number. Then do the TODOs — subtract the max before exponentiating, and subtract the mean before squaring — and rerun until both print correct, finite values. Bonus TODO at the bottom: find the smallest logit value that makes the naive softmax fail.

editor

import numpy as np

# ---------- Failure 1: softmax overflow ----------
logits = np.array([1000.0, 1001.0, 1002.0])

with np.errstate(over="ignore", invalid="ignore"):
    naive_softmax = np.exp(logits) / np.exp(logits).sum()
print("naive softmax :", naive_softmax)          # [nan nan nan]

# TODO: fix it. Subtract logits.max() from logits before exponentiating.
#       The answer should be [0.09003057 0.24472847 0.66524096].
stable_softmax = None
print("stable softmax:", stable_softmax)

# ---------- Failure 2: catastrophic cancellation in variance ----------
data = np.array([1e9 + 4.0, 1e9 + 7.0, 1e9 + 13.0, 1e9 + 16.0])
n = len(data)
mean = data.sum() / n

naive_var = (data**2).sum() / n - mean**2
print("\nnaive variance :", naive_var)           # -128.0, which is impossible

# TODO: fix it. Subtract the mean from the data FIRST, then square and average.
#       The answer should be 22.5.
stable_var = None
print("stable variance:", stable_var)

# ---------- Bonus ----------
# TODO: loop over c = 50, 80, 88, 89, 90 and find where np.exp(c) first
#       becomes inf in float32. Compare your answer to np.log(np.finfo(np.float32).max).

terminal

Press Run (⌘/Ctrl+Enter) to execute.

Ask Ada — she can read your terminal

Next: put this to work where it matters most — the numerics of low-precision training in Mixed Precision, and the conditioning story behind optimizers in Optimization.