{"id":2082,"date":"2016-10-27T18:20:07","date_gmt":"2016-10-27T23:20:07","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=2082"},"modified":"2016-10-27T18:20:07","modified_gmt":"2016-10-27T23:20:07","slug":"apologies-to-gram-schmidt","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2016\/10\/27\/apologies-to-gram-schmidt\/","title":{"rendered":"Apologies to Gram-Schmidt"},"content":{"rendered":"<div class=\"content\"><!--introduction--><p>This is a follow-up to my previous follow-up, <a href=\"https:\/\/blogs.mathworks.com\/cleve\/2016\/10\/17\/compare-gram-schmidt-and-householder-orthogonalization-algorithms-2\/\">posted several days ago.<\/a> A very careful reader, Bruno Bazzano, contributed a comment pointing out what he called \"a small typo\" in my code for the classic Gram-Schmidt algorithm.  It is more than a small typo, it is a serious blunder.  I must correct the code, then do more careful experiments and reword my conclusions.<\/p><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#20595c86-31ae-4909-8d8f-5e81055da2dd\">Classic Gram-Schmidt<\/a><\/li><li><a href=\"#24abd8c3-bb93-4561-bf11-d1ce6f870a76\">Publishing code<\/a><\/li><li><a href=\"#eefde54e-f614-4e58-8859-7e3f11200873\">Comparison<\/a><\/li><li><a href=\"#a72e206e-efac-4b80-879f-8a036854a1f8\">Well conditioned<\/a><\/li><li><a href=\"#729e0b98-7f9a-4da3-8f8a-a6f205f540e5\">Poorly conditioned Hilbert matrices<\/a><\/li><li><a href=\"#082f8bfc-0d86-49ee-ac69-0797c6564011\">Poorly conditioned stabilized magic square<\/a><\/li><li><a href=\"#8fb59b36-7dea-40b2-96e4-d079d3d027b4\">Singular<\/a><\/li><li><a href=\"#3e8d0ed0-6b0c-4e88-a1fc-f9c129b88552\">Conclusion<\/a><\/li><li><a href=\"#4f0d7f2a-fa9d-4fc9-b78f-8a4dcfe78c65\">References<\/a><\/li><\/ul><\/div><h4>Classic Gram-Schmidt<a name=\"20595c86-31ae-4909-8d8f-5e81055da2dd\"><\/a><\/h4><p>Bruno noticed that in the inner loop of my classic Gram-Schmidt function, the expression <tt>Q(:,k-1)<\/tt> should be <tt>Q(:,1:k-1)<\/tt>. I had dropped the two characters, \" <tt>1:<\/tt> \". I need to orthogonalize against all the previous vectors, not just one of them.  Here is the corrected code.<\/p><pre class=\"codeinput\">   type <span class=\"string\">gs<\/span>\r\n<\/pre><pre class=\"codeoutput\">\r\nfunction [Q,R] = gs(X)\r\n    % Classical Gram-Schmidt.  [Q,R] = gs(X);\r\n    % G. W. Stewart, \"Matrix Algorithms, Volume 1\", SIAM, 1998.\r\n    [n,p] = size(X);\r\n    Q = zeros(n,p);\r\n    R = zeros(p,p);\r\n    for k = 1:p\r\n        Q(:,k) = X(:,k);\r\n        if k ~= 1\r\n            R(1:k-1,k) = Q(:,1:k-1)'*Q(:,k);\r\n            Q(:,k) = Q(:,k) - Q(:,1:k-1)*R(1:k-1,k);\r\n        end\r\n        R(k,k) = norm(Q(:,k));\r\n        Q(:,k) = Q(:,k)\/R(k,k);\r\n     end\r\nend\r\n<\/pre><h4>Publishing code<a name=\"24abd8c3-bb93-4561-bf11-d1ce6f870a76\"><\/a><\/h4><p>By the way, this points out the importance of publishing codes.<\/p><h4>Comparison<a name=\"eefde54e-f614-4e58-8859-7e3f11200873\"><\/a><\/h4><p>For various matrices <tt>X<\/tt>, let's again check how the three algorithms perform on two tasks, accuracy and orthogonality.  How close is <tt>Q*R<\/tt> to <tt>X<\/tt>?  And, how close is <tt>Q'*Q<\/tt> to <tt>I<\/tt>.  I have modified the <tt>compare<\/tt> function in last week's post to return <tt>ortherr<\/tt>.<\/p><h4>Well conditioned<a name=\"a72e206e-efac-4b80-879f-8a036854a1f8\"><\/a><\/h4><p>First, try a well conditioned matrix, a magic square of odd order.<\/p><pre class=\"codeinput\">   n = 7;\r\n   X = magic(n)\r\n<\/pre><pre class=\"codeoutput\">X =\r\n    30    39    48     1    10    19    28\r\n    38    47     7     9    18    27    29\r\n    46     6     8    17    26    35    37\r\n     5    14    16    25    34    36    45\r\n    13    15    24    33    42    44     4\r\n    21    23    32    41    43     3    12\r\n    22    31    40    49     2    11    20\r\n<\/pre><p>Check its condition.<\/p><pre class=\"codeinput\">   kappa = condest(X)\r\n<\/pre><pre class=\"codeoutput\">kappa =\r\n    8.5681\r\n<\/pre><p>Do the comparison.<\/p><pre class=\"codeinput\">   compare(X);\r\n<\/pre><pre class=\"codeoutput\">\r\n                 Classic   Modified  Householder\r\nQR error        1.52e-16   6.09e-17   5.68e-16\r\nOrthogonality   1.87e-15   1.53e-15   1.96e-15\r\n<\/pre><p>The <tt>orthogonality<\/tt> value in the <tt>Classic<\/tt> column is now at roundoff level.  All three algorithms do well with both accuracy and orthogonality on this matrix.<\/p><h4>Poorly conditioned Hilbert matrices<a name=\"729e0b98-7f9a-4da3-8f8a-a6f205f540e5\"><\/a><\/h4><p>I want to do two experiments with matrices that are poorly conditioned, but not exactly singular.  First, Hilbert matrices of increasing order.<\/p><pre class=\"codeinput\">   kappa = zeros(8,1);\r\n   ortherr = zeros(8,3);\r\n   <span class=\"keyword\">for<\/span> n = 1:8\r\n       X = hilb(n);\r\n       kappa(n) = condest(X);\r\n       ortherr(n,:) = compare(X);\r\n   <span class=\"keyword\">end<\/span>\r\n<\/pre><p>Plot the results with a logarithmic scale.<\/p><pre class=\"codeinput\">    semilogy([ortherr eps*kappa eps*kappa.^2],<span class=\"string\">'o-'<\/span>)\r\n    axis([0 9 10^(-16) 10^5])\r\n    legend({<span class=\"string\">'gs'<\/span>,<span class=\"string\">'mgs'<\/span>,<span class=\"string\">'house'<\/span>,<span class=\"string\">'eps*kappa'<\/span>,<span class=\"string\">'eps*kappa^2'<\/span>}, <span class=\"keyword\">...<\/span>\r\n        <span class=\"string\">'location'<\/span>,<span class=\"string\">'northwest'<\/span>)\r\n    title(<span class=\"string\">'Loss of orthogonality'<\/span>)\r\n    xlabel(<span class=\"string\">'n'<\/span>)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/apologies_blog_01.png\" alt=\"\"> <p>The loss of orthogonality of GS is worse than the other two and roughly varies like the square of the condition number, <tt>kappa<\/tt>. The loss of orthogonality of MGS varies like <tt>kappa<\/tt> to the first power. Householder does well with orthogonality.<\/p><h4>Poorly conditioned stabilized magic square<a name=\"082f8bfc-0d86-49ee-ac69-0797c6564011\"><\/a><\/h4><p>Another set of nonsingular, but poorly conditioned, matrices, generated by adding a small value to the diagonal of the magic square of order 8<\/p><pre class=\"codeinput\">   M = magic(8);\r\n   I = eye(8);\r\n   <span class=\"keyword\">for<\/span> k = 1:8\r\n       s = 10^(-k);\r\n       X = M + s*I;\r\n       kappa(k) = condest(X);\r\n       ortherr(k,:) = compare(X);\r\n   <span class=\"keyword\">end<\/span>\r\n<\/pre><p>Plot the results with a logarithmic scale.<\/p><pre class=\"codeinput\">    semilogy([ortherr eps*kappa eps*kappa.^2],<span class=\"string\">'o-'<\/span>)\r\n    axis([0 9 10^(-16) 10^5])\r\n    legend({<span class=\"string\">'gs'<\/span>,<span class=\"string\">'mgs'<\/span>,<span class=\"string\">'house'<\/span>,<span class=\"string\">'eps*kappa'<\/span>,<span class=\"string\">'eps*kappa^2'<\/span>}, <span class=\"keyword\">...<\/span>\r\n        <span class=\"string\">'location'<\/span>,<span class=\"string\">'northwest'<\/span>)\r\n    title(<span class=\"string\">'Loss of orthogonality'<\/span>)\r\n    xlabel(<span class=\"string\">'log10(1\/s)'<\/span>)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/apologies_blog_02.png\" alt=\"\"> <p>Again, the orthogonality of GS varies like <tt>kappa^2<\/tt>, and the orthogonality of MGS varies like <tt>kappa<\/tt>.<\/p><h4>Singular<a name=\"8fb59b36-7dea-40b2-96e4-d079d3d027b4\"><\/a><\/h4><p>Finally, an exactly singular matrix, a magic square of even order.<\/p><pre class=\"codeinput\">   n = 8;\r\n   X = magic(n);\r\n   compare(X);\r\n<\/pre><pre class=\"codeoutput\">\r\n                 Classic   Modified  Householder\r\nQR error        4.78e-17   8.54e-17   4.85e-16\r\nOrthogonality   5.17e+00   2.16e+00   1.30e-15\r\n<\/pre><p>All three algorithms do well with accuracy. Both Gram-Schmidts fail completely with orthogonality. Householder still does well with orthogonality.<\/p><h4>Conclusion<a name=\"3e8d0ed0-6b0c-4e88-a1fc-f9c129b88552\"><\/a><\/h4><p>All three of these algorithms provide <tt>Q<\/tt> and <tt>R<\/tt> that do a good job of reproducing the data <tt>X<\/tt>, that is<\/p><div><ul><li><tt>Q<\/tt> * <tt>R<\/tt> is always close to <tt>X<\/tt> for all three algorithms.<\/li><\/ul><\/div><p>On the other hand, their behavior is very different when it comes to producing orthogonality.<\/p><div><ul><li>Classic Gram-Schmidt.  Loss of orthogonality can be much worse than   the other two when <tt>X<\/tt> is badly conditioned.<\/li><li>Modified Gram-Schmidt.  <tt>Q'*Q<\/tt> varies like the condition of <tt>X<\/tt>   and fails completely when <tt>X<\/tt> is singular.<\/li><li>Householder triangularization.  <tt>Q'*Q<\/tt> is always close to <tt>I<\/tt><\/li><\/ul><\/div><h4>References<a name=\"4f0d7f2a-fa9d-4fc9-b78f-8a4dcfe78c65\"><\/a><\/h4><p>G. W. Stewart, <i>Matrix Algorithms: Volume 1: Basic Decompositions<\/i>, SIAM, xix+458, 1998. <a href=\"http:\/\/epubs.siam.org\/doi\/book\/10.1137\/1.9781611971408\">&lt;http:\/\/epubs.siam.org\/doi\/book\/10.1137\/1.9781611971408<\/a>&gt;<\/p><p>Ake Bjorck, <i>Numerical Methods for Least Squares Problems<\/i>, SIAM, xvii+407, 1996. <a href=\"http:\/\/epubs.siam.org\/doi\/book\/10.1137\/1.9781611971484\">&lt;http:\/\/epubs.siam.org\/doi\/book\/10.1137\/1.9781611971484<\/a>&gt;<\/p><script language=\"JavaScript\"> <!-- \r\n    function grabCode_b323e99c1df1417abfff6808775e9ee3() {\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='b323e99c1df1417abfff6808775e9ee3 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' b323e99c1df1417abfff6808775e9ee3';\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 2016 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_b323e99c1df1417abfff6808775e9ee3()\"><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; R2016a<br><\/p><\/div><!--\r\nb323e99c1df1417abfff6808775e9ee3 ##### SOURCE BEGIN #####\r\n%% Apologies to Gram-Schmidt\r\n% This is a follow-up to my previous follow-up,\r\n% <https:\/\/blogs.mathworks.com\/cleve\/2016\/10\/17\/compare-gram-schmidt-and-householder-orthogonalization-algorithms-2\/\r\n% posted several days ago.>\r\n% A very careful reader, Bruno Bazzano, contributed a comment pointing\r\n% out what he called \"a small typo\" in my code for the classic\r\n% Gram-Schmidt algorithm.  It is more than a small typo, it is a\r\n% serious blunder.  I must correct the code, then do more\r\n% careful experiments and reword my conclusions.\r\n\r\n%% Classic Gram-Schmidt\r\n% Bruno noticed that in the inner loop of my classic Gram-Schmidt \r\n% function, the expression |Q(:,k-1)| should be |Q(:,1:k-1)|.\r\n% I had dropped the two characters, \" |1:| \".\r\n% I need to orthogonalize against all the previous vectors,\r\n% not just one of them.  Here is the corrected code.\r\n\r\n   type gs \r\n   \r\n%% Publishing code\r\n% By the way, this points out the importance of publishing codes.\r\n\r\n%% Comparison\r\n% For various matrices |X|, let's again check how the three algorithms\r\n% perform on two tasks, accuracy and orthogonality.  How close is |Q*R|\r\n% to |X|?  And, how close is |Q'*Q| to |I|.  I have modified the\r\n% |compare| function in last week's post to return |ortherr|.\r\n\r\n   \r\n%% Well conditioned\r\n% First, try a well conditioned matrix, a magic square of odd order.\r\n\r\n   n = 7;\r\n   X = magic(n)\r\n\r\n%%\r\n% Check its condition.\r\n\r\n   kappa = condest(X)\r\n   \r\n%%\r\n% Do the comparison.\r\n\r\n   compare(X);\r\n   \r\n%%\r\n% The |orthogonality| value in the |Classic| column is now at roundoff\r\n% level.  All three algorithms do well with both accuracy and\r\n% orthogonality on this matrix.\r\n   \r\n%% Poorly conditioned Hilbert matrices\r\n% I want to do two experiments with matrices that are poorly conditioned,\r\n% but not exactly singular.  First, Hilbert matrices of increasing order.\r\n\r\n   kappa = zeros(8,1);\r\n   ortherr = zeros(8,3);\r\n   for n = 1:8\r\n       X = hilb(n);\r\n       kappa(n) = condest(X);\r\n       ortherr(n,:) = compare(X);\r\n   end\r\n\r\n%%\r\n% Plot the results with a logarithmic scale.\r\n\r\n    semilogy([ortherr eps*kappa eps*kappa.^2],'o-')\r\n    axis([0 9 10^(-16) 10^5])\r\n    legend({'gs','mgs','house','eps*kappa','eps*kappa^2'}, ...\r\n        'location','northwest')\r\n    title('Loss of orthogonality')\r\n    xlabel('n')\r\n\r\n%%\r\n% The loss of orthogonality of GS is worse than the other two and \r\n% roughly varies like the square of the condition number, |kappa|.\r\n% The loss of orthogonality of MGS varies like |kappa| to the first power.\r\n% Householder does well with orthogonality.\r\n\r\n%% Poorly conditioned stabilized magic square\r\n% Another set of nonsingular, but poorly conditioned, matrices,\r\n% generated by adding a small value to the diagonal of the\r\n% magic square of order 8\r\n\r\n   M = magic(8);\r\n   I = eye(8);\r\n   for k = 1:8\r\n       s = 10^(-k);\r\n       X = M + s*I;\r\n       kappa(k) = condest(X);\r\n       ortherr(k,:) = compare(X);\r\n   end\r\n   \r\n%%\r\n% Plot the results with a logarithmic scale.\r\n\r\n    semilogy([ortherr eps*kappa eps*kappa.^2],'o-')\r\n    axis([0 9 10^(-16) 10^5])\r\n    legend({'gs','mgs','house','eps*kappa','eps*kappa^2'}, ...\r\n        'location','northwest')\r\n    title('Loss of orthogonality')\r\n    xlabel('log10(1\/s)')\r\n %%\r\n% Again, the orthogonality of GS varies like |kappa^2|,\r\n% and the orthogonality of MGS varies like |kappa|.\r\n\r\n%% Singular\r\n% Finally, an exactly singular matrix, a magic square of even order.\r\n\r\n   n = 8;\r\n   X = magic(n);   \r\n   compare(X);\r\n   \r\n%%\r\n% All three algorithms do well with accuracy. \r\n% Both Gram-Schmidts fail completely with orthogonality.\r\n% Householder still does well with orthogonality.\r\n   \r\n%% Conclusion\r\n% All three of these algorithms provide |Q| and |R| that do a good job\r\n% of reproducing the data |X|, that is\r\n%\r\n% * |Q| * |R| is always close to |X| for all three algorithms.\r\n\r\n%%\r\n% On the other hand, their behavior is very different when it comes to\r\n% producing orthogonality.\r\n%\r\n% * Classic Gram-Schmidt.  Loss of orthogonality can be much worse than\r\n%   the other two when |X| is badly conditioned.\r\n% * Modified Gram-Schmidt.  |Q'*Q| varies like the condition of |X|\r\n%   and fails completely when |X| is singular.\r\n% * Householder triangularization.  |Q'*Q| is always close to |I|\r\n   \r\n%% References\r\n% G. W. Stewart,\r\n% _Matrix Algorithms: Volume 1: Basic Decompositions_,\r\n% SIAM, xix+458, 1998.\r\n% <http:\/\/epubs.siam.org\/doi\/book\/10.1137\/1.9781611971408\r\n% http:\/\/epubs.siam.org\/doi\/book\/10.1137\/1.9781611971408>\r\n%\r\n% Ake Bjorck,\r\n% _Numerical Methods for Least Squares Problems_,\r\n% SIAM, xvii+407, 1996.\r\n% <http:\/\/epubs.siam.org\/doi\/book\/10.1137\/1.9781611971484\r\n% http:\/\/epubs.siam.org\/doi\/book\/10.1137\/1.9781611971484>\r\n\r\n   \r\n\r\n\r\n##### SOURCE END ##### b323e99c1df1417abfff6808775e9ee3\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img decoding=\"async\"  class=\"img-responsive\" src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/apologies_blog_01.png\" onError=\"this.style.display ='none';\" \/><\/div><!--introduction--><p>This is a follow-up to my previous follow-up, <a href=\"https:\/\/blogs.mathworks.com\/cleve\/2016\/10\/17\/compare-gram-schmidt-and-householder-orthogonalization-algorithms-2\/\">posted several days ago.<\/a> A very careful reader, Bruno Bazzano, contributed a comment pointing out what he called \"a small typo\" in my code for the classic Gram-Schmidt algorithm.  It is more than a small typo, it is a serious blunder.  I must correct the code, then do more careful experiments and reword my conclusions.... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2016\/10\/27\/apologies-to-gram-schmidt\/\">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,9,6,16],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/2082"}],"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=2082"}],"version-history":[{"count":1,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/2082\/revisions"}],"predecessor-version":[{"id":2083,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/2082\/revisions\/2083"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=2082"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=2082"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=2082"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}