≈ Numerical Computing for ML

Your formula is correct. Your loss is NaN. Both statements can be true because a computer stores a finite approximation to the numbers in the formula.

This lesson is about recognizing when the implementation changes the answer. Work through linear algebra and calculus first; information theory explains the log-probability losses used below. By the end, you should be able to stabilize a softmax, explain a precision choice, and debug the first nonfinite value instead of merely hiding it.

1. Floating point has limited range and precision

A floating-point number stores a sign, an exponent and a significand. The exponent controls the scale; the significand controls detail within that scale. Not every decimal is representable. In binary floating point, 0.1 + 0.2 commonly differs slightly from 0.3.

Format	Exponent bits	Stored fraction bits	Practical implication
float64	11	52	Useful for small reference calculations and gradient checks
float32	8	23	Common accumulation and optimizer-state format
float16	5	10	More precision near 1 than bfloat16, but much narrower range
bfloat16	8	7	Range similar to float32, with coarser precision

For normal numbers there is also an implicit leading significand bit. Spacing grows with magnitude: float32 represents every integer through 2²⁴, but not every integer above it. Adding 1 to 16,777,216 in float32 rounds back to 16,777,216 under the usual round-to-nearest rule.

Overflow exceeds the representable finite range. Underflow loses very small magnitudes, sometimes through subnormal values and sometimes by flushing them to zero depending on hardware/settings. NaN means an invalid result such as 0/0; it propagates through many operations. These are distinct from ordinary small rounding errors.

2. Algebraic equivalence is not numerical equivalence

Take logits z=[1000,1001,1002]. Directly computing exp(z) overflows even in float64. Softmax then risks inf/inf, which is NaN.

Subtract the maximum m=1002 first:

softmax(z)_i = exp(z_i−m) / Σ_j exp(z_j−m)

The shifted logits are [−2,−1,0], exponentials approximately [0.1353,0.3679,1], and probabilities approximately [0.0900,0.2447,0.6652]. A common shift cancels algebraically, so the distribution is unchanged.

For log probabilities, use log-sum-exp:

logsumexp(z) = m + ln Σ exp(z_i−m)

Here it is about 1002.4076. If the correct class is the last one, negative log-likelihood is logsumexp(z)−1002 ≈ 0.4076 nats. Prefer a framework's fused cross-entropy/log-softmax routine over taking log after a possibly underflowed softmax.

Check: would subtracting 1000 instead change the exact probabilities?

No: any common shift cancels mathematically. Subtracting the maximum is a useful numerical choice because all exponentials are at most 1 and at least one is 1. It avoids overflow in the exponentials when the input logits are finite.

An attention row with every position masked is an edge case: all logits may be −∞, so subtracting the maximum gives −∞−(−∞). Ensure each row has a valid key or define an explicit empty-row policy. Replacing NaNs with zero afterward can conceal an invalid mask.

3. Cancellation and accumulation

Subtracting nearly equal rounded numbers can destroy meaningful digits. For small x, computing ln(1+x) directly can round 1+x to 1 before taking the log. Use log1p(x). Use expm1(x) for exp(x)−1 near zero.

The variance identity E[X²]−E[X]² is exact algebraically but can be unstable when a large mean hides a small variance. For observations around one billion that differ by only a few units, both terms are enormous and nearly equal. Use a centered two-pass calculation or a stable online algorithm such as Welford's method.

Addition is not associative in floating point. With large magnitudes, (a+b)+c can differ from a+(b+c). Pairwise summation and higher-precision accumulation often reduce error. Distributed reductions can change addition order, so tiny differences across devices or runs do not automatically imply a logic bug.

Check: is exact equality a good test for two implementations of a floating-point reduction?

Usually not. Use a justified tolerance with both absolute and relative terms, and compare against a suitable reference. Absolute tolerance matters near zero; relative tolerance scales with magnitude. A very loose tolerance can also hide a real bug, so derive it from the operation and acceptable task error.

4. Conditioning belongs to the problem

Suppose a linear system is:

