Numerical Analysis/Iterative Refinement

What is Iterative Refinement?

A detailed discussion of iterative refinement can be found on the Wikipedia page.

To give a brief description, it is a technique used to improve the approximate solution ${\widehat {x}}$ to linear system $Ax=b.$ This technique is generally only used on systems that are thought or determined to be ill-conditioned.

The process involves three primary steps:

Iterative Refinement Process

Once an approximation to the solution, ${\widehat {x}}_{1}$, has been made, with either rounded-Gaussian elimination or otherwise, complete the following steps: Compute residual: $r_{i}=b-A{\widehat {x}}_{i}$ Solve $Ay_{i}=r_{i}$ Update ${\widehat {x}}_{i+1}={\widehat {x}}_{i}+y_{i}$ Continue to update step one with the newly calculated ${\widehat {x}}_{i+1}$ until a desired tolerance has been reached.

Approximating the Condition Number

In theory, the condition number of a matrix depends on the norms of the matrix and its inverse. In practice, however, calculating the inverse is subject to round-off error and the accuracy of the calculations.

If these calculations involve arithmetic with t digits of accuracy, then the approximate condition number for the matrix - call it A - is the norm of A the norm of the approximation of the inverse of A.

Assuming that the approximate solution to the linear system is being determined with t-digit arithmetic and Gaussian elimination, we can show

${\begin{Vmatrix}r\end{Vmatrix}}\approx 10^{-t}{\begin{Vmatrix}A\end{Vmatrix}}{\begin{Vmatrix}{\widehat {x}}\end{Vmatrix}},$

where r is the residual vector for the approximation ${\widehat {x}}.$

This approximation for the condition number can be obtained without having to invert matrix A.

When doing iterative refinement problems, we will have, or will calculate the values for $A,y,$ and ${\widehat {x}}$. The approximate solution, ${\widehat {y}}$ satisfies

${\widehat {y}}\approx A^{-1}r=A^{-1}(b-A{\widehat {x}})=A^{-1}b-A^{-1}A{\widehat {x}}=x-{\widehat {x}},$

so ${\widehat {y}}$ is an estimate of the error for approximate solution ${\widehat {x}}.$ Observe that

${\begin{aligned}{\begin{Vmatrix}{\widehat {y}}\end{Vmatrix}}\approx {\begin{Vmatrix}x-{\widehat {x}}\end{Vmatrix}}&={\begin{Vmatrix}A^{-1}r\end{Vmatrix}}\leq {\begin{Vmatrix}A^{-1}\end{Vmatrix}}{\begin{Vmatrix}r\end{Vmatrix}}\approx {\begin{Vmatrix}A^{-1}\end{Vmatrix}}(10^{-t}{\begin{Vmatrix}A\end{Vmatrix}}{\begin{Vmatrix}{\widehat {x}}\end{Vmatrix}})\\&=10^{-t}{\begin{Vmatrix}{\widehat {x}}\end{Vmatrix}}K(A).\end{aligned}}$

This approximates the condition number associated with solving $Ax=b$ and can be re-written and expressed as seen here:

$K(A)\approx {\frac {\begin{Vmatrix}{\widehat {y}}\end{Vmatrix}}{\begin{Vmatrix}{\widehat {x}}\end{Vmatrix}}}10^{t}.$

Example

Consider the following ill-conditioned linear system

${\begin{aligned}3.9x_{1}+1.6x_{2}&=5.5\\6.8x_{1}+2.9x_{2}&=9.7\end{aligned}}$

Part 1

Solve using Gaussian elimination and iterative refinement (and using two-digit rounding arithmetic). Continue using iterative refinement until ${\widehat {x}}_{4}$ is found. Compare with the true solution.

Solution:

First, use a calculator or computer to find the true solution $x$ of this system. We will discover that $x={\begin{bmatrix}1\\1\end{bmatrix}}.$

Set up the augmented matrix and solve using Guassian elimination, making sure to use only two-digit rounding arithmetic.

${\begin{aligned}{\begin{bmatrix}3.9&1.6&5.5\\6.8&2.9&9.7\end{bmatrix}}&={\begin{matrix}R1/3.9\\R2/6.8\end{matrix}}&\to {\begin{bmatrix}1&0.41&1.4\\1&0.43&1.4\end{bmatrix}}\\\\&={\begin{matrix}R2-R1\end{matrix}}&\to {\begin{bmatrix}1&0.41&1.4\\0&0.02&0\end{bmatrix}}\\\\&={\begin{matrix}R2/0.02\end{matrix}}&\to {\begin{bmatrix}1&0.41&1.4\\0&1&0\end{bmatrix}}\end{aligned}}$

