Iterative refinement

Template:Short description Template:Broader Template:About Template:Use dmy dates Iterative refinement is an iterative method proposed by James H. Wilkinson to improve the accuracy of numerical solutions to systems of linear equations.

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

Description

For $m=1,2,3,\dots ,$ the Template:Mvarth iteration of iterative refinement consists of three steps: Template:Ordered list

The crucial reasoning for the refinement algorithm is that although the solution for Template:Math in step (ii) may indeed be troubled by similar errors as the first solution, $\widehat{𝐱}$, the calculation of the residual Template:Math 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 (Template:Math). 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, Template:Math, closer to the desired one, ${𝐱}_{\star }.$

The iterations will stop on their own when the residual Template:Math is zero, or close enough to zero that the corresponding correction Template:Math is too small to change the solution Template:Math which produced it; alternatively, the algorithm stops when Template:Math 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 $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 $A$ is done once and is reused for the relatively inexpensive forward and back substitution to solve for Template:Math at each iteration.

Error analysis

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 Template:Math, e.g. by using quad or double extended precision IEEE 754 floating point, and if Template:Math is not too ill-conditioned (and the iteration and the rate of convergence are determined by the condition number of Template:Math).

More formally, assuming that each step (ii) can be solved reasonably accurately, i.e., in mathematical terms, for every Template:Mvar, we have $$A(I+F_m){𝐜}_m={𝐫}_m$$

where Template:Math, the relative error in the Template:Mvar-th iterate of iterative refinement satisfies $$\frac{\Vert {𝐱}_m-{𝐱}_{\star }{\Vert }_{\mathrm{\infty }}}{\Vert {𝐱}_{\star }{\Vert }_{\mathrm{\infty }}}\le (\sigma \kappa (A){\epsilon }_1{)}^m+{\mu }_1{\epsilon }_1+n\kappa (A){\mu }_2{\epsilon }_2$$

where

• Template:Math denotes the [[Uniform norm|Template:Math-norm]] of a vector,
• Template:Math is the Template:Math-condition number of Template:Math,
• Template:Mvar is the order of Template:Math,
• Template:MvarTemplate:Sub and Template:MvarTemplate:Sub are unit round-offs of floating-point arithmetic operations,
• Template:Mvar, Template:MvarTemplate:Sub and Template:MvarTemplate:Sub are constants that depend on Template:Math, Template:MvarTemplate:Sub and Template:MvarTemplate:Sub

if Template:Math is "not too badly conditioned", which in this context means

and implies that Template:MvarTemplate:Sub and Template:MvarTemplate:Sub are of order unity.

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