Condition Number, Residual, and Backward Error

A solver returns a vector, and the first thing anyone does with it is compute b - A x. If that vector is small, the answer is assumed to be good. The inference is invalid in a specific and quantifiable way, and the quantity that makes it invalid is the same one that determines how much accuracy any algorithm could have delivered.

Separating the two contributions, the algorithm’s and the problem’s, is the purpose of the backward error framework. It is the only way to answer the question “is this answer as good as it can be” without knowing the true solution, which is the situation in every real computation.

The residual is not the error

Write x for the exact solution and x_hat for the computed one. The error is x - x_hat; the residual is r = b - A x_hat. They are related by A (x - x_hat) = r, so the error is the residual passed through A^{-1}.

If A^{-1} has a large norm, a small residual maps to a large error. A matrix that is nearly singular has directions in which it barely stretches anything; a residual concentrated in such a direction produces an error far larger than the residual suggests. Conversely, a well-conditioned matrix maps a small residual to an error of comparable relative size, and the inference from residual to error is sound.

The extreme case is a matrix that is exactly singular: a nonzero residual can be orthogonal to the range, so no solution exists at all, while a residual of zero can be produced by infinitely many vectors, most of them far from any particular solution. The computed answer then has nothing to do with the question that was asked, and the residual reports success.

The inequality that splits the blame

Taking norms of x - x_hat = A^{-1} r and dividing by the norm of x gives

||x - x_hat|| / ||x||  <=  ||A|| ||A^{-1}|| * ||r|| / ||b||

because ||b|| <= ||A|| ||x||. The factor kappa(A) = ||A|| ||A^{-1}|| is the condition number, and it depends on the norm chosen; the 2-norm version is the ratio of the largest to the smallest singular value.

The inequality is sharp: for every matrix there exist a right-hand side and a perturbation that attain it. That means the condition number is not a pessimistic bound invented for the analysis; it is the exact worst-case amplification of a relative residual into a relative error.

Reading the inequality as a budget is the useful interpretation. If kappa = 1e6 and the residual is at the level of the rounding error, 1e-16 relative, then the relative error is at most about 1e-10. The solver has done its job, and the remaining eight digits of loss are a property of the matrix, not of the algorithm. No factorization, no iterative refinement, and no higher-precision arithmetic inside the same formulation can recover them, because the information was destroyed by the conditioning of the problem itself.

Backward error is what the algorithm controls

The alternative to bounding the forward error directly is to ask a different question: for what nearby problem is the computed answer exact? If (A + E) x_hat = b with ||E|| <= eta ||A||, then eta is the normwise backward error, and the answer is the exact solution of a problem within eta of the original.

This is the right frame because it separates cleanly. The backward error is a property of the algorithm and the arithmetic; the condition number is a property of the problem. Their product bounds the forward error, and each can be measured or bounded independently.

For LU with partial pivoting, the backward error is of order n^2 * rho * u * ||A||, with rho the growth factor discussed in LU Factorization with Partial Pivoting and Element Growth. For Cholesky it is n * u * ||A|| with no growth factor at all. For Householder QR it is n * u * ||A||, as described in Householder QR and Orthogonality Loss in Gram-Schmidt. In every case the bound is a modest multiple of u times the matrix norm, and none of them contains the condition number.

A backward stable algorithm is therefore one whose computed answer is exact for a nearby problem. That is the strongest guarantee available in floating point, and it is achievable for every factorization discussed in this cluster.

Estimating the condition number without an inverse

Computing kappa from its definition requires ||A^{-1}||, and forming the inverse is both expensive, at 2n^3 flops, and numerically the wrong operation: the entries of an explicit inverse are computed from the same nearly dependent columns that make the problem ill conditioned.

What is needed is a norm of A^{-1}, not the matrix. Any norm of an operator can be estimated from below by evaluating it on a few vectors, and the 1-norm has a structure that makes this effective. The 1-norm of a matrix is the maximum absolute column sum, attained by a column; the 1-norm of A^{-1} is attained by some column j of the inverse, which equals A^{-1} e_j. So the estimate can be built by solving with candidate unit vectors, examining the resulting vector, and using the sign pattern of its largest entries to construct a better candidate.

