Numerical stability: why log-sum-exp exists and how float32 silently lies

A machine learning model produces three logits: 1,000, 1,001, and 999. You apply softmax, expecting three probabilities that sum to one.

Instead, NumPy returns three NaNs.

There is no syntax error, no exception by default, and no obvious indication that your mathematical formula is wrong. Yet the calculation has failed.

Now consider a different problem. Two scores differ by exactly one point. After converting them to 32-bit floating-point numbers, they become identical. Softmax returns equal probabilities, confidently discarding the difference.

These failures reveal an uncomfortable truth about numerical computing:

A mathematically correct formula is not necessarily a numerically correct algorithm.

Floating-point arithmetic operates on a finite set of representable numbers. Intermediate results can overflow, underflow, lose precision, or accumulate rounding errors. Sometimes the failure produces NaN or infinity. More dangerously, sometimes it produces a plausible but incorrect number.

This tutorial explores how those failures happen, why log-sum-exp is essential, and how to write NumPy code that remains reliable under extreme inputs.

By the end, you will understand how to:

• Reproduce softmax overflow and catastrophic cancellation.

• Explain why float32 can silently destroy meaningful differences.

• Derive and implement the log-sum-exp trick.

• Build stable softmax, log-softmax, and cross-entropy calculations.

• Recognize cases where even a stable log-sum-exp implementation is insufficient.

• Test numerical algorithms against realistic precision failures.

Every demonstration uses NumPy and synthetic numerical data.

1. Set up a reproducible NumPy environment

The examples use the current NumPy 2.x API and work with NumPy 2.3 or later.

Use Python 3.11 or newer.

On Linux or macOS, create a virtual environment:

python -m venv .venv
source .venv/bin/activate
python -m pip install "numpy>=2.3,<3"

Verify the installation:

python -c "import numpy as np; print(np.__version__)"

The Python examples can be executed independently unless a section explicitly references a function defined earlier.

We will primarily compare two numerical formats:

• float32: 32-bit floating-point arithmetic, common in machine learning and GPU workloads.

• float64: 64-bit floating-point arithmetic, NumPy’s usual default for arrays constructed from Python floating-point literals.

Using float64 generally reduces rounding errors and expands the representable range, but it does not eliminate numerical instability.

That distinction matters throughout the tutorial.

2. Understand the limits of float32

Before diagnosing numerical failures, we need to understand what floating-point numbers actually store.

A conventional IEEE 754 binary32 value consists of:

• 1 sign bit.

• 8 exponent bits.

• 23 explicitly stored fraction bits, with 24 bits of significand precision for normal numbers.

Those 24 significant binary digits correspond to approximately 7 decimal digits of precision.

Importantly, float32 does not provide seven digits after the decimal point. It provides roughly seven significant decimal digits across the entire value.

A number near 1 has much finer absolute precision than a number near 100 million.

Inspect the actual machine limits

NumPy exposes floating-point properties through np.finfo.

import numpy as np

for dtype in (np.float32, np.float64):
    info = np.finfo(dtype)

    print(f"\n{dtype.__name__}")
    print("Machine epsilon:", info.eps)
    print("Largest finite:", info.max)
    print("Smallest positive normal:", info.smallest_normal)
    print("Smallest positive subnormal:", info.smallest_subnormal)

The important approximate values are:

Machine epsilon is the gap between 1 and the next larger representable floating-point value.

It is not a universal bound on absolute error.

As values grow, the gaps between adjacent representable numbers also grow.

When adding one does nothing

Consider the integer 16,777,216.

This is 2²⁴, the point beyond which float32 can no longer represent every consecutive integer.

import numpy as np

a = np.float32(16_777_216)
b = np.float32(16_777_217)

print("a:", a)
print("b:", b)
print("Equal:", a == b)
print("Difference:", b - a)

Output:

a: 1.6777216e+07
b: 1.6777216e+07
Equal: True
Difference: 0.0

The original integers differ by one, but the floating-point representations are identical.

Under the usual round-to-nearest, ties-to-even behavior, 16,777,217 rounds to 16,777,216 in float32.

The difference has been lost during conversion, before subtraction even begins.

This is a precision-loss problem, not an overflow problem.

Once two distinct inputs have been rounded to the same floating-point number, no subsequent numerical trick can recover their original difference.

3. Reproduce softmax overflow

Softmax converts a vector of real-valued scores into a probability distribution.

For a vector x, the probability corresponding to element i is:

