Numerical stability and floating point

Overflow, underflow and cancellation, the log-sum-exp trick, stable softmax and cross-entropy, and where float16 quietly loses accuracy.

What floating point can and cannot do

A float32 has about 7 significant decimal digits and a range roughly from 1e-38 to 3e38. Those two constraints cause almost every numerical failure in model code: overflow past the top of the range, underflow to zero below the bottom, and catastrophic cancellation when two similar large numbers are subtracted.

import numpy as np

np.finfo(np.float32).eps      # 1.1920929e-07  machine epsilon
np.finfo(np.float32).max      # 3.4028235e+38
np.finfo(np.float16).max      # 65504.0        far smaller

# cancellation: the relative error explodes when two close numbers are subtracted
a = np.float32(1.0) + np.float32(1e-4)
b = np.float32(1.0)
(a - b)                       # 9.9897385e-05, noticeably wrong
(np.float64(1.0) + 1e-4) - 1.0  # 9.9999999e-05, much closer

# accumulation: summing millions of small values in float32 loses the tail
xs = np.random.default_rng(0).random(2_000_000).astype(np.float32)
diff = abs(float(xs.sum()) - float(xs.astype(np.float64).sum()))
diff / 1_000_000              # relative error per element

• Float32 has 24 bits of mantissa, so a sum around 1e7 cannot register an addend below about 1. That is why mean loss over a huge batch should be accumulated in float64 or with a running mean.
• Subtracting two nearly equal quantities loses precision regardless of dtype. Reformulate algebraically when you can (a stable formula) rather than casting to float64.
• Huge dynamic range is float16's real problem, not precision: the same value that is fine at 6e4 overflows at 65504.
• Comparison at tolerance is the right test: np.allclose(a, b, rtol=1e-5), never a == b.

Log-sum-exp and stable softmax

def logsumexp(z, axis=-1, keepdims=False):
    m = np.max(z, axis=axis, keepdims=True)
    out = m + np.log(np.exp(z - m).sum(axis=axis, keepdims=True))
    return out if keepdims else np.squeeze(out, axis=axis)

def log_softmax(z, axis=-1):
    return z - logsumexp(z, axis=axis, keepdims=True)

def softmax(z, axis=-1):
    return np.exp(log_softmax(z, axis=axis))

big = np.array([1000.0, 1001.0, 999.0])
np.exp(big)                      # [inf, inf, inf]  -> softmax would be nan
softmax(big)                     # [0.2447, 0.6652, 0.0900]  correct
logsumexp(big)                   # 1001.4076, finite

# cross-entropy from logits, stable end to end
def cross_entropy_from_logits(logits, targets):
    ls = log_softmax(logits)
    return float(-ls[np.arange(len(targets)), targets].mean())

The trick is to subtract the maximum before exponentiating. Softmax is invariant to adding a constant to all logits, so subtracting the row maximum changes nothing mathematically but keeps every exponent at or below zero — and the largest term becomes exactly 1.

Precision choices and mixed precision

# training in mixed precision loss-scales to keep small gradients representable
scaler = 1024.0
loss = compute_loss()                 # e.g. 0.001
scaled = loss * scaler                # 1.024, comfortably above float16 denormals
scaled.backward()                     # gradients are 1024x larger in float16

# before the optimiser step, unscale and skip the update if anything is inf/nan
def step_if_finite(optimizer, params, scaler):
    ok = all(np.isfinite(p.grad).all() for p in params)
    if ok:
        for p in params:
            p.grad = p.grad / scaler
        optimizer.step()
    else:
        scaler = scaler / 2            # back off and retry next step
    return scaler

# master weights stay in float32: the optimiser updates them, then casts down
master = np.zeros(1_000_000, dtype=np.float32)
np.finfo(master.dtype).eps             # 1.19e-07

• Compute in reduced precision, accumulate in float32: reductions like softmax, layer norm and the loss should be promoted, because that is where the range is exceeded.
• Keep a float32 master copy of the weights. Directly updating a float16 copy loses small increments entirely once the value is large.
• Loss scaling multiplies the loss before backward so small gradients do not flush to zero in float16, then divides the gradients back before the step.
• Test the numerics of a suspicion cheaply: rerun a short training in float64. If the problem disappears, it is precision related, not algorithmic.

FAQ

Related

Information theory for machine learning Matrix decompositions and PCA

Last refreshed 2026-09-18.