Solving with back-substitution we find the first approximation to the system:

${\widehat {x}}_{1}={\begin{bmatrix}x_{1}\\x_{2}\end{bmatrix}}={\begin{bmatrix}1.4\\0\end{bmatrix}}$

Because we mentioned that this system is ill-conditioned, we can use the iterative refinement process to find a better approximation. Recall the steps mentioned above:

1. Compute residual: $r_{i}=b-A{\widehat {x}}_{i}$
   
2. Solve $Ay_{i}=r_{i}$
   
3. Update ${\widehat {x}}_{i+1}={\widehat {x}}_{i}+y_{i}$
   
   

STEP 1: Compute the residual $r_{1}.$

${\begin{aligned}r_{1}=b-A{\widehat {x}}_{1}&={\begin{bmatrix}5.5\\9.7\end{bmatrix}}-{\begin{bmatrix}3.9&1.6\\6.8&2.9\end{bmatrix}}{\begin{bmatrix}1.4\\0\end{bmatrix}}\\\\&={\begin{bmatrix}5.5\\9.7\end{bmatrix}}-{\begin{bmatrix}5.5\\9.5\end{bmatrix}}\\\\&={\begin{bmatrix}0\\0.2\end{bmatrix}}=r_{1}\end{aligned}}$

STEP 2: Now solve $Ay_{1}=r_{1}$ using rounded Gaussian elimination to find the $y$ vector.

${\begin{aligned}Ay_{1}=r_{1}\quad \Rightarrow {\begin{bmatrix}3.9&1.6\\6.8&2.9\end{bmatrix}}{\begin{bmatrix}y_{1}\\y_{2}\end{bmatrix}}={\begin{bmatrix}0\\0.2\end{bmatrix}}\quad &\Rightarrow {\begin{bmatrix}3.9&1.6&0\\6.8&2.9&0.2\end{bmatrix}}={\begin{bmatrix}1&0.41&0\\1&0.43&0.029\end{bmatrix}}\\\\&={\begin{bmatrix}1&0.41&0\\0&0.02&0.029\end{bmatrix}}={\begin{bmatrix}1&0.41&0\\0&1&1.5\end{bmatrix}}\\\\&={\begin{bmatrix}1&0&-0.62\\0&1&1.5\end{bmatrix}}\\\\\end{aligned}}$

$\Rightarrow y_{1}={\begin{bmatrix}-0.62\\1.5\end{bmatrix}}$

STEP 3: Update ${\widehat {x}}.$

${\begin{aligned}{\widehat {x}}_{2}&={\widehat {x}}_{1}+y_{1}\\\\&={\begin{bmatrix}1.4\\0\end{bmatrix}}+{\begin{bmatrix}-0.62\\1.5\end{bmatrix}}\\\\&={\begin{bmatrix}0.78\\1.5\end{bmatrix}}\end{aligned}}$

Note that ${\widehat {x}}_{2}$ is a better representation of $x$, than ${\widehat {x}}_{1},$ however it is still not a very accurate representation. Continue with a couple more iterations.

Repeat, as needed.

In the second iteration of Iterative Refinement, we find the following:

$r_{2}={\begin{bmatrix}0.1\\0\end{bmatrix}},\qquad y_{2}={\begin{bmatrix}0.56\\-1.3\end{bmatrix}}$

Thus,

${\begin{aligned}{\widehat {x}}_{3}&={\widehat {x}}_{2}+y_{2}={\begin{bmatrix}0.78\\1.5\end{bmatrix}}+{\begin{bmatrix}0.56\\-1.3\end{bmatrix}}={\begin{bmatrix}1.3\\0.2\end{bmatrix}}\end{aligned}}$

In the third iteration we find

$r_{3}={\begin{bmatrix}0.1\\-0.3\end{bmatrix}},\qquad y_{3}={\begin{bmatrix}-0.34\\0.9\end{bmatrix}},$

which we can use to find our final approximation, ${\widehat {x}}_{4}$ as

${\begin{aligned}{\widehat {x}}_{4}&={\widehat {x}}_{3}+y_{3}={\begin{bmatrix}1.3\\0.2\end{bmatrix}}+{\begin{bmatrix}-0.34\\0.9\end{bmatrix}}={\begin{bmatrix}0.96\\1.1\end{bmatrix}}\end{aligned}}$

