# Iterative Refinement for Solutions to Linear Systems

Iterative refinement is a technique introduced by Wilkinson for reducing the roundoff error produced during the solution of simultaneous linear equations. Higher precision arithmetic is required for the calculation of the residuals.

### Contents

#### Iterative refinement

The first research paper I ever published, in 1967, was titled "Iterative Refinement in Floating Point". It was an analysis of a technique introduced by J. H. Wilkinson almost 20 years earlier for a post processing cleanup that reduces the roundoff error generated when Gaussian elimination or a similar process is used to solve a system of simultaneous linear equations, $Ax = b$.

The iterative refinement algorithm is easily described.

- Solve $Ax = b$, saving the triangular factors.
- Compute the residuals, $r = Ax - b$.
- Use the triangular factors to solve $Ad = r$.
- Subtract the correction, $x = x - d$
- Repeat the previous three steps if desired.

#### Complexity

Almost all of the work is in the first step, which can be thought of as producing triangular factors such as $L$ and $U$ so that $A = LU$ while solving $LUx = b$. For a matrix of order $n$ the computational complexity of this step is $O(n^3)$. Saving the factorization reduces the complexity of the remaining refinement steps to something much less, only $O(n^2)$.

#### The residual

By the early 1960s we had learned from Wilkinson that if a system of simultaneous linear equations is solved by a process like Gaussian elimination or Cholesky factorization, the residual will always be order roundoff error, relative to the matrix and the computed solution, even if the system is nearly singular. This is both good news and bad news.

The good news is that $Ax$ is always close to $b$. This says that the computed solution always comes close to solving the equations, even though $x$ might not be close to the theoretical correct solution, $A^{-1}b$. The pitcher always puts the ball where the batter can hit it, even though that might not be in the strike zone.

The bad news is that it is delicate to compute the residual accurately. If the same precision arithmetic is used to compute the residual that was used to solve the system, the roundoff error involved in computing $r$ will be almost comparable to the effect of the roundoff error present in $x$, so the correction has little chance of being useful.

#### Inner product

We need to use higher precision arithmetic while computing the residual. Each component of the residual involves a sum of products and then one final subtraction. The exact product of two numbers with a certain precision has twice that precision. With the computers that Wilkinson used, and that I used early in my career, we had access to the full results of multiplications. We were able to write inner product routines that accumulated the sums with twice the working precision.

But it is not easy to write the accumulated inner product routine in modern, portable, machine independent software. It was not easy in Fortran. It is not easy in MATLAB. The original specification of the BLAS, the Basic Linear Algebra Subroutines, was deliberately silent on the matter. Subsequent proposals for extensions of the BLAS have introduced mixed precision, but these extensions have not been widely adopted. So, the key tool we need to implement iterative refinement has not been available.

In my next blog post, I will describe two MATLAB functions `residual3p` and `dot3p`. They provide enough of what I call "triple precision" arithmetic to produce an accumulated inner product. It's a hack, but it works well enough to illustrate iterative refinement.

#### Example

My example involves perhaps the world's most famous badly conditioned matrix, the Hilbert matrix. I won't begin with the Hilbert matrix itself because its elements are the fractions

$$ h_{i,j} = \frac{1}{i+j-1} $$

Many of these fractions can't be represented exactly in floating point, so I would have roundoff error before even getting started with the elimination. Fortunately, the elements of the inverse Hilbert matrix are integers that can be readily generated. There is a function `invhilb` in the MATLAB `elmat` directory. I'll choose the 8-by-8. The elements are large, so I need a custom format to display the matrix.

```
n = 8;
A = invhilb(n);
disp(sprintf('%8.0f %11.0f %11.0f %11.0f %11.0f %11.0f %11.0f %11.0f \n',A))
```

64 -2016 20160 -92400 221760 -288288 192192 -51480 -2016 84672 -952560 4656960 -11642400 15567552 -10594584 2882880 20160 -952560 11430720 -58212000 149688000 -204324120 141261120 -38918880 -92400 4656960 -58212000 304920000 -800415000 1109908800 -776936160 216216000 221760 -11642400 149688000 -800415000 2134440000 -2996753760 2118916800 -594594000 -288288 15567552 -204324120 1109908800 -2996753760 4249941696 -3030051024 856215360 192192 -10594584 141261120 -776936160 2118916800 -3030051024 2175421248 -618377760 -51480 2882880 -38918880 216216000 -594594000 856215360 -618377760 176679360

I am going to try to compute the third column of the inverse of this inverse, which is a column of the Hilbert matrix itself. The right hand side `b` is a column of the identity matrix. I am hoping to get the fractions `x = [1/3 1/4 ... 1/9 1/10]`.

b = zeros(n,1); b(3) = 1 format compact format long e x = A\b

b = 0 0 1 0 0 0 0 0 x = 3.333333289789291e-01 2.499999961540004e-01 1.999999965600743e-01 1.666666635570877e-01 1.428571400209730e-01 1.249999973935827e-01 1.111111087003338e-01 9.999999775774569e-02

Since I know what `x` is supposed to look like, I can just eyeball the output and see that I have only about half of the digits correct.

