Iterative Refinement

In the splitting A=N-P for the solution of $A\mbox{\boldmath$\space x $ } =\mbox{\boldmath$\space f $ }$, let N=LU, where L and U are the lower and upper triangular factors of Aalready found approximately by Gaussian elimination or some other method. Then P=N-A=LU-A and the iteration matrix is M=N^-1(N-A)=I-N^-1A=I-(LU)^-1A. If the factors are reasonably accurate, the elements of this matrix will be small, so that, e.g., $\left\Vert \kern.05em M \kern.05em \right\Vert _{\infty }$ will be small, and we should get good convergence. We take this slightly further by writing LU=A+E, giving

$$\begin{eqnarray*}M & = & I-(A+E)^{-1}A\\ & = & (A+E)^{-1}(A+E-A)\\ & = & (I+A^{-1}E)^{-1}A^{-1}E, \end{eqnarray*}$$

and taking norms, we find

$$\left\Vert \kern.05em M \kern.05em \right\Vert\leq \frac{\lef... ...\Vert}{1-\left\Vert \kern.05em A^{-1}E \kern.05em \right\Vert}$$

if $\left\Vert \kern.05em A^{-1}E \kern.05em \right\Vert<1$. For $\left\Vert \kern.05em M \kern.05em \right\Vert<1$, we require $\left\Vert \kern.05em A^{-1}E \kern.05em \right\Vert<\frac{1}{2}$, and this will make the iteration convergent. This will certainly be true if $\left\Vert \kern.05em A^{-1} \kern.05em \right\Vert.\left\Vert \kern.05em E \kern.05em \right\Vert<\frac{1}{2}$.

The method is implemented in the following way:

$N\mbox{\boldmath$\space x $ } _{m+1}=P\mbox{\boldmath$\space x $ } _m+\mbox{\boldmath$\space f $ }$ becomes $LU\mbox{\boldmath$\space x $ } _{m+1}=(LU-A)\mbox{\boldmath$\space x $ } _m+\mbox{\boldmath$\space f $ }$, which may be written as

$$LU(\mbox{\boldmath$ x $ } _{m+1}-\mbox{\boldmath$ x $ } _m)=\... ...$ f $ } -A\mbox{\boldmath$ x $ } _m=\mbox{\boldmath$ r $ } _m,$$

the residual at step m. This equation is easy to solve if we already have L and U available from our original solution using

$$L\mbox{\boldmath$ y $ } _m=\mbox{\boldmath$ r $ } _m,\quad U(... ...} _{m+1}-\mbox{\boldmath$ x $ } _m)=\mbox{\boldmath$ y $ } _m.$$

The solution $\mbox{\boldmath$\space x $ } _{m+1}-\mbox{\boldmath$\space x $ } _m$ is then added to $\mbox{\boldmath$\space x $ } _m$. We note that

(i)

the computation of $\mbox{\boldmath$\space r $ } _m$ must be sufficiently accurate, otherwise no improvement will result,

(ii)

as we have already seen, a small $\mbox{\boldmath$\space r $ }$ does not guarantee accuracy,

(iii)

in practice, only very ill-conditioned problems will need more than 2 iterations,

(iv)

this method should only be used as a refinement for an algorithm which computes L and U as part of its solution.