Softmax overflow NaN in training: fix with per-row max subtraction

A hand-written classifier whose loss list turns to NaN while NumPy prints RuntimeWarning: invalid value encountered in subtract at exp_a = np.exp(a-c) has two defects working together, and both hide in lines most people copy without a second look. The softmax subtracts one maximum and divides by one total for the whole batch, and the weights start at weight_init_std = 0.01 in a ReLU network deeper than two layers. This page shows a two-second test for each defect, the arithmetic behind the overflow, and the corrected code that trains the failing 784-50-50-50-10 MNIST net to a loss below 0.3.

What the warning is telling you

NumPy classifies inf - inf and NaN minus a finite number as an invalid value, not an overflow, so the warning fires inside the subtract before np.exp evaluates anything. Either the score array a or the shift c already carries a non-finite value. With c = np.max(a) taken over the whole batch, a single poisoned row makes c inf, and a - c then returns NaN for all 1,000 entries of a (100, 10) batch.

The overflow exp itself would hit sits elsewhere. float64 np.exp returns inf once its argument passes 709.78, because exp(709.78) ≈ 1.8e308 is the largest float64; float32 hits its wall at 88.72, where exp reaches 3.4e38. The shift exists to keep every exponent on the safe side of those limits:

Subtract the row's own maximum and exp overflow becomes structurally impossible: no logit can push an exponent above zero. Every NaN in this training loop therefore comes from a non-finite value upstream, from the wrong max scope, or from both.

Cause 1: the softmax subtracts and sums over the whole batch

The softmax in these training snippets is written for one sample: a 1-D vector of ten scores. Called on a batch of shape (100, 10), np.max(a) returns one number for all 1,000 values, and so does np.sum(exp_a). Two things break at once.

First, the normalization. Dividing by the batch-wide total makes each row's outputs sum to roughly 1/100 = 0.01 instead of 1, so every probability is about 100 times too small. The cross entropy floor moves from ln 10 = 2.303 to about ln 1000 = 6.91, because the correct class sits near 0.001 instead of 0.1.

Second, the shift is shared. One non-finite score anywhere in the batch turns c into inf, and the subtract then marks all 1,000 outputs NaN. That is the reported warning, and it explains why the loss list flips to NaN wholesale instead of row by row.

The test takes two lines:

import numpy as np

scores = np.random.default_rng(0).normal(size=(100, 10))  # example batch
y = softmax(scores)                            # your current implementation
print(np.allclose(y.sum(axis=1), 1.0))         # True only for a row-wise softmax

A correct softmax prints True. The batch-global version prints False, with every row summing to about 0.01 at this batch size.

Cause 2: weight_init_std = 0.01 starves a ReLU net

The second defect decides whether the net can learn at all. Each affine layer multiplies the activation standard deviation by std_W × sqrt(fan_in), and ReLU then divides it by sqrt(2), because it keeps only the positive half of the distribution:

At weight_init_std = 0.01, the first layer multiplies the input std by 0.01 × sqrt(784) = 0.28 and each 50-unit hidden layer by 0.01 × sqrt(50) = 0.071. Four layers squeeze the logits into the 1e-5 range, the softmax outputs sit at an indistinguishable 0.1, and the loss refuses to fall below ln 10 = 2.303. Gradients shrink by the same factor in reverse, so the first layer's gradient is roughly 30,000 times smaller than the last layer's.

A starved net alone gives a stuck loss, not a NaN. The NaN needs the shared shift of cause 1 to amplify the first stray non-finite value into a batch-wide failure, which is why this run diverges within its first hundred iterations. The full arithmetic follows in the worked example.

Worked example: the 784-50-50-50-10 ReLU net at weight_init_std = 0.01

Inputs, all from the failing run: batch_size = 100; layers 784-50-50-50-10; ReLU hidden activations; weight_init_std = 0.01; MNIST pixels scaled by /255, which gives a pixel std of 0.31 (the standard MNIST statistics are a mean of 0.13 and a std of 0.31).

Rule: an affine layer multiplies the activation std by std_W × sqrt(fan_in); ReLU divides it by sqrt(2) = 1.414.

Arithmetic at initialization, forward pass:

• Layer 1: 0.01 × 0.31 × sqrt(784) = 0.01 × 0.31 × 28 = 0.0868; after ReLU, 0.0868 / 1.414 = 0.0614
• Layer 2: 0.01 × sqrt(50) × 0.0614 = 0.01 × 7.071 × 0.0614 = 0.00434; after ReLU, 0.00307
• Layer 3: 0.01 × 7.071 × 0.00307 = 2.17e-4; after ReLU, 1.54e-4
• Layer 4, the logits: 0.01 × 7.071 × 1.54e-4 = 1.09e-5

Result: logits with a std near 1.1e-5, outputs pinned at 0.1 per class, and a loss floor of ln 10 = 2.303.

The same net with He initialization, std = sqrt(2/n):