(I used backslash to solve the system. My matrix happens to be symmetric and positive definite, so the elimination algorithm involves the Cholesky factorization. But I'm going to be extravagant, ignore the complexity considerations, and not save the triangular factor.)

#### Inaccurate residual

Here's my first crack at the residual. I won't do anything special about the precision this time; I'll just use an ordinary MATLAB statement.

r = A*x - b

r = -9.094947017729282e-13 1.746229827404022e-10 4.656612873077393e-10 1.862645149230957e-08 -1.490116119384766e-08 -2.980232238769531e-08 -3.725290298461914e-08 -1.862645149230957e-08

It's important to look at the size of the residual relative to the sizes of the matrix and the solution.

relative_residual = norm(r)/(norm(A)*norm(x))

relative_residual = 1.147025634044834e-17

The elements in this computed residual are the right order of magnitude, that is roundoff error, but, since I didn't use any extra precision, they are not accurate enough to provide a useful correction.

d = A\r no_help = x - d

d = -1.069920014936507e-08 -9.567761339244008e-09 -8.614990592214338e-09 -7.819389121717837e-09 -7.150997084009303e-09 -6.584022612326096e-09 -6.098163254532801e-09 -5.677765952511023e-09 no_help = 3.333333396781292e-01 2.500000057217617e-01 2.000000051750649e-01 1.666666713764768e-01 1.428571471719701e-01 1.250000039776053e-01 1.111111147984970e-01 1.000000034355116e-01

#### Accurate residual

Now I will use `residual3p`, which I intend to describe in my next blog and which employs "triple precision" accumulation of the inner products required for an accurate residual.

r = residual3p(A,x,b)

r = -4.045319634826683e-12 1.523381421009162e-10 -9.919851606809971e-10 2.429459300401504e-09 8.826383179894037e-09 -2.260851861279889e-08 -1.332933052822227e-08 -6.369845095832716e-09

Superficially, this residual looks a lot like the previous one, but it's a lot more accurate. The resulting correction works very well.

d = A\r x = x - d

d = -4.354403560053519e-09 -3.845999016894392e-09 -3.439925156187715e-09 -3.109578484769736e-09 -2.836169428940436e-09 -2.606416977917484e-09 -2.410777025186154e-09 -2.242253997222573e-09 x = 3.333333333333327e-01 2.499999999999994e-01 1.999999999999995e-01 1.666666666666662e-01 1.428571428571425e-01 1.249999999999996e-01 1.111111111111108e-01 9.999999999999969e-02

I've now got about 14 digits correct. That's almost, but not quite, full double precision accuracy.

#### Iterate

Try it again.

r = residual3p(A,x,b)

r = 3.652078639504452e-12 -1.943885052924088e-10 2.523682596233812e-09 -1.359348900109580e-08 3.645651958095186e-08 -5.142027248439263e-08 3.649529745075597e-08 -1.027348206505963e-08

Notice that the residual `r` is just about the same size as the previous one, even though the solution `x` is several orders of magnitude more accurate.

d = A\r nice_try = x - d

d = 2.733263259661321e-16 2.786131033681204e-16 2.611667424188757e-16 2.527960139656094e-16 2.492795072717761e-16 2.196895809665418e-16 2.110200076421557e-16 1.983918218604762e-16 nice_try = 3.333333333333324e-01 2.499999999999991e-01 1.999999999999992e-01 1.666666666666660e-01 1.428571428571422e-01 1.249999999999994e-01 1.111111111111106e-01 9.999999999999949e-02

The correction changed the solution, but didn't make it appreciably more accurate. I've reached the limits of my triple precision inner product.

#### More accurate residual

Bring in the big guns, the Symbolic Math Toolbox, to compute a very accurate residual. It is important to use either the 'f' or the 'd' option when converting `x` to a `sym` so that the conversion is done exactly.

% r = double(A*sym(x,'d') - b) r = double(A*sym(x,'f') - b)

r = 3.652078639504452e-12 -1.943885052924088e-10 2.523682707256114e-09 -1.359348633656055e-08 3.645651780459502e-08 -5.142027803550775e-08 3.649529478622071e-08 -1.027348206505963e-08

The correction just nudges the last two digits.

d = A\r x = x - d

d = -6.846375178532078e-16 -5.828670755614817e-16 -5.162536953886886e-16 -4.533410583885285e-16 -3.965082139594863e-16 -3.747002624289523e-16 -3.392348053271079e-16 -3.136379972458518e-16 x = 3.333333333333333e-01 2.500000000000000e-01 2.000000000000000e-01 1.666666666666667e-01 1.428571428571429e-01 1.250000000000000e-01 1.111111111111111e-01 1.000000000000000e-01

Now, with a very accurate residual, the elements I get in `x` are the floating point numbers closest to the fractions in the Hilbert matrix. That's the best I can do.

#### Further reading

There is more to say about iterative refinement. See Nick Higham's SIAM book and, especially, the 2006 TOMS paper by Demmel, Kahan and their Berkeley colleagues. A preprint is available from Kahan's web site.

#### References

Cleve Moler, Iterative Refinement in Floating Point, *Journal of the ACM* 14, pp. 316-321, 1967, <http://dl.acm.org/citation.cfm?doid=321386.321394>

Nicholas Higham, *Accuracy and Stability of Numerical Algorithms (2nd ed)*, SIAM, 2002, <http://epubs.siam.org/doi/book/10.1137/1.9780898718027>

James Demmel, Yozo Hida, William Kahan, Xiaoye S. Li, Sonil Mukherjee, and E. Jason Riedy, Error Bounds from Extra-Precise Iterative Refinement, *ACM Transactions on Mathematical Software* 32, pp. 325-351, 2006, <http://www.cs.berkeley.edu/~wkahan/p325-demmel.pdf>

**Category:**- Algorithms,
- History,
- Matrices,
- Numerical Analysis,
- Precision

## Comments

To leave a comment, please click here to sign in to your MathWorks Account or create a new one.