Iterative refinement

When solving a linear system ${\displaystyle A\mathbf{x}=\mathbf{b},}$ due to the compounded accumulation of rounding errors, the computed solution ${\displaystyle \hat{\mathbf{x}}}$ may sometimes deviate from the exact solution ${\displaystyle {\mathbf{x}}_{\star }.}$ Starting with ${\displaystyle {\mathbf{x}}_1=\hat{\mathbf{x}},}$ iterative refinement computes a sequence ${\displaystyle \{{\mathbf{x}}_1,{\mathbf{x}}_2,{\mathbf{x}}_3,\dots \}}$ which converges to ${\displaystyle {\mathbf{x}}_{\star },}$ when certain assumptions are met.

Description

For ${\displaystyle m=1,2,3,\dots ,}$ the mth iteration of iterative refinement consists of three steps:

1. Compute the residual error r_m $${\displaystyle {\mathbf{r}}_m=\mathbf{b}-A{\mathbf{x}}_m.}$$
2. Solve the system for the correction, c_m, that removes the residual error $${\displaystyle A{\mathbf{c}}_m={\mathbf{r}}_m.}$$
3. Add the correction to get the revised next solution x_m+1 $${\displaystyle {\mathbf{x}}_{m+1}={\mathbf{x}}_m+{\mathbf{c}}_m.}$$

The crucial reasoning for the refinement algorithm is that although the solution for c_m in step (ii) may indeed be troubled by similar errors as the first solution, ${\displaystyle \hat{\mathbf{x}}}$, the calculation of the residual r_m in step (i), in comparison, is numerically nearly exact: You may not know the right answer very well, but you know quite accurately just how far the solution you have in hand is from producing the correct outcome (b). If the residual is small in some sense, then the correction must also be small, and should at the very least steer the current estimate of the answer, x_m, closer to the desired one, ${\displaystyle {\mathbf{x}}_{\star }.}$

The iterations will stop on their own when the residual r_m is zero, or close enough to zero that the corresponding correction c_m is too small to change the solution x_m which produced it; alternatively, the algorithm stops when r_m is too small to convince the linear algebraist monitoring the progress that it is worth continuing with any further refinements.

Note that the matrix equation solved in step (ii) uses the same matrix ${\displaystyle A}$ for each iteration. If the matrix equation is solved using a direct method, such as Cholesky or LU decomposition, the numerically expensive factorization of ${\displaystyle A}$ is done once and is reused for the relatively inexpensive forward and back substitution to solve for c_m at each iteration.