Hager’s algorithm does exactly this, and it typically converges in a handful of iterations. The cost is a few solves with the existing factors plus a few products with A, which is negligible next to the factorization itself. LAPACK exposes it through the gecon, pocon, and trcon family, which take the factors and return an estimate of the reciprocal condition number, 1 / kappa, in a specified norm. Returning the reciprocal is deliberate: it avoids overflow when the matrix is nearly singular, and it makes the comparison against u direct.

import math

def matvec(A, x):
    return [sum(a * b for a, b in zip(row, x)) for row in A]

def norm2(v):
    return math.sqrt(sum(t * t for t in v))

def sigma_max(A, iters=100):
    """Power iteration on A^T A, which is never formed."""
    n = len(A[0])
    x = [1.0 / math.sqrt(n)] * n
    s = 0.0
    for _ in range(iters):
        y = matvec(A, x)                       # A x
        z = [sum(A[i][j] * y[i] for i in range(len(A))) for j in range(n)]  # A^T y
        s = norm2(z)
        if s == 0.0:
            return 0.0
        x = [t / math.sqrt(s) for t in z]
    return math.sqrt(s)

# sigma_min follows from the same iteration applied to (A^T A)^{-1},
# i.e. one solve with A and one with A^T per step, using the existing factors.

The estimator’s accuracy is one-sided in a useful way: it never overestimates the norm, since it evaluates the norm on a specific vector, and it is usually within a small factor of the truth. The documented worst case is a factor of n, which is why a condition estimate is a diagnostic rather than a certificate.

Where the estimate misleads

The condition number is defined relative to a norm, and different norms can differ by factors that grow with the dimension. A matrix can have a 1-norm condition number in the thousands and a 2-norm condition number in the millions, or the reverse, and neither is wrong. Comparisons between codes that report different norms are meaningless without checking which one was used.

Scaling is the larger trap. The condition number of A changes if rows or columns are rescaled, even though the solution of the scaled system is the same up to the diagonal transformation. A matrix whose rows are measured in wildly different units can have a normwise condition number of 1e15 while every entry of the solution is computed to full relative accuracy, because the sensitivity is concentrated in the scaling rather than in the underlying operator. The componentwise condition number, which measures the sensitivity of each component relative to its own scale, is the appropriate diagnostic there, and it can be far smaller than the normwise one.

Non-normal matrices add a third failure. For a matrix with a large departure from normality, the norm of a power of A can exceed the norm of A raised to that power by a large factor, and the eigenvalue-based intuition about conditioning breaks down. The singular values still describe the sensitivity of the linear system correctly, but they say nothing about the sensitivity of the eigenvalues, which for a non-normal matrix can be enormous even when the singular values are all of order one.

Digits lost, stated plainly

In double precision the unit roundoff is 2^-53, about 1.1e-16. A backward stable solve therefore delivers a relative error of roughly kappa * u, which corresponds to losing about log10(kappa) decimal digits from the sixteen that the format carries.

The arithmetic of that rule is worth internalizing. A condition number of 1e4 costs four digits, leaving twelve. A condition number of 1e12 costs twelve, leaving four. A condition number near 1e16 costs everything, and the computed answer is dominated by rounding. In single precision, with a unit roundoff near 6e-8, the corresponding threshold is a condition number near 1e7.

The rule also explains when a residual check is informative. If the computed residual is at the level of u * ||A|| * ||x||, the algorithm has achieved its backward error bound and the remaining error is the condition number’s contribution. If the residual is much larger than that, the algorithm has failed, and the cause is a bug or a breakdown rather than the problem’s conditioning. Distinguishing those two cases is the entire practical value of the backward error framework, and it is why a solver that returns a condition estimate alongside the answer is more useful than one that returns the answer alone.

Iterative refinement exploits the distinction directly: compute the residual, solve for a correction with the existing factors, and repeat. Each step reduces the forward error as long as the residual is computed accurately enough relative to the condition number, which is why the technique is paired with factors computed in a lower precision. The same quantity appears from the other direction in preconditioning. There the factorization is deliberately approximate, and the condition number of the preconditioned operator, rather than that of the original matrix, is what decides how fast an iteration converges; Incomplete LU Preconditioning and Dropped Fill describes how the dropped entries shape that spectrum.