softmax(xᵢ) = exp(xᵢ) / Σⱼ exp(xⱼ)

The output probabilities are nonnegative and sum to one in exact arithmetic.

Softmax is widely used in classification systems, attention mechanisms, and probabilistic models.

A direct NumPy implementation looks innocent:

import numpy as np

def naive_softmax(x):
    exponentials = np.exp(x)
    return exponentials / np.sum(exponentials)

scores = np.array(
    [1000.0, 1001.0, 999.0],
    dtype=np.float32,
)

with np.errstate(over="ignore", invalid="ignore"):
    probabilities = naive_softmax(scores)

print(probabilities)

Output:

[nan nan nan]

Why?

The maximum finite float32 value is approximately 3.4 × 10³⁸.

Taking the natural logarithm of that value gives approximately 88.72.

Consequently, exponentiating sufficiently large positive float32 values produces infinity.

For our logits:

• exp(999) overflows.

• exp(1000) overflows.

• exp(1001) overflows.

The denominator becomes infinity.

Every probability then involves infinity divided by infinity, which is undefined in floating-point arithmetic and produces NaN.

NumPy normally emits runtime warnings for these operations. The demonstration suppresses them so the failed result is easy to inspect.

Why switching to float64 is not a complete fix

Float64 can represent exponentially larger numbers than float32. However, exponentiating 1,000 still exceeds its finite range.

More generally, any finite floating-point format has a maximum representable exponentiation result.

Changing dtype moves the failure threshold. It does not remove the underlying problem.

The correct approach is to change the computation.

4. Derive the log-sum-exp trick

The log-sum-exp function is defined as:

LSE(x) = log(Σᵢ exp(xᵢ))

A direct implementation is:

import numpy as np

x = np.array([1000.0, 1001.0, 999.0])

with np.errstate(over="ignore"):
    result = np.log(np.sum(np.exp(x)))

print(result)

Output:

inf

Mathematically, the answer is finite.

The problem is that intermediate exponentials overflow before their logarithm can be evaluated.

The key identity

For any finite constant m:

LSE(x) = m + log(Σᵢ exp(xᵢ − m))

This follows from factoring exp(m) out of the sum:

Σᵢ exp(xᵢ) = exp(m) × Σᵢ exp(xᵢ − m)

Taking the logarithm produces the identity above.

The transformation is exact mathematically, but its numerical behavior depends critically on our choice of m.

Choose the maximum element

Set:

m = max(x)

Then every shifted value satisfies:

xᵢ − m ≤ 0

Therefore:

0 < exp(xᵢ − m) ≤ 1

At least one shifted value equals zero, so at least one exponential equals one.

For a nonempty finite vector, the shifted exponential sum lies between 1 and the number of elements.

This gives us two powerful guarantees:

1. The shifted exponentials cannot overflow.

2. The sum cannot underflow entirely to zero.

Individual very small exponentials may still underflow to zero, but they are negligible relative to the largest term. This is normally acceptable when evaluating the overall sum.

Implement the transformation

import numpy as np

x = np.array(
    [1000.0, 1001.0, 999.0],
    dtype=np.float64,
)

m = np.max(x)
shifted = x - m

result = m + np.log(np.sum(np.exp(shifted)))

print("Shifted logits:", shifted)
print("Log-sum-exp:", result)

Output:

Shifted logits: [-1.  0. -2.]
Log-sum-exp: 1001.4076059644444

The result is finite.

The computation remains stable because the exponential function only receives zero or negative inputs.

Notice an important conceptual shift: we did not approximate the original mathematical function with a different function. We rewrote the same function into a form that behaves better under floating-point arithmetic.

That is the essence of numerical stability.

5. Implement stable softmax

The same shifting transformation works for softmax.

Softmax is invariant to adding or subtracting a common constant from every input:

softmax(x) = softmax(x − m)

This follows because the common exp(m) factor cancels between the numerator and denominator.

Therefore, instead of exponentiating the original logits, exponentiate their maximum-shifted values.

import numpy as np

def stable_softmax(x):
    x = np.asarray(x, dtype=np.float64)

    shifted = x - np.max(x)
    weights = np.exp(shifted)

    return weights / np.sum(weights)

scores = np.array(
    [1000.0, 1001.0, 999.0],
    dtype=np.float32,
)

print(stable_softmax(scores))
print("Sum:", stable_softmax(scores).sum())

Output:

[0.24472847 0.66524096 0.09003057]
Sum: 0.9999999999999999

