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:

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{align} \begin{Vmatrix}\widehat{y}\end{Vmatrix} \approx \begin{Vmatrix}x - \widehat{x}\end{Vmatrix} &= \begin{Vmatrix}A^{-1}r\end{Vmatrix} \le \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{align}$$

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

Example
Consider the following ill-conditioned linear system
 * $$\begin{align}

3.9x_1 + 1.6x_2 &= 5.5\\ 6.8x_1 + 2.9x_2 &= 9.7 \end{align}$$

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{align}

\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{align}$$

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{align}

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{align}$$

STEP 2: Now solve $$Ay_1 = r_1 $$ using rounded Gaussian elimination to find the $$y$$ vector.


 * $$\begin{align}

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{align}$$


 * $$\Rightarrow

y_1 = \begin{bmatrix} -0.62\\ 1.5 \end{bmatrix}$$

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


 * $$\begin{align}

\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{align}$$

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{align}

\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{align}$$

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{align}

\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{align}$$

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{align}

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{align}$$

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
Not only does iterative refinement improve ill-conditioned matrices, it is also very help in improving well-conditioned matrices. +TRUE -FALSE
 * type=""}

{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}$$
 * type="[]"}

{ 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 = $$ { 0.96 } $$ \displaystyle \ x_2 = $$ { 1.1 }
 * type="{}"}

{ 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}. $$
 * type="{}"}

Based on this info, what is the condition number using the approximation and the typical norms-multiplication method? Using Approximation $$ \displaystyle \ K(A) = $$ { 16672 } Using Norms $$ \displaystyle \ K(A) = $$ { 15999 }