Numerical Stability and Finite Precision

The symptoms differ, the repairs differ, and only one of them announces itself. Overflow and NaN stop a run. Accumulated rounding lets the run continue and drifts instead. The drift has a schedule, and you can read it off a table.

A US government report printed that table in 1992. On 25 February 1991 a Patriot battery at Dhahran failed to intercept a Scud, and 28 Americans died. The General Accounting Office traced the failure to a conversion. The weapons control computer held time as an integer count of tenths of a second and converted it to a real number in registers that, the report states, "are only 24 bits long". Each conversion lost a little precision. The loss grew with how long the system had been switched on. The report's summary puts it in two sentences: “This problem led to an inaccurate tracking calculation that became worse the longer the system operated. At the time of the incident, the battery had been operating continuously for over 100 hours.”

Appendix II of the report is that accumulation, hour by hour. After 1 hour the calculated time is off by .0034 seconds and the range gate has shifted approximately 7 meters. After 8 hours, .0275 seconds and 55 meters. After 24 hours, .0687 seconds and 137 meters. After 48 hours, .1648 seconds and 330 meters. After 72 hours, .2472 seconds and 494 meters. After 100 hours, .3433 seconds and 687 meters. Israeli data in the same report showed "a 20 percent shift in the Patriot system's radar range gate after the system had been running for 8 consecutive hours". The Alpha Battery had been running for over a hundred.

The numerical-analysis reading of the same event ran in SIAM News a few months later. Robert Skeel reports the identical figures — "0.0275 seconds after eight hours and 0.3433 seconds after 100 hours" — and names the mechanism: "a fixed-point 24-bit representation of 0.1 in base 2". Each individual relative error is "very nearly 2-20, which is approximately 0.0001%". One ten-thousandth of one percent, per conversion. That is what 687 meters is made of.

Overflow is a magnitude exceeding the representable range and becoming infinity or an error. Underflow is tiny magnitudes rounding toward zero or into subnormals. Cancellation is subtracting close numbers and destroying the leading significant digits. Accumulated rounding is Appendix II. Ill-conditioning is different in kind from all four: the underlying problem amplifies small input perturbations no matter how you implement it.

Floating-point numbers represent a finite set of magnitudes and precisions. Most real numbers are rounded to the nearest representable value after each operation. The relative spacing between representable numbers grows with magnitude. Very large values overflow, very small values underflow, and subtracting nearby values can erase significant digits.

A formula can be algebraically correct and numerically unreliable. The tenths of a second the Patriot battery counted are exact as integers. As a binary fraction, a tenth is not exactly representable at all. The arithmetic that followed was valid at every step and wrong by 687 meters at the end. Numerical analysis asks whether the computed answer tracks the intended mathematical answer, which is a separate question from whether the derivation is correct.

Consider logits (1000, 999, 998). The overflow half of this example is standard. The underflow half has been measured, and the measurement is the more alarming of the two.

NVIDIA collected the activation gradient values of a Multibox SSD detector network during FP32 training, then asked what survives the move to half precision. In IEEE binary16 — 1 sign bit, 5 exponent bits, and 10 fractional bits — the minimum denormal is 2^-24 = ~5.96e-8 and the minimum normalized value is 2^-14 = ~6.10e-5. Below the first of those there is nothing left to represent: “When converted to FP16, 31% of these values become zeros, leaving only 5.3% as nonzeros which for this network lead to divergence during training.”

The paper behind that guide, published at ICLR in 2018, states the floor plainly: "any value whose magnitude is smaller than 2^-24 becomes zero in FP16". It reports that roughly 5% of weight-gradient values fall below that exponent, and confirms of the detector that "This particular network diverges when gradients are not scaled". The remedy is arithmetic, not architecture: "we propose scaling the loss appropriately to handle the loss of information with half-precision gradients", together with an FP32 master copy of the weights. The scale factor has a ceiling set by the other end of the format: "Choose a value so that its product with the maximum absolute gradient value is below 65,504 (the maximum value representable in FP16)." PyTorch ships the same remedy as torch.amp gradient scaling.

Separating these concepts makes numerical diagnosis far clearer, and each box has a worked instance elsewhere in this lesson.

Well-conditioned and stable is the ideal setting. Small input and rounding errors produce small output errors, results agree across reasonable implementations, higher precision refines predictably, and the remaining error is controllable.

Well-conditioned and unstable is the benign problem with an amplifying algorithm, and reformulation fixes it. The textbook variance formula on NIST's NumAcc4 is this box: the certified standard deviation is exactly 0.1, the data are unremarkable, and a shipped spreadsheet still returned zero correct significant digits. Explicit matrix inversion is the other common culprit, cancellation usually dominates, and alternative factorizations help.

