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.

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>

Get the MATLAB code

Published with MATLAB® R2014b