Numerical stability and floating point

Computers do not compute exactly with decimal numbers. 0.1 + 0.2 becomes 0.30000000000000004, because 0.1 cannot be written exactly in base 2 — just as 1/3 cannot be written exactly in base 10.

Four traps:

The range of the common formats:

float16 overflows already at 65 504. That is why mixed-precision training needs loss scaling.

The log-sum-exp trick is the single most important device. You want to compute

$\log\sum_i e^{x_i}$

With $x = [1000, 1001]$, $e^{1000}$ overflows. But since

$\log\sum_i e^{x_i} = c + \log\sum_i e^{x_i - c}$

for any $c$, you choose $c = \max_i x_i$. Then the largest exponent is exactly 0, and $e^0 = 1$. No overflow, the same answer.

The same thing in softmax:

$\frac{e^{x_i}}{\sum_j e^{x_j}} = \frac{e^{x_i - c}}{\sum_j e^{x_j - c}}$

This is not a detail — it is built into every serious implementation of softmax and cross-entropy. It is also why cross_entropy in PyTorch takes logits, not probabilities: it combines log-softmax and NLL in one numerically stable operation. Running softmax first and then log is both slower and unstable.

Cancellation arises when two nearly equal numbers are subtracted: the shared digits cancel and only the junk at the end remains.

Never == on floating point. Use math.isclose or np.allclose with a tolerance. And in tests: give both an absolute and a relative tolerance, since one suffices for small numbers and the other for large ones.

A practical rule in ML: work in log space as long as you can, convert back as late as possible, and trust the libraries' logsumexp, log_softmax and cross_entropy instead of writing them yourself.

import math, numpy as np
from scipy.special import logsumexp

# 1. Overflow
x = np.array([1000.0, 1001.0, 1002.0])
try:
    print(np.log(np.exp(x).sum()))
except Exception:
    pass
print(np.log(np.exp(x).sum()))          # inf  (with a warning)
print(logsumexp(x))                     # 1002.4076059644443  ← stable

# The same thing by hand
c = x.max()
print(c + np.log(np.exp(x - c).sum()))  # 1002.4076059644443

# 2. Underflow
print(math.prod([0.1] * 400))                       # 0.0  — all the information gone
print(sum(math.log(0.1) for _ in range(400)))       # -921.03

# 3. Cancellation
print((1 + 1e-16) - 1)                   # 0.0
print(math.log(1 + 1e-16), math.log1p(1e-16))       # 0.0  1e-16

# 4. Comparison
print(0.1 + 0.2 == 0.3)                             # False
print(0.1 + 0.2)                                    # 0.30000000000000004
print(math.isclose(0.1 + 0.2, 0.3))                 # True

# 5. Softmax: naive against stable
def softmax_naive(z):
    e = np.exp(z); return e / e.sum()

def softmax(z):
    e = np.exp(z - np.max(z)); return e / e.sum()

z = np.array([1000.0, 1001.0, 1002.0])
print(softmax_naive(z))    # [nan nan nan]
print(softmax(z))          # [0.09003057 0.24472847 0.66524096]

# 6. Why cross_entropy takes logits, not probabilities
import torch
logits = torch.tensor([[1000.0, 1001.0, 1002.0]])
target = torch.tensor([2])
print(torch.nn.functional.cross_entropy(logits, target))         # tensor(0.4076)
p = torch.softmax(logits, -1)
print(torch.nn.functional.nll_loss(torch.log(p), target))        # tensor(0.4076) but by an unstable route