192.16.184.89 at 16:22, 2 February 2024

2024-02-02T16:22:38Z

New page

{{Short description|Method to improve accuracy of numerical solutions to systems of linear equations}}
{{broader|Iterative method}}
{{about|iterative refinement in mathematics|iterative refinement in software development|Iterative and incremental development}}
{{use dmy dates|date=July 2020}}
'''Iterative refinement''' is an [[iterative method]] proposed by [[James H. Wilkinson]] to improve the accuracy of numerical solutions to [[systems of linear equations]].<ref name=Wilkinson1963>{{cite book |first=James H. |last=Wilkinson |year=1963 |title=Rounding Errors in Algebraic Processes |publisher=[[Prentice Hall]] |place=Englewood Cliffs, NJ |url=https://archive.org/details/roundingerrorsin0000wilk |url-access=registration}}</ref><ref name=Moler1967>{{cite journal |first=Cleve B. |last=Moler |date=April 1967 |title=Iterative refinement in floating point |journal=[[Journal of the ACM]] |volume=14 |issue=2 |pages=316–321 |publisher=[[Association for Computing Machinery]] |place=New York, NY|doi=10.1145/321386.321394|doi-access=free }}</ref>

When solving a linear system <math>A \mathbf{x} = \mathbf{b} \,,</math> due to the compounded accumulation of [[Round-off error|rounding error]]s, the computed solution <math>\hat{\mathbf{x}}</math> may sometimes deviate from the exact solution <math>\mathbf{x}_\star\,.</math> Starting with <math>\mathbf{x}_1 = \hat{\mathbf{x}}\,,</math> iterative refinement computes a sequence <math>\{\mathbf{x}_1,\, \mathbf{x}_2,\, \mathbf{x}_3,\dots\}</math> which converges to <math>\mathbf{x}_\star\,,</math> when certain assumptions are met.

==Description==
For <math>m = 1, 2, 3, \dots\,,</math> the {{mvar|m}}th iteration of iterative refinement consists of three steps:
{{ordered list | list-style-type = lower-roman
| Compute the residual error {{math|'''r'''{{sub|''m''}}}}
<math display="block">\mathbf{r}_m = \mathbf{b} - A \mathbf{x}_m\,.</math>
| Solve the system for the correction, {{math|'''c'''{{sub|''m''}}}}, that removes the residual error
<math display="block">A \mathbf{c}_m = \mathbf{r}_m\,.</math>
| Add the correction to get the revised next solution {{math|'''x'''{{sub|''m''+1}}}}
<math display="block">\mathbf{x}_{m+1} = \mathbf{x}_m + \mathbf{c}_m\,.</math>
}}

The crucial reasoning for the refinement algorithm is that although the solution for {{math|'''c'''{{sub|''m''}}}} in step (ii) may indeed be troubled by similar errors as the first solution, <math>\hat\mathbf{x}</math>, the calculation of the residual {{math|'''r'''{{sub|''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 ({{math|'''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, {{math|'''x'''{{sub|''m''}}}}, closer to the desired one, <math>\mathbf{x}_\star\,.</math>

The iterations will stop on their own when the residual {{math|'''r'''{{sub|''m''}}}} is zero, or close enough to zero that the corresponding correction {{math|'''c'''{{sub|''m''}}}} is too small to change the solution {{math|'''x'''{{sub|''m''}}}} which produced it; alternatively, the algorithm stops when {{math|'''r'''{{sub|''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 <math>A</math> for each iteration. If the matrix equation is solved using a direct method, such as
[[Cholesky decomposition|Cholesky]] or [[LU decomposition|LU]] [[Matrix decomposition|decomposition]], the numerically expensive factorization of <math>A</math> is done once and is reused for the relatively inexpensive [[forward substitution|forward]] and [[back substitution]] to solve for {{math|'''c'''{{sub|''m''}}}} at each iteration.<ref name=Moler1967/>

==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 {{math|'''r'''}}, e.g. by using [[Quadruple-precision floating-point format|quad]] or [[extended precision|double extended]] precision [[IEEE 754]] [[floating point]], and if {{math|A}} is not too ill-conditioned (and the iteration and the rate of convergence are determined by the condition number of {{math|A}}).<ref>{{cite book |first=Nicholas |last=Higham |title=Accuracy and Stability of Numerical Algorithms |edition=2 |publisher=SIAM |year=2002 |pages=232}}</ref>

More formally, assuming that each step (ii) can be solved reasonably accurately, i.e., in mathematical terms, for every {{mvar|m}}, we have
<math display="block">A \left( I + F_m \right) \mathbf{c}_m = \mathbf{r}_m</math>

where {{math|‖F{{sub|''m''}}‖{{sub|∞}} < 1}}, the [[Approximation error|relative error]] in the {{mvar|m}}-th iterate of iterative refinement satisfies
<math display="block">\frac{ \lVert \mathbf{x}_m - \mathbf{x}_\star \rVert_\infty }{ \lVert \mathbf{x}_\star \rVert_\infty} \leq \bigl( \sigma\,\kappa(A)\,\varepsilon_1\bigr)^m + \mu_1\, \varepsilon_1 + n\, \kappa(A)\, \mu_2\, \varepsilon_2</math>

where
* {{math|‖·‖<sub>∞</sub>}} denotes the [[Uniform norm|{{math|∞}}-norm]] of a vector,
* {{math|''κ''(A)}} is the {{math|∞}}-[[condition number]] of {{math|A}},
* {{mvar|n}} is the order of {{math|A}},
* {{mvar|ε}}{{sub|1}} and {{mvar|ε}}{{sub|2}} are [[Machine epsilon|unit round-offs]] of [[floating-point arithmetic]] operations,
* {{mvar|σ}}, {{mvar|μ}}{{sub|1}} and {{mvar|μ}}{{sub|2}} are constants that depend on {{math|A}}, {{mvar|ε}}{{sub|1}} and {{mvar|ε}}{{sub|2}}
if {{math|A}} is "not too badly conditioned", which in this context means

{{block indent | em = 1.5 | text = {{math|0 < ''σ κ''(A) ''ε''{{sub|1}} ≪ 1}}}}

and implies that {{mvar|μ}}{{sub|1}} and {{mvar|μ}}{{sub|2}} are of order unity.

The distinction of {{mvar|ε}}{{sub|1}} and {{mvar|ε}}{{sub|2}} is intended to allow mixed-precision evaluation of {{math|'''r'''{{sub|''m''}}}} where intermediate results are computed with unit round-off {{mvar|ε}}{{sub|2}} before the final result is rounded (or truncated) with unit round-off {{mvar|ε}}{{sub|1}}. All other computations are assumed to be carried out with unit round-off {{mvar|ε}}{{sub|1}}.

==References==
{{reflist|25em}}

{{Numerical linear algebra}}

[[Category:Numerical linear algebra]]

Iterative refinement - Revision history

192.16.184.89 at 16:22, 2 February 2024