12.3 Condition Numbers and Numerical Stability: Preventing Error Amplification

Your model's loss suddenly becomes NaN. The math on paper doesn't involve division by zero, and the input data isn't infinite. The issue lies in the computation path: an exponential term overflows, two nearly equal large numbers undergo catastrophic cancellation, and a tiny input error is ultimately amplified into a completely different result.

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. Condition Description Problems, Stability of Algorithms

If the input $x$ is perturbed slightly by $\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 given matrix norm is:

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

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.

The numerical stability of an algorithm asks 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 $\Delta x$ such that $\hat y = f(x + \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.

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.

4. Stable softmax and log-sum-exp

Direct computation of

$$ \operatorname{softmax}(x_i)=\frac{e^{x_i}}{\sum_j e^{x_j}} $$

may suffer from overflow due to large positive exponents. Leveraging the translation invariance of softmax, let $m = \max_i x_i$:

$$ \operatorname{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 May not be empty")
    if not all(math.isfinite(value) for value in values):
        raise ValueError("The example accepts only limited input")

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

print(softmax([1000.0, 1001.0, 1002.0]))

The corresponding log-sum-exp is:

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

Subtracting the maximum value ensures that the largest exponent becomes 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.

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 high-precision 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.

import math

values = [1e16, 1.0, -1e16]
print(sum(values))        # Common results:0.0
print(math.fsum(values))  # 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\times2$ 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.

The Developer Workshop's core content concludes here: logic helps you articulate propositions, discrete structures help you organize relationships, probability enables you to handle uncertainty, linear algebra and calculus empower you to model and optimize, while information theory and numerical analysis reveal the fundamental limits of representation and computation. The story of the workshop and its mentors is not yet over; after full technical validation across the system, these insights will be reconnected to the characters' decisions, failures, and next steps.