{"id":2528,"date":"2017-06-07T15:55:24","date_gmt":"2017-06-07T20:55:24","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=2528"},"modified":"2017-06-07T15:55:24","modified_gmt":"2017-06-07T20:55:24","slug":"hilbert-matrices-2","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2017\/06\/07\/hilbert-matrices-2\/","title":{"rendered":"Hilbert Matrices"},"content":{"rendered":"<div class=\"content\"><!--introduction--><p>I first encountered the Hilbert matrix when I was doing individual studies under <a href=\"https:\/\/blogs.mathworks.com\/cleve\/2015\/12\/07\">Professor John Todd<\/a> at Caltech in 1960. It has been part of my professional life ever since.<\/p><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#03acbff0-1993-452a-897d-08f1b30fe0eb\">David Hilbert<\/a><\/li><li><a href=\"#817234d7-fc77-4212-9104-d7096909f19f\">Least squares<\/a><\/li><li><a href=\"#449edfd1-6cb5-4d34-b135-38c0f7f1a302\">Properties<\/a><\/li><li><a href=\"#f8cfa34c-baaf-4046-94a5-b32e0a20598f\">hilb<\/a><\/li><li><a href=\"#92415b66-c11f-4218-b2e2-0b8b9692b045\">Inverse<\/a><\/li><li><a href=\"#994b20ac-6682-48db-9147-85a328bdf9b9\">invhilv<\/a><\/li><li><a href=\"#b73a7b29-b17f-49ef-b5a7-6f45efd954a2\">Rookie experiment<\/a><\/li><li><a href=\"#5f919fa3-707a-4589-ae6e-eea641220b58\">Deeper explanation<\/a><\/li><\/ul><\/div><h4>David Hilbert<a name=\"03acbff0-1993-452a-897d-08f1b30fe0eb\"><\/a><\/h4><p>Around the turn of the 20th century, David Hilbert was the world's most famous mathematician.  He introduced the matrix that now bears his name in a paper in 1895.  The elements of the matrix, which are reciprocals of consecutive positive integers, are constant along the antidiagonals.<\/p><p>$$ h_{i,j} = \\frac{1}{i+j-1}, \\ \\ i,j = 1:n $$<\/p><pre class=\"codeinput\">    format <span class=\"string\">rat<\/span>\r\n    H5 = hilb(5)\r\n<\/pre><pre class=\"codeoutput\">H5 =\r\n       1              1\/2            1\/3            1\/4            1\/5     \r\n       1\/2            1\/3            1\/4            1\/5            1\/6     \r\n       1\/3            1\/4            1\/5            1\/6            1\/7     \r\n       1\/4            1\/5            1\/6            1\/7            1\/8     \r\n       1\/5            1\/6            1\/7            1\/8            1\/9     \r\n<\/pre><pre class=\"codeinput\">    latex(sym(H5));\r\n<\/pre><p>$$ H_5 = \\left(\\begin{array}{ccccc}\r\n1 &amp; \\frac{1}{2} &amp; \\frac{1}{3} &amp; \\frac{1}{4} &amp; \\frac{1}{5}\\\\\r\n\\frac{1}{2} &amp; \\frac{1}{3} &amp; \\frac{1}{4} &amp; \\frac{1}{5} &amp; \\frac{1}{6}\\\\\r\n\\frac{1}{3} &amp; \\frac{1}{4} &amp; \\frac{1}{5} &amp; \\frac{1}{6} &amp; \\frac{1}{7}\\\\\r\n\\frac{1}{4} &amp; \\frac{1}{5} &amp; \\frac{1}{6} &amp; \\frac{1}{7} &amp; \\frac{1}{8}\\\\\r\n\\frac{1}{5} &amp; \\frac{1}{6} &amp; \\frac{1}{7} &amp; \\frac{1}{8} &amp; \\frac{1}{9}\r\n\\end{array}\\right) $$<\/p><p>Here's a picture.  The continuous surface generated is very smooth.<\/p><pre class=\"codeinput\">    H12 = hilb(12);\r\n    surf(log(H12))\r\n    view(60,60)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/hilbert_blog_01.png\" alt=\"\"> <h4>Least squares<a name=\"817234d7-fc77-4212-9104-d7096909f19f\"><\/a><\/h4><p>Hilbert was interested in this matrix because it comes up in the least squares approximation of a continuous function on the unit interval by polynomials, expressed in the conventional basis as linear combinations of monomials.<\/p><p>$$ p(x) = \\sum_{j=1}^n c_j x^{j-1} $$<\/p><p>The coefficient matrix for the normal equations has elements<\/p><p>$$ \\int_0^1 x^{i+j-2} dx \\ = \\ \\frac{1}{i+j-1} $$<\/p><h4>Properties<a name=\"449edfd1-6cb5-4d34-b135-38c0f7f1a302\"><\/a><\/h4><p>A Hilbert matrix has many useful properties.<\/p><div><ul><li>Symmetric.<\/li><li>Positive definite.<\/li><li>Hankel, $a_{i,j}$ is a function of $i+j$.<\/li><li>Cauchy, $a_{i,j} = 1\/(x_i - y_j)$.<\/li><li>Nearly singular.<\/li><\/ul><\/div><p>Each column of a Hilbert matrix is nearly a multiple of the other columns.  So the columns are nearly linearly dependent and the matrix is close to, but not exactly, singular.<\/p><h4>hilb<a name=\"f8cfa34c-baaf-4046-94a5-b32e0a20598f\"><\/a><\/h4><p>MATLAB has always had functions <tt>hilb<\/tt> and <tt>invhilb<\/tt> that compute the Hilbert matrix and its inverse.  The body of <tt>hilb<\/tt> is now only two lines.<\/p><pre class=\"language-matlab\">J = 1:n;\r\nH = 1.\/(J'+J-1);\r\n<\/pre><p>We often cite this as a good example of <i>singleton expansion<\/i>. A column vector is added to a row vector to produce a matrix, then a scalar is subtracted from that matrix, and finally the reciprocals of the elements produce the result.<\/p><h4>Inverse<a name=\"92415b66-c11f-4218-b2e2-0b8b9692b045\"><\/a><\/h4><p>It is possible to express the elements of the inverse of the Hilbert matrix in terms of binomial coefficients.  For reasons that I've now forgotten, I always use $T$ for $H^{-1}$.<\/p><p>$$ t_{i,j} = (-1)^{i+j} (i+j-1) {{n+i-1} \\choose {n-j}}\r\n {{n+j-1} \\choose {n-i}} {{i+j-2} \\choose {i-1}}^2 $$<\/p><p>The elements of the inverse Hilbert matrix are integers, but they are <i>large<\/i> integers.  The largest ones are on the diagonal.  For order 13, the largest element is<\/p><p>$$ (T_{13})_{9,9} \\ = \\ 100863567447142500 $$<\/p><p>This is over $10^{17}$ and is larger than double precision <tt>flintmax<\/tt>.<\/p><pre class=\"codeinput\">    format <span class=\"string\">longe<\/span>\r\n    flintmax = flintmax\r\n<\/pre><pre class=\"codeoutput\">flintmax =\r\n     9.007199254740992e+15\r\n<\/pre><p>So, it is possible to represent the largest elements exactly only if $n \\le 12$.<\/p><p>The HELP entry for <tt>invhilb<\/tt> includes a sentence inherited from my original Fortran MATLAB.<\/p><pre>  The result is exact for  N  less than about 15.<\/pre><p>Now that's misleading.  It should say<\/p><pre>  The result is exact for  N  &lt;= 12.<\/pre><p>(I'm filing a bug report.)<\/p><h4>invhilv<a name=\"994b20ac-6682-48db-9147-85a328bdf9b9\"><\/a><\/h4><p>The body of <tt>invhilb<\/tt> begins by setting <tt>p<\/tt> to the order <tt>n<\/tt>. The doubly nested <tt>for<\/tt> loops then evaluate the binomial coefficient formula recursively, avoiding unnecessary integer overflow.<\/p><pre class=\"codeinput\">    dbtype <span class=\"string\">invhilb<\/span> <span class=\"string\">18:28<\/span>\r\n<\/pre><pre class=\"codeoutput\">\r\n18    p = n;\r\n19    for i = 1:n\r\n20        r = p*p;\r\n21        H(i,i) = r\/(2*i-1);\r\n22        for j = i+1:n\r\n23            r = -((n-j+1)*r*(n+j-1))\/(j-1)^2;\r\n24            H(i,j) = r\/(i+j-1);\r\n25            H(j,i) = r\/(i+j-1);\r\n26        end\r\n27        p = ((n-i)*p*(n+i))\/(i^2);\r\n28    end\r\n<\/pre><p>I first programmed this algorithm in machine language for the Burroughs 205 Datatron at Caltech almost 60 years ago.  I was barely out of my teens.<\/p><p>Here's the result for <tt>n = 6<\/tt>.<\/p><pre class=\"codeinput\">    format <span class=\"string\">short<\/span>\r\n    T6 = invhilb(6)\r\n<\/pre><pre class=\"codeoutput\">T6 =\r\n          36        -630        3360       -7560        7560       -2772\r\n        -630       14700      -88200      211680     -220500       83160\r\n        3360      -88200      564480    -1411200     1512000     -582120\r\n       -7560      211680    -1411200     3628800    -3969000     1552320\r\n        7560     -220500     1512000    -3969000     4410000    -1746360\r\n       -2772       83160     -582120     1552320    -1746360      698544\r\n<\/pre><p>A checkerboard sign pattern with large integers in the inverse cancels the smooth surface of the Hilbert matrix itself.<\/p><pre class=\"codeinput\">    T12 = invhilb(12);\r\n    spy(T12 &gt; 0)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/hilbert_blog_02.png\" alt=\"\"> <p>A log scale mitigates the jaggedness.<\/p><pre class=\"codeinput\">    surf(sign(T12).*log(abs(T12)))\r\n    view(60,60)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/hilbert_blog_03.png\" alt=\"\"> <h4>Rookie experiment<a name=\"b73a7b29-b17f-49ef-b5a7-6f45efd954a2\"><\/a><\/h4><p>Now using MATLAB, I am going to repeat the experiment that I did on the Burroughs 205 when I was still a rookie.  I had just written my first program that used Gaussian elimination to invert matrices. I proceeded to test it by inverting Hilbert matrices and comparing the results with the exact inverses.  (Today's code uses this utility function that picks out the largest element in a matrix.)<\/p><pre class=\"codeinput\">    maxabs = @(X) max(double(abs(X(:))));\r\n<\/pre><p>Here is <tt>n = 10<\/tt>.<\/p><pre class=\"codeinput\">    n = 10\r\n    H = hilb(n);\r\n    X = inv(H);        <span class=\"comment\">% Computed inverse<\/span>\r\n    T = invhilb(n);    <span class=\"comment\">% Theoretical inverse<\/span>\r\n    E = X - T;         <span class=\"comment\">% The error<\/span>\r\n    err = maxabs(E)\r\n<\/pre><pre class=\"codeoutput\">n =\r\n    10\r\nerr =\r\n   5.0259e+08\r\n<\/pre><p>At first I might have said:<\/p><p><i>Wow!  The error is $10^8$.  That's a pretty big error. Can I trust my matrix inversion code?   What went wrong?<\/i><\/p><p>But I knew the elements of the inverse are huge. We should be looking at <i>relative<\/i> error.<\/p><pre class=\"codeinput\">   relerr = maxabs(E)\/maxabs(T)\r\n<\/pre><pre class=\"codeoutput\">relerr =\r\n   1.4439e-04\r\n<\/pre><p><i>OK.  The relative error is $10^{-4}$.  That still seems like a lot.<\/i><\/p><p>I knew that the Hilbert matrix is nearly singular. That's why I was using it.  John Todd was one of the first people to write about condition numbers.  An error estimate that involves nearness to singularity and the floating point accuracy would be expressed today by<\/p><pre class=\"codeinput\">   esterr = cond(H)*eps\r\n<\/pre><pre class=\"codeoutput\">esterr =\r\n    0.0036\r\n<\/pre><p>That was about all I understood at the time.  The roundoff error in the inversion process is magnified by the condition number of the matrix.  And, the error I observe is less than the estimate that this simple analysis provides.  So my inversion code passes this test.<\/p><h4>Deeper explanation<a name=\"5f919fa3-707a-4589-ae6e-eea641220b58\"><\/a><\/h4><p>I met <a href=\"https:\/\/blogs.mathworks.com\/cleve\/2013\/02\/18\">Jim Wilkinson<\/a> a few years later and came to realize that there is more to the story. I'm not actually inverting the Hilbert matrix.  There are roundoff errors involved in computing <tt>H<\/tt> even before it is passed to the inversion routine.<\/p><p>Today, the Symbolic Math Toolbox helps provide a deeper explanation. The <tt>'f'<\/tt> flag on the <tt>sym<\/tt> constructor says to convert double precision floating point arguments exactly to their rational representation. Here's how it works in this situation.  To save space, I'll print just the first column.<\/p><pre class=\"codeinput\">    H = hilb(n);\r\n    F = sym(H,<span class=\"string\">'f'<\/span>);\r\n    F(:,1)\r\n<\/pre><pre class=\"codeoutput\">ans =\r\n                                  1\r\n                                1\/2\r\n 6004799503160661\/18014398509481984\r\n                                1\/4\r\n 3602879701896397\/18014398509481984\r\n 6004799503160661\/36028797018963968\r\n 2573485501354569\/18014398509481984\r\n                                1\/8\r\n 2001599834386887\/18014398509481984\r\n 3602879701896397\/36028797018963968\r\n<\/pre><p>The elements of <tt>H<\/tt> that are not exact binary fractions become ratios of large integers.  The denominators are powers of two; in this case $2^{54}$ and $2^{55}$.  The numerators are these denominators divided by $3$, $5$, etc. and then rounded to the nearest integer. The elements of <tt>F<\/tt> are as close to the exact elements of a Hilbert matrix as we can get in binary floating point.<\/p><p>Let's invert <tt>F<\/tt>, using the exact rational arithmetic provided by the Symbolic Toolbox.  (I couldn't do this in 1960.)<\/p><pre class=\"codeinput\">    S = inv(F);\r\n<\/pre><p>We now have three inverse Hilbert matrices, <tt>X<\/tt>, <tt>S<\/tt>, and <tt>T<\/tt>.<\/p><div><ul><li><tt>X<\/tt> is the approximate inverse computed with floating point arithmetic by the routine I was testing years ago, or by MATLAB <tt>inv<\/tt> function today.<\/li><li><tt>S<\/tt> is the exact inverse of the floating point matrix that was actually passed to the inversion routine.<\/li><\/ul><\/div><div><ul><li><tt>T<\/tt> is the exact Hilbert inverse, obtained from the binomial coefficient formula.<\/li><\/ul><\/div><p>Let's print the first columns alongside each other.<\/p><pre class=\"codeinput\">   fprintf(<span class=\"string\">'%12s %16s %15s\\n'<\/span>,<span class=\"string\">'X'<\/span>,<span class=\"string\">'S'<\/span>,<span class=\"string\">'T'<\/span>)\r\n   fprintf(<span class=\"string\">'%16.4f %16.4f %12.0f\\n'<\/span>,[X(:,1) S(:,1) T(:,1)]')\r\n<\/pre><pre class=\"codeoutput\">           X                S               T\r\n         99.9961          99.9976          100\r\n      -4949.6667       -4949.7926        -4950\r\n      79192.8929       79195.5727        79200\r\n    -600535.3362     -600559.6914      -600600\r\n    2522211.3665     2522327.5182      2522520\r\n   -6305451.1288    -6305770.4041     -6306300\r\n    9608206.6797     9608730.4926      9609600\r\n   -8750253.0592    -8750759.2546     -8751600\r\n    4375092.6697     4375358.4162      4375800\r\n    -923624.4113     -923682.8529      -923780\r\n<\/pre><p>It looks like <tt>X<\/tt> is closer to <tt>S<\/tt> than <tt>S<\/tt> is to <tt>T<\/tt>.  Let's confirm by computing two relative errors, the difference between <tt>X<\/tt> and <tt>S<\/tt>, and the difference between <tt>S<\/tt> and <tt>T<\/tt>.<\/p><pre class=\"codeinput\">    format <span class=\"string\">shorte<\/span>\r\n    relerr(1) = maxabs(X - S)\/maxabs(T);\r\n    relerr(2) = maxabs(S - T)\/maxabs(T)\r\n    relerrtotal = sum(relerr)\r\n<\/pre><pre class=\"codeoutput\">relerr =\r\n   5.4143e-05   9.0252e-05\r\nrelerrtotal =\r\n   1.4439e-04\r\n<\/pre><p>The error in the computed inverse comes from two sources -- generating the matrix in the first place and then computing the inverse. The first of these is actually larger than the second, although the two are comparable.<\/p><script language=\"JavaScript\"> <!-- \r\n    function grabCode_b5980438b46c4cb5acd20123a1a034fd() {\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='b5980438b46c4cb5acd20123a1a034fd ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' b5980438b46c4cb5acd20123a1a034fd';\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 2017 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_b5980438b46c4cb5acd20123a1a034fd()\"><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; R2017a<br><\/p><\/div><!--\r\nb5980438b46c4cb5acd20123a1a034fd ##### SOURCE BEGIN #####\r\n%% Hilbert Matrices\r\n% I first encountered the Hilbert matrix when I was doing\r\n% individual studies under\r\n% <https:\/\/blogs.mathworks.com\/cleve\/2015\/12\/07\r\n% Professor John Todd> at Caltech in 1960.\r\n% It has been part of my professional life ever since.\r\n\r\n%% David Hilbert\r\n% Around the turn of the 20th century, David Hilbert was the world's\r\n% most famous mathematician.  He introduced the matrix that now bears\r\n% his name in a paper in 1895.  The elements of the matrix, which are\r\n% reciprocals of consecutive positive integers, are constant along the\r\n% antidiagonals.\r\n%\r\n% $$ h_{i,j} = \\frac{1}{i+j-1}, \\ \\ i,j = 1:n $$\r\n\r\n%%\r\n    format rat\r\n    H5 = hilb(5)\r\n    \r\n%%\r\n    latex(sym(H5));\r\n \r\n%%\r\n% $$ H_5 = \\left(\\begin{array}{ccccc} \r\n% 1 & \\frac{1}{2} & \\frac{1}{3} & \\frac{1}{4} & \\frac{1}{5}\\\\ \r\n% \\frac{1}{2} & \\frac{1}{3} & \\frac{1}{4} & \\frac{1}{5} & \\frac{1}{6}\\\\ \r\n% \\frac{1}{3} & \\frac{1}{4} & \\frac{1}{5} & \\frac{1}{6} & \\frac{1}{7}\\\\\r\n% \\frac{1}{4} & \\frac{1}{5} & \\frac{1}{6} & \\frac{1}{7} & \\frac{1}{8}\\\\\r\n% \\frac{1}{5} & \\frac{1}{6} & \\frac{1}{7} & \\frac{1}{8} & \\frac{1}{9}\r\n% \\end{array}\\right) $$\r\n\r\n\r\n%%\r\n% Here's a picture.  The continuous surface generated is very smooth.\r\n\r\n    H12 = hilb(12);\r\n    surf(log(H12))\r\n    view(60,60)\r\n\r\n%% Least squares\r\n% Hilbert was interested in this matrix because it comes up in the\r\n% least squares approximation of a continuous function on the unit\r\n% interval by polynomials, expressed in the conventional basis\r\n% as linear combinations of monomials.\r\n%\r\n% $$ p(x) = \\sum_{j=1}^n c_j x^{j-1} $$\r\n% \r\n% The coefficient matrix for the normal equations has elements\r\n%\r\n% $$ \\int_0^1 x^{i+j-2} dx \\ = \\ \\frac{1}{i+j-1} $$\r\n\r\n%% Properties\r\n% A Hilbert matrix has many useful properties.\r\n%\r\n% * Symmetric.\r\n% * Positive definite.\r\n% * Hankel, $a_{i,j}$ is a function of $i+j$.\r\n% * Cauchy, $a_{i,j} = 1\/(x_i - y_j)$.\r\n% * Nearly singular.\r\n\r\n%%\r\n% Each column of a Hilbert matrix is nearly a multiple of\r\n% the other columns.  So the columns are nearly linearly\r\n% dependent and the matrix is close to, but not exactly, singular.\r\n\r\n%% hilb\r\n% MATLAB has always had functions |hilb| and |invhilb| that compute\r\n% the Hilbert matrix and its inverse.  The body of |hilb| is now only\r\n% two lines.\r\n%\r\n%   J = 1:n;\r\n%   H = 1.\/(J'+J-1);\r\n%\r\n% We often cite this as a good example of _singleton expansion_.\r\n% A column vector is added to a row vector to produce a matrix,\r\n% then a scalar is subtracted from that matrix, and finally the\r\n% reciprocals of the elements produce the result.\r\n   \r\n%% Inverse\r\n% It is possible to express the elements of the inverse of the Hilbert\r\n% matrix in terms of binomial coefficients.  For reasons that I've now\r\n% forgotten, I always use $T$ for $H^{-1}$.\r\n%\r\n% $$ t_{i,j} = (-1)^{i+j} (i+j-1) {{n+i-1} \\choose {n-j}}\r\n%  {{n+j-1} \\choose {n-i}} {{i+j-2} \\choose {i-1}}^2 $$\r\n\r\n%%\r\n% The elements of the inverse Hilbert matrix are integers, but they\r\n% are _large_ integers.  The largest ones are on the diagonal.  For\r\n% order 13, the largest element is \r\n%\r\n% $$ (T_{13})_{9,9} \\ = \\ 100863567447142500 $$\r\n%\r\n% This is over $10^{17}$ and is larger than double precision |flintmax|.\r\n\r\n    format longe\r\n    flintmax = flintmax\r\n    \r\n%%\r\n% So, it is possible to represent the largest\r\n% elements exactly only if $n \\le 12$.\r\n\r\n%%\r\n% The HELP entry for |invhilb| includes a sentence inherited from\r\n% my original Fortran MATLAB.\r\n%\r\n%    The result is exact for  N  less than about 15.\r\n%\r\n% Now that's misleading.  It should say \r\n%\r\n%    The result is exact for  N  <= 12.\r\n% \r\n% (I'm filing a bug report.)\r\n\r\n%% invhilv\r\n% The body of |invhilb| begins by setting |p| to the order |n|.\r\n% The doubly nested |for| loops then evaluate the binomial coefficient\r\n% formula recursively, avoiding unnecessary integer overflow.\r\n\r\n    dbtype invhilb 18:28\r\n\r\n%%\r\n% I first programmed this algorithm in machine language for the Burroughs\r\n% 205 Datatron at Caltech almost 60 years ago.  I was barely out of my\r\n% teens.\r\n\r\n%%\r\n% Here's the result for |n = 6|.\r\n\r\n    format short\r\n    T6 = invhilb(6)\r\n    \r\n%%\r\n% A checkerboard sign pattern with large integers in the inverse cancels\r\n% the smooth surface of the Hilbert matrix itself.\r\n\r\n    T12 = invhilb(12);\r\n    spy(T12 > 0)\r\n\r\n%%\r\n% A log scale mitigates the jaggedness.\r\n\r\n    surf(sign(T12).*log(abs(T12)))\r\n    view(60,60)\r\n    \r\n%% Rookie experiment\r\n% Now using MATLAB, I am going to repeat the experiment that I did on\r\n% the Burroughs 205 when I was still a rookie.  I had just written my\r\n% first program that used Gaussian elimination to invert matrices.\r\n% I proceeded to test it by inverting Hilbert matrices and comparing\r\n% the results with the exact inverses.  (Today's code uses\r\n% this utility function that picks out the largest element in\r\n% a matrix.)\r\n\r\n    maxabs = @(X) max(double(abs(X(:))));\r\n\r\n%%\r\n% Here is |n = 10|.\r\n\r\n    n = 10\r\n    H = hilb(n);\r\n    X = inv(H);        % Computed inverse\r\n    T = invhilb(n);    % Theoretical inverse\r\n    E = X - T;         % The error\r\n    err = maxabs(E)\r\n    \r\n%%\r\n% At first I might have said:\r\n%\r\n% _Wow!  The error is $10^8$.  That's a pretty big error.\r\n% Can I trust my matrix inversion code?   What went wrong?_\r\n\r\n%%\r\n% But I knew the elements of the inverse are huge.\r\n% We should be looking at _relative_ error.\r\n\r\n   relerr = maxabs(E)\/maxabs(T)\r\n   \r\n%%\r\n% _OK.  The relative error is $10^{-4}$.  That still seems like a lot._\r\n\r\n%%\r\n% I knew that the Hilbert matrix is nearly singular. That's why I\r\n% was using it.  John Todd was one of the first people to write about\r\n% condition numbers.  An error estimate that involves nearness to\r\n% singularity and the floating point accuracy would be expressed today\r\n% by\r\n\r\n   esterr = cond(H)*eps\r\n   \r\n%%\r\n% That was about all I understood at the time.  The roundoff error\r\n% in the inversion process is magnified by the condition number of\r\n% the matrix.  And, the error I observe is less than the estimate that\r\n% this simple analysis provides.  So my inversion code passes this test.\r\n\r\n%% Deeper explanation\r\n% I met <https:\/\/blogs.mathworks.com\/cleve\/2013\/02\/18 Jim Wilkinson>\r\n% a few years later and came to realize that there is more to the story.\r\n% I'm not actually inverting the Hilbert matrix.  There are roundoff\r\n% errors involved in computing |H| even before it is passed to the\r\n% inversion routine.\r\n\r\n%%\r\n% Today, the Symbolic Math Toolbox helps provide a deeper explanation.\r\n% The |'f'| flag on the |sym| constructor says to convert double precision\r\n% floating point arguments exactly to their rational representation.\r\n% Here's how it works in this situation.  To save space, I'll print just\r\n% the first column.\r\n\r\n    H = hilb(n);\r\n    F = sym(H,'f');\r\n    F(:,1)\r\n    \r\n%%\r\n% The elements of |H| that are not exact binary fractions become ratios\r\n% of large integers.  The denominators are powers of two; in this case\r\n% $2^{54}$ and $2^{55}$.  The numerators are these denominators divided by\r\n% $3$, $5$, etc. and then rounded to the nearest integer.\r\n% The elements of |F| are as close to the exact elements of a Hilbert\r\n% matrix as we can get in binary floating point.\r\n\r\n%%\r\n% Let's invert |F|, using the exact rational arithmetic provided\r\n% by the Symbolic Toolbox.  (I couldn't do this in 1960.)\r\n\r\n    S = inv(F);\r\n    \r\n%%\r\n% We now have three inverse Hilbert matrices, |X|, |S|, and |T|.\r\n%\r\n% * |X| is the approximate inverse computed with floating point arithmetic\r\n% by the routine I was testing years ago, or by MATLAB |inv| function\r\n% today.\r\n% * |S| is the exact inverse of the floating point matrix \r\n% that was actually passed to the inversion routine.  \r\n%\r\n% * |T| is the exact Hilbert inverse, obtained from the binomial \r\n% coefficient formula.\r\n\r\n%%\r\n% Let's print the first columns alongside each other.\r\n\r\n   fprintf('%12s %16s %15s\\n','X','S','T')\r\n   fprintf('%16.4f %16.4f %12.0f\\n',[X(:,1) S(:,1) T(:,1)]')\r\n\r\n%%\r\n% It looks like |X| is closer to |S| than |S| is to |T|.  Let's confirm\r\n% by computing two relative errors, the difference between\r\n% |X| and |S|, and the difference between |S| and |T|.\r\n\r\n    format shorte\r\n    relerr(1) = maxabs(X - S)\/maxabs(T);\r\n    relerr(2) = maxabs(S - T)\/maxabs(T)\r\n    relerrtotal = sum(relerr)\r\n    \r\n%%\r\n% The error in the computed inverse comes from two sources REPLACE_WITH_DASH_DASH generating\r\n% the matrix in the first place and then computing the inverse.\r\n% The first of these is actually larger than the second, although\r\n% the two are comparable.\r\n\r\n##### SOURCE END ##### b5980438b46c4cb5acd20123a1a034fd\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img decoding=\"async\"  class=\"img-responsive\" src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/hilbert_blog_01.png\" onError=\"this.style.display ='none';\" \/><\/div><!--introduction--><p>I first encountered the Hilbert matrix when I was doing individual studies under <a href=\"https:\/\/blogs.mathworks.com\/cleve\/2015\/12\/07\">Professor John Todd<\/a> at Caltech in 1960. It has been part of my professional life ever since.... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2017\/06\/07\/hilbert-matrices-2\/\">read more >><\/a><\/p>","protected":false},"author":78,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[4,6,7,20],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/2528"}],"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=2528"}],"version-history":[{"count":1,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/2528\/revisions"}],"predecessor-version":[{"id":2529,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/2528\/revisions\/2529"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=2528"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=2528"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=2528"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}