{"id":1166,"date":"2015-02-16T12:00:28","date_gmt":"2015-02-16T17:00:28","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=1166"},"modified":"2015-02-14T22:14:28","modified_gmt":"2015-02-15T03:14:28","slug":"iterative-refinement-for-solutions-to-linear-systems","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2015\/02\/16\/iterative-refinement-for-solutions-to-linear-systems\/","title":{"rendered":"Iterative Refinement for Solutions to Linear Systems"},"content":{"rendered":"<div class=\"content\"><!--introduction--><p>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.<\/p><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#e5240c85-349c-4b35-8fc7-824f006c03b0\">Iterative refinement<\/a><\/li><li><a href=\"#f226db51-9352-4ff9-a4d6-15b5b07b081f\">Complexity<\/a><\/li><li><a href=\"#eab96bb1-5e7a-4630-9999-7b6dc518ce81\">The residual<\/a><\/li><li><a href=\"#858cef39-84a1-4bb1-9504-c943a2fbecbc\">Inner product<\/a><\/li><li><a href=\"#fd399589-901b-403f-9edd-ee9b4b532dd5\">Example<\/a><\/li><li><a href=\"#2efbd771-6b31-4228-9c0d-22707b6a418d\">Inaccurate residual<\/a><\/li><li><a href=\"#6741b5c8-a579-44e6-a709-f449b72ee33c\">Accurate residual<\/a><\/li><li><a href=\"#eeb3d3f2-1a50-4e8b-908e-07ffb324bce8\">Iterate<\/a><\/li><li><a href=\"#6d6841f1-6dd2-4c75-b305-0432bf0e0b75\">More accurate residual<\/a><\/li><li><a href=\"#a767050c-a946-4f48-9da2-39766b07f7fa\">Further reading<\/a><\/li><li><a href=\"#ddf4faf8-9088-41ce-b4c9-088227de4a5d\">References<\/a><\/li><\/ul><\/div><h4>Iterative refinement<a name=\"e5240c85-349c-4b35-8fc7-824f006c03b0\"><\/a><\/h4><p>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$.<\/p><p>The iterative refinement algorithm is easily described.<\/p><div><ul><li>Solve $Ax = b$, saving the triangular factors.<\/li><li>Compute the residuals, $r = Ax - b$.<\/li><li>Use the triangular factors to solve $Ad = r$.<\/li><li>Subtract the correction, $x = x - d$<\/li><li>Repeat the previous three steps if desired.<\/li><\/ul><\/div><h4>Complexity<a name=\"f226db51-9352-4ff9-a4d6-15b5b07b081f\"><\/a><\/h4><p>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)$.<\/p><h4>The residual<a name=\"eab96bb1-5e7a-4630-9999-7b6dc518ce81\"><\/a><\/h4><p>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.<\/p><p>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.<\/p><p>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.<\/p><h4>Inner product<a name=\"858cef39-84a1-4bb1-9504-c943a2fbecbc\"><\/a><\/h4><p>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.<\/p><p>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.<\/p><p>In my next blog post, I will describe two MATLAB functions <tt>residual3p<\/tt> and <tt>dot3p<\/tt>.  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.<\/p><h4>Example<a name=\"fd399589-901b-403f-9edd-ee9b4b532dd5\"><\/a><\/h4><p>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<\/p><p>$$ h_{i,j} = \\frac{1}{i+j-1} $$<\/p><p>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 <tt>invhilb<\/tt> in the MATLAB <tt>elmat<\/tt> directory.  I'll choose the 8-by-8.  The elements are large, so I need a custom format to display the matrix.<\/p><pre class=\"codeinput\">   n = 8;\r\n   A = invhilb(n);\r\n   disp(sprintf(<span class=\"string\">'%8.0f %11.0f %11.0f %11.0f %11.0f %11.0f %11.0f %11.0f \\n'<\/span>,A))\r\n<\/pre><pre class=\"codeoutput\">      64       -2016       20160      -92400      221760     -288288      192192      -51480 \r\n   -2016       84672     -952560     4656960   -11642400    15567552   -10594584     2882880 \r\n   20160     -952560    11430720   -58212000   149688000  -204324120   141261120   -38918880 \r\n  -92400     4656960   -58212000   304920000  -800415000  1109908800  -776936160   216216000 \r\n  221760   -11642400   149688000  -800415000  2134440000 -2996753760  2118916800  -594594000 \r\n -288288    15567552  -204324120  1109908800 -2996753760  4249941696 -3030051024   856215360 \r\n  192192   -10594584   141261120  -776936160  2118916800 -3030051024  2175421248  -618377760 \r\n  -51480     2882880   -38918880   216216000  -594594000   856215360  -618377760   176679360 \r\n\r\n<\/pre><p>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 <tt>b<\/tt> is a column of the identity matrix.  I am hoping to get the fractions <tt>x = [1\/3 1\/4 ... 1\/9 1\/10]<\/tt>.<\/p><pre class=\"codeinput\">   b = zeros(n,1);\r\n   b(3) = 1\r\n   format <span class=\"string\">compact<\/span>\r\n   format <span class=\"string\">long<\/span> <span class=\"string\">e<\/span>\r\n   x = A\\b\r\n<\/pre><pre class=\"codeoutput\">b =\r\n     0\r\n     0\r\n     1\r\n     0\r\n     0\r\n     0\r\n     0\r\n     0\r\nx =\r\n     3.333333289789291e-01\r\n     2.499999961540004e-01\r\n     1.999999965600743e-01\r\n     1.666666635570877e-01\r\n     1.428571400209730e-01\r\n     1.249999973935827e-01\r\n     1.111111087003338e-01\r\n     9.999999775774569e-02\r\n<\/pre><p>Since I know what <tt>x<\/tt> is supposed to look like, I can just eyeball the output and see that I have only about half of the digits correct.<\/p><p>(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.)<\/p><h4>Inaccurate residual<a name=\"2efbd771-6b31-4228-9c0d-22707b6a418d\"><\/a><\/h4><p>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.<\/p><pre class=\"codeinput\">   r = A*x - b\r\n<\/pre><pre class=\"codeoutput\">r =\r\n    -9.094947017729282e-13\r\n     1.746229827404022e-10\r\n     4.656612873077393e-10\r\n     1.862645149230957e-08\r\n    -1.490116119384766e-08\r\n    -2.980232238769531e-08\r\n    -3.725290298461914e-08\r\n    -1.862645149230957e-08\r\n<\/pre><p>It's important to look at the size of the residual relative to the sizes of the matrix and the solution.<\/p><pre class=\"codeinput\">   relative_residual = norm(r)\/(norm(A)*norm(x))\r\n<\/pre><pre class=\"codeoutput\">relative_residual =\r\n     1.147025634044834e-17\r\n<\/pre><p>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.<\/p><pre class=\"codeinput\">   d = A\\r\r\n   no_help = x - d\r\n<\/pre><pre class=\"codeoutput\">d =\r\n    -1.069920014936507e-08\r\n    -9.567761339244008e-09\r\n    -8.614990592214338e-09\r\n    -7.819389121717837e-09\r\n    -7.150997084009303e-09\r\n    -6.584022612326096e-09\r\n    -6.098163254532801e-09\r\n    -5.677765952511023e-09\r\nno_help =\r\n     3.333333396781292e-01\r\n     2.500000057217617e-01\r\n     2.000000051750649e-01\r\n     1.666666713764768e-01\r\n     1.428571471719701e-01\r\n     1.250000039776053e-01\r\n     1.111111147984970e-01\r\n     1.000000034355116e-01\r\n<\/pre><h4>Accurate residual<a name=\"6741b5c8-a579-44e6-a709-f449b72ee33c\"><\/a><\/h4><p>Now I will use <tt>residual3p<\/tt>, which I intend to describe in my next blog and which employs \"triple precision\" accumulation of the inner products required for an accurate residual.<\/p><pre class=\"codeinput\">   r = residual3p(A,x,b)\r\n<\/pre><pre class=\"codeoutput\">r =\r\n    -4.045319634826683e-12\r\n     1.523381421009162e-10\r\n    -9.919851606809971e-10\r\n     2.429459300401504e-09\r\n     8.826383179894037e-09\r\n    -2.260851861279889e-08\r\n    -1.332933052822227e-08\r\n    -6.369845095832716e-09\r\n<\/pre><p>Superficially, this residual looks a lot like the previous one, but it's a lot more accurate.  The resulting correction works very well.<\/p><pre class=\"codeinput\">   d = A\\r\r\n   x = x - d\r\n<\/pre><pre class=\"codeoutput\">d =\r\n    -4.354403560053519e-09\r\n    -3.845999016894392e-09\r\n    -3.439925156187715e-09\r\n    -3.109578484769736e-09\r\n    -2.836169428940436e-09\r\n    -2.606416977917484e-09\r\n    -2.410777025186154e-09\r\n    -2.242253997222573e-09\r\nx =\r\n     3.333333333333327e-01\r\n     2.499999999999994e-01\r\n     1.999999999999995e-01\r\n     1.666666666666662e-01\r\n     1.428571428571425e-01\r\n     1.249999999999996e-01\r\n     1.111111111111108e-01\r\n     9.999999999999969e-02\r\n<\/pre><p>I've now got about 14 digits correct.  That's almost, but not quite, full double precision accuracy.<\/p><h4>Iterate<a name=\"eeb3d3f2-1a50-4e8b-908e-07ffb324bce8\"><\/a><\/h4><p>Try it again.<\/p><pre class=\"codeinput\">   r = residual3p(A,x,b)\r\n<\/pre><pre class=\"codeoutput\">r =\r\n     3.652078639504452e-12\r\n    -1.943885052924088e-10\r\n     2.523682596233812e-09\r\n    -1.359348900109580e-08\r\n     3.645651958095186e-08\r\n    -5.142027248439263e-08\r\n     3.649529745075597e-08\r\n    -1.027348206505963e-08\r\n<\/pre><p>Notice that the residual <tt>r<\/tt> is just about the same size as the previous one, even though the solution <tt>x<\/tt> is several orders of magnitude more accurate.<\/p><pre class=\"codeinput\">   d = A\\r\r\n   nice_try = x - d\r\n<\/pre><pre class=\"codeoutput\">d =\r\n     2.733263259661321e-16\r\n     2.786131033681204e-16\r\n     2.611667424188757e-16\r\n     2.527960139656094e-16\r\n     2.492795072717761e-16\r\n     2.196895809665418e-16\r\n     2.110200076421557e-16\r\n     1.983918218604762e-16\r\nnice_try =\r\n     3.333333333333324e-01\r\n     2.499999999999991e-01\r\n     1.999999999999992e-01\r\n     1.666666666666660e-01\r\n     1.428571428571422e-01\r\n     1.249999999999994e-01\r\n     1.111111111111106e-01\r\n     9.999999999999949e-02\r\n<\/pre><p>The correction changed the solution, but didn't make it appreciably more accurate.  I've reached the limits of my triple precision inner product.<\/p><h4>More accurate residual<a name=\"6d6841f1-6dd2-4c75-b305-0432bf0e0b75\"><\/a><\/h4><p>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 <tt>x<\/tt> to a <tt>sym<\/tt> so that the conversion is done exactly.<\/p><pre class=\"codeinput\"><span class=\"comment\">%  r = double(A*sym(x,'d') - b)<\/span>\r\n   r = double(A*sym(x,<span class=\"string\">'f'<\/span>) - b)\r\n<\/pre><pre class=\"codeoutput\">r =\r\n     3.652078639504452e-12\r\n    -1.943885052924088e-10\r\n     2.523682707256114e-09\r\n    -1.359348633656055e-08\r\n     3.645651780459502e-08\r\n    -5.142027803550775e-08\r\n     3.649529478622071e-08\r\n    -1.027348206505963e-08\r\n<\/pre><p>The correction just nudges the last two digits.<\/p><pre class=\"codeinput\">   d = A\\r\r\n   x = x - d\r\n<\/pre><pre class=\"codeoutput\">d =\r\n    -6.846375178532078e-16\r\n    -5.828670755614817e-16\r\n    -5.162536953886886e-16\r\n    -4.533410583885285e-16\r\n    -3.965082139594863e-16\r\n    -3.747002624289523e-16\r\n    -3.392348053271079e-16\r\n    -3.136379972458518e-16\r\nx =\r\n     3.333333333333333e-01\r\n     2.500000000000000e-01\r\n     2.000000000000000e-01\r\n     1.666666666666667e-01\r\n     1.428571428571429e-01\r\n     1.250000000000000e-01\r\n     1.111111111111111e-01\r\n     1.000000000000000e-01\r\n<\/pre><p>Now, with a very accurate residual, the elements I get in <tt>x<\/tt> are the floating point numbers closest to the fractions in the Hilbert matrix. That's the best I can do.<\/p><h4>Further reading<a name=\"a767050c-a946-4f48-9da2-39766b07f7fa\"><\/a><\/h4><p>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.<\/p><h4>References<a name=\"ddf4faf8-9088-41ce-b4c9-088227de4a5d\"><\/a><\/h4><p>Cleve Moler, Iterative Refinement in Floating Point, <i>Journal of the ACM<\/i> 14, pp. 316-321, 1967, <a href=\"http:\/\/dl.acm.org\/citation.cfm?doid=321386.321394\">&lt;http:\/\/dl.acm.org\/citation.cfm?doid=321386.321394<\/a>&gt;<\/p><p>Nicholas Higham, <i>Accuracy and Stability of Numerical Algorithms (2nd ed)<\/i>, SIAM, 2002, <a href=\"http:\/\/epubs.siam.org\/doi\/book\/10.1137\/1.9780898718027\">&lt;http:\/\/epubs.siam.org\/doi\/book\/10.1137\/1.9780898718027<\/a>&gt;<\/p><p>James Demmel, Yozo Hida, William Kahan, Xiaoye S. Li, Sonil Mukherjee, and E. Jason Riedy, Error Bounds from Extra-Precise Iterative Refinement, <i>ACM Transactions on Mathematical Software<\/i> 32, pp. 325-351, 2006, <a href=\"http:\/\/www.cs.berkeley.edu\/~wkahan\/p325-demmel.pdf\">&lt;http:\/\/www.cs.berkeley.edu\/~wkahan\/p325-demmel.pdf<\/a>&gt;<\/p><script language=\"JavaScript\"> <!-- \r\n    function grabCode_156eab9b18d8415cbed32b8c5f1d2594() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='156eab9b18d8415cbed32b8c5f1d2594 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 156eab9b18d8415cbed32b8c5f1d2594';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        copyright = 'Copyright 2015 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add copyright line at the bottom if specified.\r\n        if (copyright.length > 0) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\\n');\r\n\r\n        d.title = title + ' (MATLAB code)';\r\n        d.close();\r\n    }   \r\n     --> <\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_156eab9b18d8415cbed32b8c5f1d2594()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n      the MATLAB code <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; R2014b<br><\/p><\/div><!--\r\n156eab9b18d8415cbed32b8c5f1d2594 ##### SOURCE BEGIN #####\r\n%% Iterative Refinement for Solutions to Linear Systems\r\n% Iterative refinement is a technique introduced by Wilkinson for reducing\r\n% the roundoff error produced during the solution of simultaneous linear\r\n% equations.  Higher precision arithmetic is required for the calculation\r\n% of the residuals.\r\n\r\n%% Iterative refinement\r\n% The first research paper I ever published, in 1967, was titled \"Iterative\r\n% Refinement in Floating Point\".  It was an analysis of a technique introduced\r\n% by J. H. Wilkinson almost 20 years earlier for a post processing cleanup\r\n% that reduces the roundoff error generated when Gaussian elimination or a\r\n% similar process is used to solve a system of simultaneous linear equations,\r\n% $Ax = b$.\r\n\r\n%%\r\n% The iterative refinement algorithm is easily described.\r\n%\r\n% * Solve $Ax = b$, saving the triangular factors.\r\n% * Compute the residuals, $r = Ax - b$.\r\n% * Use the triangular factors to solve $Ad = r$.\r\n% * Subtract the correction, $x = x - d$\r\n% * Repeat the previous three steps if desired.\r\n\r\n%% Complexity\r\n% Almost all of the work is in the first step, which can be thought of as\r\n% producing triangular factors such as $L$ and $U$ so that $A = LU$ while\r\n% solving $LUx = b$.  For a matrix of order $n$ the computational complexity\r\n% of this step is $O(n^3)$.  Saving the factorization reduces the complexity\r\n% of the remaining refinement steps to something much less, only $O(n^2)$.  \r\n\r\n%% The residual\r\n% By the early 1960s we had learned from Wilkinson that if a system of\r\n% simultaneous linear equations is solved by a process like Gaussian\r\n% elimination or Cholesky factorization, the residual will always be order\r\n% roundoff error, relative to the matrix and the computed solution, even if\r\n% the system is nearly singular.  This is both good news and bad news. \r\n\r\n%%\r\n% The good news is that $Ax$ is always close to $b$.  This says that the\r\n% computed solution always comes close to solving the equations, even\r\n% though $x$ might not be close to the theoretical correct solution,\r\n% $A^{-1}b$.  The pitcher always puts the ball where the batter can hit it,\r\n% even though that might not be in the strike zone.\r\n\r\n%%\r\n% The bad news is that it is delicate to compute the residual accurately.\r\n% If the same precision arithmetic is used to compute the residual that was\r\n% used to solve the system, the roundoff error involved in computing $r$\r\n% will be almost comparable to the effect of the roundoff error present in\r\n% $x$, so the correction has little chance of being useful.\r\n\r\n%% Inner product\r\n% We need to use higher precision arithmetic while computing the residual.\r\n% Each component of the residual involves a sum of products and then one\r\n% final subtraction.  The exact product of two numbers with a certain\r\n% precision has twice that precision.  With the computers that Wilkinson used,\r\n% and that I used early in my career, we had access to the full results of\r\n% multiplications.  We were able to write inner product routines that\r\n% accumulated the sums with twice the working precision.\r\n\r\n%%\r\n% But it is not easy to write the accumulated inner product routine in\r\n% modern, portable, machine independent software.  It was not easy in Fortran.\r\n% It is not easy in MATLAB.  The original specification of the BLAS, the\r\n% Basic Linear Algebra Subroutines, was deliberately silent on the matter.\r\n% Subsequent proposals for extensions of the BLAS have introduced mixed\r\n% precision, but these extensions have not been widely adopted.  So, the\r\n% key tool we need to implement iterative refinement has not been available.\r\n\r\n%%\r\n% In my next blog post, I will describe two MATLAB functions |residual3p|\r\n% and |dot3p|.  They provide enough of what I call \"triple precision\"\r\n% arithmetic to produce an accumulated inner product.  It's a hack, but\r\n% it works well enough to illustrate iterative refinement.\r\n\r\n%% Example\r\n% My example involves perhaps the world's most famous badly conditioned matrix,\r\n% the Hilbert matrix.  I won't begin with the Hilbert matrix itself because\r\n% its elements are the fractions\r\n%\r\n% $$ h_{i,j} = \\frac{1}{i+j-1} $$\r\n%\r\n% Many of these fractions can't be represented exactly in floating point,\r\n% so I would have roundoff error before even getting started with the\r\n% elimination.  Fortunately, the elements of the inverse Hilbert matrix are\r\n% integers that can be readily generated.  There is a function |invhilb| in\r\n% the MATLAB |elmat| directory.  I'll choose the 8-by-8.  The elements are\r\n% large, so I need a custom format to display the matrix.\r\n\r\n   n = 8;\r\n   A = invhilb(n);\r\n   disp(sprintf('%8.0f %11.0f %11.0f %11.0f %11.0f %11.0f %11.0f %11.0f \\n',A))\r\n\r\n%%\r\n% I am going to try to compute the third column of the inverse of this inverse,\r\n% which is a column of the Hilbert matrix itself.  The right hand side |b|\r\n% is a column of the identity matrix.  I am hoping to get the fractions\r\n% |x = [1\/3 1\/4 ... 1\/9 1\/10]|. \r\n\r\n   b = zeros(n,1);\r\n   b(3) = 1\r\n   format compact\r\n   format long e\r\n   x = A\\b\r\n\r\n%%\r\n% Since I know what |x| is supposed to look like, I can just eyeball the\r\n% output and see that I have only about half of the digits correct.\r\n\r\n%%\r\n% (I used backslash to solve the system.  My matrix happens to be symmetric\r\n% and positive definite, so the elimination algorithm involves the Cholesky\r\n% factorization.  But I'm going to be extravagant, ignore the complexity\r\n% considerations, and not save the triangular factor.)\r\n\r\n%% Inaccurate residual\r\n% Here's my first crack at the residual.  I won't do anything special about\r\n% the precision this time; I'll just use an ordinary MATLAB statement.\r\n\r\n   r = A*x - b\r\n\r\n%%\r\n% It's important to look at the size of the residual relative to the sizes\r\n% of the matrix and the solution.\r\n\r\n   relative_residual = norm(r)\/(norm(A)*norm(x))\r\n\r\n%%\r\n% The elements in this computed residual are the right order of magnitude,\r\n% that is roundoff error,\r\n% but, since I didn't use any extra precision, they are not accurate enough\r\n% to provide a useful correction.\r\n\r\n   d = A\\r\r\n   no_help = x - d\r\n\r\n%% Accurate residual\r\n% Now I will use |residual3p|, which I intend to describe in my next blog\r\n% and which employs \"triple precision\" accumulation of the inner products\r\n% required for an accurate residual.\r\n\r\n   r = residual3p(A,x,b)\r\n\r\n%%\r\n% Superficially, this residual looks a lot like the previous one,\r\n% but it's a lot more accurate.  The resulting correction works very well.\r\n\r\n   d = A\\r\r\n   x = x - d\r\n\r\n%%\r\n% I've now got about 14 digits correct.  That's almost, but not quite, full\r\n% double precision accuracy.\r\n\r\n%% Iterate\r\n% Try it again.\r\n\r\n   r = residual3p(A,x,b)\r\n\r\n%%\r\n% Notice that the residual |r| is just about the same size as the previous one,\r\n% even though the solution |x| is several orders of magnitude more accurate.\r\n\r\n   d = A\\r\r\n   nice_try = x - d\r\n%%\r\n% The correction changed the solution, but didn't make it appreciably more\r\n% accurate.  I've reached the limits of my triple precision inner product.\r\n\r\n%% More accurate residual\r\n% Bring in the big guns, the Symbolic Math Toolbox, to compute a very accurate\r\n% residual.  It is important to use either the 'f' or the 'd' option when\r\n% converting |x| to a |sym| so that the conversion is done exactly.\r\n\r\n%  r = double(A*sym(x,'d') - b)\r\n   r = double(A*sym(x,'f') - b)\r\n\r\n%%\r\n% The correction just nudges the last two digits.\r\n\r\n   d = A\\r\r\n   x = x - d\r\n\r\n%%\r\n% Now, with a very accurate residual, the elements I get in |x| are the\r\n% floating point numbers closest to the fractions in the Hilbert matrix.\r\n% That's the best I can do.\r\n\r\n%% Further reading\r\n% There is more to say about iterative refinement.  See Nick Higham's SIAM\r\n% book and, especially, the 2006 TOMS paper by Demmel, Kahan and their\r\n% Berkeley colleagues.  A preprint is available from Kahan's web site.\r\n\r\n%% References\r\n% Cleve Moler, Iterative Refinement in Floating Point,\r\n% _Journal of the ACM_ 14, pp. 316-321, 1967,\r\n% <http:\/\/dl.acm.org\/citation.cfm?doid=321386.321394\r\n% http:\/\/dl.acm.org\/citation.cfm?doid=321386.321394>\r\n\r\n%%\r\n% Nicholas Higham, _Accuracy and Stability of Numerical Algorithms (2nd ed)_,\r\n% SIAM, 2002, \r\n% <http:\/\/epubs.siam.org\/doi\/book\/10.1137\/1.9780898718027\r\n% http:\/\/epubs.siam.org\/doi\/book\/10.1137\/1.9780898718027>\r\n\r\n%%\r\n% James Demmel, Yozo Hida, William Kahan, Xiaoye S. Li, Sonil Mukherjee,\r\n% and E. Jason Riedy, Error Bounds from Extra-Precise Iterative Refinement,\r\n% _ACM Transactions on Mathematical Software_ 32, pp. 325-351, 2006,\r\n% <http:\/\/www.cs.berkeley.edu\/~wkahan\/p325-demmel.pdf\r\n% http:\/\/www.cs.berkeley.edu\/~wkahan\/p325-demmel.pdf>\r\n\r\n##### SOURCE END ##### 156eab9b18d8415cbed32b8c5f1d2594\r\n-->","protected":false},"excerpt":{"rendered":"<!--introduction--><p>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.... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2015\/02\/16\/iterative-refinement-for-solutions-to-linear-systems\/\">read more >><\/a><\/p>","protected":false},"author":78,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[11,4,6,16,7],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/1166"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/users\/78"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/comments?post=1166"}],"version-history":[{"count":4,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/1166\/revisions"}],"predecessor-version":[{"id":1168,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/1166\/revisions\/1168"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=1166"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=1166"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=1166"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}