Iterative refinement

Remarks

Let $e_m = x - x_m$, so that $r_m = A e_m$. Then $e_{m+1} = e_m - (A + E_m)^{-1} A e_m = (A + E_m)^{-1} E_m e_m$. Since $\norm{A^{-1} E_m} \le \eta < 1$, the matrix $A + E_m = A(I + A^{-1} E_m)$ is invertible with $\norm{(A + E_m)^{-1}} \le \norm{A^{-1}} / (1 - \eta)$, which gives the bound. The matrix $E_m$ models the error of an inexact solver, such as a computed LU decomposition reused at each step, or a solve in lower precision. If the solver is too inaccurate, so that $\eta \ge 1/2$, the bound gives no contraction and the process can diverge. This happens when the condition number of $A$ times the solver's relative accuracy is not small. In floating-point practice, the residual is the critical step: computing $r_m$ in higher precision lets refinement reach a solution accurate to the working precision, even for fairly ill-conditioned $A$. In exact arithmetic with $E_m = 0$, one step gives the exact solution.

Going deeper

For a student

Iterative refinement turns a cheap or low-precision solver into an accurate one. It reuses one factorization of $A$, so each extra round costs only $O(n^2)$ operations after the $O(n^3)$ factorization. Modern mixed-precision solvers factor in low precision and refine in higher precision to get both speed and accuracy.