Day 62 · 2026.08.23

Numerical Analysis & the Floating-Point World

There are no real numbers inside a computer — only a substitute whose ruler stretches, and a discipline that keeps its books
"If rounding errors vanished, 95% of numerical analysis would remain." — Nick Trefethen

Floating Point & Machine Epsilon

Scattering points on a logarithmic scale
Representation
Intuition

The real line is continuous; 64 bits hold only about $1.8\times10^{19}$ numbers. The design question is not how many you store but where you put them. Fixed point spreads them evenly — the same millimetre ruler for atoms and for galaxies. Floating point instead spreads them evenly on a logarithmic scale: every binary octave $[2^e,2^{e+1})$ gets the same number of points. The ruler stretches with the value, so absolute precision keeps changing while relative precision stays constant.

That is where $0.1+0.2\neq0.3$ comes from: the denominator of $1/10$ has a factor of 5, so in binary it repeats forever, exactly as $1/3$ never terminates in decimal. This is not insufficient precision — decimal 0.1 simply does not exist in the binary world.

Floating point: same count per octave → spacing doubles with magnitude 1 2 4 8 16 Fixed point: constant spacing → wasted when small, too coarse when large Same 25 points: floating point buys relative precision, fixed point buys absolute
Formal definition
$$\mathrm{fl}(x)=x\,(1+\delta),\qquad |\delta|\le u=2^{-53}\approx1.11\times10^{-16}$$

$\mathrm{fl}$ is the map "round to the nearest representable number", $\delta$ is the relative error it introduces, and $u$ is called the unit roundoff. The whole identity says one thing: rounding replaces your input by a neighbour less than $u$ away. A double is $\pm(1.b_1\cdots b_{52})_2\times2^{e}$ — 53 significant bits (the leading 1 costs no storage) and an 11-bit exponent covering $10^{\pm308}$.

Why it is beautiful

One 64-bit pattern spans 600-odd orders of magnitude, and the relative error is the same everywhere in that range — floating point's uniformity lives in the multiplicative group, not the additive one. The consequence is sharp: multiplication and division are nearly lossless (exponents add, mantissas multiply, each taking one rounding), so addition and subtraction are the dangerous half — floating point is a number system designed for multiplication. The edges are elegant too: subnormals, $\pm0$, $\infty$ and NaN close the arithmetic over every input. Kahan's philosophy was to make the exception a number rather than a crash.

Applications

Every precision decision in AI hardware is a reading of this table. fp16 is 5 exponent bits plus 10 mantissa bits; bf16 keeps fp32's full 8-bit exponent and cuts the mantissa to 7 — training fears gradients overflowing to $\infty$, not a couple of lost decimals. Same 16 bits, a different cut, a different failure mode. The counterexample is finance, which uses decimal fixed point because 0.1 of a unit of currency must be exactly one tenth.

Essence + question
Floating point is not an approximation to the reals; it is its own number system whose ruler stretches with the value — it promises relative precision and never absolute.
Question: if your quantities span 20 orders of magnitude in one computation, did floating point help you, or just hide the problem?

The Condition Number

A sickness of the problem, not a fault of the algorithm
Conditioning
Intuition

Split "we can't compute it accurately" into two independent questions: is the problem itself sensitive to perturbed input? and does the algorithm manufacture extra error? The condition number measures only the first, and is indifferent to your language, hardware and algorithm.

Solving $Ax=b$ is finding where two lines cross. When they meet at a healthy angle, nudging a line nudges the crossing; when they are nearly parallel, the same nudge throws the crossing far away — $\kappa(A)$ measures exactly "how parallel". $\kappa=10^{8}$ means double precision's 16 significant digits go in and 8 come out; those 8 were not lost by the algorithm, the problem never carried them through in the first place.

Two nearly parallel lines → a wildly unstable crossing the crossing runs away Dashed = the second line nudged slightly (a tiny data perturbation) κ(A) = worst-case ratio of relative output change to relative input change
Formal definition
$$\frac{\|\delta x\|}{\|x\|}\;\le\;\kappa(A)\,\frac{\|\delta b\|}{\|b\|},\qquad \kappa(A)=\|A\|\,\|A^{-1}\|=\frac{\sigma_{\max}}{\sigma_{\min}}$$

The left side is the relative change in the solution, the right the relative change in the data, and $\kappa$ is the worst-case amplification between them. The $\sigma$ are singular values: $A$ stretches the unit ball into an ellipsoid whose longest and shortest axes are $\sigma_{\max}$ and $\sigma_{\min}$, so the condition number is how flat that ellipsoid is. Rule of thumb: one solve costs you about $\log_{10}\kappa$ significant digits.

Why it is beautiful

