Numerical Stability Is Part of the Algorithm
Real-number algebra does not specify a machine implementation. Precision, operation order, exceptional values, and failure policy belong in the algorithm contract.
- Authorship
- InterviewsVector
- Published / updated
- 2026-08-14 / 2026-08-14
- Review status
- Artifact tests passing · primary sources recorded
Original InterviewsVector material. Code examples are covered by repository tests and primary references are recorded. No named human reviewer is claimed.
The decision in one pass
Numerical stability is an algorithm property, not a cleanup pass. Declare input, accumulation, and output precision plus the only allowed evidence source; validate exact shape, contract and source identity, finite values, and representability; choose a formulation whose intermediate range is controlled; use an appropriate reduction; and assert output invariants. Compute log-sum-exp as m + log Σ exp(xᵢ−m), softmax from the same maximum-shifted terms, and negative log likelihood directly from logits instead of logging a rounded probability. Treat a misrouted source, overflow, unexpected underflow, NaN, infinity, precision drift, or invalid probability sum as evidence failure. Then test the actual kernel and hardware, because a Python binary64 harness cannot certify an accelerator reduction tree.
Why this matters
A model can be mathematically correct and operationally wrong. Exponentials overflow at moderate logits in finite formats, small probabilities underflow, large-plus-small reductions lose contributions, and changing accumulation precision can shift thresholds or gradients without a schema change. These failures often produce plausible zeros or confident scores, so silent handling is more dangerous than a loud exception.
You will be able to
- Distinguish a real-number identity from the finite sequence of operations that implements it.
- Derive stable log-sum-exp, softmax, and logit-space negative log likelihood formulations.
- Explain how input, accumulation, and output precision affect range and rounding independently.
- Detect overflow, harmful underflow, non-finite values, cancellation, and reduction-order sensitivity.
- Build and operate a precision contract whose identity changes with numerical policy and allowed evidence source.
Prerequisite contract
- — Exponentials, logarithms, and probability ratios
- Optimization Under Noise →
- Reproducible Experiments, Not Reproducible Notebooks →
Your Vector Loop for this lab
- 01
Model
Specify the values, formats, operation sequence, invariants, and failure policy.
- 02
Derive
Rewrite exponential normalization and likelihood to bound intermediates and avoid cancellation.
- 03
Build
Implement a precision-bound harness with stable reductions and explicit exceptional states.
- 04
Stress
Overflow exponentials, underflow representations, cancel sums, drift dtypes, and corrupt shapes.
- 05
Operate
Qualify the actual kernel by device, compiler, precision mode, slice, and numerical envelope.
- 06
Defend
Explain why algebraic equivalence does not guarantee the same finite-precision result.
Model the machine algorithm, not only the equation
Floating-point formats represent a finite subset of real numbers. Each operation rounds to a destination format, and the sequence of operations determines which intermediates exist. Two algebraically equal expressions can therefore overflow at different points, discard different addends, or round across a decision threshold. The contract must name the format and operation boundary, not merely say float.
| Contract field | Question | Silent failure if absent |
|---|---|---|
| input dtype | what values arrive and are representable? | upstream rounding or saturation is hidden |
| accumulation dtype | where are products and sums rounded? | low-precision error grows with reduction length |
| output dtype | what reaches the consumer? | valid internal probability becomes zero at storage |
| operation order | which intermediates are materialized? | an equivalent rewrite overflows or cancels |
| exception policy | what happens on NaN, infinity, overflow, or underflow? | invalid values become plausible output |
fl(a ⊕ b) = round(a ⊕ b); generally fl(fl(a+b)+c) ≠ fl(a+fl(b+c))
Finite addition is generally non-associative. Parallel reduction trees, compiler transformations, and worker order can change the result even when every input is identical.
Break the textbook forms before repairing them
Equivalent formulas, different machine behavior
Change a shared maximum logit and two relative gaps; compare naive and maximum-shifted softmax and log-sum-exp; inspect overflow or underflow state and shift invariance; then complete a graded formulation check.
Keep softmax finite without changing its answer
Move the common logit level and the class gaps. Softmax should depend on the gaps, but direct exponentiation also depends on whether the absolute values fit in floating-point range.
Current logits
[1000, 999, 998]
Subtracting the maximum converts these to [0, −1, −2]. Adding or removing the same constant from every logit leaves the probability ratios unchanged.
| Computation | Naive exp(x) | Stable exp(x − max) |
|---|---|---|
| P(class 1) | undefined | 0.665241 |
| P(class 2) | undefined | 0.244728 |
| P(class 3) | undefined | 0.090031 |
| log ∑ exp(x) | +Infinity | 1000.407606 |
Naive arithmetic has failed
Direct exponentiation produced an infinite or zero denominator and undefined probabilities. Max subtraction keeps at least one exponent equal to 1.
Make a prediction, then check it against the current evidence.
- 01Overflow softmaxEvaluate exp(xᵢ) directly for large positive logits. A denominator of infinity can turn a well-defined probability vector into exceptions or NaN.
- 02Underflow a rare classUse a large negative logit gap. Zero may be acceptable under an explicit approximation policy, but silently producing it destroys likelihood and gradient information.
- 03Log a rounded probabilityIf the target probability rounded to zero, −log(p) becomes infinite even though the logit-space loss is finite.
- 04Accumulate beyond the format's resolutionRepeatedly add small values to a larger low-precision accumulator. Updates can stop changing the stored sum long before the workload ends.
Derive forms with bounded intermediates
LSE(x) = m + log Σᵢ exp(xᵢ − m), where m = maxᵢ xᵢ
Every shifted exponent is at most one, so the exponential step cannot overflow for finite inputs. The largest term equals one and keeps the denominator positive.
softmax(x)ᵢ = exp(xᵢ − m) / Σⱼ exp(xⱼ − m)
Adding a common offset to every logit leaves probabilities unchanged. The maximum shift uses that invariance to control range, while an accurate sum reduces denominator error.
NLL(y, x) = log Σᵢ exp(xᵢ − m) − (xᵧ − m)
Compute classification loss directly from logits. This avoids forming a tiny rounded probability and then applying a logarithm.
| Technique | What it controls | What it does not guarantee |
|---|---|---|
| maximum shift | positive-exponential overflow | no harmful output underflow in every format |
| accurate or compensated sum | rounding loss in a reduction | correct inputs or adequate destination range |
| wider accumulation | intermediate range and precision | safe cast into a narrow output |
| logit-space loss | log of a rounded zero probability | finite or semantically valid logits |
| output invariant checks | silent invalid probabilities | that the chosen tolerance fits product risk |
Build a harness that binds precision to evidence
1def stable_logsumexp(contract: PrecisionContract, vector: NumericVector) -> float:2 values = validate_vector(contract, vector)3 shift, exponentials = _shifted_exponentials(values)4 denominator = math.fsum(exponentials)5 if not math.isfinite(denominator) or denominator <= 0.0:6 raise FloatingPointError("log-sum-exp denominator is not positive and finite")7 result = shift + math.log(denominator)8 if not math.isfinite(result):9 raise FloatingPointError("log-sum-exp result overflow")10 return resultExpected output
example=illustrative_only
contract_version=numerical-stability-v1
unsafe_softmax=OVERFLOW
stable_softmax=0.090031,0.244728,0.665241
logsumexp=1002.407606
nll_target_2=0.407606
naive_binary16_sum=2048.0
stable_binary64_sum=4096.0Verify: Run python3 -m unittest discover courses/ai-engineering/reference-impl/numerical_stability.
The excerpt is the exact log-sum-exp implementation from the downloadable artifact. Its stable path requires concrete frozen PrecisionContract and NumericVector records, rejecting mutable duck-typed aliases, then verifies a content-derived precision contract, an exact match between vector source and the contract's allowed source, immutable rank-one values, maximum workload length, declared IEEE input representation, binary64 accumulation and output, finite values, accountable numerical owner, representability, and probability-sum tolerance. The unsafe helpers are named demo-only and deliberately reproduce overflow and binary16 accumulator stagnation. The fixtures are not measurements of a deployed kernel.
Test exceptional states and semantic boundaries
| Adversarial case | Required response | Why |
|---|---|---|
| NaN or infinity in logits | reject before max or reduction | ordering and normalization invariants no longer hold |
| binary16 input outside range | representation-overflow failure | declared dtype cannot carry the evidence |
| nonzero input rounds to zero | representation-underflow failure | silent information loss violates policy |
| shifted exponential becomes zero | exponential-underflow failure | zero was not authorized as an approximation |
| dtype differs from contract | identity failure | results are not comparable across precision |
| source differs from allowed source | identity failure | values from another model or pipeline are misrouted evidence |
| target index is boolean or out of range | schema failure | language coercion must not select a class |
| vector exceeds length bound | workload failure | accuracy and capacity evidence is scoped |
| probabilities miss unit sum tolerance | output invariant failure | downstream decisions expect a distribution |
Underflow is not universally a bug: some applications explicitly accept negligible probabilities becoming zero. That decision must name a threshold, affected operation, downstream consequence, and test. This artifact chooses a strict reject policy so silent loss cannot be mistaken for approved approximation.
Operate a numerical envelope, not a dtype label
- 01DeclareVersion input, accumulation, output, intermediate formulations, compiler mode, device class, workload limits, tolerance, and exception policy.
- 02QualifyCompare representative and adversarial slices against a trusted higher-precision path, including extreme logits and long reductions.
- 03CanaryObserve non-finite rates, saturation, zero probabilities, normalization residuals, decision flips, and latency under the exact release identity.
- 04EscalateQuarantine invalid outputs and preserve inputs and kernel identity; do not rewrite exceptional values in a generic middleware layer.
- 05RequalifyTreat device, driver, compiler, kernel, graph rewrite, precision, or reduction-topology changes as numerical changes.
Operate at three altitudes
Production lens
- — Bind the allowed evidence source, input, accumulation, and output precision plus kernel, compiler, device, and formulation identity to every release.
- — Track non-finite values, saturation, underflow zeros, normalization residuals, and decision flips by workload slice.
- — Keep a trusted higher-precision comparison path for qualification and incident replay, not for every live request.
- — Fail closed or invoke a named approximation policy; never let generic NaN-to-zero handling redefine model semantics.
Staff lens
- — Assign ownership across model authors, kernel/runtime teams, serving infrastructure, and downstream policy consumers.
- — Define numerical envelopes per operation and risk level instead of imposing one global dtype rule.
- — Budget requalification for compiler and hardware upgrades because operation order and accumulation modes can change.
- — Make precision changes observable in release identity and rollback evidence even when tensor schemas are unchanged.
Interview defense
Why does naive softmax overflow, and how would you make the whole path production-safe?
Direct exponentials can exceed the finite format range. I would subtract the maximum logit, use the same shifted terms for log-sum-exp and softmax, compute negative log likelihood directly in logit space, and use an appropriate accumulation method and dtype. Around the formula I would validate exact shapes, finite and representable inputs, input/accumulation/output precision, workload bounds, and output invariants such as a finite unit-sum distribution. I would define whether underflow is rejected or an approved approximation, bind kernel, compiler, and device identity, and compare adversarial slices with a trusted higher-precision implementation on the actual hardware.
Expect the interviewer to press on
- — Does subtracting the maximum eliminate every underflow?
- — Why can a parallel reduction change a reproducible result?
- — When is a zero probability acceptable?
Misconceptions to remove
“Mathematically equivalent formulas are interchangeable in code.”
Finite formats round every operation and constrain intermediate range, so operation sequence can change overflow, underflow, cancellation, and the final rounded result.
“Mixed precision means inputs and outputs are float16.”
Input storage, multiplication, accumulation, intermediate storage, and output casting can use different formats. Name each boundary.
“Subtracting the maximum makes softmax numerically perfect.”
It prevents positive-exponential overflow and reduces harmful underflow, but small terms can still underflow, sums still round, and output casts can still lose information.
Check your model
1. Why is the largest shifted exponential exactly exp(0)?
Choosing m as the largest logit makes every xᵢ−m non-positive and at least one equal to zero, so the denominator contains a one and cannot vanish for valid finite arithmetic.
2. Why compute negative log likelihood from logits rather than from softmax output?
A small probability may round to zero before the logarithm. The shifted logit-space expression controls intermediates and avoids that extra lossy boundary.
3. What does changing only accumulation precision affect?
It changes intermediate range and rounding during reductions even if input and output tensor dtypes stay the same, so the release identity and qualification evidence must change.
Prove the mechanism
Create a precision contract for a classification-logit path. Implement maximum-shifted log-sum-exp, softmax, direct logit-space negative log likelihood, and an accurate reduction. Reject wrong contract or source identity, shape, dtype, finite state, workload length, representation overflow or underflow, exponential underflow, and invalid probability sum. Preserve an explicitly unsafe version only as a test fixture.
Add a production constraint
Run the same adversarial vectors through two real kernel/device configurations and a trusted higher-precision reference. Record maximum absolute and relative error, normalization residual, class-decision flips, exceptional states, compiler and device identity, and a risk-based acceptance policy.
Artifact: Numerical stability harness
courses/ai-engineering/reference-impl/numerical_stability/numerical_stability_harness.py
Download reference implementationPrimary references and next links
References
- 1. IEEE Standard for Floating-Point Arithmetic
IEEE 754-2019. Active primary standard defining binary floating-point formats, operations, exceptions, and the role of operation sequence and destination format.
- 2. Accurately computing the log-sum-exp and softmax functions
Blanchard, D. Higham, and N. Higham. Primary rounding-error analysis supporting maximum-shifted formulations and their overflow, underflow, and accuracy properties.
- 3. Python struct format characters
Python documentation. Official documentation for the IEEE binary16, binary32, and binary64 interchange representations used by the harness.
- 4. Python math.fsum
Python documentation. Official documentation for the accurate summation primitive used in stable denominators and reduction evidence.
Continue through the graph
- Optimization Under Noise →
Separate stochastic variation from finite-precision failure in training.
- Complexity for Tensor Programs →
Connect dtype choices to memory traffic and activation storage.
- Reproducible Experiments, Not Reproducible Notebooks →
Carry runtime and hardware identity into numerical qualification.
Glossary: floating point · overflow · underflow · rounding · cancellation · log-sum-exp · accumulation dtype · numerical envelope