Forward Error, Backward Error, and Condition Numbers

A numerical result can be wrong for two entirely different reasons: the algorithm lost accuracy that the problem did not force it to lose, or the problem itself turns small changes in the input into large changes in the output. Conflating the two leads to the wrong remedy: adding precision when the algorithm is unstable, or rewriting the algorithm when the problem is simply ill-conditioned. The pair of concepts that separates them is forward error against backward error, and the quantity that connects them is the condition number.

Two questions about the same discrepancy

Let y = f(x) be the exact value and ŷ the computed one. The forward error is the difference between them, usually measured relatively: |ŷ − y| / |y|. It answers the question a user asks, how far the answer is from the truth, and it is often impossible to measure, because the true value is what the computation was for.

The backward error asks a different question. Find the smallest perturbation Δx such that ŷ = f(x + Δx) exactly. If such a Δx exists with |Δx| / |x| on the order of the unit roundoff, then the computed result is the exact answer to a problem very close to the one that was posed. An algorithm with that property for every input is backward stable.

The two are related by a first-order identity that is the reason both concepts exist:

$$ \text{forward error} ;\lesssim; \kappa(x) \cdot \text{backward error} $$

where κ is the condition number of f at x. For linear problems in a norm, the inequality is exact up to the definition of the norms; for general nonlinear functions it is a first-order statement that holds when the perturbation is small, which is the regime the unit roundoff lives in. The product form is the useful part: it says the forward error is the product of a property of the algorithm and a property of the problem, and the two can be improved independently, or not at all.

Backward error comes free from the rounding model

The rounding model supplies a backward error without any analysis of the algorithm. For round-to-nearest in the normal range, fl(x + y) = (x + y)(1 + δ) with |δ| ≤ u. Rewrite that as

$$ \operatorname{fl}(x + y) = x(1 + \delta) + y(1 + \delta) $$

and the computed sum is exactly the sum of two perturbed inputs, each perturbed by a relative amount no larger than u. The backward error of floating-point addition is therefore u, with no dependence on the magnitudes involved. Every one of the four basic operations has the same property, and so does the square root.

Composition preserves it, with accumulation. A sum of n terms computed left to right is the exact sum of perturbed terms with relative perturbations bounded by γ_{n−1}, the same quantity that appeared in the forward bound of Compensated Summation and How Summation Order Bounds Error. A dot product is the exact dot product of perturbed vectors. A Horner evaluation of a polynomial is the exact evaluation of a polynomial with perturbed coefficients. In each case the statement is about the input data, not about the result, and it holds regardless of whether the terms cancel.

That last point is what makes backward error the right diagnostic for cancellation. When a sum of large cancelling terms produces a small result with no correct digits, the backward error is still γ_{n−1}: the algorithm did exactly what it promised. The forward error is large because κ for that sum is Σ|xᵢ| / |Σxᵢ|, which is large precisely when the terms cancel. The algorithm is stable; the problem is ill-conditioned; and the only cures are a reformulation of the problem or a wider format.

The condition number is a property of the problem

For a scalar function of one variable, the relative condition number is

$$ \kappa(x) = \left| \frac{x f’(x)}{f(x)} \right| $$

which is the ratio of the relative change in the output to a relative change in the input, in the limit of small changes. It depends only on f and x, not on how f is computed, which is what makes it a statement about the problem.

Three examples show how differently it behaves.

For f(x) = x − y with y fixed, κ = |x| / |x − y|, which grows without bound as x → y. This is the cancellation condition number from Catastrophic Cancellation and the Cases Where Subtraction Is Exact, derived here rather than asserted.

For f(x) = ln x, κ = 1/|ln x|. Near x = 1 the logarithm is close to zero, so a relative perturbation of the input produces a much larger relative perturbation of the output. The function is ill-conditioned there in the relative sense even though it is perfectly smooth, which is why log1p exists: it does not change the conditioning, but it prevents the input 1 + x from being rounded before the function sees it.

For polynomial evaluation, κ = Σ|aᵢ||x|^i / |p(x)|. A polynomial with large coefficients that nearly cancel at x has a large condition number at that point, and the relative accuracy of p(x) is limited accordingly, however the evaluation is arranged.

