{"id":2007,"date":"2016-10-03T12:00:40","date_gmt":"2016-10-03T17:00:40","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=2007"},"modified":"2016-10-03T18:02:12","modified_gmt":"2016-10-03T23:02:12","slug":"householder-reflections-and-the-qr-decomposition","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2016\/10\/03\/householder-reflections-and-the-qr-decomposition\/","title":{"rendered":"Householder Reflections and the QR Decomposition"},"content":{"rendered":"\r\n<div class=\"content\"><!--introduction--><p>The QR decomposition is often the first step in algorithms for solving many different matrix problems, including linear systems, eigenvalues, and singular values.  Householder reflections are the preferred tool for computing the QR decomposition.<\/p><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#39ee045d-7f6e-41f2-b85a-698db74253a6\">Alston Householder<\/a><\/li><li><a href=\"#d19172f6-591a-4323-b3d6-8db9600cc655\">Pete Stewart<\/a><\/li><li><a href=\"#742ad311-5d1b-4b5d-8768-11799d40f72d\">QR Decomposition<\/a><\/li><li><a href=\"#3c155450-18cf-4e5d-a095-c5c9e12e2e28\">Householder reflections<\/a><\/li><li><a href=\"#61f5929e-c83d-41d4-837c-ec9dfb34ec1c\">house_gen<\/a><\/li><li><a href=\"#b5d7aac6-ed0c-4c04-8d18-fde02305f992\">Householder matrix<\/a><\/li><li><a href=\"#3dc3b553-761c-4907-a307-acacce2e07c4\">house_qr<\/a><\/li><li><a href=\"#842737ae-39e8-40c5-a230-612a4b08c1b4\">Magic square example<\/a><\/li><li><a href=\"#807fd7a0-e67e-4659-824e-1c9cb7c8cddb\">house_apply<\/a><\/li><li><a href=\"#23ae622a-619d-4821-8803-0d1321a78b17\">Q at last<\/a><\/li><li><a href=\"#2475446b-1abb-47f1-b3c9-6ee4c614ad73\">References<\/a><\/li><\/ul><\/div><h4>Alston Householder<a name=\"39ee045d-7f6e-41f2-b85a-698db74253a6\"><\/a><\/h4><p>Alston Householder (1904-1993) is one of the pioneers of modern numerical linear algebra.  He was a member of the mathematics division of Oak Ridge National Laboratory for over 20 years, from 1946 until 1969, and was also a Professor at the University of Tennessee.<\/p><p>Alston served as President of both SIAM and ACM. He introduced what he called <i>elementary Hermitian matrices<\/i> in a paper in the Journal of the ACM in 1958.<\/p><p>Alston was head of the organizing committee for the Gatlinburg Conferences on Numerical Algebra.  A photo of the 1964 committee is available in your MATLAB <tt>demos<\/tt> directory.<\/p><pre class=\"language-matlab\">load <span class=\"string\">gatlin<\/span>\r\nimage(X)\r\ncolormap(map)\r\n<\/pre><p>The Gatlinburg Conferences are now called the Householder Conferences. I wrote about them in <a href=\"https:\/\/www.mathworks.com\/company\/newsletters\/articles\/the-gatlinburg-and-householder-symposia.html\">MATLAB News &amp; Notes<\/a>.<\/p><h4>Pete Stewart<a name=\"d19172f6-591a-4323-b3d6-8db9600cc655\"><\/a><\/h4><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 \"Pete\".  He has never been able to satisfactorily explain the origins of \"Pete\" to me.  It somehow goes back through his father to his grandfather and maybe great grandfather, who were also nicknamed \"Pete\".<\/p><p>Pete was Householder's Ph.D. student at the University of Tennessee.<\/p><p>Pete has written several books on numerical linear algebra.  The reference for my blog today is his book \"Matrix Algorithms, Volume I: Basic Decompositions\", published by SIAM. His pseudocode is MATLAB ready.<\/p><h4>QR Decomposition<a name=\"742ad311-5d1b-4b5d-8768-11799d40f72d\"><\/a><\/h4><p>The QR decomposition expresses a matrix as the product of an orthogonal matrix and an upper triangular matrix.  The letter Q is a substitute for the letter O from \"orthogonal\" and the letter R is from \"right\", an alternative for \"upper\".<\/p><p>The decomposition is available explicitly from the MATLAB function <tt>qr<\/tt>. But, more importantly, the decomposition is part of the fundamental MATLAB linear equation solver denoted by backslash, \" <tt>\\<\/tt> \", as well as both the <tt>eig<\/tt> and <tt>svd<\/tt> functions for dense matrices.<\/p><p>For any matrix $A$ we write<\/p><p>$$ A = QR $$<\/p><p>Operations based upon orthogonal matrices are very desirable in numeric computation because they do not magnify errors, either those inherited from the underlying data, or those introduced by floating point arithmetic.<\/p><p>A nifty example, taken from the <i>Wikipedia<\/i> page on \"QR decomposition\", is unusual because $A$, $R$, and a renormalization of $Q$ all have integer entries.<\/p><pre class=\"codeinput\">   A = [12 -51 4; 6 167 -68; -4 24 -41]\r\n   [Q,R] = qr(A)\r\n<\/pre><pre class=\"codeoutput\">A =\r\n    12   -51     4\r\n     6   167   -68\r\n    -4    24   -41\r\nQ =\r\n   -0.8571    0.3943    0.3314\r\n   -0.4286   -0.9029   -0.0343\r\n    0.2857   -0.1714    0.9429\r\nR =\r\n  -14.0000  -21.0000   14.0000\r\n         0 -175.0000   70.0000\r\n         0         0  -35.0000\r\n<\/pre><p>Scale $Q$ to get integers.<\/p><pre class=\"codeinput\">   cQ = 175*Q\r\n<\/pre><pre class=\"codeoutput\">cQ =\r\n -150.0000   69.0000   58.0000\r\n  -75.0000 -158.0000   -6.0000\r\n   50.0000  -30.0000  165.0000\r\n<\/pre><p>Let's check that $Q$ and $R$ reproduce $A$.<\/p><pre class=\"codeinput\">   QR = Q*R\r\n<\/pre><pre class=\"codeoutput\">QR =\r\n   12.0000  -51.0000    4.0000\r\n    6.0000  167.0000  -68.0000\r\n   -4.0000   24.0000  -41.0000\r\n<\/pre><h4>Householder reflections<a name=\"3c155450-18cf-4e5d-a095-c5c9e12e2e28\"><\/a><\/h4><p>A Householder reflection is characterized by a vector $u$, which, following Pete's convention, is normalized to have<\/p><p>$$ ||u|| = \\sqrt{2} $$<\/p><p>It is usual to define a Householder reflection by a matrix.  But I want to define it in a different way, by a MATLAB anonymous function.<\/p><pre class=\"codeinput\">   H = @(u,x) x - u*(u'*x);\r\n<\/pre><p>This emphasizes the way in which the reflection should be computed. For any vector $x$, find its projection onto $u$ and then subtract that multiple of $u$ from $x$.  In the following picture $u_p$ is a line perpendicular to $u$.  In more dimensions $u_p$ would be the plane through the origin perpendicular to $u$.  With the $\\sqrt{2}$ normalization of $u$, the effect of $H$ is to transform any vector $x$ to its mirror image $Hx$ on the other side of this plane.  (Vectors in the plane are left there by $H$.)<\/p><pre class=\"codeinput\">   house_pic\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/house_blog_01.png\" alt=\"\"> <h4>house_gen<a name=\"61f5929e-c83d-41d4-837c-ec9dfb34ec1c\"><\/a><\/h4><p>Where do we get the vector <tt>u<\/tt> that characterizes the reflection? We want to operate on the columns of a matrix to introduce zeros below the diagonal.  Let's begin with the first column, call it <tt>x<\/tt>. We want to find <tt>u<\/tt> so that <tt>Hx = H(u,x)<\/tt> is zero below the first element.  And, since the reflection is to preserve length, the only nonzero element in <tt>Hx<\/tt> should have <tt>abs(Hx(1)) == norm(x)<\/tt>.<\/p><p>Start with <tt>x<\/tt> normalized to have length one.<\/p><pre class=\"language-matlab\">u = x\/norm(x)\r\n<\/pre><p>Now add <tt>+1<\/tt> or <tt>-1<\/tt> to <tt>u(1)<\/tt>, choosing the sign to match. The other choice involves cancellation and is less desirable numerically.<\/p><pre class=\"language-matlab\">u(1) = u(1) + sign(u(1))\r\n<\/pre><p>Finally, scale <tt>u<\/tt> by <tt>sqrt(abs(u(1)))<\/tt>.<\/p><pre class=\"language-matlab\">u = u\/sqrt(abs(u(1)))\r\n<\/pre><p>It turns out that this makes <tt>norm(u)<\/tt> equal to <tt>sqrt(2)<\/tt>.  And a bit more algebra shows that the resulting reflection zeros out all but the first element of <tt>x<\/tt>.<\/p><p>The figure above illustrates the situation because <tt>Hx<\/tt> has only one nonzero component.<\/p><p>Here is the code, including the case where <tt>x<\/tt> is already all zero.<\/p><pre class=\"codeinput\">   type <span class=\"string\">house_gen<\/span>\r\n<\/pre><pre class=\"codeoutput\">\r\nfunction u = house_gen(x)\r\n    % u = house_gen(x)\r\n    % Generate Householder reflection.\r\n    % u = house_gen(x) returns u with norm(u) = sqrt(2), and\r\n    % H(u,x) = x - u*(u'*x) = -+ norm(x)*e_1.\r\n    \r\n    % Modify the sign function so that sign(0) = 1.\r\n    sig = @(u) sign(u) + (u==0);\r\n    \r\n    nu = norm(x);\r\n    if nu ~= 0\r\n        u = x\/nu;\r\n        u(1) = u(1) + sig(u(1));\r\n        u = u\/sqrt(abs(u(1)));\r\n    else\r\n        u = x;\r\n        u(1) = sqrt(2);\r\n    end\r\nend\r\n<\/pre><p>Let's try this on the first column of the <i>Wikipedia<\/i> example<\/p><pre class=\"codeinput\">   x = [12 6 -4]'\r\n<\/pre><pre class=\"codeoutput\">x =\r\n    12\r\n     6\r\n    -4\r\n<\/pre><p>The vector <tt>u<\/tt> generated is<\/p><pre class=\"codeinput\">   u  = house_gen(x)\r\n<\/pre><pre class=\"codeoutput\">u =\r\n    1.3628\r\n    0.3145\r\n   -0.2097\r\n<\/pre><p>The resulting reflection has the desired effect.<\/p><pre class=\"codeinput\">   Hx = H(u,x)\r\n<\/pre><pre class=\"codeoutput\">Hx =\r\n  -14.0000\r\n   -0.0000\r\n    0.0000\r\n<\/pre><p><tt>Hx(1)<\/tt> is equal to <tt>-norm(x)<\/tt> and the other elements of <tt>Hx<\/tt> are zero.<\/p><h4>Householder matrix<a name=\"b5d7aac6-ed0c-4c04-8d18-fde02305f992\"><\/a><\/h4><p>Our anonymous function can be represented by a matrix.  This is the usual way of defining Householder reflections<\/p><pre class=\"codeinput\">   I = eye(3);\r\n   M = H(u,I)\r\n<\/pre><pre class=\"codeoutput\">M =\r\n   -0.8571   -0.4286    0.2857\r\n   -0.4286    0.9011    0.0659\r\n    0.2857    0.0659    0.9560\r\n<\/pre><p>In general<\/p><p>$$M = I - uu'$$<\/p><p>Multiplication by this matrix produces the same reflection as the anonymous function, but requires $n^2$ operations, instead of $2n$.<\/p><p>Recall that Gaussian elimination can be described by a sequence of matrix multiplications, but it is actually carried out by subtracting multiples of rows from other rows.  There is an anonymous function for this elimination operation, but it is not so easy to express the pivoting.<\/p><h4>house_qr<a name=\"3dc3b553-761c-4907-a307-acacce2e07c4\"><\/a><\/h4><p>We are now ready to compute the QR decomposition.   Pass over the columns of the input matrix, using Householder reflections to introduce zeros below the diagonal.  This produces the R matrix. It is inefficient and usually not necessary to actually compute Q.  Just save the <tt>u<\/tt> 's.<\/p><p>Here is where the way we have expressed the anonymous function is important.  When <tt>x<\/tt> is a vector, <tt>u'*x<\/tt> is a scalar and we could have written it in front of <tt>u<\/tt>, like this.<\/p><pre class=\"language-matlab\">H = @(u,x) x - (u'*x)*u;\r\n<\/pre><p>But we want to apply <tt>H<\/tt> to several columns of a matrix at once, so we have written <tt>(u'*x)<\/tt> after <tt>u<\/tt>, like this,<\/p><pre class=\"language-matlab\">H = @(u,x) x - u*(u'*x);\r\n<\/pre><pre class=\"codeinput\">   type <span class=\"string\">house_qr<\/span>\r\n<\/pre><pre class=\"codeoutput\">\r\nfunction [R,U] = house_qr(A)\r\n    % Householder reflections for QR decomposition.\r\n    % [R,U] = house_qr(A) returns\r\n    % R, the upper triangular factor, and\r\n    % U, the reflector generators for use by house_apply.    \r\n    H = @(u,x) x - u*(u'*x);\r\n    [m,n] = size(A);\r\n    U = zeros(m,n);\r\n    R = A;\r\n    for j = 1:min(m,n)\r\n        u = house_gen(R(j:m,j));\r\n        U(j:m,j) = u;\r\n        R(j:m,j:n) = H(u,R(j:m,j:n));\r\n        R(j+1:m,j) = 0;\r\n    end\r\nend\r\n<\/pre><h4>Magic square example<a name=\"842737ae-39e8-40c5-a230-612a4b08c1b4\"><\/a><\/h4><p>Use a magic square as an example.<\/p><pre class=\"codeinput\">   A = magic(6)\r\n<\/pre><pre class=\"codeoutput\">A =\r\n    35     1     6    26    19    24\r\n     3    32     7    21    23    25\r\n    31     9     2    22    27    20\r\n     8    28    33    17    10    15\r\n    30     5    34    12    14    16\r\n     4    36    29    13    18    11\r\n<\/pre><p>Compute its decomposition.<\/p><pre class=\"codeinput\">   [R,U] = house_qr(A)\r\n<\/pre><pre class=\"codeoutput\">R =\r\n  -56.3471  -16.4693  -30.0459  -39.0969  -38.0321  -38.6710\r\n         0  -54.2196  -34.8797  -23.1669  -25.2609  -23.2963\r\n         0         0   32.4907   -8.9182  -11.2895   -7.9245\r\n         0         0         0   -7.6283    3.9114   -7.4339\r\n         0         0         0         0   -3.4197   -6.8393\r\n         0         0         0         0         0   -0.0000\r\nU =\r\n    1.2732         0         0         0         0         0\r\n    0.0418    1.2568         0         0         0         0\r\n    0.4321    0.0451   -1.1661         0         0         0\r\n    0.1115    0.3884    0.4557    1.0739         0         0\r\n    0.4182   -0.0108    0.5942   -0.6455    1.0796         0\r\n    0.0558    0.5171    0.2819   -0.6558   -0.9135    1.4142\r\n<\/pre><p>The fact that <tt>R(6,6)<\/tt> is equal to zero tells us that the magic square of order 6 has rank less than 6 and so is singular<\/p><h4>house_apply<a name=\"807fd7a0-e67e-4659-824e-1c9cb7c8cddb\"><\/a><\/h4><pre class=\"codeinput\">   type <span class=\"string\">house_apply<\/span>\r\n   type <span class=\"string\">house_apply_transpose<\/span>\r\n<\/pre><pre class=\"codeoutput\">\r\nfunction Z = house_apply(U,X)\r\n    % Apply Householder reflections.\r\n    % Z = house_apply(U,X), with U from house_qr\r\n    % computes Q*X without actually computing Q.\r\n    H = @(u,x) x - u*(u'*x);\r\n    Z = X;\r\n    [~,n] = size(U);\r\n    for j = n:-1:1\r\n        Z = H(U(:,j),Z);\r\n    end\r\nend\r\n\r\nfunction Z = house_apply_transpose(U,X)\r\n    % Apply Householder transposed reflections.\r\n    % Z = house_apply(U,X), with U from house_qr\r\n    % computes Q'*X without actually computing Q'.\r\n    H = @(u,x) x - u*(u'*x);\r\n    Z = X;\r\n    [~,n] = size(U);\r\n    for j = 1:n\r\n        Z = H(U(:,j),Z);\r\n    end\r\nend\r\n<\/pre><h4>Q at last<a name=\"23ae622a-619d-4821-8803-0d1321a78b17\"><\/a><\/h4><p>We can use <tt>house_apply<\/tt> to get the matrix $Q$ of the QR decomposition by applying the transformations to the identity matrix.<\/p><pre class=\"codeinput\">   I = eye(size(U));\r\n   Q = house_apply(U,I)\r\n<\/pre><pre class=\"codeoutput\">Q =\r\n   -0.6211    0.1702   -0.2070   -0.4998    0.2062    0.5000\r\n   -0.0532   -0.5740   -0.4500   -0.2106   -0.6487   -0.0000\r\n   -0.5502    0.0011   -0.4460    0.4537    0.2062   -0.5000\r\n   -0.1420   -0.4733    0.3763   -0.5034    0.3329   -0.5000\r\n   -0.5324    0.0695    0.6287    0.2096   -0.5220   -0.0000\r\n   -0.0710   -0.6424    0.1373    0.4501    0.3329    0.5000\r\n<\/pre><p>Check that <tt>Q<\/tt> is orthogonal.<\/p><pre class=\"codeinput\">   QQ = Q'*Q\r\n<\/pre><pre class=\"codeoutput\">QQ =\r\n    1.0000   -0.0000    0.0000    0.0000   -0.0000   -0.0000\r\n   -0.0000    1.0000   -0.0000   -0.0000   -0.0000    0.0000\r\n    0.0000   -0.0000    1.0000   -0.0000   -0.0000    0.0000\r\n    0.0000   -0.0000   -0.0000    1.0000   -0.0000   -0.0000\r\n   -0.0000   -0.0000   -0.0000   -0.0000    1.0000    0.0000\r\n   -0.0000    0.0000    0.0000   -0.0000    0.0000    1.0000\r\n<\/pre><p>And that <tt>Q*R<\/tt> regenerates our magic square.<\/p><pre class=\"codeinput\">   QR = Q*R\r\n<\/pre><pre class=\"codeoutput\">QR =\r\n   35.0000    1.0000    6.0000   26.0000   19.0000   24.0000\r\n    3.0000   32.0000    7.0000   21.0000   23.0000   25.0000\r\n   31.0000    9.0000    2.0000   22.0000   27.0000   20.0000\r\n    8.0000   28.0000   33.0000   17.0000   10.0000   15.0000\r\n   30.0000    5.0000   34.0000   12.0000   14.0000   16.0000\r\n    4.0000   36.0000   29.0000   13.0000   18.0000   11.0000\r\n<\/pre><h4>References<a name=\"2475446b-1abb-47f1-b3c9-6ee4c614ad73\"><\/a><\/h4><p>Alston S. Householder, \"Unitary Triangularization of a Nonsymmetric Matrix\", Journal of the ACM 5, 339-242, 1958. <a href=\"http:\/\/dl.acm.org\/citation.cfm?doid=320941.320947\">&lt;http:\/\/dl.acm.org\/citation.cfm?doid=320941.320947<\/a>&gt;<\/p><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><script language=\"JavaScript\"> <!-- \r\n    function grabCode_c61743dea388452e869eb5a578b169b1() {\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='c61743dea388452e869eb5a578b169b1 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' c61743dea388452e869eb5a578b169b1';\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_c61743dea388452e869eb5a578b169b1()\"><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\nc61743dea388452e869eb5a578b169b1 ##### SOURCE BEGIN #####\r\n%% Householder Reflections and the QR Decomposition\r\n% The QR decomposition is often the first step in algorithms for solving\r\n% many different matrix problems, including linear systems,\r\n% eigenvalues, and singular values.  Householder reflections are\r\n% the preferred tool for computing the QR decomposition.\r\n\r\n%% Alston Householder\r\n% Alston Householder (1904-1993) is one of the pioneers of modern\r\n% numerical linear algebra.  He was a member of the mathematics division\r\n% of Oak Ridge National Laboratory for over 20 years, from 1946 until 1969,\r\n% and was also a Professor at the University of Tennessee.\r\n\r\n%%\r\n% Alston served as President of both SIAM and ACM.\r\n% He introduced what he called _elementary Hermitian matrices_ in a paper\r\n% in the Journal of the ACM in 1958.\r\n\r\n%%\r\n% Alston was head of the organizing committee for the Gatlinburg\r\n% Conferences on Numerical Algebra.  A photo of the 1964 committee is\r\n% available in your MATLAB |demos| directory.\r\n%\r\n%   load gatlin\r\n%   image(X)\r\n%   colormap(map)\r\n%\r\n% The Gatlinburg Conferences are now called the Householder Conferences.\r\n% I wrote about them in\r\n% <https:\/\/www.mathworks.com\/company\/newsletters\/articles\/the-gatlinburg-and-householder-symposia.html\r\n% MATLAB News & Notes>.\r\n\r\n%% Pete Stewart\r\n% My colleague and friend G. W. Stewart is a Distinguished University\r\n% Professor Emeritus at the Department of Computer Science, University\r\n% of Maryland.  Everybody knows him as \"Pete\".  He has never been able\r\n% to satisfactorily explain the origins of \"Pete\" to me.  It somehow goes\r\n% back through his father to his grandfather and maybe great grandfather,\r\n% who were also nicknamed \"Pete\".\r\n\r\n%%\r\n% Pete was Householder's Ph.D. student at the University of Tennessee.\r\n\r\n%%\r\n% Pete has written several books on numerical linear algebra.  The\r\n% reference for my blog today is his book \"Matrix Algorithms,\r\n% Volume I: Basic Decompositions\", published by SIAM.\r\n% His pseudocode is MATLAB ready.\r\n\r\n%% QR Decomposition\r\n% The QR decomposition expresses a matrix as the product of an orthogonal\r\n% matrix and an upper triangular matrix.  The letter Q is a substitute\r\n% for the letter O from \"orthogonal\" and the letter R is from \"right\",\r\n% an alternative for \"upper\".\r\n\r\n%%\r\n% The decomposition is available explicitly from the MATLAB function |qr|.\r\n% But, more importantly, the decomposition is part of the fundamental\r\n% MATLAB linear equation solver denoted by backslash, \" |\\| \", as well as\r\n% both the |eig| and |svd| functions for dense matrices.\r\n\r\n%%\r\n% For any matrix $A$ we write\r\n%\r\n% $$ A = QR $$\r\n%\r\n% Operations based upon orthogonal matrices are very desirable in\r\n% numeric computation because they do not magnify errors, either those\r\n% inherited from the underlying data, or those introduced by floating\r\n% point arithmetic.\r\n\r\n%%\r\n% A nifty example, taken from the _Wikipedia_ page on \"QR decomposition\",\r\n% is unusual because $A$, $R$, and a renormalization of $Q$ all have\r\n% integer entries.\r\n\r\n   A = [12 -51 4; 6 167 -68; -4 24 -41]\r\n   [Q,R] = qr(A)\r\n   \r\n%%\r\n% Scale $Q$ to get integers.\r\n\r\n   cQ = 175*Q\r\n   \r\n%%\r\n% Let's check that $Q$ and $R$ reproduce $A$.\r\n   QR = Q*R\r\n   \r\n%% Householder reflections\r\n% A Householder reflection is characterized by a vector $u$,\r\n% which, following Pete's convention, is normalized to have\r\n%\r\n% $$ ||u|| = \\sqrt{2} $$\r\n% \r\n% It is usual to define a Householder reflection by a matrix.  But I want\r\n% to define it in a different way, by a MATLAB anonymous function.\r\n\r\n   H = @(u,x) x - u*(u'*x);\r\n \r\n%%\r\n% This emphasizes the way in which the reflection should be computed.\r\n% For any vector $x$, find its projection onto $u$ and then subtract\r\n% that multiple of $u$ from $x$.  In the following picture $u_p$ is a line\r\n% perpendicular to $u$.  In more dimensions $u_p$ would be the plane\r\n% through the origin perpendicular to $u$.  With the $\\sqrt{2}$\r\n% normalization of $u$, the effect of $H$ is to transform any vector $x$ \r\n% to its mirror image $Hx$ on the other side of this plane.  (Vectors in\r\n% the plane are left there by $H$.)\r\n\r\n   house_pic\r\n   \r\n%% house_gen\r\n% Where do we get the vector |u| that characterizes the reflection?\r\n% We want to operate on the columns of a matrix to introduce zeros\r\n% below the diagonal.  Let's begin with the first column, call it |x|.\r\n% We want to find |u| so that |Hx = H(u,x)| is zero below the first\r\n% element.  And, since the reflection is to preserve length, the only\r\n% nonzero element in |Hx| should have |abs(Hx(1)) == norm(x)|.\r\n\r\n%%\r\n% Start with |x| normalized to have length one.\r\n%\r\n%   u = x\/norm(x)\r\n%\r\n% Now add |+1| or |-1| to |u(1)|, choosing the sign to match.\r\n% The other choice involves cancellation and is less desirable numerically.\r\n%\r\n%   u(1) = u(1) + sign(u(1))\r\n%\r\n% Finally, scale |u| by |sqrt(abs(u(1)))|.\r\n%\r\n%   u = u\/sqrt(abs(u(1)))\r\n% \r\n% It turns out that this makes |norm(u)| equal to |sqrt(2)|.  And a bit\r\n% more algebra shows that the resulting reflection zeros out all but the\r\n% first element of |x|.\r\n\r\n%%\r\n% The figure above illustrates the situation because |Hx| has only one\r\n% nonzero component.\r\n\r\n%%\r\n% Here is the code, including the case where |x| is already all zero.\r\n\r\n   type house_gen\r\n   \r\n%%\r\n% Let's try this on the first column of the _Wikipedia_ example \r\n\r\n   x = [12 6 -4]'\r\n   \r\n%%\r\n% The vector |u| generated is\r\n\r\n   u  = house_gen(x)\r\n   \r\n%%\r\n% The resulting reflection has the desired effect.\r\n\r\n   Hx = H(u,x)\r\n   \r\n%% \r\n% |Hx(1)| is equal to |-norm(x)| and the other elements of |Hx| are zero.\r\n   \r\n%% Householder matrix\r\n% Our anonymous function can be represented by a matrix.  This is the\r\n% usual way of defining Householder reflections\r\n\r\n   I = eye(3);\r\n   M = H(u,I)\r\n   \r\n%%\r\n% In general\r\n%\r\n% $$M = I - uu'$$\r\n\r\n%%\r\n% Multiplication by this matrix produces the same reflection as the\r\n% anonymous function, but requires $n^2$ operations, instead of $2n$.\r\n\r\n%%\r\n% Recall that Gaussian elimination can be described by a sequence of\r\n% matrix multiplications, but it is actually carried out by subtracting\r\n% multiples of rows from other rows.  There is an anonymous function\r\n% for this elimination operation, but it is not so easy to express\r\n% the pivoting.\r\n\r\n%% house_qr\r\n% We are now ready to compute the QR decomposition.   Pass over the\r\n% columns of the input matrix, using Householder reflections to introduce\r\n% zeros below the diagonal.  This produces the R matrix. It is\r\n% inefficient and usually not necessary to actually compute Q.  Just save\r\n% the |u| 's.\r\n\r\n%%\r\n% Here is where the way we have expressed the anonymous function is\r\n% important.  When |x| is a vector, |u'*x| is a scalar and we could\r\n% have written it in front of |u|, like this.\r\n%\r\n%   H = @(u,x) x - (u'*x)*u;\r\n%\r\n% But we want to apply |H| to several columns of a matrix at once, \r\n% so we have written |(u'*x)| after |u|, like this,\r\n%\r\n%   H = @(u,x) x - u*(u'*x);\r\n%\r\n \r\n   type house_qr\r\n   \r\n%% Magic square example\r\n% Use a magic square as an example.\r\n\r\n   A = magic(6)\r\n   \r\n%%\r\n% Compute its decomposition.\r\n\r\n   [R,U] = house_qr(A)\r\n   \r\n%%\r\n% The fact that |R(6,6)| is equal to zero tells us that the magic square\r\n% of order 6 has rank less than 6 and so is singular\r\n   \r\n%% house_apply\r\n\r\n   type house_apply\r\n   type house_apply_transpose\r\n   \r\n%% Q at last\r\n% We can use |house_apply| to get the matrix $Q$ of the QR decomposition\r\n% by applying the transformations to the identity matrix.\r\n\r\n   I = eye(size(U));\r\n   Q = house_apply(U,I)\r\n   \r\n%%\r\n% Check that |Q| is orthogonal.\r\n\r\n   QQ = Q'*Q\r\n\r\n%%\r\n% And that |Q*R| regenerates our magic square.\r\n\r\n   QR = Q*R\r\n   \r\n%% References\r\n% Alston S. Householder,\r\n% \"Unitary Triangularization of a Nonsymmetric Matrix\",\r\n% Journal of the ACM 5, 339-242, 1958.\r\n% <http:\/\/dl.acm.org\/citation.cfm?doid=320941.320947\r\n% http:\/\/dl.acm.org\/citation.cfm?doid=320941.320947>\r\n%\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##### SOURCE END ##### c61743dea388452e869eb5a578b169b1\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img decoding=\"async\"  class=\"img-responsive\" src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/house_blog_01.png\" onError=\"this.style.display ='none';\" \/><\/div><!--introduction--><p>The QR decomposition is often the first step in algorithms for solving many different matrix problems, including linear systems, eigenvalues, and singular values.  Householder reflections are the preferred tool for computing the QR decomposition.... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2016\/10\/03\/householder-reflections-and-the-qr-decomposition\/\">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,4,6,8],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/2007"}],"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=2007"}],"version-history":[{"count":5,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/2007\/revisions"}],"predecessor-version":[{"id":2026,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/2007\/revisions\/2026"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=2007"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=2007"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=2007"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}