x+y=2; x+1.0001y=2.0001

The solution is x=1,y=1. Change the second right-hand side to 2.0002—a tiny input perturbation—and the solution becomes x=0,y=2. The equations are almost redundant, so recovering x and y separately is sensitive.

A condition number measures sensitivity to perturbations; for an invertible matrix in the 2-norm, κ(A)=σ_max/σ_min. Large κ means some input errors can be amplified. A stable algorithm avoids introducing much more error than small input perturbations would explain. Stability concerns the algorithm; conditioning concerns the problem. Even a stable solver cannot make an ill-conditioned problem insensitive.

For least squares, forming XᵀX squares the 2-norm condition number when X has full column rank. Solve using QR or SVD rather than explicitly computing (XᵀX)⁻¹. Feature scaling can improve geometry; regularization can improve conditioning but changes the estimation problem. Adding an arbitrary ε is not a free accuracy fix.

5. Gradient checks need a sensible step

The central-difference approximation for parameter i is:

[L(w+h e_i)−L(w−h e_i)]/(2h)

Its truncation error shrinks with h for a sufficiently smooth loss, but rounding and subtraction error grow in relative importance as h gets tiny. Try a range of h, use float64 on a small deterministic case, and compare relative as well as absolute error. Avoid ReLU kinks, random dropout and changing batches.

Automatic differentiation differentiates the implemented computation graph. It avoids finite-difference truncation, but can still suffer underflow, overflow, unstable operations and inappropriate gradient conventions. Check gradient shapes and scales as well as finiteness.

6. Mixed precision: choose which quantities can be coarse

Mixed-precision training commonly computes selected operations in a lower precision while retaining sensitive accumulations or optimizer state in a higher precision. The exact layout depends on the framework, optimizer, kernels and distributed strategy; “16-bit training” does not mean every tensor uses two bytes.

FP16 has a maximum finite value of 65,504 and a narrower exponent range than BF16. Small FP16 gradients may underflow. Loss scaling multiplies the loss by a positive scale before backpropagation, then divides gradients by that scale before the optimizer step. Dynamic scaling reduces the scale after nonfinite gradients and may skip that update. Unscale before norm clipping so the threshold has the intended meaning.

BF16's wider range often avoids the need for FP16-style loss scaling, but its coarse significand still creates rounding error. Higher-precision accumulation helps many matrix operations; actual kernel behavior must be checked. Low-precision formats do not guarantee a speedup on hardware without efficient support.

Check: if you scale the loss by 1,024, what must happen before updating weights?

The gradients must be divided by 1,024, checked for nonfinite values, and then clipped if clipping is used. Otherwise the optimizer sees a different scale. The scaler/framework may manage the order, so follow its API rather than applying the factor twice.

7. A practical debugging order

When a run first produces a nonfinite loss, reproduce the smallest failing batch. Inspect the inputs and labels, then intermediate activations, loss terms, gradients and optimizer state to find the first invalid value. Check division by zero, logarithm domains, invalid masks, large logits, overflowed gradients and empty reductions.

Compare a small float32/float64 reference with the lower-precision path. If both fail, investigate the formula or data. If only one fails, investigate range, accumulation and precision policy. Record dtype, tensor shape, scale, seed and kernel settings. Do not “fix” the run by globally replacing NaNs; first explain why they appeared.

Lower precision also underlies quantization, where values map to discrete levels with scales and possible zero-points. That introduces a different approximation from merely choosing BF16. Distributed training adds reduction-order and state-layout choices; model serving adds throughput and latency constraints.

Continue to ML Fundamentals for the full train/validation workflow. Keep this lesson nearby when a correct derivation fails in code.

Sources and further practice

• What Every Computer Scientist Should Know About Floating-Point Arithmetic — Goldberg's explanation of rounding, cancellation and stability.
• PyTorch: numerical accuracy — precision, reductions and backend-dependent behavior.
• PyTorch: automatic mixed precision examples — scaling, unscaling and clipping in the correct order.
• SciPy: logsumexp — stable log-sum-exp API and examples.