#include <math.h>

/* p(x) by the direct formula, plus the relative condition number of the
   evaluation at x. The condition number is a property of the polynomial. */
static double poly_with_cond(const double *a, int n, double x, double *kappa)
{
    double p = 0.0, mag = 0.0, xp = 1.0;
    for (int i = 0; i <= n; i++) {
        p   += a[i] * xp;
        mag += fabs(a[i]) * xp;   /* the same sum with absolute values */
        xp  *= x;
    }
    *kappa = (p != 0.0) ? mag / fabs(p) : INFINITY;
    return p;
}

The second accumulator is the whole computation: the condition number is the magnitude of the terms divided by the magnitude of the result, which is why it can be estimated during an evaluation at the cost of one extra sum.

For linear systems the condition number is normwise, κ(A) = ‖A‖‖A⁻¹‖ in whatever norm is in use, and it bounds the amplification of a relative perturbation in b or in A when solving Ax = b. It is always at least 1, equals 1 for a scalar multiple of an orthogonal matrix in the 2-norm, and is expensive to compute exactly, so the usual practice is to estimate it with a few matrix-vector products rather than to form A⁻¹.

Residuals, and the trap of a small one

The residual of a computed solution x̂ to Ax = b is r = b − Ax̂. It is cheap to compute and it is not the error. The relationship is

$$ \frac{|\hat{x} - x|}{|x|} \le \kappa(A) \frac{|r|}{|b|} $$

so a residual of order u‖b‖ guarantees an error of order κ(A)·u. When κ(A) is 10^10 in binary64, the residual can look like perfect accuracy while the solution has about six correct digits. When κ(A) approaches 1/u, the computed solution carries no information at all, and the matrix is numerically singular even if it is exactly nonsingular.

This is the single most common misreading of numerical output. A solver that reports a residual norm has reported the backward error, not the forward error, and the two differ by the condition number. The honest report is either the condition number estimate alongside the residual, or a forward error bound obtained by a different route.

The residual also has a backward-error interpretation that is exact: if r is the computed residual, then x̂ is the exact solution of (A + ΔA)x̂ = b with ‖ΔA‖ of order ‖r‖/‖x̂‖. A small residual means the computed solution solves a nearby system exactly, which is precisely the definition of backward stability for a linear solver. Gaussian elimination with partial pivoting is backward stable in this sense, with a bound on ΔA of order n²u‖A‖ multiplied by a growth factor that depends on the pivoting and the matrix.

What the product of the two buys

The decomposition turns a vague question into two measurable ones, and the measurements point at different remedies.

A large forward error with a small backward error and a large condition number is not an algorithm defect. The algorithm solved a nearby problem exactly; the problem amplifies the neighbourhood. The remedies are to reformulate the problem so that the condition number is smaller, to compute in a wider format so that u is smaller, or to accept the loss and report it. Iterative refinement sits in the third category with a twist: it reduces the backward error further by solving for the correction in working precision, which helps when the condition number is moderate but does nothing when it is at the level of 1/u.

A large forward error with a large backward error is an algorithm defect, and the condition number is irrelevant to the diagnosis. An unstable recurrence, a summation that loses the compensation because the compiler removed it, a matrix factorisation without pivoting on a matrix that needs it — all produce results that are not the exact answer to any nearby problem. Adding precision to such an algorithm reduces the backward error and can reduce the forward error with it, which is why wider formats sometimes appear to fix an instability; they are shrinking the perturbation, not repairing the method.

The case that resists both diagnoses is a problem whose condition number cannot be reduced and whose data is already only known to a few digits. There, the forward error is dominated by the uncertainty in the input, and the floating-point arithmetic is not the limiting factor. Distinguishing that case from an unstable algorithm requires the condition number, which is why it is worth computing even when the algorithm is known to be stable. The machinery for computing with the residuals that backward error analysis produces — exact products, exact sums, and the extended-precision formats built from them — is the subject of Error-Free Transformations and Double-Double Arithmetic.