Bezen

NumPy · Module 11 of 13 · 7 concepts · about 9 min

Numerical Stability

Why correct-looking formulas produce inf and nan on real data, and the standard rewrites that stop them.

At work: Keeping calculations like probabilities and risk scores from quietly turning into infinity or “not a number” in production.

Floating point, briefly

What it is

float64 holds about 15-16 significant decimal digits, float32 about 7. Values are spaced logarithmically, so adding a small number to a large one can change nothing at all.

Why it matters

It is why 0.1 + 0.2 != 0.3, why summing a million small floats loses accuracy, and why a float32 accumulator of a long sum can stall entirely.

In the real world: Why 0.1 + 0.2 isn't exactly 0.3, and why money is often stored in whole paise or cents in payment systems.

0.1 + 0.2 == 0.3                          # False
np.isclose(0.1 + 0.2, 0.3)                # True — the right comparison
np.float32(1e8) + np.float32(1.0) == np.float32(1e8)   # True: 1.0 vanished
np.arange(1e6, dtype=np.float32).sum()     # noticeably off; pass dtype=np.float64

Comparing floats with ==, including against a value you just computed. Use np.isclose or np.allclose, with atol set deliberately when the numbers can be near zero.

Practise it: Running Share of a Total

inf, nan, and how they propagate

What it is

np.inf comes from overflow or division of a nonzero by zero; np.nan comes from 0/0, inf - inf, sqrt of a negative, and any arithmetic involving an existing nan. nan is not equal to itself.

Why it matters

One nan in a batch makes every gradient computed from it nan, so the training loss goes to nan and stays there. Finding the first one is a debugging skill.

In the real world: One “nan” in a feature column can make a model's predictions all “nan”, which then fills a dashboard with blanks.

np.array([1.0]) / 0.0        # inf, with a RuntimeWarning
np.array([0.0]) / 0.0        # nan
np.nan == np.nan             # False — always use np.isnan
np.isfinite(X).all()         # the assertion worth putting in a pipeline
np.argwhere(~np.isfinite(X))[:5]   # where the first offenders are

Testing for nan with x == np.nan, which is always False, or sorting an array containing nan — every comparison with it is False, so nan sorts to the end and argmax can return it.

Overflow in exp, and the log-sum-exp trick

What it is

exp(800) is already inf in float64. Subtracting a constant from every input before exponentiating leaves ratios unchanged, so logsumexp(x) = max(x) + log(sum(exp(x - max(x)))) is exact and finite.

Why it matters

Every softmax, every log-likelihood and every cross-entropy is built on it. Without the shift, real logits — which routinely reach the hundreds — produce inf / inf, that is, nan.

In the real world: Keeping a recommendation model's scores from overflowing to infinity when some items score very high.

x = np.array([1000.0, 1001.0])
np.log(np.exp(x).sum())                    # inf: overflowed first
m = x.max()
m + np.log(np.exp(x - m).sum())            # 1001.31... — correct
from scipy.special import logsumexp        # the library version, same trick

Worrying about underflow instead. exp(-800) becoming 0.0 is harmless — it was negligible — while overflow to inf destroys the whole result. Shift by the maximum, never by the minimum.

Practise it: Softmax That Does Not Overflow

Softmax, done properly

What it is

Subtract the per-axis maximum with keepdims=True, exponentiate, then divide by the per-axis sum, also with keepdims=True. Two reductions, both keeping their axis.

Why it matters

It is the most frequently written numerically sensitive expression in machine learning, and both of its failure modes — overflow, and broadcasting along the wrong axis — are silent.

In the real world: Turning a classifier's raw scores into probabilities like “92% spam” that always add up to 100%.

Visual example

Softmax

Subtract the row max first: the largest exponent becomes exp(0) = 1, so nothing can overflow and the ratios are unchanged.

logits — shape (2, 3)
012
a10001001999
b01-1
shift, exp, normalise
softmax(axis=1)
012
a0.2450.6650.09
b0.2450.6650.09

Row a, step by step

a: max = 1001

a: shifted = [-1, 0, -2]

a: exp = [0.37, 1, 0.14]

a: / 1.51 → [0.245, 0.665, 0.09]

Why this matters

The two rows differ by 1000 and give identical answers: softmax depends only on the differences. The shift exploits that to keep every exponent at or below 1.

Common misconception

Wrong idea: the shift is an approximation that trades accuracy for safety.

