Iterative refinement

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

(i)	Compute the residual error rm	${\displaystyle \mathbf {r} _{m}=\mathbf {b} -\mathrm {A} \mathbf {x} _{m}\,.}$
(ii)	Solve the system for the correction, cm,that removes the residual error	${\displaystyle \mathrm {A} \,\mathbf {c} _{m}=\mathbf {r} _{m}\,.}$
(iii)	Add the correction to get the revised next solution xm+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.

As a rule of thumb, iterative refinement for Gaussian elimination produces a solution correct to working precision if double the working precision is used in the computation of r, e.g. by using quad or double extended precision IEEE 754 floating point, and if A is not too ill-conditioned (and the iteration and the rate of convergence are determined by the condition number of A).[3]

More formally, assuming that each step (ii) can be solved reasonably accurately, i.e., in mathematical terms, for every m, we have

${\displaystyle \mathrm {A} \,\left(\mathrm {I+F} _{m}\right)\mathbf {c} _{m}=\mathbf {r} _{m}}$

where ‖F_m‖_∞ < 1, the relative error in the mth iterate of iterative refinement satisfies

${\displaystyle {\frac {\lVert \mathbf {x} _{m}-\mathbf {x} _{\star }\rVert _{\infty }}{\lVert \mathbf {x} _{\star }\rVert _{\infty }}}\leq {\bigl (}\sigma \,\kappa (\mathrm {A} )\,\varepsilon _{1}{\bigr )}^{m}+\mu _{1}\,\varepsilon _{1}+n\,\kappa (\mathrm {A} )\,\mu _{2}\,\varepsilon _{2}}$

where

• ‖·‖_∞ denotes the ∞-norm of a vector,
• κ(A) is the ∞-condition number of A,
• n is the order of A,
• ε₁ and ε₂ are unit round-offs of floating-point arithmetic operations,
• σ, μ₁ and μ₂ are constants that depend on A, ε₁ and ε₂

if A is "not too badly conditioned", which in this context means

0 < σ κ(A) ε₁ ≪ 1

and implies that μ₁ and μ₂ are of order unity.

The distinction of ε₁ and ε₂ is intended to allow mixed-precision evaluation of r_m where intermediate results are computed with unit round-off ε₂ before the final result is rounded (or truncated) with unit round-off ε₁. All other computations are assumed to be carried out with unit round-off ε₁.