The sum differs from 1.0 only in the last binary digit. That is ordinary rounding, not a bug. Output above is from NumPy 2.5.4; older releases may print the scalar values in a different format.

We obtain valid probabilities despite the extremely large original logits.

The probabilities reflect relative differences rather than absolute magnitudes.

This is exactly what softmax is supposed to do.

One caveat: the displayed sum is one for this example. Floating-point rounding means that probability sums are generally only approximately one, not mathematically guaranteed to be exactly one.

Why the largest probability survives

If the scores are extremely different, some smaller probabilities may underflow to zero.

For example, a logit that is 2,000 below the maximum will have a shifted exponential too small to represent in float64.

That does not necessarily make the normalization wrong.

A tiny contribution can be safely discarded when it is below the representable precision of the desired result.

The important distinction is between:

• Losing a negligible contribution whose effect on the answer is insignificant.

• Losing meaningful information and producing an incorrect result.

Stable softmax prevents unnecessary overflow and avoids a zero denominator for finite, nonempty inputs.

It does not guarantee that every mathematically positive probability remains representable as a positive floating-point number.

6. Underflow: the opposite failure

Overflow is only half the story.

Underflow occurs when a nonzero mathematical result is too small to represent in the chosen floating-point format.

Consider two extremely small quantities represented by their natural logarithms.

import numpy as np

a = -800.0
b = -801.0

with np.errstate(divide="ignore"):
    naive = np.log(np.exp(a) + np.exp(b))

stable = np.logaddexp(a, b)

print("Naive:", naive)
print("Stable:", stable)

Output:

Naive: -inf
Stable: -799.6867383124818

The direct calculation fails because both exponentials underflow to zero in float64.

The logarithm of their sum becomes the logarithm of zero, producing negative infinity.

However, the logarithm of the true sum is finite.

NumPy already provides a stable primitive

The np.logaddexp function computes:

logaddexp(a, b) = log(exp(a) + exp(b))

It does so without naively exponentiating both original arguments.

For two finite values, the key transformation is:

logaddexp(a, b) = max(a, b) + log(1 + exp(−|a − b|))

Notice the similarity to the maximum-shifting strategy.

NumPy also supports reducing multiple elements with np.logaddexp.reduce:

import numpy as np

log_values = np.array(
    [-800.0, -801.0, -805.0, -900.0],
    dtype=np.float64,
)

combined = np.logaddexp.reduce(log_values)

print("Combined log value:", combined)

This is useful for combining extremely small probabilities while retaining them in logarithmic form.

For larger arrays, a vectorized maximum-shifted implementation can be preferable for performance and may have different rounding behavior. Reduction order matters, so the two approaches are not guaranteed to produce identical final bits.

Think in the log domain

Suppose you are calculating the likelihood of a sequence of independent events.

Multiplying hundreds of small probabilities can quickly underflow, even in float64.

Instead of multiplying probabilities, add their logarithms:

log(p₁ × p₂ × … × pₙ) = log(p₁) + log(p₂) + … + log(pₙ)

And when you need to add probabilities stored in logarithmic form, use log-sum-exp.

This pattern is fundamental to numerically robust probabilistic computing.

7. How float32 silently changes softmax predictions

Overflow often announces itself through infinity or NaN.

Precision loss is more subtle because the computation may finish successfully and return perfectly reasonable-looking probabilities.

Imagine an e-commerce ranking model assigning scores to two products.

The scores are:

• Product A: 100,000,000

• Product B: 100,000,001

The second product should have a higher softmax probability because its score is one point larger.

But what happens if those values are stored as float32?

import numpy as np

def softmax(x):
    shifted = x - np.max(x)
    weights = np.exp(shifted)
    return weights / np.sum(weights)

raw_scores = np.array(
    [100_000_000, 100_000_001],
    dtype=np.int64,
)

scores32 = raw_scores.astype(np.float32)
scores64 = raw_scores.astype(np.float64)

print("float32 scores:", scores32)
print("float32 spacing:", np.spacing(scores32[0]))

print("float32 probabilities:", softmax(scores32))
print("float64 probabilities:", softmax(scores64))

Output:

float32 scores: [1.e+08 1.e+08]
float32 spacing: 8.0
float32 probabilities: [0.5 0.5]
float64 probabilities: [0.26894142 0.73105858]

This example is particularly important.

The float32 computation does not return NaN. It does not overflow. The probabilities even sum to one.

Yet the result is wrong relative to the original integer scores.

What happened?

At approximately 100 million, adjacent float32 values are separated by 8.