Correct idea: it is exact. Multiplying numerator and denominator by exp(-c) cannot change a ratio — the naive version is the one that is wrong, because exp(1001) is inf.

  1. 1ShiftSubtract x.max(axis=axis, keepdims=True), so the largest input becomes 0.
  2. 2ExponentiateEvery value is now in (0, 1]: no overflow possible.
  3. 3NormaliseDivide by the sum along the same axis, again with keepdims=True.
def softmax(x, axis=-1):
    x = np.asarray(x, dtype=float)
    shifted = x - x.max(axis=axis, keepdims=True)
    e = np.exp(shifted)
    return e / e.sum(axis=axis, keepdims=True)

Dropping keepdims=True. On a square array of logits the division still broadcasts — across the wrong axis — and returns rows that do not sum to 1 without raising anything.

Practise it: Softmax That Does Not Overflow

Safe division

What it is

np.divide(a, b, out=..., where=b != 0) computes the quotient only where the denominator is usable and leaves your chosen default elsewhere. Flooring the denominator with np.maximum(b, eps) is the other standard form.

Why it matters

Every per-group mean, rate and normalisation has a denominator that can be zero — an empty group, a constant column, a zero vector — and the choice of what that case returns is yours to make rather than the arithmetic's.

In the real world: Computing conversion rate (orders ÷ visits) for every page without crashing on pages that had zero visits.

Visual example

Safe division

A guarded division never evaluates 0 / 0 — the empty group gets the value you chose, not a nan.

sums and counts
sumcount
0302
100
251
divide(where=count > 0)
means
mean
015
10
25

Group by group

0: 30 / 2 = 15

1: 0 / 0 → skipped, 0.0

2: 5 / 1 = 5

Why this matters

An empty group is normal, not exceptional: a user with no clicks, a class absent from a batch, a segment with no rows. The pipeline has to keep running.

Common misconception

Wrong idea: np.where(count > 0, sums / counts, 0) avoids the division by zero.

Correct idea: it computes every quotient first, warning and all, then discards some. where= on the ufunc skips the work itself.

  1. 1Provide the defaultout=np.zeros_like(sums) decides what an unusable group returns.
  2. 2Mask the workwhere=counts > 0 computes only the legal positions.
  3. 3Nothing to clean upNo warning is raised and no nan ever exists.
sums = np.array([30.0, 0.0, 5.0])
counts = np.array([2.0, 0.0, 1.0])
np.divide(sums, counts, out=np.zeros_like(sums), where=counts > 0)
# array([15.,  0.,  5.]) — no warning, no nan
rates = clicks / np.maximum(impressions, 1)    # the floor form

np.where(counts > 0, sums / counts, 0.0). The division is still evaluated everywhere, so the warning is still raised and the nan is still created before being thrown away — and in a float32 pipeline the discarded inf can still trip an error handler.

Practise it: Group Sums Without a GroupBy, Intersection Over Union for Every Box Pair

clip, and clamping to a valid range

What it is

np.clip(x, lo, hi) bounds every element; either bound can be None. np.maximum(x, 0) is the one-sided form, and the standard way to say 'a negative overlap is no overlap'.

Why it matters

It is the guard that turns an impossible intermediate into a correct one: a negative intersection width, a probability slightly above 1, a squared distance of -1e-16, an index that ran off the end of a row.

In the real world: Keeping a predicted discount between 0% and 50% even if the model suggests something wild.

w = np.maximum(x2_min - x1_max, 0.0)        # a negative width means no overlap
p = np.clip(p, 1e-7, 1 - 1e-7)              # before log, for cross-entropy
np.sqrt(np.maximum(d2, 0.0))                 # cancellation can go negative
idx = np.clip(idx, 0, x.shape[1] - 1)        # keep a gather in bounds

Clipping only the upper bound of the intersection in an IoU calculation. A negative width multiplied by a negative height is a *positive* area, so two boxes that do not overlap at all report an overlap — the one bug everybody writes once.

Practise it: Clip Outliers to a Percentile Range

errstate and seterr

What it is

NumPy's floating-point warnings are configurable: np.seterr(divide='raise', invalid='raise') turns them into exceptions, and with np.errstate(divide='ignore'): silences them for one block.

Why it matters

Raising is the right default in a pipeline — it turns a silent nan into a stack trace at the line that made it — and the context manager is how you allow the one place where the inf is intentional.

In the real world: Making a risk-calculation job stop loudly on a division by zero, instead of quietly producing wrong numbers overnight.

np.seterr(divide='raise', invalid='raise')     # fail loudly during development
with np.errstate(divide='ignore'):
    log_p = np.log(p)                          # -inf here is deliberate

Leaving warnings at their default 'warn' setting, which prints once per location and then stays quiet. The second million nans are silent, so the log gives no sense of how widespread the problem is.