{"id":1246,"date":"2015-08-24T12:00:05","date_gmt":"2015-08-24T17:00:05","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=1246"},"modified":"2016-12-05T14:06:20","modified_gmt":"2016-12-05T19:06:20","slug":"golub-matrices-deceptively-ill-conditioned","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2015\/08\/24\/golub-matrices-deceptively-ill-conditioned\/","title":{"rendered":"Golub Matrices, Deceptively Ill Conditioned"},"content":{"rendered":"<div class=\"content\">\r\n\r\n<!--introduction-->I got the idea for these test matrices from Gene Golub years ago, but the mathematical foundation comes from a paper by Divakar Viswanath and Nick Trefethen.\r\n\r\n<!--\/introduction-->\r\n<h3>Contents<\/h3>\r\n<div>\r\n<ul>\r\n \t<li><a href=\"#fc5195e1-df62-42a9-aa8a-6f9634724cf9\">Two programs from NCM<\/a><\/li>\r\n \t<li><a href=\"#e367873d-22d1-4ff9-81c1-4bdb8cdafb96\">A 6-by-6 example<\/a><\/li>\r\n \t<li><a href=\"#509e139e-0e70-4b6f-b86b-b5b4d6036233\">Product of triangular factors<\/a><\/li>\r\n \t<li><a href=\"#e64f34b8-5c00-489b-bce8-8cb8195b7f5a\">Experiment with pivoting<\/a><\/li>\r\n \t<li><a href=\"#e335aab0-f3e2-4b06-ac1d-5cf4cf788b43\">Resurrecting ancestors<\/a><\/li>\r\n \t<li><a href=\"#1efecaec-2cf6-46c5-a964-5194450d0480\">Diagonal of U<\/a><\/li>\r\n \t<li><a href=\"#63db1065-2ccb-49b3-8edb-9a0a9b7e7123\">Orthogonal factorization<\/a><\/li>\r\n \t<li><a href=\"#df687528-dcf1-430e-88f4-a23dc7092eb0\">References<\/a><\/li>\r\n<\/ul>\r\n<\/div>\r\n<h4>Two programs from NCM<a name=\"fc5195e1-df62-42a9-aa8a-6f9634724cf9\"><\/a><\/h4>\r\nThis post is about two MATLAB programs included in the collection from <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/37976-numerical-computing-with-matlab\"><i>Numerical Computing with MATLAB<\/i><\/a>. The program <a href=\"https:\/\/www.mathworks.com\/content\/dam\/mathworks\/mathworks-dot-com\/moler\/ncm\/golub.m\"><tt>golub<\/tt><\/a> generates badly conditioned integer test matrices. The program <a href=\"https:\/\/www.mathworks.com\/content\/dam\/mathworks\/mathworks-dot-com\/moler\/ncm\/lugui.m\"><tt>lugui<\/tt><\/a> is an interactive graphical interface that allows you to experiment with pivot strategies in Gaussian elimination.\r\n<h4>A 6-by-6 example<a name=\"e367873d-22d1-4ff9-81c1-4bdb8cdafb96\"><\/a><\/h4>\r\nSet the random number generator to the state it has when MATLAB starts and generate a 6-by-6 example.\r\n<pre class=\"codeinput\">   rng(<span class=\"string\">'default'<\/span>)\r\n   A = golub(6)\r\n<\/pre>\r\n<pre class=\"codeoutput\">A =\r\n     1     3    11     0   -11   -15\r\n    18    55   209    15  -198  -277\r\n   -23   -33   144   532   259    82\r\n     9    55   405   437  -100  -285\r\n     3    -4  -111  -180    39   219\r\n   -13    -9   202   346   401   253\r\n<\/pre>\r\nIt's not obvious that this matrix is badly conditioned, but it is.\r\n<pre class=\"codeinput\">   cond(A)\r\n<\/pre>\r\n<pre class=\"codeoutput\">ans =\r\n   3.0115e+12\r\n<\/pre>\r\nIn case you needed another demonstration that the determinant does not indicate near singularity.\r\n<pre class=\"codeinput\">   format <span class=\"string\">short<\/span>\r\n   det(A)\r\n<\/pre>\r\n<pre class=\"codeoutput\">ans =\r\n    1.0000\r\n<\/pre>\r\nThe elements of the inverse vary over nine orders of magnitude and show perverse scaling that is a direct consequence of the way this matrix was created.\r\n<pre class=\"codeinput\">   format <span class=\"string\">short<\/span> <span class=\"string\">e<\/span>\r\n   X = inv(A)\r\n<\/pre>\r\n<pre class=\"codeoutput\">X =\r\n   2.8508e+09  -1.5513e+08   3.1110e+06   1.4854e+06  -1.6559e+05  -2.0380e+04\r\n  -1.3362e+09   7.2711e+07  -1.4582e+06  -6.9624e+05   7.7615e+04   9.5520e+03\r\n   1.0422e+08  -5.6713e+06   1.1373e+05   5.4306e+04  -6.0540e+03  -7.4500e+02\r\n   1.2586e+07  -6.8489e+05   1.3735e+04   6.5580e+03  -7.3100e+02  -9.0000e+01\r\n  -8.4225e+05   4.5833e+04  -9.1900e+02  -4.3900e+02   4.9000e+01   6.0000e+00\r\n  -1.3830e+05   7.5260e+03  -1.5100e+02  -7.2000e+01   8.0000e+00   1.0000e+00\r\n<\/pre>\r\n<h4>Product of triangular factors<a name=\"509e139e-0e70-4b6f-b86b-b5b4d6036233\"><\/a><\/h4>\r\nIn a 1998 paper Divakar Viswanath and Nick Trefethen demonstrate that the 2-norm condition number of an <i>n<\/i> -by- <i>n<\/i> triangular matrix with normally distributed entries below the diagonal and ones on the diagonal grows exponentially almost surely at a rate <i>c^n<\/i> with <i>c<\/i> = 1.30568...\r\n\r\nI want my matrices to have integer entries so I use randomly distributed integers, but the exponential behavior must be nearly the same.\r\n\r\nHere is the code for <tt>golub<\/tt>.\r\n<pre class=\"codeinput\">   type <span class=\"string\">golub<\/span>\r\n<\/pre>\r\n<pre class=\"codeoutput\">function A = golub(n)\r\n%GOLUB  Badly conditioned integer test matrices.\r\n%   GOLUB(n) is the product of two random integer n-by-n matrices,\r\n%   one of them unit lower triangular and one unit upper triangular.\r\n%   LU factorization without pivoting fails to reveal that such\r\n%   matrices are badly conditioned.\r\n%   See also LUGUI.\r\n\r\n%   Copyright 2014 Cleve Moler\r\n%   Copyright 2014 The MathWorks, Inc.\r\n\r\ns = 10;\r\nL = tril(round(s*randn(n)),-1)+eye(n);\r\nU = triu(round(s*randn(n)),1)+eye(n);\r\nA = L*U;\r\n<\/pre>\r\n<h4>Experiment with pivoting<a name=\"e64f34b8-5c00-489b-bce8-8cb8195b7f5a\"><\/a><\/h4>\r\n<tt>Lugui<\/tt> lets you use the mouse to select the pivot elements for Gaussian elimination. Or, you can choose diagonal pivoting, partial pivoting or complete pivoting. Here are animated gifs of these three choices.\r\n\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/lugui_diagonal.gif\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/lugui_partial.gif\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/lugui_complete.gif\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n<h4>Resurrecting ancestors<a name=\"e335aab0-f3e2-4b06-ac1d-5cf4cf788b43\"><\/a><\/h4>\r\nDiagonal pivoting actually defies the conventional wisdom and picks the <i>smallest<\/i> possible element at each stage as the pivot, namely 1.0. But this choice means that division by the pivot involves no roundoff. Moreover the choice leads to the reconstruction of the original triangular factors.\r\n<pre class=\"codeinput\">   [L,U] = lugui_noshow(A,<span class=\"string\">'diagonal'<\/span>)\r\n   cond_L = cond(L)\r\n   cond_U = cond(U)\r\n<\/pre>\r\n<pre class=\"codeoutput\">L =\r\n     1     0     0     0     0     0\r\n    18     1     0     0     0     0\r\n   -23    36     1     0     0     0\r\n     9    28    -2     1     0     0\r\n     3   -13    -1     7     1     0\r\n   -13    30    15    16    -8     1\r\nU =\r\n     1     3    11     0   -11   -15\r\n     0     1    11    15     0    -7\r\n     0     0     1    -8     6   -11\r\n     0     0     0     1    11    24\r\n     0     0     0     0     1    -6\r\n     0     0     0     0     0     1\r\ncond_L =\r\n   8.4806e+06\r\ncond_U =\r\n   7.9889e+05\r\n<\/pre>\r\n<h4>Diagonal of U<a name=\"1efecaec-2cf6-46c5-a964-5194450d0480\"><\/a><\/h4>\r\nPartial pivoting gives only a weak indication that this matrix is badly conditioned.\r\n<pre class=\"codeinput\">   [L,U,p] = lugui_noshow(A,<span class=\"string\">'partial'<\/span>);\r\n   condition_estimate_partial = max(abs(diag(U)))\/min(abs(diag(U)))\r\n<\/pre>\r\n<pre class=\"codeoutput\">condition_estimate_partial =\r\n   5.8208e+06\r\n<\/pre>\r\nComplete pivoting does a much better job.\r\n<pre class=\"codeinput\">   [L,U,p,q] = lugui_noshow(A,<span class=\"string\">'complete'<\/span>);\r\n   condition_estimate_complete = max(abs(diag(U)))\/min(abs(diag(U)))\r\n<\/pre>\r\n<pre class=\"codeoutput\">condition_estimate_complete =\r\n   1.5166e+12\r\n<\/pre>\r\n<h4>Orthogonal factorization<a name=\"63db1065-2ccb-49b3-8edb-9a0a9b7e7123\"><\/a><\/h4>\r\nThe QR decomposition needs column pivoting to reveal the ill conditioning.\r\n<pre class=\"codeinput\">   [Q,R] = qr(A);\r\n   condition_estimate_without = max(abs(diag(R)))\/min(abs(diag(R)))\r\n\r\n   [Q,R,E] = qr(A);\r\n   condition_estimate_with = max(abs(diag(R)))\/min(abs(diag(R)))\r\n<\/pre>\r\n<pre class=\"codeoutput\">condition_estimate_without =\r\n   6.6062e+06\r\ncondition_estimate_with =\r\n   2.2595e+12\r\n<\/pre>\r\n<h4>References<a name=\"df687528-dcf1-430e-88f4-a23dc7092eb0\"><\/a><\/h4>\r\nDivakar Viswanath and Nick Trefethen, \"Condition Numbers of Random Triangular Matrices,\" <i>SIAM J. Matrix Anal. Appl.<\/i> 19, 564-581, 1998, <a href=\"http:\/\/epubs.siam.org\/doi\/pdf\/10.1137\/S0895479896312869\">&lt;http:\/\/epubs.siam.org\/doi\/pdf\/10.1137\/S0895479896312869<\/a>&gt;. Also available at <a href=\"https:\/\/people.maths.ox.ac.uk\/trefethen\/publication\/PDF\/1998_77.pdf\">https:\/\/people.maths.ox.ac.uk\/trefethen\/publication\/PDF\/1998_77.pdf<\/a>\r\n\r\nCleve Moler, <i>Numerical Computing with MATLAB<\/i>, 2004. <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/37976-numerical-computing-with-matlab\">https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/37976-numerical-computing-with-matlab<\/a>.\r\n\r\n<script>\/\/ <![CDATA[\r\nfunction grabCode_f64f73f9f8ae4eac9b63af2e6ab39f89() {\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='f64f73f9f8ae4eac9b63af2e6ab39f89 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' f64f73f9f8ae4eac9b63af2e6ab39f89';\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('\r\n\r\n\r\n\r\n\r\n\r\n\r\n\r\n\r\n\r\n\r\n\r\n<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>\r\n\r\n\r\n\r\n\r\n\r\n\r\n\r\n\r\n\r\n\r\n\r\n\r\n\\n');\r\n\r\n        d.title = title + ' (MATLAB code)';\r\n        d.close();\r\n    }\r\n\/\/ ]]><\/script>\r\n<p style=\"text-align: right; font-size: xx-small; font-weight: lighter; font-style: italic; color: gray;\"><a><span style=\"font-size: x-small; font-style: italic;\">Get\r\nthe MATLAB code<noscript>(requires JavaScript)<\/noscript><\/span><\/a><\/p>\r\nPublished with MATLAB\u00ae R2015a\r\n\r\n<\/div>\r\n<!--\r\nf64f73f9f8ae4eac9b63af2e6ab39f89 ##### SOURCE BEGIN #####\r\n%% Golub Matrices, Deceptively Ill Conditioned\r\n% I got the idea for these test matrices from Gene Golub years ago,\r\n% but the mathematical foundation comes from a paper by Divakar Viswanath\r\n% and Nick Trefethen.\r\n\r\n%% Two programs from NCM\r\n% This post is about two MATLAB programs included in the collection from\r\n% <https:\/\/www.mathworks.com\/moler.htmlncmfilelist.html % _Numerical Computing with MATLAB_>.  The program\r\n% <https:\/\/www.mathworks.com\/moler.htmlncm\/golub.m |golub|> generates\r\n% badly conditioned integer test matrices.\r\n% The program <https:\/\/www.mathworks.com\/moler.htmlncm\/lugui.m |lugui|> is\r\n% an interactive graphical interface that allows you to experiment with\r\n% pivot strategies in Gaussian elimination.\r\n\r\n%% A 6-by-6 example\r\n% Set the random number generator to the state it has when MATLAB starts\r\n% and generate a 6-by-6 example.\r\n\r\nrng('default')\r\nA = golub(6)\r\n\r\n%%\r\n% It's not obvious that this matrix is badly conditioned, but it is.\r\n\r\ncond(A)\r\n\r\n%%\r\n% In case you needed another demonstration that the determinant does not\r\n% indicate near singularity.\r\n\r\nformat short\r\ndet(A)\r\n\r\n%%\r\n% The elements of the inverse vary over nine orders of magnitude and\r\n% show perverse scaling that is a direct consequence of the way this matrix\r\n% was created.\r\n\r\nformat short e\r\nX = inv(A)\r\n\r\n%% Product of triangular factors\r\n% In a 1998 paper Divakar Viswanath and Nick Trefethen demonstrate that\r\n% the 2-norm condition number of an _n_ -by- _n_ triangular matrix with\r\n% normally distributed entries below the diagonal and ones on the diagonal\r\n% grows exponentially almost surely at a rate _c^n_ with _c_ = 1.30568...\r\n\r\n%%\r\n% I want my matrices to have integer entries so I use randomly distributed\r\n% integers, but the exponential behavior must be nearly the same.\r\n\r\n%%\r\n% Here is the code for |golub|.\r\n\r\ntype golub\r\n\r\n%% Experiment with pivoting\r\n% |Lugui| lets you use the mouse to select the pivot elements for Gaussian\r\n% elimination.  Or, you can choose diagonal pivoting, partial pivoting or\r\n% complete pivoting.  Here are animated gifs of these three choices.\r\n%\r\n% <<lugui_diagonal.gif>>\r\n%\r\n% <<lugui_partial.gif>>\r\n%\r\n% <<lugui_complete.gif>>\r\n\r\n%% Resurrecting ancestors\r\n% Diagonal pivoting actually defies the conventional wisdom and picks the\r\n% _smallest_ possible element at each stage as the pivot, namely 1.0.\r\n% But this choice means that division by the pivot involves no roundoff.\r\n% Moreover the choice leads to the reconstruction of the original triangular\r\n% factors.\r\n\r\n[L,U] = lugui_noshow(A,'diagonal')\r\ncond_L = cond(L)\r\ncond_U = cond(U)\r\n\r\n%% Diagonal of U\r\n% Partial pivoting gives only a weak indication that this matrix is\r\n% badly conditioned.\r\n\r\n[L,U,p] = lugui_noshow(A,'partial');\r\ncondition_estimate_partial = max(abs(diag(U)))\/min(abs(diag(U)))\r\n\r\n%%\r\n% Complete pivoting does a much better job.\r\n\r\n[L,U,p,q] = lugui_noshow(A,'complete');\r\ncondition_estimate_complete = max(abs(diag(U)))\/min(abs(diag(U)))\r\n\r\n%% Orthogonal factorization\r\n% The QR decomposition needs column pivoting to reveal the ill conditioning.\r\n\r\n[Q,R] = qr(A);\r\ncondition_estimate_without = max(abs(diag(R)))\/min(abs(diag(R)))\r\n\r\n[Q,R,E] = qr(A);\r\ncondition_estimate_with = max(abs(diag(R)))\/min(abs(diag(R)))\r\n\r\n%% References\r\n% Divakar Viswanath and Nick Trefethen,\r\n% \"Condition Numbers of Random Triangular Matrices,\"\r\n% _SIAM J. Matrix Anal. Appl._ 19, 564-581, 1998,\r\n% <http:\/\/epubs.siam.org\/doi\/pdf\/10.1137\/S0895479896312869 % http:\/\/epubs.siam.org\/doi\/pdf\/10.1137\/S0895479896312869>.\r\n% Also available at\r\n% <https:\/\/people.maths.ox.ac.uk\/trefethen\/publication\/PDF\/1998_77.pdf % https:\/\/people.maths.ox.ac.uk\/trefethen\/publication\/PDF\/1998_77.pdf>\r\n\r\n%%\r\n% Cleve Moler, _Numerical Computing with MATLAB_, 2004.\r\n% <https:\/\/www.mathworks.com\/moler.htmlncmfilelist.html % https:\/\/www.mathworks.com\/moler.htmlncmfilelist.html>.\r\n\r\n##### SOURCE END ##### f64f73f9f8ae4eac9b63af2e6ab39f89\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img decoding=\"async\"  class=\"img-responsive\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/lugui_diagonal.gif\" onError=\"this.style.display ='none';\" \/><\/div><!--introduction-->I got the idea for these test matrices from Gene Golub years ago, but the mathematical foundation comes from a paper by Divakar Viswanath and Nick Trefethen.\r\n\r\n<!--\/introduction-->... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2015\/08\/24\/golub-matrices-deceptively-ill-conditioned\/\">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,16],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/1246"}],"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=1246"}],"version-history":[{"count":7,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/1246\/revisions"}],"predecessor-version":[{"id":2193,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/1246\/revisions\/2193"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=1246"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=1246"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=1246"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}