The input difference of 1 is smaller than the local spacing between representable floating-point numbers.

Both scores round to the same value.

Softmax cannot distinguish them anymore.

The stable softmax implementation works correctly on the numbers it receives, but those numbers no longer contain the information needed to reproduce the original answer.

Numerical stability cannot recover information destroyed by earlier rounding.

This distinction is critical in machine learning pipelines.

A stable algorithm can still produce an incorrect result when its inputs were quantized too aggressively.

Practical implications

Pay particular attention when processing:

• Large timestamps or counters converted to floating-point values.

• Large accumulated scores with small differences.

• High-magnitude coordinates and sensor measurements.

• Model logits containing a large common offset.

• Financial amounts where exact decimal precision matters.

If the meaningful information lives in small differences between large values, float32 may be inappropriate.

Possible solutions include retaining float64, subtracting a common reference value before converting to float32, or representing inherently exact quantities as integers or decimals where appropriate.

Centering before conversion works only if the differences still exist in the higher-precision source.

If the source values have already been rounded to the same float32 value, that information is permanently lost.

8. Catastrophic cancellation: subtracting nearly equal numbers

Another major source of numerical error is catastrophic cancellation.

It occurs when subtracting two nearly equal floating-point values.

The exact difference may be very small, while the input values are much larger.

Rounding errors in those inputs can dominate the resulting difference.

A small but revealing example

Consider:

f(x) = √(1 + x²) − 1

For x = 0.0001, the true result is approximately 5 × 10⁻⁹.

Let’s evaluate it using float32.

import numpy as np

x = np.float32(1e-4)

direct = np.sqrt(np.float32(1) + x * x) - np.float32(1)

print("Direct result:", direct)

Output:

Direct result: 0.0

The mathematical answer is nonzero.

Why did NumPy return zero?

The term x² is approximately 10⁻⁸.

When added to 1 in float32, this contribution is too small to change the rounded result.

The square root therefore evaluates to 1.

Subtracting 1 gives exactly zero.

The final subtraction exposes the precision lost in the preceding operation.

Fix the algebra, not just the dtype

We can rationalize the expression:

√(1 + x²) − 1 = x² / (√(1 + x²) + 1)

This identity is mathematically exact.

However, the revised formula avoids subtracting two nearly equal numbers.

import numpy as np

x = np.float32(1e-4)
one = np.float32(1)

direct = np.sqrt(one + x * x) - one

stable = (x * x) / (np.sqrt(one + x * x) + one)

print("Direct:", direct)
print("Stable:", stable)

Output:

Direct: 0.0
Stable: 4.9999995e-09

The stable version preserves the meaningful result.

The sum inside the square root can still round to 1, but that rounding no longer eliminates the small numerator.

This is an important general technique:

When a calculation subtracts nearly equal expressions, look for an algebraically equivalent form that avoids the subtraction.

Other approaches include specialized numerical functions, compensated algorithms, and higher-precision arithmetic.

Merely changing float32 to float64 is not always enough.

9. Use specialized functions for values near zero

NumPy provides functions specifically designed to prevent precision loss in common expressions.

Two especially useful functions are:

• np.log1p(x), which computes log(1 + x).

• np.expm1(x), which computes exp(x) − 1.

Why do these functions exist when NumPy already provides logarithm and exponential functions?

Because their naive equivalents become inaccurate when x is very close to zero.

Reproduce the failure

import numpy as np

x = np.float64(1e-16)

print("Naive log:", np.log(1 + x))
print("Stable log:", np.log1p(x))

print("Naive exponential difference:", np.exp(x) - 1)
print("Stable exponential difference:", np.expm1(x))

Output:

Naive log: 0.0
Stable log: 1e-16
Naive exponential difference: 0.0
Stable exponential difference: 1e-16

The problem is the same general phenomenon as catastrophic cancellation.

In the logarithm example, adding a tiny number to one loses the tiny increment.

In the exponential example, exp(x) rounds to one, so subtracting one eliminates the meaningful result.

The specialized functions use numerical methods designed to retain accuracy near zero.

They are also useful when implementing other stable transformations.

For example, the two-argument log-sum-exp identity can use log1p to accurately evaluate its small correction term.

The broader lesson is that numerically robust software often relies on specialized mathematical primitives rather than mechanically translating textbook formulas.

10. Build production-oriented log-sum-exp and softmax functions

So far, our implementations have operated on one-dimensional vectors.

Real applications frequently process entire batches of logits.

