12.3 Condition Numbers and Numerical Stability: Preventing Error Amplification

Finite inputs can still produce a NaN loss: direct softmax can overflow its exponentials, then form inf/inf. Cancellation between nearby values is a separate failure mode, often producing an inaccurate finite result first. Inspect the computation path as well as the formula.

Numerical analysis must distinguish two key questions: whether the problem itself is sensitive to input changes, and whether the algorithm further amplifies that sensitivity.

Learning Objectives

• Distinguish forward error, backward error, and problem conditioning;
• Identify catastrophic cancellation, overflow, underflow, and accumulated error;
• Implement numerically stable softmax, log-sum-exp, and summation operations;
• Design boundary tests and error checks for numerical code.

1. Conditioning Describes the Problem; Stability Describes the Algorithm

If the input $x$ is perturbed slightly by $\mathrm{\Delta }x$, the output $f(x)$ may change only slightly or dramatically. The condition number quantifies how sensitive a problem is to input perturbations.

For a linear system $A\mathbf{x}=\mathbf{b}$ with an invertible matrix, the condition number under a corresponding induced matrix norm is:

$$\kappa (A)=\Vert A\Vert \Vert A^{-1}\Vert .$$

A large condition number means that small relative errors in the input or due to rounding can be significantly amplified in the solution. A condition number of infinity corresponds to a singular matrix. The condition number depends on the chosen norm and does not provide an exact measure of the actual error in any individual computation.

Here we assess backward stability by asking whether the computed result is close to the exact solution of a slightly perturbed input. An ill-conditioned problem cannot be made insensitive by using a stable algorithm, but an unstable algorithm can turn a well-conditioned problem into one that produces poor results.

2. Forward and Backward Error

Let the exact solution be $y=f(x)$, and the computed result be $\hat{y}$:

• The forward error compares $\hat{y}$ with $y$;
• The backward error seeks a small perturbation $\mathrm{\Delta }x$ such that $\hat{y}=f(x+\mathrm{\Delta }x)$.

Backward stability means that the algorithm's output can be interpreted as the exact result of a slightly perturbed input. If the problem is well-conditioned, small backward errors typically lead to small forward errors; if the problem is ill-conditioned, even backward stability may result in large forward errors.

This explains why simply printing a few decimal digits of the result is insufficient. You must also examine residuals, input scales, and condition numbers. For example, in linear systems, one can check:

$$\mathbf{r}=\mathbf{b}-A\hat{\mathbf{x}}.$$

A small residual is a significant indicator, but in ill-conditioned systems, a small residual does not necessarily imply that the solution is close to the true parameters.

This NumPy 2.x diagonal example changes the right-hand side by only about 1e-8 in relative norm. The second solution component changes from 1 to 2, yet its residual against the original equation is only 1e-8:

import numpy as np

A = np.diag([1.0, 1e-8])
b = np.array([1.0, 1e-8])
perturbed_b = b + np.array([0.0, 1e-8])
x = np.linalg.solve(A, b)
changed_x = np.linalg.solve(A, perturbed_b)
relative_input = np.linalg.norm(perturbed_b - b) / np.linalg.norm(b)
residual = np.linalg.norm(b - A @ changed_x)
print(f"{np.linalg.cond(A):.0e} {relative_input:.0e} {residual:.0e}")
print(x, changed_x)
# 1e+08 1e-08 1e-08
# [1. 1.] [1. 2.]

The second equation itself has a large relative perturbation; it looks small only against the whole vector norm. Normwise and componentwise error models answer different questions, and unit scaling can affect the assessment. A bare residual is insufficient evidence of parameter accuracy.

3. Catastrophic Cancellation

When subtracting two nearly identical floating-point numbers, the significant digits cancel out, and the existing rounding errors in the input are amplified relative to the result.

A typical example is evaluating the expression:

$$\sqrt{1+x}-1$$

for small values of $x$. By rationalizing the numerator:

$$\sqrt{1+x}-1=\frac{x}{\sqrt{1+x}+1}.$$

we obtain a mathematically equivalent form that avoids subtracting two numbers close to 1.

import math

def unstable(x: float) -> float:
    return math.sqrt(1.0 + x) - 1.0

def stable(x: float) -> float:
    return x / (math.sqrt(1.0 + x) + 1.0)

for x in [1e-8, 1e-12, 1e-16]:
    print(x, unstable(x), stable(x))

"The rule of avoiding large minus small" is not accurate; the real danger lies in subtracting two nearly equal numbers, especially when the result is much smaller than the operands.

Both real-valued formulas require $x\ge -1$. At x=1e-16, binary64 rounds 1+x to 1, so direct evaluation returns zero. The reformulation retains x in the numerator and gives about 5e-17. Subtraction itself can sometimes be exact; catastrophic cancellation exposes errors already present in its operands.

4. Stable softmax and log-sum-exp

Direct computation of

$$\mathrm{softmax}(x_i)=\frac{e^{x_i}}{\sum _j{e}^{x_j}}$$

may suffer from overflow due to large positive exponents. Using translation invariance of softmax, let $m=\underset{i}{max}{x}_i$:

$$\mathrm{softmax}(x_i)=\frac{e^{x_i-m}}{\sum _j{e}^{x_j-m}}.$$

import math

def softmax(values: list[float]) -> list[float]:
    if not values:
        raise ValueError("values must not be empty")
    if not all(math.isfinite(value) for value in values):
        raise ValueError("only finite inputs are accepted")

    maximum = max(values)
    shifted = [math.exp(value - maximum) for value in values]
    denominator = math.fsum(shifted)
    return [value / denominator for value in shifted]

