Standard Error Analysis

We illustrate standard error analysis with the simple example of evaluating the scalar function y=f(z). Let the output of the subroutine which implements f(z) be denoted ${\rm alg}(z)$; this includes the effects of roundoff. If ${\rm alg}(z) = f(z+\delta)$ where $\delta$ is small, then we say ${\rm alg}$ is a backward stable algorithm for f, or that the backward error $\delta$ is small. In other words, ${\rm alg}(z)$ is the exact value of f at a slightly perturbed input $z+\delta$.^4.5

Suppose now that f is a smooth function, so that we may approximate it near z by a straight line: $f(z+\delta) \approx f(z) + f'(z) \cdot \delta$. Then we have the simple error estimate

$${\rm alg}(z)-f(z) = f(z+\delta) - f(z) \approx f'(z) \cdot \delta .$$

Thus, if $\delta$ is small, and the derivative f'(z) is moderate, the error ${\rm alg}(z)-f(z)$ will be small^4.6. This is often written in the similar form

$$\left\vert \frac{{\rm alg}(z)-f(z)}{f(z)} \right\vert \appro... ...iv \kappa (f,z) \cdot \left\vert \frac{\delta}{z} \right\vert .$$

This approximately bounds the relative error $\left\vert \frac{{\rm alg}(z)-f(z)}{f(z)} \right\vert$ by the product of the condition number of f at z, $\kappa (f,z)$, and the relative backward error $\vert\frac{\delta}{z}\vert$. Thus we get an error bound by multiplying a condition number and a backward error (or bounds for these quantities). We call a problem ill-conditioned if its condition number is large, and ill-posed if its condition number is infinite (or does not exist)^4.7.

If f and z are vector quantities, then f'(z) is a matrix (the Jacobian). So instead of using absolute values as before, we now measure $\delta$ by a vector norm $\Vert \delta \Vert$ and f'(z) by a matrix norm |f'(z)|. The conventional (and coarsest) error analysis uses a norm such as the infinity norm. We therefore call this normwise backward stability. For example, a normwise stable method for solving a system of linear equations Ax=b will produce a solution $\hat{x}$ satisfying $(A+E)\hat{x}=b+f$ where $\Vert E\Vert _{\infty}/ \Vert A\Vert _{\infty}$ and $\Vert f\Vert _{\infty}/ \Vert b\Vert _{\infty}$ are both small (close to machine epsilon). In this case the condition number is $\kappa_{\infty}(A) = \Vert A\Vert _{\infty}\cdot \Vert A^{-1}\Vert _{\infty}$ (see section 4.4 below).

Almost all of the algorithms in LAPACK (as well as LINPACK and EISPACK) are stable in the sense just described^4.8: when applied to a matrix A they produce the exact result for a slightly different matrix A+E, where $\Vert E\Vert _{\infty}/ \Vert A\Vert _{\infty}$ is of order $\epsilon$. In a certain sense, a user can hardly ask for more, provided the data is at all uncertain.

It is often possible to compute the norm |E| of the actual backward error by computing a residual r, such as r=Ax-b or $r=Ax - \lambda x$, and suitably scaling its norm |r|. The expert driver routines for solving Ax=b do this, for example. For details see [55,67,85,95].

Condition numbers may be expensive to compute exactly. For example, it costs about $\frac{2}{3} n^3$ operations to solve Ax=b for a general matrix A, and computing $\kappa_{\infty}(A)$ exactly costs an additional $\frac{4}{3} n^3$ operations, or twice as much. But $\kappa_{\infty}(A)$ can be estimated in only O(n²) operations beyond those $\frac{2}{3} n^3$ necessary for solution, a tiny extra cost. Therefore, most of LAPACK's condition numbers and error bounds are based on estimated condition numbers, using the method of [59,62,63]. The price one pays for using an estimated rather than an exact condition number is occasional (but very rare) underestimates of the true error; years of experience attest to the reliability of our estimators, although examples where they badly underestimate the error can be constructed [65]. Note that once a condition estimate is large enough, (usually $O( 1/ \epsilon )$), it confirms that the computed answer may be completely inaccurate, and so the exact magnitude of the condition estimate conveys little information.