A classification model may produce an array shaped like:

• Rows: examples.

• Columns: candidate classes.

We need to normalize each row independently.

We also need predictable dtype handling and explicit behavior for invalid inputs.

Implement reusable NumPy functions

The following implementation supports arbitrary reduction axes, promotes inputs to float64, and rejects nonfinite logits.

It deliberately requires finite inputs. Applications using negative infinity for masking need additional handling, discussed below.

import numpy as np


def _checked_logits(x, axis):
    z = np.asarray(x, dtype=np.float64)

    if z.ndim == 0 or z.size == 0:
        raise ValueError("Expected a nonempty array")

    if not np.isfinite(z).all():
        raise ValueError("Expected finite logits")

    return z


def logsumexp(x, axis=-1, keepdims=False):
    z = _checked_logits(x, axis)

    maximum = np.max(z, axis=axis, keepdims=True)

    shifted = z - maximum

    result = maximum + np.log(
        np.sum(
            np.exp(shifted),
            axis=axis,
            keepdims=True,
        )
    )

    if keepdims:
        return result

    return np.squeeze(result, axis=axis)


def softmax(x, axis=-1):
    z = _checked_logits(x, axis)

    shifted = z - np.max(
        z,
        axis=axis,
        keepdims=True,
    )

    weights = np.exp(shifted)

    denominator = np.sum(
        weights,
        axis=axis,
        keepdims=True,
    )

    return weights / denominator


def log_softmax(x, axis=-1):
    z = _checked_logits(x, axis)

    shifted = z - np.max(
        z,
        axis=axis,
        keepdims=True,
    )

    log_denominator = np.log(
        np.sum(
            np.exp(shifted),
            axis=axis,
            keepdims=True,
        )
    )

    return shifted - log_denominator


scores = np.array(
    [
        [1000.0, 1001.0, 999.0],
        [-1000.0, -1001.0, -999.0],
    ],
    dtype=np.float32,
)

print("Log-sum-exp:")
print(logsumexp(scores, axis=1))

print("\nSoftmax:")
print(softmax(scores, axis=1))

print("\nLog-softmax:")
print(log_softmax(scores, axis=1))

All three functions produce finite results for the demonstrated inputs.

Why we promote to float64

The conversion to float64 provides greater precision for intermediate calculations.

In CPU-based NumPy workloads, this is often a reasonable default when correctness matters more than memory or throughput.

However, it is not universally optimal.

For large models, float32 may be substantially more memory-efficient, especially on accelerators.

Production systems should choose accumulation and output dtypes deliberately, based on numerical requirements and performance constraints.

Also remember that converting an existing float32 array to float64 cannot restore information that was lost when the array was originally created.

Why we use keepdims=True

The maximum and denominator are calculated along an axis while retaining that axis as a dimension of size one.

This allows NumPy broadcasting to subtract or divide the correct value across each row.

Without proper axis handling, a softmax implementation can accidentally normalize an entire batch rather than each example.

That error may produce plausible numbers while completely changing the meaning of the result.

Handling infinity and masked elements

Many attention implementations use negative infinity to mark positions that must receive zero probability.

Our reference implementation rejects these values intentionally, making its contract simple and explicit.

A masked softmax needs additional rules:

• Individual negative-infinity logits can represent excluded positions.

• An entirely masked row has no valid probability distribution and requires an explicit policy.

• NaN inputs should not pass unnoticed.

• Positive infinity requires a defined convention rather than blindly subtracting the maximum.

For an all-masked row, returning uniform probabilities would generally be misleading.

Whether the correct behavior is an exception, a special mask result, or a framework-specific convention depends on the application.

The important principle is to define these cases rather than accidentally inheriting NaNs from undefined arithmetic.

11. A subtle trap: stable log-sum-exp can still produce unstable normalization

Here is a surprising edge case.

We have learned to compute log-sum-exp with maximum shifting.

It is tempting to write softmax as:

softmax(xᵢ) = exp(xᵢ − LSE(x))

That identity is mathematically correct.

Unfortunately, implementing it literally can still cause a precision problem.

Consider two enormous, identical logits

import numpy as np

scores = np.array(
    [1e20, 1e20],
    dtype=np.float64,
)

maximum = np.max(scores)

lse = maximum + np.log(
    np.sum(np.exp(scores - maximum))
)

via_lse = np.exp(scores - lse)

shifted_weights = np.exp(scores - maximum)
direct_normalization = shifted_weights / shifted_weights.sum()