Recall that the true solution of the linear system is $x={\begin{bmatrix}1\\1\end{bmatrix}}.$ After completing three iterations of iterative refinement we have approximate the solution ${\widehat {x}}_{4}={\begin{bmatrix}0.96\\1.1\end{bmatrix}},$ which is significantly closer than our original approximation ${\widehat {x}}_{1}={\begin{bmatrix}1.4\\0\end{bmatrix}}.$

Part 2

Estimate the condition number $K(A)$ in Part 1 using the method discussed above, and compare this with the condition number calculated using norms.

Solution:

We can find the true condition number by calculating it using norms - in the way it is described on the condition number Wikipedia . (Keep in mind, we will still use 2 significant digits in this example as well.)

$A={\begin{bmatrix}3.9&1.6\\6.8&2.9\end{bmatrix}},\quad$ and $A^{-1}={\begin{bmatrix}6.7&-3.7\\-16&9.1\end{bmatrix}}.$

We can then solve for the condition number $K(A).$

$K(A)={\begin{Vmatrix}A\end{Vmatrix}}{\begin{Vmatrix}A^{-1}\end{Vmatrix}}={\begin{Vmatrix}{\begin{bmatrix}3.9&1.6\\6.8&2.9\end{bmatrix}}\end{Vmatrix}}{\begin{Vmatrix}{\begin{bmatrix}6.7&-3.7\\-16&9.1\end{bmatrix}}\end{Vmatrix}}={\begin{Vmatrix}9.7\end{Vmatrix}}{\begin{Vmatrix}25\end{Vmatrix}}=240.$

We can use our approximations, ${\widehat {x}}_{1}$ and $y_{1}$, found in the first iteration of Example 1 to use in the method we learned above, on Approximating the Condition Number.

${\begin{aligned}K(A)\approx {\frac {\begin{Vmatrix}{\widehat {y}}_{1}\end{Vmatrix}}{\begin{Vmatrix}{\widehat {x}}_{1}\end{Vmatrix}}}10^{t}\quad =\quad {\frac {\begin{Vmatrix}{\begin{bmatrix}-0.62\\1.5\end{bmatrix}}\end{Vmatrix}}{\begin{Vmatrix}{\begin{bmatrix}1.4\\0\end{bmatrix}}\end{Vmatrix}}}10^{2}={\frac {1.5}{1.4}}10^{2}=1.07*10^{2}=107.\end{aligned}}$

So we have found the true condition number to be $240$ and the approximation to be $107$. These numbers are both in the realm of $10^{2}$, meaning both for the true and approximate condition number, we can see there may be a loss of two digits of accuracy.

Quiz

1 Not only does iterative refinement improve ill-conditioned matrices, it is also very help in improving well-conditioned matrices.

	TRUE
	FALSE

2 Which of the following matrices might benefit from using iterative refinement?

	${\begin{bmatrix}1&2\\1.0001&2\end{bmatrix}}$
	${\begin{bmatrix}1&2\\10,001&2\end{bmatrix}}$
	${\begin{bmatrix}1&2\\10,001&20,002\end{bmatrix}}$
	${\begin{bmatrix}1&20,002\\10,001&2\end{bmatrix}}$

3

If you calculate another step of iterative refinement on Example 1, what will be the new approximation, ${\widehat {x}}_{5}={\begin{bmatrix}x_{1}\\x_{2}\end{bmatrix}}$  ?

$\displaystyle \ x_{1}=$

$\displaystyle \ x_{2}=$

4

A linear system is given by

${\begin{bmatrix}3.3330&15920&-10.333\\2.2220&16.710&9.6120\\1.5611&5.1791&1.6852\end{bmatrix}}$ ${\begin{bmatrix}x_{1}\\x_{2}\\x_{3}\end{bmatrix}}$ = ${\begin{bmatrix}15913\\28.544\\8.4254\end{bmatrix}},$ with ${\widehat {x}}_{1}$ = ${\begin{bmatrix}1.2001\\0.99991\\0.92538\end{bmatrix}},$ and $y_{1}$ = ${\begin{bmatrix}-0.20008\\8.9987\times 10^{-5}\\0.074607\end{bmatrix}}.$

Based on this info, what is the condition number using the approximation and the typical norms-multiplication method?

Using Approximation $\displaystyle \ K(A)=$

Using Norms $\displaystyle \ K(A)=$

References

Burden and Faires. Numerical Analysis, 3rd Edition. ISBN 0-87150-857-5