Numerical Stability and Conditioning

The previous chapter showed that every floating-point operation carries a tiny error. This chapter asks the deeper question: what happens when you chain millions of them together?

The answer depends on two independent properties:

• Conditioning is a property of the problem. A problem is well-conditioned if small input changes cause small output changes. Ill-conditioned problems amplify error no matter how cleverly you compute.
• Stability is a property of the algorithm. A stable algorithm produces an answer that is the correct answer to a slightly perturbed problem. An unstable algorithm produces garbage even on well-conditioned problems.

You always need to ask both questions.

Conditioning, intuitively

If $f$ is a function and you evaluate $f(x + \delta x)$ for small $\delta x$, how big is $f(x + \delta x) - f(x)$? The ratio of relative output change to relative input change is the condition number of $f$ at $x$.

For evaluating $f$:

$\kappa_f(x) \approx \left| \frac{x f'(x)}{f(x)} \right|$

For a linear system $A x = b$:

$\kappa(A) = \|A\| \, \|A^{-1}\|$

A condition number near $1$ is excellent. A condition number near $1/\epsilon \approx 10^{16}$ means no algorithm can give you more than a couple of correct bits.

Near $x = \pi$, $\sin(x)$ is nearly zero, so the relative output change for a tiny input change is enormous. Evaluating $\sin$ near a root is intrinsically ill-conditioned — there is no clever algorithm that fixes this.

Stability, intuitively

An algorithm is backward stable if the result it returns is the exact answer to a slightly different problem. The standard example is Gaussian elimination with partial pivoting: the answer you get to $Ax = b$ is the exact answer to $(A + E)x = b$, with $\|E\| / \|A\| = O(\epsilon)$. That is the best you could possibly hope for at finite precision.

An algorithm is unstable if it can produce arbitrarily wrong answers on well-conditioned problems. The classic example: solving a $2\times2$ linear system with Cramer's rule on a matrix whose determinant is near zero but well-conditioned.

Let us watch both happen.

The condition number $\kappa(A) \approx 4 \times 10^{13}$ tells us we have lost roughly 13 digits of precision no matter what. solve recovers the remaining 3 digits cleanly; inv is acceptable but worse; cramer is the worst of the three because it does extra cancellation. None of them can do better than $\kappa$ allows.

A canonical instability: forward vs backward recursion

Consider the integrals

$I_n = \int_0^1 \frac{x^n}{x + 10} \, dx$

These satisfy a tidy recurrence $I_n = \frac{1}{n} - 10 \, I_{n-1}$ with $I_0 = \ln(1.1)$. The forward recurrence is numerically unstable — every step multiplies the error by $-10$ — even though the integrals themselves are perfectly well-defined small positive numbers.

The recursion produces negative values, then enormous oscillations — all wrong. The integrand is positive and the integral has to be less than $1/(10n)$. The error introduced at step 1 is multiplied by 10 at every later step.

The cure: run the recurrence backwards, where the multiplier is $1/10$ instead of $10$. Start from a guess at large $n$ — even a wildly wrong guess — and iterate downward. The error shrinks by $10$ per step.

Same equation. Same data. Two algorithms; one is useless, the other is correct to all 16 digits. This is the entire content of "numerical stability."

Norms, error bounds, and "how wrong can it be?"

In numerical linear algebra, we measure error with norms. The most useful identities are:

For $A x = b$, the relative error in $x$ is bounded by:

$$\frac{\lVert \hat x - x \rVert}{\lVert x \rVert} \;\le\; \kappa(A) \left( \frac{\lVert \delta A \rVert}{\lVert A \rVert} + \frac{\lVert \delta b \rVert}{\lVert b \rVert} \right)$$

A condition number of $10^{14}$ and a perturbation of $10^{-16}$ gives a relative output error of $10^{-2}$ — two correct digits. That number is your truth ceiling, set by the problem, not the algorithm.

Watch how the relative error scales with $\kappa$. Once $\kappa(H) > 10^{16}$, you have lost all your precision.

Designing stable algorithms

A handful of rules turn most numerical code from fragile to robust.

1. Avoid subtracting nearly equal numbers. Rewrite the formula if you can (the 1 - cos x example), or compute in higher precision.
2. Sum in increasing order of magnitude, or use a compensated sum like Kahan summation, when adding many tiny contributions.
3. Prefer factorizations over inverses. np.linalg.solve(A, b) beats np.linalg.inv(A) @ b every time.
4. Scale your inputs. If one column of your matrix has values in $[10^{-9}]$ and another in $[10^{9}]$, condition number explodes for no real reason; rescale them.
5. Pick algorithms with proven backward error bounds. LAPACK, SciPy, Eigen, and PETSc routines all come with these. Custom code usually does not.

Here is a beautiful illustration of point 2 — Kahan summation recovers most of the digits a naive sum throws away.

NumPy's .sum() uses pairwise summation, which is much better than a naive left-to-right loop and almost as good as Kahan, for a fraction of the cost. Trust the library; reach for Kahan only when you really need it.

Check your understanding

QuestionSelect one

A problem has condition number $\kappa \approx 10^{12}$. You compute the answer with a backward-stable algorithm at 64-bit precision. Roughly how many correct decimal digits should you expect?

All 16

About 12

About 4

Zero — the problem is unsolvable

QuestionSelect one

Which statement most cleanly distinguishes conditioning from stability?

They mean the same thing

Conditioning is a property of the problem; stability is a property of the algorithm. An ill-conditioned problem cannot be saved by a clever algorithm; an unstable algorithm can ruin a perfectly well-conditioned problem

Conditioning is about speed and stability is about accuracy

Both are properties of the algorithm