Floating point and numerical precision

Numerical bugs are silent: understanding float behaviour is what makes results trustworthy.

float64 stores about 15-16 significant decimal digits; float32 about 7. Most decimal fractions, like 0.1, have no exact binary representation, so tiny rounding errors appear: 0.1 + 0.2 == 0.3 is False. np.finfo(np.float64) reports the limits; its eps (about 2.2e-16) is the gap between 1.0 and the next float.

Never compare computed floats with ==. Use np.isclose(a, b) or np.allclose (or np.testing.assert_allclose in tests). They check |a - b| <= atol + rtol * |b| with defaults rtol=1e-05, atol=1e-08. The absolute tolerance is a trap for tiny values: np.isclose(1e-10, 2e-10) is True. Set atol=0 when your values are naturally small.

Classic hazards:

Catastrophic cancellation: subtracting nearly equal numbers wipes out significant digits. Use stable formulas: np.log1p(x) instead of np.log(1 + x), np.expm1, np.hypot.
Overflow: np.exp(1000) is inf. Compute softmax as exp(x - x.max()).
Accumulation: summing millions of values loses precision. np.sum uses pairwise summation, much better than a naive loop, while np.cumsum is sequential. math.fsum gives a correctly rounded sum when it really matters.
Lower precision: float32 and float16 save memory but errors grow fast. Accumulate in float64 with a.sum(dtype=np.float64).

Special values follow IEEE 754: inf, -inf, and nan, which is not equal to anything including itself, so test with np.isnan. np.errstate(over="raise", invalid="raise") turns silent infs and NaNs into exceptions while debugging.

Example

import numpy as np
import math

print(0.1 + 0.2 == 0.3, np.isclose(0.1 + 0.2, 0.3))   # False True
print(np.finfo(np.float64).eps)                      # 2.220446049250313e-16
print(np.isclose(1e-10, 2e-10), np.isclose(1e-10, 2e-10, atol=0))  # True False

x = 1e-12
print(np.log(1 + x), np.log1p(x))   # ~1.00009e-12 vs 1e-12 (accurate)

z = np.array([1000.0, 1001.0])
soft = np.exp(z - z.max()); soft /= soft.sum()
print(soft.round(3))                 # [0.269 0.731], no overflow

v = np.full(10_000_000, 0.1, dtype=np.float32)
print(v.sum(dtype=np.float64))       # 1000000.0149... (float32 0.1 != 0.1)
print(math.fsum([0.1] * 10))         # 1.0

print(np.nan == np.nan, np.isnan(np.nan))   # False True
with np.errstate(over="raise"):
    try:
        np.exp(np.array([1000.0]))
    except FloatingPointError as e:
        print("caught:", e)

Remember

Compare floats with isclose/allclose, and think about atol for small values
Use stable functions (log1p, expm1, max-shift softmax) to avoid cancellation and overflow
np.sum is pairwise and fairly accurate; accumulate low-precision data in float64
NaN != NaN: test with np.isnan; errstate can make float errors raise

Watch out for

Using np.isclose defaults on values around 1e-9 and getting True for numbers that differ by 100%.

Check yourself

Why does `np.isclose(1e-9, 5e-9)` return True?
Which is the numerically better way to compute log(1 + x) for very small x?

Official docs: numpy.org/doc/stable/reference/generated/numpy.finfo.html

Example

Remember

Watch out for

Check yourself

Builds on

Unlocks