• Layer 1: sqrt(2/784) = 0.0505, giving a pre-activation std of 0.0505 × 0.31 × 28 = 0.438 and 0.310 after ReLU
• Hidden layers: sqrt(2/50) = 0.200, and each returns 0.200 × 7.071 × 0.310 = 0.438 before ReLU, 0.310 after
• Logits: std 0.438 — some 40,000 times larger than the 0.01 run and far inside the float64 exp range of ±709.78

Decision: a ReLU net deeper than two hidden layers gets std = sqrt(2/n), never 0.01. At 0.01 the run sits at the ln 10 floor and the first hundred iterations end in NaN through cause 1. With He init and the corrected softmax of the next section, this exact net trains below a loss of 0.3 within its first thousand batches at Momentum, lr = 0.01.

The fix, line by line

Two changes: the softmax works per row, and each weight matrix gets its own He standard deviation.

import numpy as np

def softmax(a: np.ndarray) -> np.ndarray:
    """Row-wise softmax for a (batch, classes) array."""
    c = np.max(a, axis=1, keepdims=True)      # one shift per row
    exp_a = np.exp(a - c)                     # every exponent <= 0
    return exp_a / np.sum(exp_a, axis=1, keepdims=True)

def he_std(fan_in: int) -> float:
    return float(np.sqrt(2.0 / fan_in))

# W1 = np.random.randn(784, 50) * he_std(784)   # std 0.0505
# W2 = np.random.randn(50, 50) * he_std(50)     # std 0.200
# W4 = np.random.randn(50, 10) * np.sqrt(1/50)  # softmax-fed layer, std 0.141

axis=1 with keepdims=True keeps the (100, 1) shape, so the subtract and the divide broadcast down the rows instead of collapsing into one number. The proof is the same test at a hostile scale:

scores = np.random.default_rng(1).normal(loc=0.0, scale=1000.0, size=(100, 10))
y = softmax(scores)
assert np.isfinite(y).all()              # no inf, no NaN
assert np.allclose(y.sum(axis=1), 1.0)   # every row a proper probability vector

Unshifted, logits of several thousand overflow float64 exp past 709.78; shifted, both asserts pass. Keep the delta = 1e-7 guard in the cross entropy — it is harmless once the softmax is correct.

Verify before the long run

Four checks, in order:

• Rows sum to 1: np.allclose(softmax(scores).sum(axis=1), 1.0) on a random (100, 10) batch.
• The loss stays finite for the first 100 iterations and is below 1.0 by then; below 0.3 by iteration 1000 at the failing run's settings (Momentum, lr = 0.01).
• Nothing raises while np.seterr(invalid='raise') is set.
• Held-out accuracy clears 95% after a few thousand batches; if it stalls near 10%, the labels are not one-hot or the loss is disconnected from the update.

A loss that is high but finite and stops falling is a different diagnosis — anomaly detection early stopping with high training loss covers that case.

When the same code runs in a production vision pipeline

The numerics follow the code onto the line, and three things change:

• Deployment pipelines often run float32, where the exp wall drops from 709.78 to 88.72. The per-row shift holds exponents at or below zero in both precisions, so the same fix carries over; only an unshifted softmax cares where the wall is.
• Check np.isfinite(logits).all() on a validation batch before export. A preprocessing mismatch between training and deployment shifts the score scale and can push borderline logits past the limit — training-versus-deployment confidence mismatches are diagnosed that way.
• Keep the loop vectorized. A Python for-loop over rows reopens the slow path that per-element loops cause throughout vision code; the per-pixel loop fix applies to training arrays just the same.

Where a job like this trains or infers inside a station, it runs on one of the vision controllers; filter their product pages by LAN port count, USB 3 ports and built-in light channels to match your camera count. Training throughput is a separate question: see why GPUs beat CPUs for deep learning training.

If it still goes NaN

Work forward one layer at a time and find the first non-finite value:

1. Print np.isfinite(a).all() immediately before the softmax each iteration; the first False places the defect upstream of the output layer.
2. Check the input batch: np.isfinite(x_batch).all() catches a NaN pixel from a divide-by-zero in preprocessing, which poisons every layer above it.
3. Check the targets: t.sum(axis=1) must equal 1 on every row. Labels that are not one-hot corrupt the (y - t) gradient rather than the softmax itself.
4. Halve the learning rate. The corrected softmax tolerates large logits, but momentum with an oversized step can still grow weights until something other than exp overflows.

More loop-level fixes sit on the systems topic page.

The decision

Three symptoms, three answers:

• RuntimeWarning at np.exp(a - c) with a batch-wide NaN: subtract and sum per row with axis=1, keepdims=True. The warning disappears because no exponent can exceed zero.
• Loss stuck near 2.303 on a ReLU net deeper than two hidden layers: replace weight_init_std = 0.01 with sqrt(2/fan_in) — 0.0505 at 784 inputs, 0.200 at 50.
• The same net moving to float32: keep the per-row shift and the 88.72 wall never matters; drop the shift and it does.

If this training runs beside a camera line rather than on a desk, match the machine to your camera count and frame rates through the contact page, and we will point you at the right vision controller.