It promotes "numerically hard" from an engineering complaint to an intrinsic invariant of the problem, and recognises four unrelated-looking phenomena as one: ill-conditioned linear systems; polynomial roots that are hypersensitive to coefficients (Wilkinson's degree-20 polynomial, where perturbing one coefficient by $2^{-23}$ sends integer roots off into the complex plane); the Lyapunov exponent of chaos; and exploding or vanishing gradients in deep networks, which are the layerwise Jacobian condition numbers multiplied together. One quantity, four disciplines, four different names.

Applications

Once you know $\kappa$, a whole family of "numerical tricks" reveals a single motive: shrink it. Ridge regression's $\lambda I$ lifts $\sigma_{\min}$; conjugate gradient's iteration count scales as $\sqrt{\kappa}$, so a preconditioner $M^{-1}A$ exists purely to swap in an equivalent problem with a smaller condition number; spectral normalisation and orthogonal initialisation push the Jacobian's spectrum toward 1; and Adam is a diagonal preconditioner — scaling gradients coordinatewise is exactly rounding out that ellipsoid.

Essence + question
The condition number is the information decay rate built into the problem: it fixes how many digits you can possibly get right before any algorithm gets a say.
Question: last time you lowered the learning rate until training "finally converged" — did you fix the optimisation, or merely sidestep a parametrisation with a huge $\kappa$?

Catastrophic Cancellation

Subtraction does not create error; it exposes it
Cancellation
Intuition

Subtract two nearby numbers: $1.2345678-1.2345600=0.0000078$. Each input carried 8 significant digits; the result has 2 — and the subtraction got every digit right. All it did was promote the uncertainty hiding in the last place to the first. The noise was always there, screened by the leading digits. Relative error merely accumulates through multiplication and division; only subtracting nearby numbers can amplify it without bound. This is the one genuinely catastrophic step in the floating-point world.

The textbook quadratic formula walks straight into it: when $b>0$ and $4ac\ll b^2$, the numerator of $x=\frac{-b+\sqrt{b^2-4ac}}{2a}$ becomes a difference of two nearly equal positive numbers. The fix is one line: compute the stable root $x_1=\frac{-b-\sqrt{b^2-4ac}}{2a}$ first, then recover the other from Vieta's relation, $x_2=c/(ax_1)$. One identity, two computational paths, worlds apart numerically.

a b a−b first 14 digits cancel to zero 2 digits left, one of them already noise red = rounding noise blue = trustworthy digits the noise did not grow — it was promoted from last place to first
Formal definition
$$\frac{\big|(\hat a-\hat b)-(a-b)\big|}{|a-b|}\;\le\;\frac{|a|+|b|}{|a-b|}\cdot u$$

Here $\hat a,\hat b$ are the stored values, already carrying rounding error, and $u$ is machine precision. The factor $\frac{|a|+|b|}{|a-b|}$ on the right is the amplification, and it blows up as $a\to b$. Note that it is precisely the condition number of subtraction — cancellation is not a new phenomenon, it is conditioning showing its face in the simplest possible operation.

Why it is beautiful

It punctures an assumption we take for granted: an identity that holds over the reals is not the same program in floating point. $a^2-b^2$ and $(a-b)(a+b)$ are one number with two fates — abstract equivalence torn open by finite word length, and the whole of numerical analysis grew out of that tear. The repair is prettier still: Kahan summation carries a compensation variable that recovers the low-order bits each addition discards, pulling the error bound from $O(nu)$ down to something independent of $n$. Those "lost" bits never left the register; you only have to stoop and pick them up.

Applications

The one-pass variance formula $E[X^2]-(E[X])^2$ returns negative numbers when the mean dwarfs the standard deviation, which is why Welford's online algorithm exists; expm1 and log1p exist for the single reason that $e^x-1$ and $\log(1+x)$ cancel away every significant bit for tiny $x$; and on GPUs the summation order inside an all-reduce is not fixed, which is a leading reason large-model training is hard to reproduce bit for bit.

Essence + question
Subtracting nearby numbers creates no error; it only removes the cover — your significant digits were eaten somewhere else long before.
Question: is there a place in your code computing $A-B$ where $A$ and $B$ are each long accumulated sums?

Backward Error & Low Precision

Pushing the error back to the input
Stability
Intuition

In the 1960s Wilkinson inverted the question. The traditional one — how far is my answer from the true answer (the forward error) — is hard to answer, because it blends the problem's sickness with the algorithm's mistakes. Wilkinson asked instead: "the answer I computed — which problem is it the exact solution of?"

If $\hat x$ exactly satisfies $(A+\Delta A)\hat x=b$ with $\|\Delta A\|/\|A\|$ only $10^{-15}$, the algorithm is beyond reproach: it handed you the exact solution of a neighbouring problem, and $A$ came from measurement anyway, already carrying far larger uncertainty. A good algorithm is thus defined as: one that returns the exact solution of a nearby problem.

data A, b true solution x nearby data A+ΔA computed x̂ exact solve exact solve backward error η forward error floating-point algorithm forward error ≤ condition number × backward error Wilkinson diagram: the computed answer exactly solves a nearby problem
Formal definition
$$\underbrace{\frac{\|\hat x-x\|}{\|x\|}}_{\text{forward error}}\;\lesssim\;\underbrace{\kappa(A)}_{\text{problem's sickness}}\;\times\;\underbrace{\eta}_{\text{backward error}}$$

This is the fundamental equation of numerical analysis, and it splits the blame cleanly in two: $\kappa$ cannot be repaired (only the formulation of the problem can change), while $\eta$ can — and that is where algorithm design happens. Gaussian elimination with partial pivoting satisfies $\eta\le c_n\rho u$, where $\rho$ is the growth factor of the entries during elimination, and pivoting exists to hold it down: the worst-case bound is still $2^{n-1}$, yet in practice $\rho$ sits near $O(n)$ — a gap that still has no satisfying explanation.

Why it is beautiful

The framework adapts to precision by itself: swap $u$ from $2^{-53}$ to bf16's $2^{-8}$ and not a symbol of the inequality changes. So "how many bits should I use" turns from folklore into arithmetic — and mixed-precision iterative refinement is the answer: do the $O(n^3)$ factorisation in low precision, compute the $O(n^2)$ residual in high precision and correct, and a few rounds buy back double-precision accuracy. Cheap compute produces the answer; expensive compute audits it.

Applications

Every LAPACK and cuSOLVER routine documents its backward stability — not academic courtesy but an interface contract. Mixed-precision training keeps the same books: fp16 forward and backward, an fp32 master copy of the weights, and loss scaling to lift gradients clear of fp16's subnormal range. Quantising an LLM to int8 or int4 changes the question in the same way: not "how much accuracy did we lose" but "the quantised model is the exact execution of which slightly different set of weights" — so what matters is the direction of the perturbation, and GPTQ and AWQ are essentially about steering that error into subspaces the activations do not notice.

Essence + question
Rather than asking how wrong the answer is, ask which question it answers exactly — pushing the error back to the input is what separates a sick problem from a faulty algorithm.
Question: when quantisation costs a model only 0.3% accuracy, is that a good quantiser, or a model that was well conditioned to begin with?

Going Deeper

Floating-point error is entirely deterministic — so why does the same code give different results on a GPU twice?
The cause is not randomness but the fact that floating-point addition is not associative: $(a+b)+c$ and $a+(b+c)$ can differ in the last place. "Summation" is therefore not a well-defined operation in floating point — it depends on a particular reduction tree, and multithreaded reductions, atomicAdd, and size-dependent kernel selection all change that tree. Bit-exact reproducibility means fixing the reduction order, and you pay for it in throughput. The deeper point: we assume that mathematically equal means the same computation, whereas in floating point the algorithm is the semantics.
Low-precision training works — does that mean neural networks are inherently well conditioned?
More precisely: the gradient noise SGD already carries is far larger than quantisation noise — a gradient is a minibatch estimate to begin with, and one more layer of rounding perturbation leaves its statistics intact. But conditioning is layered. The forward pass is well conditioned; long accumulation is not. Optimiser state and weight updates are the classic "large value absorbs a small increment" cancellation scenario, which is exactly why an fp32 master copy of the weights is mandatory.
Could interval arithmetic or exact rational arithmetic abolish floating-point error once and for all?
They give rigorous conclusions at a cost that is usually unacceptable. Interval arithmetic guarantees the true value lies inside the interval, but the width grows exponentially with iteration count and often becomes uninformative; exact rationals blow up the digit counts of numerator and denominator. The real answer is this: error is not to be abolished, it is to be accounted for. Accepting it and proving a bound is far more tractable than eliminating it.
The condition number says you lose $\log_{10}\kappa$ digits — can "information" in computation be measured like Shannon entropy?
It is more than an analogy. Significant digits are bits, and $\log_2\kappa$ is the bit loss of a single solve: an ill-conditioned linear system is a lossy channel and $\kappa$ is its attenuation (echoing Day 13). The mechanism differs — a channel adds random noise, while here the collapse is deterministic: $A$ compresses certain directions down toward $\sigma_{\min}$, and information along them is simply truncated at finite word length. A practical criterion follows: if $\kappa>10^{16}$, more precision only postpones the loss by one digit, and the right move is to reformulate.