print("Log-sum-exp:", lse)
print("Via log-sum-exp:", via_lse)
print("Direct shifted normalization:", direct_normalization)

Output:

Log-sum-exp: 1e+20
Via log-sum-exp: [1. 1.]
Direct shifted normalization: [0.5 0.5]

The first probability calculation sums to 2, not 1.

How is that possible when log-sum-exp itself was computed using a stable transformation?

The missing correction term

The exact mathematical answer is:

LSE([10²⁰, 10²⁰]) = 10²⁰ + log(2)

But at 10²⁰, float64 does not have enough absolute precision to retain an addition as small as log(2).

The computed log-sum-exp rounds back to 10²⁰.

Subtracting that rounded result from each original logit gives zero.

Exponentiating produces two ones.

The maximum-shifted normalization avoids reconstructing and subsequently subtracting this enormous log-sum-exp result.

It keeps the small but essential normalization correction separate from the large common offset.

The stable way to compute an intermediate value is not automatically the stable way to use that value in every downstream expression.

This is why the log_softmax implementation above computes its logarithmic denominator in shifted coordinates.

It avoids subtracting a rounded, large-magnitude log-sum-exp value from the original logits.

12. Stable cross-entropy without taking the log of zero

A common classification loss is negative log-likelihood.

For the correct class y, the loss is:

Loss = −log(softmax(x)ᵧ)

A direct implementation first calculates probabilities, selects the correct-class probability, and takes its logarithm.

But what happens when that probability underflows to zero?

The result becomes infinity, even if the true negative log-likelihood is finite.

Reproduce the problem

The following example uses a three-class product-category classifier.

The first example is intentionally assigned an extremely low score for its correct category.

This block reuses the softmax and log_softmax functions from Section 10.

import numpy as np

logits = np.array(
    [
        [1000.0, 0.0, -1000.0],
        [2.0, 1.0, 0.0],
    ],
    dtype=np.float64,
)

targets = np.array([2, 1])

row_indices = np.arange(targets.size)

with np.errstate(divide="ignore"):
    naive_losses = -np.log(
        softmax(logits)[row_indices, targets]
    )

stable_losses = -log_softmax(logits)[
    row_indices,
    targets,
]

print("Naive losses:", naive_losses)
print("Stable losses:", stable_losses)
print("Mean stable loss:", stable_losses.mean())

The naive calculation returns infinity for the first example.

The stable version returns a finite loss of 2,000.

Why?

The correct class is 2,000 logit units below the highest-scoring class.

Its probability is too small to represent as a positive float64 number.

But its logarithm remains representable.

Computing the logarithmic probability directly avoids unnecessary underflow.

Why this matters during model training

A nonfinite loss can contaminate downstream calculations.

For example:

• Loss monitoring may report infinity.

• Custom gradient calculations may become nonfinite.

• Model comparisons may fail.

• Training-control logic may react incorrectly.

• Invalid values may propagate through later computations.

A stable loss implementation cannot fix every possible training instability, but it prevents avoidable failures caused by materializing tiny probabilities.

For classification from logits, the right approach is to use a numerically stable log-softmax or cross-entropy implementation instead of manually chaining naive softmax and logarithm operations.

13. Reduction order is another hidden source of error

Even operations as simple as summation can produce different results depending on evaluation order.

In exact arithmetic:

(a + b) + c = a + (b + c)

Floating-point arithmetic does not generally preserve this associativity.

Demonstrate with three numbers

import numpy as np

values = np.array(
    [1e8, 1.0, -1e8],
    dtype=np.float32,
)

sum32 = np.sum(values, dtype=np.float32)
sum64 = np.sum(values, dtype=np.float64)

print("float32 accumulation:", sum32)
print("float64 accumulation:", sum64)

Output:

float32 accumulation: 0.0
float64 accumulation: 1.0

The mathematically exact sum of the represented inputs is one.

The float32 accumulation loses the small contribution when it is combined with a much larger value.

Using float64 for the accumulator preserves the contribution in this example.

This is an important distinction between storage precision and accumulation precision.

You may store data in float32 but perform selected reductions in float64.

That can improve accuracy without requiring every array in the application to use twice as much memory.

However, there is a limit: wider accumulation cannot repair rounding errors that already occurred when individual input values were created.

Why parallel processing complicates this

Modern numerical libraries often use parallel reductions.

Instead of summing elements sequentially, an implementation may divide the array into smaller blocks and combine partial sums.

Different execution orders produce different patterns of rounding.

