{"id":15,"date":"2012-06-03T04:01:06","date_gmt":"2012-06-03T09:01:06","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=15"},"modified":"2013-05-02T10:06:45","modified_gmt":"2013-05-02T15:06:45","slug":"fibonacci-matrices","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2012\/06\/03\/fibonacci-matrices\/","title":{"rendered":"Fibonacci Matrices"},"content":{"rendered":"<!DOCTYPE html\r\n  PUBLIC \"-\/\/W3C\/\/DTD HTML 4.01 Transitional\/\/EN\">\r\n<style type=\"text\/css\">\r\n\r\nh1 { font-size:18pt; }\r\nh2.titlebg { font-size:13pt; }\r\nh3 { color:#4A4F55; padding:0px; margin:5px 0px 5px; font-family:Arial, Helvetica, sans-serif; font-size:11pt; font-weight:bold; line-height:140%; border-bottom:1px solid #d6d4d4; display:block; }\r\nh4 { color:#4A4F55; padding:0px; margin:0px 0px 5px; font-family:Arial, Helvetica, sans-serif; font-size:10pt; font-weight:bold; line-height:140%; border-bottom:1px solid #d6d4d4; display:block; }\r\n   \r\np { padding:0px; margin:0px 0px 20px; }\r\nimg { padding:0px; margin:0px 0px 20px; border:none; }\r\np img, pre img, tt img, li img { margin-bottom:0px; } \r\n\r\nul { padding:0px; margin:0px 0px 20px 23px; list-style:square; }\r\nul li { padding:0px; margin:0px 0px 7px 0px; background:none; }\r\nul li ul { padding:5px 0px 0px; margin:0px 0px 7px 23px; }\r\nul li ol li { list-style:decimal; }\r\nol { padding:0px; margin:0px 0px 20px 0px; list-style:decimal; }\r\nol li { padding:0px; margin:0px 0px 7px 23px; list-style-type:decimal; }\r\nol li ol { padding:5px 0px 0px; margin:0px 0px 7px 0px; }\r\nol li ol li { list-style-type:lower-alpha; }\r\nol li ul { padding-top:7px; }\r\nol li ul li { list-style:square; }\r\n\r\npre, tt, code { font-size:12px; }\r\npre { margin:0px 0px 20px; }\r\npre.error { color:red; }\r\npre.codeinput { padding:10px; border:1px solid #d3d3d3; background:#f7f7f7; }\r\npre.codeoutput { padding:10px 11px; margin:0px 0px 20px; color:#4c4c4c; }\r\n\r\n@media print { pre.codeinput, pre.codeoutput { word-wrap:break-word; width:100%; } }\r\n\r\nspan.keyword { color:#0000FF }\r\nspan.comment { color:#228B22 }\r\nspan.string { color:#A020F0 }\r\nspan.untermstring { color:#B20000 }\r\nspan.syscmd { color:#B28C00 }\r\n\r\n.footer { width:auto; padding:10px 0px; margin:25px 0px 0px; border-top:1px dotted #878787; font-size:0.8em; line-height:140%; font-style:italic; color:#878787; text-align:left; float:none; }\r\n.footer p { margin:0px; }\r\n\r\n  <\/style><div class=\"content\"><h3>Contents<\/h3><div><ul><li><a href=\"#3090c586-3dc6-48b5-a37c-3fb8ab591d01\">A bug report<\/a><\/li><li><a href=\"#c88c211f-704a-441b-b352-91d13f4037b3\">Powers of a Fibonacci matrix<\/a><\/li><li><a href=\"#f3ad9510-d468-4b52-8146-ba05f2ddb23a\">Determinants<\/a><\/li><li><a href=\"#c903db35-0b1f-4282-850f-a01391424783\">What happened?<\/a><\/li><li><a href=\"#1a5747bf-9ec1-4989-90b0-6f5fe655c6f3\">A log plot of the error<\/a><\/li><li><a href=\"#b5aef0f6-ee24-4c98-b4df-cfa449f45f4c\">Computing determinants<\/a><\/li><li><a href=\"#10d2ca58-3e6c-4461-ac98-afac011eef65\">Golden Ratio<\/a><\/li><li><a href=\"#b2248602-2c3b-4042-a4c0-72a116c325c7\">Condition<\/a><\/li><li><a href=\"#1dfc9840-8e09-4c5f-8c4c-36f39ec9e494\">Backward error<\/a><\/li><li><a href=\"#0dcf96a9-d2d5-45d4-a590-f9b0700aa81d\">The high school formula<\/a><\/li><li><a href=\"#2ca1eece-064a-4c83-a6ad-51cf9b0e88df\">Historical irony<\/a><\/li><li><a href=\"#a3a47226-7b2b-4a90-b918-2c0aee6ddf10\"><tt>eig<\/tt> and <tt>svd<\/tt><\/a><\/li><\/ul><\/div><h4>A bug report<a name=\"3090c586-3dc6-48b5-a37c-3fb8ab591d01\"><\/a><\/h4><p>MathWorks recently received a bug report involving the matrix<\/p><pre class=\"codeinput\">X = [  63245986, 102334155\r\n      102334155, 165580141]\r\n<\/pre><pre class=\"codeoutput\">\r\nX =\r\n\r\n    63245986   102334155\r\n   102334155   165580141\r\n\r\n<\/pre><p>The report claimed that MATLAB computed an inaccurate result for the determinant of <tt>X<\/tt>.<\/p><pre class=\"codeinput\">format <span class=\"string\">long<\/span>\r\ndetX = det(X)\r\n<\/pre><pre class=\"codeoutput\">\r\ndetX =\r\n\r\n   1.524897739291191\r\n\r\n<\/pre><p>The familiar high school formula gives a significantly different result.<\/p><pre class=\"codeinput\">detX = X(1,1)*X(2,2) - X(1,2)*X(2,1)\r\n<\/pre><pre class=\"codeoutput\">\r\ndetX =\r\n\r\n     2\r\n\r\n<\/pre><p>While checking out the report, I computed the elementwise ratio of the rows of <tt>X<\/tt> and discovered an old friend.<\/p><pre class=\"codeinput\">ratio = X(2,:).\/X(1,:)\r\n<\/pre><pre class=\"codeoutput\">\r\nratio =\r\n\r\n   1.618033988749895   1.618033988749895\r\n\r\n<\/pre><p>This is $\\phi$, the Golden ratio.<\/p><pre class=\"codeinput\">phi = (1 + sqrt(5))\/2\r\n<\/pre><pre class=\"codeoutput\">\r\nphi =\r\n\r\n   1.618033988749895\r\n\r\n<\/pre><p>So, I decided to investigate further.<\/p><h4>Powers of a Fibonacci matrix<a name=\"c88c211f-704a-441b-b352-91d13f4037b3\"><\/a><\/h4><p>Let's call the following 2-by-2 matrix the Fibonacci matrix.<\/p><p>$$F = \\pmatrix{0 &amp; 1 \\cr 1 &amp; 1}$$<\/p><p>Generate the matrix in MATLAB<\/p><pre class=\"codeinput\">F = [0 1; 1 1]\r\n<\/pre><pre class=\"codeoutput\">\r\nF =\r\n\r\n     0     1\r\n     1     1\r\n\r\n<\/pre><p>The elements of powers of $F$ are Fibonacci numbers. Let $f_n$ be the $n$-th Fibanacci number.  Then the $n$-th power of $F$ contains three successive $f_n$.<\/p><p>$$F^n = \\pmatrix{f_{n-1} &amp; f_n \\cr f_n &amp; f_{n+1}}$$<\/p><p>For example<\/p><pre class=\"codeinput\">F2 = F*F\r\nF3 = F*F2\r\nF4 = F*F3\r\nF5 = F*F4\r\n<\/pre><pre class=\"codeoutput\">\r\nF2 =\r\n\r\n     1     1\r\n     1     2\r\n\r\n\r\nF3 =\r\n\r\n     1     2\r\n     2     3\r\n\r\n\r\nF4 =\r\n\r\n     2     3\r\n     3     5\r\n\r\n\r\nF5 =\r\n\r\n     3     5\r\n     5     8\r\n\r\n<\/pre><p>The matrix in the bug report is $F^{40}$.<\/p><pre class=\"codeinput\">X = F^40\r\n<\/pre><pre class=\"codeoutput\">\r\nX =\r\n\r\n    63245986   102334155\r\n   102334155   165580141\r\n\r\n<\/pre><p>These matrix powers are computed without any roundoff error. Their elements are \"flints\", floating point numbers whose values are integers.<\/p><h4>Determinants<a name=\"f3ad9510-d468-4b52-8146-ba05f2ddb23a\"><\/a><\/h4><p>The determinant of $F$ is clearly $-1$.  Hence<\/p><p>$$\\mbox{det}(F^n) = (-1)^n$$<\/p><p>So <tt>det(F^40)<\/tt> should be <tt>1<\/tt>, not <tt>1.5249<\/tt>, and not <tt>2<\/tt>.<\/p><p>Let's examine the computation of determinants using floating point arithmetic. We expect these to be +1 or -1.<\/p><pre class=\"codeinput\">det1 = det(F)\r\ndet2 = det(F^2)\r\ndet3 = det(F^3)\r\n<\/pre><pre class=\"codeoutput\">\r\ndet1 =\r\n\r\n    -1\r\n\r\n\r\ndet2 =\r\n\r\n     1\r\n\r\n\r\ndet3 =\r\n\r\n    -1\r\n\r\n<\/pre><p>So far, so good.  However,<\/p><pre class=\"codeinput\">det40 = det(F^40)\r\n<\/pre><pre class=\"codeoutput\">\r\ndet40 =\r\n\r\n   1.524897739291191\r\n\r\n<\/pre><p>This has barely one digit of accuracy.<\/p><h4>What happened?<a name=\"c903db35-0b1f-4282-850f-a01391424783\"><\/a><\/h4><p>It is instructive to look at the accuracy of all the determinants.<\/p><pre class=\"codeinput\">d = zeros(40,1);\r\n<span class=\"keyword\">for<\/span> n = 1:40\r\n   d(n) = det(F^n);\r\n<span class=\"keyword\">end<\/span>\r\nformat <span class=\"string\">long<\/span>\r\nd\r\n<\/pre><pre class=\"codeoutput\">\r\nd =\r\n\r\n  -1.000000000000000\r\n   1.000000000000000\r\n  -1.000000000000000\r\n   0.999999999999999\r\n  -0.999999999999996\r\n   1.000000000000000\r\n  -0.999999999999996\r\n   0.999999999999996\r\n  -1.000000000000057\r\n   0.999999999999872\r\n  -0.999999999999943\r\n   0.999999999997726\r\n  -0.999999999998209\r\n   1.000000000000227\r\n  -1.000000000029786\r\n   1.000000000138243\r\n  -0.999999999990905\r\n   1.000000000261935\r\n  -1.000000000591172\r\n   0.999999999999091\r\n  -0.999999996060069\r\n   1.000000003768946\r\n  -0.999999926301825\r\n   0.999999793712050\r\n  -1.000000685962732\r\n   0.999998350307578\r\n  -0.999995890248101\r\n   0.999983831774443\r\n  -0.999970503151417\r\n   1.000005397945643\r\n  -0.999914042185992\r\n   0.999138860497624\r\n  -0.997885793447495\r\n   0.998510751873255\r\n  -0.996874589473009\r\n   0.973348170518875\r\n  -0.989943694323301\r\n   0.873688660562038\r\n  -0.471219316124916\r\n   1.524897739291191\r\n\r\n<\/pre><p>We see that the accuracy deteriorates as <tt>n<\/tt> increases. In fact, the number of correct digits is roughly a linear function of <tt>n<\/tt>, reaching zero around n = 40.<\/p><h4>A log plot of the error<a name=\"1a5747bf-9ec1-4989-90b0-6f5fe655c6f3\"><\/a><\/h4><pre class=\"codeinput\">semilogy(abs(abs(d)-1),<span class=\"string\">'.'<\/span>)\r\nset(gca,<span class=\"string\">'ydir'<\/span>,<span class=\"string\">'rev'<\/span>)\r\ntitle(<span class=\"string\">'error in det(F^n)'<\/span>)\r\nxlabel(<span class=\"string\">'n'<\/span>)\r\nylabel(<span class=\"string\">'error'<\/span>)\r\naxis([0 41 eps 1])\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/fibmat_01.png\" alt=\"\"> <h4>Computing determinants<a name=\"b5aef0f6-ee24-4c98-b4df-cfa449f45f4c\"><\/a><\/h4><p>MATLAB computes determinants as the product of the diagonal elements of the triangular factors resulting from Gaussian elimination.<\/p><p>Let's look at $F^{12}$.<\/p><pre class=\"codeinput\">format <span class=\"string\">long<\/span> <span class=\"string\">e<\/span>\r\nF12 = F^12\r\n[L,U] = lu(F12)\r\n<\/pre><pre class=\"codeoutput\">\r\nF12 =\r\n\r\n    89   144\r\n   144   233\r\n\r\n\r\nL =\r\n\r\n     6.180555555555555e-01     1.000000000000000e+00\r\n     1.000000000000000e+00                         0\r\n\r\n\r\nU =\r\n\r\n     1.440000000000000e+02     2.330000000000000e+02\r\n                         0    -6.944444444428655e-03\r\n\r\n<\/pre><p>Since $det(L) = -1$,<\/p><pre class=\"codeinput\">det(L)\r\n<\/pre><pre class=\"codeoutput\">\r\nans =\r\n\r\n    -1\r\n\r\n<\/pre><p>we have<\/p><pre class=\"codeinput\">det12 = -prod(diag(U))\r\n<\/pre><pre class=\"codeoutput\">\r\ndet12 =\r\n\r\n     9.999999999977263e-01\r\n\r\n<\/pre><p>We can see that, for $F^{12}$, the crucial element <tt>U(2,2)<\/tt> has lost about 5 significant digits in the elimination and that this is reflected in the accuracy of the determinant.<\/p><p>How about $F^{40}$?<\/p><pre class=\"codeinput\">[L,U] = lu(F^40)\r\ndet40 = -prod(diag(U))\r\n<\/pre><pre class=\"codeoutput\">\r\nL =\r\n\r\n     6.180339887498949e-01     1.000000000000000e+00\r\n     1.000000000000000e+00                         0\r\n\r\n\r\nU =\r\n\r\n     1.023341550000000e+08     1.655801410000000e+08\r\n                         0    -1.490116119384766e-08\r\n\r\n\r\ndet40 =\r\n\r\n     1.524897739291191e+00\r\n\r\n<\/pre><p>The troubling value of <tt>1.5249<\/tt> from the bug report is a direct result of the subtractive cancellation involved in computing <tt>U(2,2)<\/tt>. In order for the computed determinant to be <tt>1<\/tt>, the value of <tt>U(2,2)<\/tt> should have been<\/p><pre class=\"codeinput\">U22 = -1\/U(1,1)\r\n<\/pre><pre class=\"codeoutput\">\r\nU22 =\r\n\r\n    -9.771908508943080e-09\r\n\r\n<\/pre><p>Compare <tt>U22<\/tt> with <tt>U(2,2)<\/tt>.  The value of <tt>U(2,2)<\/tt> resulting from elimination has lost almost all its accuracy.<\/p><p>This prompts us to check the condition of $F^{40}$.<\/p><pre class=\"codeinput\">cond40 = cond(F^40)\r\n<\/pre><pre class=\"codeoutput\">\r\ncond40 =\r\n\r\n   Inf\r\n\r\n<\/pre><p>Well, that's not much help. Something is going wrong in the computation of <tt>cond(F^40)<\/tt>.<\/p><h4>Golden Ratio<a name=\"10d2ca58-3e6c-4461-ac98-afac011eef65\"><\/a><\/h4><p>There is an intimate relationship between Fibonacci numbers and the Golden Ratio. The eigenvalues of $F$ are $\\phi$ and its negative reciprocal.<\/p><p>$$\\phi = (1+\\sqrt{5})\/2$$<\/p><p>$$\\bar{\\phi} = (1-\\sqrt{5})\/2$$<\/p><pre class=\"codeinput\">format <span class=\"string\">long<\/span>\r\nphibar = (1-sqrt(5))\/2\r\nphi = (1+sqrt(5))\/2\r\neigF = eig(F)\r\n<\/pre><pre class=\"codeoutput\">\r\nphibar =\r\n\r\n  -0.618033988749895\r\n\r\n\r\nphi =\r\n\r\n   1.618033988749895\r\n\r\n\r\neigF =\r\n\r\n  -0.618033988749895\r\n   1.618033988749895\r\n\r\n<\/pre><p>Because powers of $\\phi$ dominate powers of $\\bar{\\phi}$, it is possible to generate $F^n$ by rounding a scaled matrix of powers of $\\phi$ to the nearest integers.<\/p><p>$$F^n = \\mbox{round}\\left(\\pmatrix{\\phi^{n-1} &amp; \\phi^n \\cr \\phi^n &amp; \\phi^{n+1}}\/\\sqrt{5}\\right)$$<\/p><pre class=\"codeinput\">n = 40\r\nF40 = round( [phi^(n-1)  phi^n;  phi^n  phi^(n+1)]\/sqrt(5) )\r\n<\/pre><pre class=\"codeoutput\">\r\nn =\r\n\r\n    40\r\n\r\n\r\nF40 =\r\n\r\n    63245986   102334155\r\n   102334155   165580141\r\n\r\n<\/pre><p>Before rounding, the matrix of powers of $\\phi$ is clearly singular, so the rounded matrix $F^n$ must be regarded as \"close to singular\", even though its determinant is $+1$.  This is yet another example of the fact that the size of the determinant cannot be a reliable indication of nearness to singularity.<\/p><h4>Condition<a name=\"b2248602-2c3b-4042-a4c0-72a116c325c7\"><\/a><\/h4><p>The condition number of $F$ is the ratio of its singular values, and because $F$ is symmetric, its singular values are the absolute value of its eigenvalues.  So<\/p><p>$$\\mbox{cond}(F) = \\phi\/|\\bar{\\phi}| = \\phi^2 = 1+\\phi$$<\/p><pre class=\"codeinput\">condF = 1+phi\r\n<\/pre><pre class=\"codeoutput\">\r\ncondF =\r\n\r\n   2.618033988749895\r\n\r\n<\/pre><p>This agrees with<\/p><pre class=\"codeinput\">condF = cond(F)\r\n<\/pre><pre class=\"codeoutput\">\r\ncondF =\r\n\r\n   2.618033988749895\r\n\r\n<\/pre><p>With the help of the Symbolic Toolbox, we can get an exact expression for $\\mbox{cond}(F^{40}) = \\mbox{cond}(F)^{40}$<\/p><pre class=\"codeinput\">Phi = sym(<span class=\"string\">'(1+sqrt(5))\/2'<\/span>);\r\ncond40 = expand((1+Phi)^40)\r\n<\/pre><pre class=\"codeoutput\"> \r\ncond40 =\r\n \r\n(23416728348467685*5^(1\/2))\/2 + 52361396397820127\/2\r\n \r\n<\/pre><p>The numerical value is<\/p><pre class=\"codeinput\">cond40 = double(cond40)\r\n<\/pre><pre class=\"codeoutput\">\r\ncond40 =\r\n\r\n     5.236139639782013e+16\r\n\r\n<\/pre><p>This is better than the <tt>Inf<\/tt> we obtain from <tt>cond(F^40)<\/tt>. Note that it is more than an order of magnitude larger than <tt>1\/eps<\/tt>.<\/p><h4>Backward error<a name=\"1dfc9840-8e09-4c5f-8c4c-36f39ec9e494\"><\/a><\/h4><p>J. H. Wilkinson showed us that the <tt>L<\/tt> and <tt>U<\/tt> factors computing by Gaussian elimination are the exact factors of some matrix within roundoff error of the given matrix.  With hindsight and the Symbolic Toolbox, we can find that matrix.<\/p><pre class=\"codeinput\">digits(25)\r\nX = vpa(L)*vpa(U)\r\n<\/pre><pre class=\"codeoutput\"> \r\nX =\r\n \r\n[ 63245986.00000000118877885, 102334154.9999999967942319]\r\n[                102334155.0,                165580141.0]\r\n \r\n<\/pre><p>Our computed <tt>L<\/tt> and <tt>U<\/tt> are the actual triangular decomposition of this extended precision matrix.  And the determinant of this <tt>X<\/tt> is the shaky value that prompted this investigation.<\/p><pre class=\"codeinput\">det(X)\r\n<\/pre><pre class=\"codeoutput\"> \r\nans =\r\n \r\n1.524897739291191101074219\r\n \r\n<\/pre><p>The double precision matrix closest to <tt>X<\/tt> is precisely $F^{40}$.<\/p><pre class=\"codeinput\">double(X)\r\n<\/pre><pre class=\"codeoutput\">\r\nans =\r\n\r\n    63245986   102334155\r\n   102334155   165580141\r\n\r\n<\/pre><p>Retrospective backward error analysis confirms Wilkinson's theory in this example.<\/p><h4>The high school formula<a name=\"0dcf96a9-d2d5-45d4-a590-f9b0700aa81d\"><\/a><\/h4><p>The familiar formula for 2-by-2 determinants<\/p><p>$$\\mbox{det}(X) = x_{1,1} x_{2,2} - x_{1,2} x_{2,1}$$<\/p><p>gives +1 or -1 with no roundoff error for all $X = F^n$ with $n &lt; 40$. It fails completely for $n &gt; 40$. The behavior for exactly $n = 40$ is interesting. Set the output format to<\/p><pre class=\"codeinput\">format <span class=\"string\">bank<\/span>\r\n<\/pre><p>This format allows us to see all the digits generated when converting binary floating point numbers to decimal for printing. Again, let<\/p><pre class=\"codeinput\">X = F^40\r\n<\/pre><pre class=\"codeoutput\">\r\nX =\r\n\r\n   63245986.00  102334155.00\r\n  102334155.00  165580141.00\r\n\r\n<\/pre><p>Look at the last digits of these two products.<\/p><pre class=\"codeinput\">p1 = X(1,1)*X(2,2)\r\np2 = X(2,1)*X(1,2)\r\n<\/pre><pre class=\"codeoutput\">\r\np1 =\r\n\r\n 10472279279564026.00\r\n\r\n\r\np2 =\r\n\r\n 10472279279564024.00\r\n\r\n<\/pre><p>The <tt>6<\/tt> at the end of <tt>p1<\/tt> is correct because <tt>X(1,1)<\/tt> ends in a <tt>6<\/tt> and <tt>X(2,2)<\/tt> ends in a <tt>1<\/tt>.  But the <tt>4<\/tt> at the end of <tt>p2<\/tt> is incorrect. It should be a <tt>5<\/tt> because both <tt>X(2,1)<\/tt> and <tt>X(1,2)<\/tt> end in <tt>5<\/tt>. However, the spacing between floating point numbers of this magnitude is<\/p><pre class=\"codeinput\">format <span class=\"string\">short<\/span>\r\ndelta = eps(p2)\r\n<\/pre><pre class=\"codeoutput\">\r\ndelta =\r\n\r\n     2\r\n\r\n<\/pre><p>So near <tt>p2<\/tt> only even flints can be represented. In fact, <tt>p1<\/tt> is the next floating point number after <tt>p2<\/tt>. The true product of <tt>X(2,1)<\/tt> and <tt>X(1,2)<\/tt> falls halfway between <tt>p1<\/tt> and <tt>p2<\/tt> and must be rounded to one or the other. The familiar formula cannot possibly produce the correct result.<\/p><h4>Historical irony<a name=\"2ca1eece-064a-4c83-a6ad-51cf9b0e88df\"><\/a><\/h4><p>For many years, the <tt>det<\/tt> function in MATLAB would apply the <tt>round<\/tt> function to the computed value if all the matrix entries were integers. So, these old versions would have returned exactly <tt>+1<\/tt> or <tt>-1<\/tt> for <tt>det(F^n)<\/tt> with <tt>n &lt; 40<\/tt>.  But they would round <tt>det(F^40)<\/tt> to <tt>2<\/tt>. This kind of behavior was the reason we got rid of the rounding.<\/p><h4><tt>eig<\/tt> and <tt>svd<\/tt><a name=\"a3a47226-7b2b-4a90-b918-2c0aee6ddf10\"><\/a><\/h4><p>Let's try to compute eigenvalues and singular values of $F^{40}$.<\/p><pre class=\"codeinput\">format <span class=\"string\">long<\/span> <span class=\"string\">e<\/span>\r\neig(F^40)\r\nsvd(F^40)\r\n<\/pre><pre class=\"codeoutput\">\r\nans =\r\n\r\n           0\r\n   228826127\r\n\r\n\r\nans =\r\n\r\n     2.288261270000000e+08\r\n                         0\r\n\r\n<\/pre><p>Unfortunately, the small eigenvalue and small singular value are completely lost.  We know that the small value should be<\/p><pre class=\"codeinput\">phibar^40\r\n<\/pre><pre class=\"codeoutput\">\r\nans =\r\n\r\n     4.370130339181083e-09\r\n\r\n<\/pre><p>But this is smaller than roundoff error in the large value<\/p><pre class=\"codeinput\">eps(phi^40)\r\n<\/pre><pre class=\"codeoutput\">\r\nans =\r\n\r\n     2.980232238769531e-08\r\n\r\n<\/pre><p>Log plots of the accuracy loss in the computed small eigenvalue and the small singular value are similar to our plot for the determinant. I'm not sure how LAPACK computes eigenvalues and singular values of 2-by-2 matrices.  Perhaps that can be the subject of a future blog.<\/p><script language=\"JavaScript\"> <!-- \r\n    function grabCode_a3a22f3c6e84406ca3c3172c788adfcd() {\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='a3a22f3c6e84406ca3c3172c788adfcd ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' a3a22f3c6e84406ca3c3172c788adfcd';\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 2012 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_a3a22f3c6e84406ca3c3172c788adfcd()\"><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; 7.14<br><\/p><p class=\"footer\"><br>\r\n      Published with MATLAB&reg; 7.14<br><\/p><\/div><!--\r\na3a22f3c6e84406ca3c3172c788adfcd ##### SOURCE BEGIN #####\r\n%% A bug report\r\n% MathWorks recently received a bug report involving the matrix\r\nX = [  63245986, 102334155\r\n      102334155, 165580141]\r\n\r\n%%\r\n% The report claimed that MATLAB computed an inaccurate result for the\r\n% determinant of |X|.\r\nformat long\r\ndetX = det(X)\r\n\r\n%%\r\n% The familiar high school formula gives a significantly different result.\r\ndetX = X(1,1)*X(2,2) - X(1,2)*X(2,1)\r\n\r\n%%\r\n% While checking out the report, I computed the elementwise ratio of the\r\n% rows of |X| and discovered an old friend.\r\nratio = X(2,:).\/X(1,:)\r\n\r\n%%\r\n% This is $\\phi$, the Golden ratio.\r\nphi = (1 + sqrt(5))\/2\r\n\r\n%%\r\n% So, I decided to investigate further.\r\n\r\n%% Powers of a Fibonacci matrix\r\n% Let's call the following 2-by-2 matrix the Fibonacci matrix.\r\n%\r\n% $$F = \\pmatrix{0 & 1 \\cr 1 & 1}$$\r\n%\r\n% Generate the matrix in MATLAB\r\nF = [0 1; 1 1]\r\n\r\n%%\r\n% The elements of powers of $F$ are Fibonacci numbers.\r\n% Let $f_n$ be the $n$-th Fibanacci number.  Then the\r\n% $n$-th power of $F$ contains three successive $f_n$.\r\n%\r\n% $$F^n = \\pmatrix{f_{n-1} & f_n \\cr f_n & f_{n+1}}$$\r\n%\r\n% For example\r\nF2 = F*F\r\nF3 = F*F2\r\nF4 = F*F3\r\nF5 = F*F4\r\n\r\n%%\r\n% The matrix in the bug report is $F^{40}$.\r\nX = F^40\r\n\r\n%%\r\n% These matrix powers are computed without any roundoff error.\r\n% Their elements are \"flints\", floating point numbers whose values are integers.\r\n\r\n%% Determinants\r\n% The determinant of $F$ is clearly $-1$.  Hence\r\n%\r\n% $$\\mbox{det}(F^n) = (-1)^n$$\r\n%\r\n% So |det(F^40)| should be |1|, not |1.5249|, and not |2|.\r\n\r\n%%\r\n% Let's examine the computation of determinants using floating point arithmetic.\r\n% We expect these to be +1 or -1.\r\n\r\ndet1 = det(F)\r\ndet2 = det(F^2)\r\ndet3 = det(F^3)\r\n\r\n%%\r\n% So far, so good.  However,\r\ndet40 = det(F^40)\r\n\r\n%%\r\n% This has barely one digit of accuracy.\r\n\r\n%% What happened?\r\n% It is instructive to look at the accuracy of all the determinants.\r\n\r\nd = zeros(40,1);\r\nfor n = 1:40\r\n   d(n) = det(F^n);\r\nend\r\nformat long\r\nd\r\n\r\n%%\r\n% We see that the accuracy deteriorates as |n| increases.\r\n% In fact, the number of correct digits is roughly a linear function of |n|,\r\n% reaching zero around n = 40.\r\n\r\n%% A log plot of the error\r\nsemilogy(abs(abs(d)-1),'.')\r\nset(gca,'ydir','rev')\r\ntitle('error in det(F^n)')\r\nxlabel('n')\r\nylabel('error')\r\naxis([0 41 eps 1])\r\n\r\n%% Computing determinants\r\n% MATLAB computes determinants as the product of the diagonal elements\r\n% of the triangular factors resulting from Gaussian elimination.\r\n\r\n%%\r\n% Let's look at $F^{12}$.\r\n\r\nformat long e\r\nF12 = F^12\r\n[L,U] = lu(F12)\r\n\r\n%%\r\n% Since $det(L) = -1$,\r\n\r\ndet(L)\r\n\r\n%%\r\n% we have\r\n\r\ndet12 = -prod(diag(U))\r\n\r\n%%\r\n% We can see that, for $F^{12}$, the crucial element |U(2,2)| has lost about 5\r\n% significant digits in the elimination and that this is reflected in the accuracy\r\n% of the determinant.\r\n\r\n%%\r\n% How about $F^{40}$?\r\n\r\n[L,U] = lu(F^40)\r\ndet40 = -prod(diag(U))\r\n\r\n%%\r\n% The troubling value of |1.5249| from the bug report is a direct result of the\r\n% subtractive cancellation involved in computing |U(2,2)|.\r\n% In order for the computed determinant to be |1|, the value of |U(2,2)|\r\n% should have been\r\n\r\nU22 = -1\/U(1,1)\r\n\r\n%%\r\n% Compare |U22| with |U(2,2)|.  The value of |U(2,2)| resulting from elimination\r\n% has lost almost all its accuracy.\r\n\r\n%%\r\n% This prompts us to check the condition of $F^{40}$.\r\ncond40 = cond(F^40)\r\n\r\n%%\r\n% Well, that's not much help.\r\n% Something is going wrong in the computation of |cond(F^40)|.\r\n\r\n%% Golden Ratio\r\n% There is an intimate relationship between Fibonacci numbers and\r\n% the Golden Ratio.\r\n% The eigenvalues of $F$ are $\\phi$ and its negative reciprocal.\r\n%\r\n% $$\\phi = (1+\\sqrt{5})\/2$$\r\n%\r\n% $$\\bar{\\phi} = (1-\\sqrt{5})\/2$$\r\n\r\nformat long\r\nphibar = (1-sqrt(5))\/2\r\nphi = (1+sqrt(5))\/2\r\neigF = eig(F)\r\n\r\n%%\r\n% Because powers of $\\phi$ dominate powers of $\\bar{\\phi}$, it is possible\r\n% to generate $F^n$ by rounding a scaled matrix of powers of $\\phi$ to the\r\n% nearest integers.\r\n%\r\n% $$F^n = \\mbox{round}\\left(\\pmatrix{\\phi^{n-1} & \\phi^n \\cr \\phi^n & \\phi^{n+1}}\/\\sqrt{5}\\right)$$\r\n\r\nn = 40\r\nF40 = round( [phi^(n-1)  phi^n;  phi^n  phi^(n+1)]\/sqrt(5) )\r\n\r\n%%\r\n% Before rounding, the matrix of powers of $\\phi$ is clearly singular, so the\r\n% rounded matrix $F^n$ must be regarded as \"close to singular\", even though its\r\n% determinant is $+1$.  This is yet another example of the fact that the size\r\n% of the determinant cannot be a reliable indication of nearness to singularity.\r\n\r\n%% Condition\r\n% The condition number of $F$ is the ratio of its singular values, and\r\n% because $F$ is symmetric, its singular values are the absolute value of\r\n% its eigenvalues.  So\r\n%\r\n% $$\\mbox{cond}(F) = \\phi\/|\\bar{\\phi}| = \\phi^2 = 1+\\phi$$\r\n\r\ncondF = 1+phi\r\n\r\n%%\r\n% This agrees with\r\n\r\ncondF = cond(F)\r\n\r\n%%\r\n% With the help of the Symbolic Toolbox, we can get an exact expression\r\n% for $\\mbox{cond}(F^{40}) = \\mbox{cond}(F)^{40}$\r\n\r\nPhi = sym('(1+sqrt(5))\/2');\r\ncond40 = expand((1+Phi)^40)\r\n\r\n%%\r\n% The numerical value is\r\n\r\ncond40 = double(cond40)\r\n\r\n%%\r\n% This is better than the |Inf| we obtain from |cond(F^40)|.\r\n% Note that it is more than an order of magnitude larger than |1\/eps|.\r\n\r\n%% Backward error\r\n% J. H. Wilkinson showed us that the |L| and |U| factors computing by Gaussian\r\n% elimination are the exact factors of some matrix within roundoff error\r\n% of the given matrix.  With hindsight and the Symbolic Toolbox, we can\r\n% find that matrix.\r\n\r\ndigits(25)\r\nX = vpa(L)*vpa(U)\r\n\r\n%%\r\n% Our computed |L| and |U| are the actual triangular decomposition of this\r\n% extended precision matrix.  And the determinant of this |X| is the shaky value\r\n% that prompted this investigation.\r\n\r\ndet(X)\r\n\r\n%%\r\n% The double precision matrix closest to |X| is precisely $F^{40}$.\r\n\r\ndouble(X)\r\n\r\n%%\r\n% Retrospective backward error analysis confirms Wilkinson's theory\r\n% in this example.\r\n\r\n%% The high school formula\r\n% The familiar formula for 2-by-2 determinants\r\n%\r\n% $$\\mbox{det}(X) = x_{1,1} x_{2,2} - x_{1,2} x_{2,1}$$\r\n%\r\n% gives +1 or -1 with no roundoff error for all $X = F^n$ with $n < 40$.\r\n% It fails completely for $n > 40$.\r\n% The behavior for exactly $n = 40$ is interesting.\r\n% Set the output format to\r\n\r\nformat bank\r\n\r\n%%\r\n% This format allows us to see all the digits generated when converting\r\n% binary floating point numbers to decimal for printing.\r\n% Again, let\r\n\r\nX = F^40\r\n\r\n%%\r\n% Look at the last digits of these two products.\r\n\r\np1 = X(1,1)*X(2,2)\r\np2 = X(2,1)*X(1,2)\r\n\r\n%%\r\n% The |6| at the end of |p1| is correct because |X(1,1)| ends in a |6|\r\n% and |X(2,2)| ends in a |1|.  But the |4| at the end of |p2| is incorrect.\r\n% It should be a |5| because both |X(2,1)| and |X(1,2)| end in |5|.\r\n% However, the spacing between floating point numbers of this magnitude is\r\n\r\nformat short\r\ndelta = eps(p2)\r\n\r\n%%\r\n% So near |p2| only even flints can be represented.\r\n% In fact, |p1| is the next floating point number after |p2|.\r\n% The true product of |X(2,1)| and |X(1,2)| falls halfway between |p1|\r\n% and |p2| and must be rounded to one or the other.\r\n% The familiar formula cannot possibly produce the correct result.\r\n\r\n%% Historical irony\r\n% For many years, the |det| function in MATLAB would apply the |round|\r\n% function to the computed value if all the matrix entries were integers.\r\n% So, these old versions would have returned exactly |+1| or |-1|\r\n% for |det(F^n)| with |n < 40|.  But they would round |det(F^40)| to |2|.\r\n% This kind of behavior was the reason we got rid of the rounding.\r\n\r\n%% |eig| and |svd|\r\n% Let's try to compute eigenvalues and singular values of $F^{40}$.\r\n\r\nformat long e\r\neig(F^40)\r\nsvd(F^40)\r\n\r\n%%\r\n% Unfortunately, the small eigenvalue and small singular value are\r\n% completely lost.  We know that the small value should be\r\n\r\nphibar^40 \r\n  \r\n%%\r\n% But this is smaller than roundoff error in the large value\r\n\r\neps(phi^40)\r\n\r\n%%\r\n% Log plots of the accuracy loss in the computed small eigenvalue and\r\n% the small singular value are similar to our plot for the determinant.\r\n% I'm not sure how LAPACK computes eigenvalues and singular values\r\n% of 2-by-2 matrices.  Perhaps that can be the subject of a future blog.\r\n\r\n##### SOURCE END ##### a3a22f3c6e84406ca3c3172c788adfcd\r\n-->","protected":false},"excerpt":{"rendered":"<p>\r\n\r\n\r\nh1 { font-size:18pt; }\r\nh2.titlebg { font-size:13pt; }\r\nh3 { color:#4A4F55; padding:0px; margin:5px 0px 5px; font-family:Arial, Helvetica, sans-serif; font-size:11pt; font-weight:bold;... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2012\/06\/03\/fibonacci-matrices\/\">read more >><\/a><\/p>","protected":false},"author":78,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[6],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/15"}],"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=15"}],"version-history":[{"count":55,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/15\/revisions"}],"predecessor-version":[{"id":45,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/15\/revisions\/45"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=15"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=15"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=15"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}