The a priori bounds in equations (4.4.4) and (4.4.6) are generally not at all tight and are of limited usefulness in monitoring the stability. We therefore now look at a posteriori methods of monitoring stability. If instability is detected, further work will be necessary to get a solution with as much accuracy as the data warrants. This might involve iterative refinement (Section 4.13), choosing the pivotal sequence afresh, or even working with greater precision.

Two forms of stability monitoring are suggested by the discussion of Section 4.3. After the factorization has been found we may compute the factorization-error matrix H = LU — A, or after the solution has been found we may compute the residual vector r = b — Ax.

Computing the residual r is easier, and in practice this measure for the stability of the solution is usually employed. Indeed, if ||r|| is small compared with ||b||, we have obviously solved a nearby problem since Ax = b — r; if ||r|| is small compared with ||A|| ||x||, then x is the exact solution of the equation (A + H)x = b where ||H|| is small compared with ||A|| (see Exercise 4.4). Thus, if ||r|| is small compared with ||b||, ||Ax||, or ||A|| ||x||, we have done a good job in solving the equation.

For the converse, it is important to compare ||r|| with ||A|| ||x||, rather than with | b| since | r| may be large compared with | b| in spite of an accurately computed approximate solution. The set of equations

has as exact solution

yet the approximate solution

which would be considered very accurate on a 3-decimal computer, has residual

and certainly ||r|| is large relative to ||b|| (although not with respect to ||A|| ||x||). This is an example of ill-conditioning, a subject to be explored later in this chapter.

Moreover, r tells us nothing about the behaviour for other vectors b. For example, the linear system with matrix (4.3.1) and right-hand side vector

has the exact solution

and this would be produced exactly by our 3-digit computer with the algorithm of Section 4.3. This does not imply that the factorization is stable, but rather that the instability does not affect the solution for this particular b. We cannot assume that the same will be true for other vectors b, and we know that the algorithm of Section 4.3 is unstable.

The alternative of stability monitoring by checking the size of

is expensive. Even computing p = maxj_{jjk} a^{(k)} is expensive, since it involves

testing over all i,j,k, which is an O(n^{3}) operation for a dense matrix.

Attempts have been made to find bounds for p (for example, Erisman and Reid (1974)) but, for general matrices, the bounds are too loose to be useful. We describe a practical application where this has been used effectively in Section 4.8.

The best way of monitoring stability is to compute a posteriori the scaled residual

Note that this will not tell you how accurate the solution is, but only how close the problem that was actually solved is to the original problem.