Consequently, numerically valid implementations can differ in their least significant bits.

This is one reason exact bitwise equality is often too strict when testing floating-point reductions across implementations or execution environments.

Use explicit error tolerances and evaluate accuracy against the numerical requirements of the application.

14. Cherry on the cake: a real 2026 database precision bug

Numerical precision problems are not limited to educational examples or machine learning experiments.

In July 2026, the ClickHouse project merged a fix for a floating-point parsing inconsistency that could cause the same decimal input to be stored as different binary floating-point values depending on how it entered the database.

The underlying issue was documented in ClickHouse issue #60146.

A reported example involved the decimal value represented by 1.234567891209 × 10¹².

Under the affected behavior, one parsing path produced:

1234567891209

Another produced:

1234567891209.0002

The results differed by approximately one unit in the last place, or ULP, at that magnitude.

That distinction is small in relative terms but can matter for exact comparisons, repeated data processing, or consistency checks.

What caused the discrepancy?

ClickHouse supported a precise floating-point parsing option, but historically that option did not apply consistently to every parsing path.

Some conversions used the precise parser, while certain input formats and numeric literals still used a faster algorithm that could produce a slightly different approximation.

The issue was not that IEEE 754 arithmetic itself was malfunctioning.

The problem was inconsistent conversion of decimal text into binary floating-point numbers.

The merged fix, pull request #109086, completed on July 6, 2026, made precise parsing the default and extended its use across input formats including CSV, TSV, JSON, and VALUES.

The change was included in ClickHouse 26.7.

This was a documented correctness bug in production-grade database software, not a claimed outage or an invented incident.

What should engineers learn from it?

The subtle part is that the bug occurred at the boundary between text and numbers.

Even before running calculations, different input paths could produce different floating-point values.

This matters because many systems assume that parsing the same decimal representation yields identical stored data.

For high-reliability numerical pipelines:

• Test decimal-to-binary conversions across ingestion paths.

• Validate round-trip behavior for serialized numerical data.

• Avoid assuming display formatting proves bitwise equality.

• Define whether values require approximate numerical equality or exact representation.

• Test values close to representational boundaries.

A numerical pipeline begins at data ingestion, not at the first mathematical operation.

15. Build a numerical stability test suite

A numerically robust implementation needs more than a few tests on small random values.

Ordinary inputs often conceal numerical weaknesses.

Testing should deliberately target the boundary conditions where floating-point arithmetic behaves differently from ideal real-number arithmetic.

Important test categories

Test probability normalization

After defining the stable functions from Section 10, run:

import numpy as np

rng = np.random.default_rng(2026)

scores = rng.normal(size=(64, 10))

probabilities = softmax(scores, axis=1)

np.testing.assert_allclose(
    probabilities.sum(axis=1),
    np.ones(scores.shape[0]),
    rtol=1e-14,
    atol=1e-14,
)

assert np.isfinite(probabilities).all()
assert np.all(probabilities >= 0)

print("Normalization tests passed")

These tests verify that:

• Every row approximately sums to one.

• Every output is finite.

• No probability is negative.

The tolerances are intentionally appropriate for float64 results on this well-conditioned test input.

They should not be copied blindly into unrelated applications.

Test shift invariance

Softmax should produce the same probabilities when a common constant is added to the logits, provided the transformed inputs retain sufficient precision.

import numpy as np

rng = np.random.default_rng(2026)

scores = rng.normal(size=(64, 10))

original = softmax(scores, axis=1)
shifted = softmax(scores + 1000.0, axis=1)

np.testing.assert_allclose(
    original,
    shifted,
    rtol=1e-12,
    atol=1e-12,
)

print("Shift invariance test passed")

This is a useful metamorphic test.

Rather than requiring a separately calculated expected probability for every input, it checks a mathematical property that the function should preserve.

Such tests can detect implementation errors that ordinary example-based tests miss.

Note that shift invariance is a mathematical property. If the shift is so enormous that adding it destroys differences between finite-precision logits, the input arrays are no longer equivalent at the machine level.

A test should not demand precision the chosen dtype cannot represent.

Test the conditions that actually break algorithms

Randomized tests are useful, but random sampling rarely reaches the most dangerous floating-point boundaries.

Add explicit tests for extreme magnitudes, underflow regions, large common offsets, cancellation-prone expressions, and dtype conversions.

Use high-precision references or analytically known answers where practical.

The test suite should answer two separate questions:

1. Does the implementation produce finite, structurally valid outputs?