print([round(p, 6) for p in softmax([1000.0, 1001.0, 1002.0])])
# [0.090031, 0.244728, 0.665241]

The corresponding log-sum-exp is:

$$\log \sum _i{e}^{x_i}=m+\log \sum _i{e}^{x_i-m}.$$

Subtracting the maximum makes the largest exponential argument zero, giving an exponential value of 1, preventing positive overflow. Very small terms may still underflow to zero, but their relative contribution to the sum is inherently negligible given the current precision; whether this is acceptable depends on the required accuracy.

Continuing in Python 3.9+, we can implement log-sum-exp as a complete function. This version rejects empty and nonfinite inputs explicitly; other libraries can define different special-value semantics:

def logsumexp(values: list[float]) -> float:
    if not values or not all(math.isfinite(x) for x in values):
        raise ValueError("require a nonempty list of finite values")
    maximum = max(values)
    return maximum + math.log(math.fsum(math.exp(x - maximum) for x in values))

print(f"{logsumexp([1000.0, 1001.0, 1002.0]):.6f}")  # 1002.407606
try:
    math.log(math.fsum(math.exp(x) for x in [1000.0, 1001.0, 1002.0]))
except OverflowError:
    print("direct exponentiation overflowed")

For nonempty finite input, at least one term is 1, so the softmax denominator cannot underflow entirely to zero. Individual tiny probabilities can still become zero. For log probabilities, retain the shifted calculation (x_i-m) - log(sum_j exp(x_j-m)). This avoids both log(0) and subtracting a very large common offset after it has already swallowed the small logarithmic correction.

5. Cumulative Sum Also Requires Algorithms

When summing a large number of values spanning several orders of magnitude, small numbers may be lost due to rounding. Improved approaches include:

• Starting the summation with the smallest absolute values;
• Pairwise summation, which slows the growth of error with the number of terms;
• Kahan-style compensated summation, which retains some of the lost low-order bits;
• Using accurate summation implementations provided by the language or library, such as Python math.fsum.

These techniques improve numerical accuracy, but do not guarantee exact results for all inputs. In parallel systems, trade-offs must be made between speed, determinism, and reproducibility.

Python 3.12 improved built-in floating-point sum, so the older sum(values)==0 result is not a portable expectation. The example uses explicit sequential addition to demonstrate loss and separately prints the current built-in result. math.fsum tracks multiple partial sums for accuracy; it is not arbitrary-precision arithmetic. See the Python3.12 change and fsum documentation.

import math
import sys

values = [1e16, 1.0, -1e16]
sequential = 0.0
for value in values:
    sequential += value
print(sequential)         # 0.0
print(math.fsum(values))  # 1.0
print(sys.version_info[:2], sum(values))  # Observed with CPython 3.12: (3, 12) 1.0

6. When Choosing Formulas, First Consider Boundaries

Numerically stable implementations often rely on specialized functions:

• For small $x$, compute $\log (1+x)$ using log1p(x);
• For small $x$, compute $e^x-1$ using expm1(x);
• Compute $\sqrt{x^2+y^2}$ using hypot(x, y) to reduce risks of intermediate overflow or underflow;
• Solve linear systems via decomposition and solve, avoiding explicit inversion;
• For probability losses, compute directly from logits rather than first rounding to probabilities.

Mathematical equivalence does not guarantee floating-point execution equivalence. When designing algorithms, consider the range of intermediate values, cancellation, continuity of branches, and the number of rounding operations.

7. Checklist for Validating Numerical Code

1. Establish small-scale benchmarks using analytical solutions or high-precision results.
2. Test edge cases such as zero, very small or very large values, values close to each other, and sign changes.
3. Clearly define whether the code accepts NaN, infinity, or subnormal numbers.
4. Evaluate absolute error, relative error, residual, or conservation quantities.
5. Vary computation order, batch size, and parallelism to assess sensitivity.
6. Record convergence history for iterative algorithms, rather than the final output alone.
7. Log data types, rounding rules, library versions, and hardware assumptions.

Error thresholds should be tied to the scale of the problem. Asserting that "values agree to six decimal places" may be too lenient near zero and too strict for large numbers.

Common Misconceptions

• Correct formulas imply correct implementation: Intermediate computations may still overflow, underflow, or cancel out.
• High precision can fix any problem: It delays the exposure of errors but does not correct ill-posed models or flawed algorithms.
• Small residuals mean accurate parameters: In ill-conditioned problems, residuals and parameter accuracy can be entirely decoupled.
• A stable algorithm produces exact answers: Stability controls error propagation, but does not eliminate input or rounding errors.

Exercise

1. Compare the relative error of the two formulas sqrt(1+x)-1 at different values of $x$.
2. Implement a stable version of logsumexp and compare its output against the direct formula for large input values.
3. Construct a $2\times 2$ linear system with a large condition number, perturb the vector $b$ slightly, and compare the resulting solutions.
4. Design two modes for a batch aggregation service: "fast non-deterministic sum" and "reproducible sum".

Summary

The condition number describes a problem's inherent sensitivity to perturbations, while stability indicates whether the algorithm itself amplifies errors. To transform mathematical formulas into reliable software, essential steps include identifying error-reducing techniques, controlling intermediate value ranges, selecting stable functions, and validating responses to extreme inputs.

• LAPACK: Error Bounds for Linear Equation Solving