{"id":1887,"date":"2016-07-25T12:00:03","date_gmt":"2016-07-25T17:00:03","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=1887"},"modified":"2016-07-26T20:23:14","modified_gmt":"2016-07-27T01:23:14","slug":"compare-gram-schmidt-and-householder-orthogonalization-algorithms","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2016\/07\/25\/compare-gram-schmidt-and-householder-orthogonalization-algorithms\/","title":{"rendered":"Compare Gram-Schmidt and Householder Orthogonalization Algorithms"},"content":{"rendered":"<div class=\"content\"><!--introduction--><\/p>\n<p>Classical Gram-Schmidt and Modified Gram-Schmidt are two algorithms for orthogonalizing a set of vectors.  Householder elementary reflectors can be used for the same task.  The three algorithms have very different roundoff error properties.<\/p>\n<p><!--\/introduction--><\/p>\n<h3>Contents<\/h3>\n<div>\n<ul>\n<li><a href=\"#81c7eda2-ceb0-4d77-ad28-4b34e14b15c2\">G. W. (Pete) Stewart<\/a><\/li>\n<li><a href=\"#fce6f086-0555-4bf6-9f01-a4cc9b1eacf5\">Classic Gram-Schmidt<\/a><\/li>\n<li><a href=\"#56bc9f27-5693-4279-a435-f45848422d17\">X = Q R<\/a><\/li>\n<li><a href=\"#1ad80787-3697-40ef-9e94-2fa47f17100c\">Modified Gram-Schmidt<\/a><\/li>\n<li><a href=\"#be7846a9-920d-44d9-9482-43b278f38c4b\">Householder Reflections<\/a><\/li>\n<li><a href=\"#1e45f288-685a-4ff0-8d1e-1f04a6199cbc\">Householder QR factorization<\/a><\/li>\n<li><a href=\"#c67235bd-03f7-489d-a938-52902903e596\">Comparison<\/a><\/li>\n<li><a href=\"#ff8ccf57-9fbe-438b-9628-1918c9eaeef0\">Reference<\/a><\/li>\n<\/ul>\n<\/div>\n<h4>G. W. (Pete) Stewart<a name=\"81c7eda2-ceb0-4d77-ad28-4b34e14b15c2\"><\/a><\/h4>\n<p>My colleague and friend G. W. Stewart is a Distinguished University Professor Emeritus at the Department of Computer Science, University of Maryland.  Everybody knows him as &#8220;Pete&#8221;.  He has never been able to satisfactorily explain the origins of &#8220;Pete&#8221; to me.  It somehow goes back through his father to his grandfather and maybe great grandfather, who were also nicknamed &#8220;Pete&#8221;.<\/p>\n<p>Pete has written several books on numerical linear algebra.  For my blog today I am going to rely on the descriptions and pseudocode from his book &#8220;Matrix Algorithms, Volume I: Basic Decompositions&#8221;. His pseudocode is MATLAB ready.<\/p>\n<h4>Classic Gram-Schmidt<a name=\"fce6f086-0555-4bf6-9f01-a4cc9b1eacf5\"><\/a><\/h4>\n<p>The classic Gram-Schmidt algorithm is the first thing you might think of for producing an orthogonal set of vectors.  For each vector in your data set, remove its projection onto the data set, normalize what is left, and add it to the orthogonal set.  Here is the code. <tt>X<\/tt> is the original set of vectors, <tt>Q<\/tt> is the resulting set of orthogonal vectors, and <tt>R<\/tt> is the set of coefficients, organized into an upper triangular matrix.<\/p>\n<pre class=\"codeinput\">   type <span class=\"string\">gs<\/span>\r\n<\/pre>\n<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(:,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>\n<h4>X = Q R<a name=\"56bc9f27-5693-4279-a435-f45848422d17\"><\/a><\/h4>\n<p>The entire process can be expressed by matrix multiplication, with orthogonal $Q$ and upper triangular $R$.<\/p>\n<p>$$X = QR$$<\/p>\n<h4>Modified Gram-Schmidt<a name=\"1ad80787-3697-40ef-9e94-2fa47f17100c\"><\/a><\/h4>\n<p>This is a rather different algorithm, not just a simple modification of classical Gram-Schmidt.  The idea is to orthogonalize against the emerging set of vectors instead of against the original set.  There are two variants, a column-oriented one and a row-oriented one.  They produce the same results, in different order.  Here is the column version,<\/p>\n<pre class=\"codeinput\">   type <span class=\"string\">mgs<\/span>\r\n<\/pre>\n<pre class=\"codeoutput\">\r\nfunction [Q,R] =  mgs(X)\r\n    % Modified Gram-Schmidt.  [Q,R] = mgs(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        for i = 1:k-1\r\n            R(i,k) = Q(:,i)'*Q(:,k);\r\n            Q(:,k) = Q(:,k) - R(i,k)*Q(:,i);\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>\n<h4>Householder Reflections<a name=\"be7846a9-920d-44d9-9482-43b278f38c4b\"><\/a><\/h4>\n<p>A Householder reflection $H$ transforms a given vector $x$ into a multiple of a unit vector $e_k$ while preserving length, so $Hx = \\pm \\sigma e_k$ where $\\sigma = ||x||$ .<\/p>\n<p>The matrix $H$ is never formed.  The reflection is expressed as a rank-one modification of the identity<\/p>\n<p>$$H = I &#8211; uu^T \\ \\mbox{where} \\ ||u|| = \\sqrt{2}$$<\/p>\n<p>(The use of $\\sqrt{2}$ here is Pete&#8217;s signature normalization.)<\/p>\n<pre class=\"codeinput\">   type <span class=\"string\">housegen<\/span>\r\n<\/pre>\n<pre class=\"codeoutput\">\r\nfunction [u,nu] = housegen(x)\r\n    % [u,nu] = housegen(x)\r\n    % Generate Householder reflection.\r\n    % G. W. Stewart, \"Matrix Algorithms, Volume 1\", SIAM, 1998.\r\n    % [u,nu] = housegen(x).\r\n    % H = I - uu' with Hx = -+ nu e_1\r\n    %    returns nu = norm(x).\r\n    u = x;\r\n    nu = norm(x);\r\n    if nu == 0\r\n        u(1) = sqrt(2);\r\n        return\r\n    end\r\n    u = x\/nu;\r\n    if u(1) &gt;= 0\r\n        u(1) = u(1) + 1;\r\n        nu = -nu;\r\n    else\r\n        u(1) = u(1) - 1;\r\n    end\r\n    u = u\/sqrt(abs(u(1)));\r\nend\r\n<\/pre>\n<p>Now to apply $H$ to any other $y$, compute the inner product<\/p>\n<p>$$\\tau = u^Ty$$<\/p>\n<p>and then simply subtract<\/p>\n<p>$$Hy = x &#8211; \\tau u$$<\/p>\n<h4>Householder QR factorization<a name=\"1e45f288-685a-4ff0-8d1e-1f04a6199cbc\"><\/a><\/h4>\n<p>This program does not actually compute the QR orthogonalization, but rather computes <tt>R<\/tt> and a matrix <tt>U<\/tt> containing vectors that generate the Householder reflectors whose product is Q.<\/p>\n<pre class=\"codeinput\">   type <span class=\"string\">hqrd<\/span>\r\n<\/pre>\n<pre class=\"codeoutput\">\r\nfunction [U,R] = hqrd(X)  \r\n    % Householder triangularization.  [U,R] = hqrd(X);\r\n    % Generators of Householder reflections stored in U.\r\n    % H_k = I - U(:,k)*U(:,k)'.\r\n    % prod(H_m ... H_1)X = [R; 0]\r\n    % where m = min(size(X))\r\n    % G. W. Stewart, \"Matrix Algorithms, Volume 1\", SIAM, 1998.\r\n    [n,p] = size(X);\r\n    U = zeros(size(X));\r\n    m = min(n,p);\r\n    R = zeros(m,m);\r\n    for k = 1:min(n,p)\r\n        [U(k:n,k),R(k,k)] = housegen(X(k:n,k));\r\n        v = U(k:n,k)'*X(k:n,k+1:p);\r\n        X(k:n,k+1:p) = X(k:n,k+1:p) - U(k:n,k)*v;\r\n        R(k,k+1:p) = X(k,k+1:p);\r\n    end\r\nend\r\n<\/pre>\n<h4>Comparison<a name=\"c67235bd-03f7-489d-a938-52902903e596\"><\/a><\/h4>\n<p>All three of these algorithms compute <tt>Q<\/tt> and <tt>R<\/tt> that do a good job of reproducing the data <tt>X<\/tt>, that is<\/p>\n<div>\n<ul>\n<li><tt>Q<\/tt> * <tt>R<\/tt> is always close to <tt>X<\/tt> for all three algorithms.<\/li>\n<\/ul>\n<\/div>\n<p>On the other hand, their behavior is very different when it comes to producing orthogonality.  Is <tt>Q'*Q<\/tt> close to the identity?<\/p>\n<div>\n<ul>\n<li>Classic Gram-Schmidt.  Usually very poor orthogonality.<\/li>\n<li>Modified Gram-Schmidt.  Depends upon condition of <tt>X<\/tt>.   Fails completely when <tt>X<\/tt> is singular.<\/li>\n<li>Householder triangularization.  Always good orthogonality.<\/li>\n<\/ul>\n<\/div>\n<pre class=\"codeinput\">   type <span class=\"string\">compare.m<\/span>\r\n<\/pre>\n<pre class=\"codeoutput\">\r\nfunction compare(X);\r\n% compare(X)\r\n% Compare three QR decompositions,\r\n%\r\nI = eye(size(X));\r\n\r\n%% Classic Gram Schmidt\r\n[Q,R] = gs(X);\r\nqrerr_gs = norm(Q*R-X,inf)\/norm(X,inf);\r\northerr_gs = norm(Q'*Q-I,inf);\r\n\r\n%% Modified Gram Schmidt\r\n[Q,R] = mgs(X);\r\nqrerr_mgs = norm(Q*R-X,inf)\/norm(X,inf);\r\northerr_mgs = norm(Q'*Q-I,inf);\r\n\r\n%% Householder QR Decomposition\r\n[U,R] = hqrd(X);\r\nQR = R;\r\nE = I;\r\nfor k = size(X,2):-1:1\r\n    uk = U(:,k);\r\n    QR = QR - uk*(uk'*QR);\r\n    E = E - uk*(uk'*E) - (E*uk)*uk' + uk*(uk'*E*uk)*uk';\r\nend\r\nqrerr_h = norm(QR-X,inf)\/norm(X,inf);\r\northerr_h = norm(E-I,inf);\r\n\r\n%% Report results \r\nfprintf('QR error\\n')\r\nfprintf('Classic:     %10.3e\\n',qrerr_gs)\r\nfprintf('Modified:    %10.3e\\n',qrerr_mgs)\r\nfprintf('Householder: %10.3e\\n',qrerr_h)\r\nfprintf('\\n')\r\nfprintf('Orthogonality error\\n')\r\nfprintf('Classic:     %10.3e\\n',ortherr_gs)\r\nfprintf('Modified:    %10.3e\\n',ortherr_mgs)\r\nfprintf('Householder: %10.3e\\n',ortherr_h)\r\n<\/pre>\n<p>Well conditioned.<\/p>\n<pre class=\"codeinput\">   compare(magic(7))\r\n<\/pre>\n<pre class=\"codeoutput\">QR error\r\nClassic:      1.726e-16\r\nModified:     6.090e-17\r\nHouseholder:  3.654e-16\r\n\r\nOrthogonality error\r\nClassic:      3.201e+00\r\nModified:     1.534e-15\r\nHouseholder:  1.069e-15\r\n<\/pre>\n<p>Poorly conditioned, nonsingular.<\/p>\n<pre class=\"codeinput\">   compare(hilb(7))\r\n<\/pre>\n<pre class=\"codeoutput\">QR error\r\nClassic:      5.352e-17\r\nModified:     5.352e-17\r\nHouseholder:  7.172e-16\r\n\r\nOrthogonality error\r\nClassic:      5.215e+00\r\nModified:     1.219e-08\r\nHouseholder:  1.686e-15\r\n<\/pre>\n<p>Singular.<\/p>\n<pre class=\"codeinput\">   compare(magic(8))\r\n<\/pre>\n<pre class=\"codeoutput\">QR error\r\nClassic:      1.435e-16\r\nModified:     8.540e-17\r\nHouseholder:  2.460e-16\r\n\r\nOrthogonality error\r\nClassic:      5.413e+00\r\nModified:     2.162e+00\r\nHouseholder:  2.356e-15\r\n<\/pre>\n<h4>Reference<a name=\"ff8ccf57-9fbe-438b-9628-1918c9eaeef0\"><\/a><\/h4>\n<p>G. W. Stewart, &#8220;Matrix Algorithms, Volume 1: Basic Decompositions&#8221;, SIAM, 1998.<\/p>\n<p><a href=\"https:\/\/www.amazon.com\/Matrix-Algorithms-1-Basic-Decompositions\/dp\/0898714141\">https:\/\/www.amazon.com\/Matrix-Algorithms-1-Basic-Decompositions\/dp\/0898714141<\/a><\/p>\n<p><script language=\"JavaScript\"> <!-- \n    function grabCode_82a105f1559648d4a953ae06677dd8b1() {\n        \/\/ Remember the title so we can use it in the new page\n        title = document.title;\n\n        \/\/ Break up these strings so that their presence\n        \/\/ in the Javascript doesn't mess up the search for\n        \/\/ the MATLAB code.\n        t1='82a105f1559648d4a953ae06677dd8b1 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 82a105f1559648d4a953ae06677dd8b1';\n    \n        b=document.getElementsByTagName('body')[0];\n        i1=b.innerHTML.indexOf(t1)+t1.length;\n        i2=b.innerHTML.indexOf(t2);\n \n        code_string = b.innerHTML.substring(i1, i2);\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\n\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \n        \/\/ in the XML parser.\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\n        \/\/ doesn't go ahead and substitute the less-than character. \n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\n\n        copyright = 'Copyright 2016 The MathWorks, Inc.';\n\n        w = window.open();\n        d = w.document;\n        d.write('\n\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>\n\n\\n');\n\n        d.title = title + ' (MATLAB code)';\n        d.close();\n    }   \n     --> <\/script><\/p>\n<p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><a href=\"javascript:grabCode_82a105f1559648d4a953ae06677dd8b1()\"><span style=\"font-size: x-small;        font-style: italic;\">Get<br \/>\n      the MATLAB code <noscript>(requires JavaScript)<\/noscript><\/span><\/a><\/p>\n<p>      Published with MATLAB&reg; R2016a<\/p>\n<\/div>\n<p><!--\n82a105f1559648d4a953ae06677dd8b1 ##### SOURCE BEGIN #####\n%% Compare Gram-Schmidt and Householder Orthogonalization Algorithms\n% Classical Gram-Schmidt and Modified Gram-Schmidt are two algorithms for\n% orthogonalizing a set of vectors.  Householder elementary reflectors can\n% be used for the same task.  The three algorithms have very\n% different roundoff error properties.\n\n%% G. W. (Pete) Stewart\n% My colleague and friend G. W. Stewart is a Distinguished University\n% Professor Emeritus at the Department of Computer Science, University\n% of Maryland.  Everybody knows him as \"Pete\".  He has never been able\n% to satisfactorily explain the origins of \"Pete\" to me.  It somehow goes\n% back through his father to his grandfather and maybe great grandfather,\n% who were also nicknamed \"Pete\".\n\n%%\n% Pete has written several books on numerical linear algebra.  For my\n% blog today I am going to rely on the descriptions and pseudocode from\n% his book \"Matrix Algorithms, Volume I: Basic Decompositions\".\n% His pseudocode is MATLAB ready.\n\n%% Classic Gram-Schmidt\n% The classic Gram-Schmidt algorithm is the first thing you might thinke\n% of for producing an orthogonal set of vectors.  For each vector in\n% your data set, remove its projection onto the data set, normalize what\n% is left, and add it to the orthogonal set.  Here is the code.\n% |X| is the original set of vectors, |Q| is the resulting set of\n% orthogonal vectors, and |R| is the set of coefficients, organized into\n% an upper triangular matrix.\n\n   type gs\n   \n\n%% X = Q R\n% The entire process can be expressed by matrix multiplication, with\n% orthogonal $Q$ and upper triangular $R$.\n%\n% $$X = QR$$\n\n%% Modified Gram-Schmidt\n% This is a rather different algorithm, not just a simple modification\n% of classical Gram-Schmidt.  The idea is to orthogonalize against the\n% emerging set of vectors instead of against the original set.  There\n% are two variants, a column-oriented one and a row-oriented one.  They\n% produce the same results, in different order.  Here is the column\n% version,\n\n   type mgs\n   \n%% Householder Reflections\n% A Householder reflection $H$ transforms a given vector $x$ into a multiple\n% of a unit vector $e_k$ while preserving length, so $Hx = \\pm \\sigma e_k$\n% where $\\sigma = ||x||$ .\n\n%%\n% The matrix $H$ is never formed.  The reflection is expressed as a \n% rank-one modification of the identity\n%\n% $$H = I - uu^T \\ \\mbox{where} \\ ||u|| = \\sqrt{2}$$\n%\n% (The use of $\\sqrt{2}$ here is Pete's signature normalization.)\n\n   type housegen\n\n%%\n% Now to apply $H$ to any other $y$, compute the inner product\n%\n% $$\\tau = u^Ty$$\n%\n% and then simply subtract\n%\n% $$Hy = x - \\tau u$$\n\n%% Householder QR factorization\n% This program does not actually compute the QR orthogonalization,\n% but rather computes |R| and a matrix |U| containing vectors\n% that generate the Householder reflectors whose product is Q.\n\n   type hqrd\n   \n%% Comparison\n% All three of these algorithms compute |Q| and |R| that do a good job\n% of reproducing the data |X|, that is\n%\n% * |Q| * |R| is always close to |X| for all three algorithms.\n\n%%\n% On the other hand, their behavior is very different when it comes to\n% producing orthogonality.  Is |Q'*Q| close to the identity?\n%\n% * Classic Gram-Schmidt.  Usually very poor orthogonality.\n% * Modified Gram-Schmidt.  Depends upon condition of |X|.\n%   Fails completely when |X| is singular.\n% * Householder triangularization.  Always good orthogonality.\n\n   type compare.m\n   \n%%\n% Well conditioned.\n\n   compare(magic(7))\n   \n%%\n% Poorly conditioned, nonsingular.\n\n   compare(hilb(7))\n   \n%%\n% Singular.\n\n   compare(magic(8))\n   \n%% Reference\n% G. W. Stewart, \"Matrix Algorithms, Volume 1: Basic Decompositions\",\n% SIAM, 1998.\n%\n% <https:\/\/www.amazon.com\/Matrix-Algorithms-1-Basic-Decompositions\/dp\/0898714141\n% https:\/\/www.amazon.com\/Matrix-Algorithms-1-Basic-Decompositions\/dp\/0898714141>\n\n   \n\n\n##### SOURCE END ##### 82a105f1559648d4a953ae06677dd8b1\n--><\/p>\n","protected":false},"excerpt":{"rendered":"<p><!--introduction--><\/p>\n<p>Classical Gram-Schmidt and Modified Gram-Schmidt are two algorithms for orthogonalizing a set of vectors.  Householder elementary reflectors can be used for the same task.  The three algorithms have very different roundoff error properties&#8230;. <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2016\/07\/25\/compare-gram-schmidt-and-householder-orthogonalization-algorithms\/\">read more >><\/a><\/p>\n","protected":false},"author":78,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[11,4,6,16],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/1887"}],"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=1887"}],"version-history":[{"count":2,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/1887\/revisions"}],"predecessor-version":[{"id":1911,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/1887\/revisions\/1911"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=1887"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=1887"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=1887"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}