2. Are those outputs sufficiently accurate for the intended task?

A probability vector summing to one does not prove that its individual probabilities are correct.

The earlier float32 score-collision example demonstrates exactly why.

16. Practical rules for stable numerical computing

The examples lead to several principles worth applying across machine learning, data engineering, and scientific computing.

Avoid dangerous intermediate expressions

If a formula contains exponentials, logarithms, subtraction of nearby values, or division by very small quantities, investigate whether a stable reformulation exists.

For example:

• Replace naive log-sum-exp with a maximum-shifted computation.

• Replace naive softmax with shifted exponentials followed by normalization.

• Use log-softmax for log-likelihood calculations.

• Use np.log1p for log(1 + x) when x is small.

• Use np.expm1 for exp(x) − 1 when x is small.

• Consider algebraic rationalization when subtraction causes cancellation.

Choose dtype based on the information that must survive

Float32 is an excellent practical format for many machine learning workloads.

But suitability depends on the magnitude and required resolution of the values being represented.

Use float64 where the greater precision is justified, particularly for sensitive reductions or computations involving small differences between large values.

For exact decimal accounting, an appropriate decimal or scaled-integer representation is generally more suitable than binary floating-point arithmetic.

Track precision across the entire pipeline

The final computation may be numerically stable while the data has already been damaged upstream.

Examine:

• Initial parsing and ingestion.

• Numeric dtype conversions.

• Intermediate storage and serialization.

• Arithmetic and reduction operations.

• Output formatting and comparison semantics.

A stable numerical kernel is only one part of a reliable system.

Treat warnings as useful diagnostics

NumPy runtime warnings can reveal overflow, underflow, invalid operations, and division by zero.

During development, make suspicious operations explicit and investigate warnings rather than suppressing them globally.

For example, NumPy’s np.errstate can temporarily configure floating-point error behavior for a particular block of code.

Warnings do not catch ordinary finite rounding errors, however. Precision loss can occur without producing any warning.

Do not confuse stability with exactness

A stable algorithm is designed to avoid unnecessary amplification of numerical errors.

It is not an arbitrary-precision arithmetic system.

Stable softmax can still return zero for a probability too small to represent.

Stable log-sum-exp can still round away a small correction when the final result is extremely large.

Float64 still has finite precision.

The correct objective is not perfect arithmetic, which ordinary floating-point formats cannot provide. It is sufficient numerical accuracy under a clearly defined operating range.

17. Final mental model: four ways a number can betray you

When diagnosing numerical problems, distinguish these four failure mechanisms.

These mechanisms can also interact.

For example, computing a log-likelihood using naive softmax might first overflow, then produce NaN, and finally contaminate a downstream loss calculation.

A numerical problem often appears far from the operation that originally caused it.

That is why diagnosing the first precision-losing transformation is more useful than simply checking where NaN eventually appears.

Further reading

The following current references provide implementation details and numerical contracts:

• NumPy: Floating-point machine limits — inspect dtype precision, spacing, and representable ranges.

• NumPy: logaddexp — stable addition of exponentials in logarithmic form.

• NumPy: log1p — accurate computation of log(1 + x).

• NumPy: expm1 — accurate computation of exp(x) − 1.

• ClickHouse precision parsing fix, July 2026 — a concrete example of how numerical correctness can depend on data-ingestion behavior.

Conclusion: make numerical correctness part of code quality

Floating-point arithmetic is not broken. It is finite, approximate arithmetic with precisely defined representational limits.

The real problem is assuming that mathematical identities remain equally accurate under every computational arrangement.

Log-sum-exp exists because computing the logarithm of a sum of exponentials directly is often unnecessarily dangerous. Maximum shifting converts the same mathematical expression into a substantially more reliable computation.

But the deeper lesson extends beyond exponential functions.

A model can produce incorrect probabilities without emitting a warning. A subtraction can erase a meaningful result. A database can parse the same numerical text differently through separate input paths. A correctly implemented algorithm can receive inputs that have already lost critical precision.

Reliable numerical software requires deliberate choices about representation, arithmetic, and validation.

Your next step: take one existing NumPy workflow that contains exponentials, logarithms, probability normalization, or large-value reductions. Test it using extreme positive and negative values, nearly equal inputs, and both float32 and float64. Replace unstable formulas, add assertions around the expected numerical properties, and document the precision guarantees your application actually requires.

The goal is not merely to make NaNs disappear. It is to ensure that a plausible-looking numerical result is also a trustworthy one.