Ill-conditioned and stable is a sound algorithm on a sensitive problem. Input uncertainty caps achievable accuracy, regularization or new information may be needed, more precision helps only partly, and condition numbers communicate the limit.

Ill-conditioned and unstable is both at once. Outputs may be meaningless, debugging by seeds is insufficient, the estimand itself should be reconsidered, and what gets reported is a sensitivity range rather than a point result. Note which box the Patriot failure is not in. The tracking calculation was not sensitive in principle, and no seed would have revealed it, because the error was a function of uptime alone.

Softmax should subtract the maximum logit before exponentiation. Log-sum-exp should use a shifted form. Probabilities should often be accumulated in log space. The expression log(1+x) is inaccurate for tiny x if computed by ordinary addition and logarithm; specialized log1p routines preserve precision. Similar functions exist for exp(x)−1. Linear systems should use stable solvers and factorizations rather than explicit inverses. Covariance and least-squares calculations need algorithms matched to conditioning.

Shipped software gets this wrong, and an agency built a dataset to prove it. NIST's Statistical Reference Datasets include NumAcc4, whose file header states: “This is a constructed/fabricated data set to test accuracy in summary statistic calculations. The numbers are 9-digit floating point values and differ only in the last decimal place.” It is 1001 observations: one value of 10000000.2 plus 500 pairs of 10000000.1 and 10000000.3. Its certified sample mean is 10000000.2, its certified sample standard deviation is exactly 0.1, and its certified lag-1 autocorrelation is -0.999. It carries the "Higher Level of Difficulty" label. The project's stated purpose is "to improve the accuracy of statistical software by providing reference datasets with certified computational results that enable the objective evaluation of statistical software."

The benchmark was run against shipped spreadsheets in 2010, in the Journal of Statistical Software. That study notes that the certified values were "calculated using multiple precision arithmetic to obtain 500 digits answers that were rounded to fifteen significant digits", and that "NumAcc4 is the only high difficulty dataset for univariate summary statistics". Its table of correct significant digits for the sample standard deviation gives NumAcc4 as Calc 8.2, Excel 2007 8.2, Excel 2008 0, Gnumeric 8.2, NeoOffice 8.2, and Oleo no correct digits at all. The true answer is 0.1. One product returned none of its digits, and the best of them returned 8.2 of 15 — on a formula every one of them implements and every one of their users believes.

The standard reference here is still a 1991 survey in ACM Computing Surveys: David Goldberg's What Every Computer Scientist Should Know About Floating-Point Arithmetic.

The expression log Σ exp(xᵢ) can overflow even when the final logarithm is moderate. In binary16 a single exponential above 65,504 is enough to lose the whole sum.

Run these checks before attributing instability to the model or data. Inspect ranges first — log minima, maxima, norms, NaNs, infinities and subnormal-sensitive quantities. Then prefer stable primitives: logsumexp, log1p, factorization solvers and fused operations. Then vary precision across float16, bfloat16, float32 and float64 on reduced tests. Then perturb and reorder, testing sensitivity to input noise, summation order, batch order and seeds. Then set tolerances from impact, judging numerical differences by predictions, metrics and decisions rather than exact bytes alone.

The third step is the one usually run without knowing what is being compared. float16 and bfloat16 are both sixteen bits, and they are not the same trade. bfloat16 keeps float32's 8-bit exponent and takes the saved bits out of the significand, buying range with precision. A 2019 study from Intel and Facebook states the consequence: “BFLOAT16 is attractive for Deep Learning training for two reasons: the range of values it can represent is the same as that of IEEE 754 floating-point format (FP32) and conversion to/from FP32 is simple.” Google's own documentation for the format says the same thing independently: "The dynamic range of bfloat16 and float32 are equivalent."

The consequence is a hyper-parameter that does or does not exist. That study covered image classification, speech recognition, language modelling, GANs and industrial recommendation, and reported state-of-the-art results in the same number of iterations with no hyper-parameter changes. It says why: "Maintaining the same range as FP32 is important to ensure that no hyper-parameter tuning is required for convergence; e.g., IEEE 754 compliant half-precision floating point (FP16) requires hyper-parameter tuning." The tuning it means is loss scaling — the 2^-24 floor from the softmax section. So step 3 is not a sweep for the fastest option. It is a question about the exponent field, and its answer decides whether step 4 has an extra knob in it.