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.

Contents

Alston Householder

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.

Alston served as President of both SIAM and ACM. He introduced what he called elementary Hermitian matrices in a paper in the Journal of the ACM in 1958.

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 demos directory.

load gatlin
image(X)
colormap(map)

The Gatlinburg Conferences are now called the Householder Conferences. I wrote about them in MATLAB News & Notes.

Pete Stewart

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".

Pete was Householder's Ph.D. student at the University of Tennessee.

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.

QR Decomposition

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".

The decomposition is available explicitly from the MATLAB function qr. But, more importantly, the decomposition is part of the fundamental MATLAB linear equation solver denoted by backslash, " \ ", as well as both the eig and svd functions for dense matrices.

For any matrix $A$ we write

$$ A = QR $$

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.

A nifty example, taken from the Wikipedia page on "QR decomposition", is unusual because $A$, $R$, and a renormalization of $Q$ all have integer entries.

   A = [12 -51 4; 6 167 -68; -4 24 -41]
   [Q,R] = qr(A)

A =
    12   -51     4
     6   167   -68
    -4    24   -41
Q =
   -0.8571    0.3943    0.3314
   -0.4286   -0.9029   -0.0343
    0.2857   -0.1714    0.9429
R =
  -14.0000  -21.0000   14.0000
         0 -175.0000   70.0000
         0         0  -35.0000

Scale $Q$ to get integers.

   cQ = 175*Q

cQ =
 -150.0000   69.0000   58.0000
  -75.0000 -158.0000   -6.0000
   50.0000  -30.0000  165.0000

Let's check that $Q$ and $R$ reproduce $A$.

   QR = Q*R

QR =
   12.0000  -51.0000    4.0000
    6.0000  167.0000  -68.0000
   -4.0000   24.0000  -41.0000

Householder reflections

A Householder reflection is characterized by a vector $u$, which, following Pete's convention, is normalized to have

$$ ||u|| = \sqrt{2} $$

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.

   H = @(u,x) x - u*(u'*x);

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$.)

   house_pic

house_gen

Where do we get the vector u 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 x. We want to find u so that Hx = H(u,x) is zero below the first element. And, since the reflection is to preserve length, the only nonzero element in Hx should have abs(Hx(1)) == norm(x).

Start with x normalized to have length one.

u = x/norm(x)

Now add +1 or -1 to u(1), choosing the sign to match. The other choice involves cancellation and is less desirable numerically.

u(1) = u(1) + sign(u(1))

Finally, scale u by sqrt(abs(u(1))).

u = u/sqrt(abs(u(1)))

It turns out that this makes norm(u) equal to sqrt(2). And a bit more algebra shows that the resulting reflection zeros out all but the first element of x.

The figure above illustrates the situation because Hx has only one nonzero component.

Here is the code, including the case where x is already all zero.

   type house_gen

function u = house_gen(x)
    % u = house_gen(x)
    % Generate Householder reflection.
    % u = house_gen(x) returns u with norm(u) = sqrt(2), and
    % H(u,x) = x - u*(u'*x) = -+ norm(x)*e_1.
    
    % Modify the sign function so that sign(0) = 1.
    sig = @(u) sign(u) + (u==0);
    
    nu = norm(x);
    if nu ~= 0
        u = x/nu;
        u(1) = u(1) + sig(u(1));
        u = u/sqrt(abs(u(1)));
    else
        u = x;
        u(1) = sqrt(2);
    end
end

Let's try this on the first column of the Wikipedia example

   x = [12 6 -4]'

x =
    12
     6
    -4

The vector u generated is

   u  = house_gen(x)

u =
    1.3628
    0.3145
   -0.2097

The resulting reflection has the desired effect.

   Hx = H(u,x)

Hx =
  -14.0000
   -0.0000
    0.0000

Hx(1) is equal to -norm(x) and the other elements of Hx are zero.

Householder matrix

Our anonymous function can be represented by a matrix. This is the usual way of defining Householder reflections

   I = eye(3);
   M = H(u,I)

M =
   -0.8571   -0.4286    0.2857
   -0.4286    0.9011    0.0659
    0.2857    0.0659    0.9560

In general

$$M = I - uu'$$

Multiplication by this matrix produces the same reflection as the anonymous function, but requires $n^2$ operations, instead of $2n$.

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.

house_qr

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 u 's.

Here is where the way we have expressed the anonymous function is important. When x is a vector, u'*x is a scalar and we could have written it in front of u, like this.

H = @(u,x) x - (u'*x)*u;

But we want to apply H to several columns of a matrix at once, so we have written (u'*x) after u, like this,

H = @(u,x) x - u*(u'*x);

   type house_qr

function [R,U] = house_qr(A)
    % Householder reflections for QR decomposition.
    % [R,U] = house_qr(A) returns
    % R, the upper triangular factor, and
    % U, the reflector generators for use by house_apply.    
    H = @(u,x) x - u*(u'*x);
    [m,n] = size(A);
    U = zeros(m,n);
    R = A;
    for j = 1:min(m,n)
        u = house_gen(R(j:m,j));
        U(j:m,j) = u;
        R(j:m,j:n) = H(u,R(j:m,j:n));
        R(j+1:m,j) = 0;
    end
end

Magic square example

Use a magic square as an example.

   A = magic(6)

A =
    35     1     6    26    19    24
     3    32     7    21    23    25
    31     9     2    22    27    20
     8    28    33    17    10    15
    30     5    34    12    14    16
     4    36    29    13    18    11

Compute its decomposition.

   [R,U] = house_qr(A)

R =
  -56.3471  -16.4693  -30.0459  -39.0969  -38.0321  -38.6710
         0  -54.2196  -34.8797  -23.1669  -25.2609  -23.2963
         0         0   32.4907   -8.9182  -11.2895   -7.9245
         0         0         0   -7.6283    3.9114   -7.4339
         0         0         0         0   -3.4197   -6.8393
         0         0         0         0         0   -0.0000
U =
    1.2732         0         0         0         0         0
    0.0418    1.2568         0         0         0         0
    0.4321    0.0451   -1.1661         0         0         0
    0.1115    0.3884    0.4557    1.0739         0         0
    0.4182   -0.0108    0.5942   -0.6455    1.0796         0
    0.0558    0.5171    0.2819   -0.6558   -0.9135    1.4142

The fact that R(6,6) is equal to zero tells us that the magic square of order 6 has rank less than 6 and so is singular

house_apply

   type house_apply
   type house_apply_transpose

function Z = house_apply(U,X)
    % Apply Householder reflections.
    % Z = house_apply(U,X), with U from house_qr
    % computes Q*X without actually computing Q.
    H = @(u,x) x - u*(u'*x);
    Z = X;
    [~,n] = size(U);
    for j = n:-1:1
        Z = H(U(:,j),Z);
    end
end

function Z = house_apply_transpose(U,X)
    % Apply Householder transposed reflections.
    % Z = house_apply(U,X), with U from house_qr
    % computes Q'*X without actually computing Q'.
    H = @(u,x) x - u*(u'*x);
    Z = X;
    [~,n] = size(U);
    for j = 1:n
        Z = H(U(:,j),Z);
    end
end

Q at last

We can use house_apply to get the matrix $Q$ of the QR decomposition by applying the transformations to the identity matrix.

   I = eye(size(U));
   Q = house_apply(U,I)

Q =
   -0.6211    0.1702   -0.2070   -0.4998    0.2062    0.5000
   -0.0532   -0.5740   -0.4500   -0.2106   -0.6487   -0.0000
   -0.5502    0.0011   -0.4460    0.4537    0.2062   -0.5000
   -0.1420   -0.4733    0.3763   -0.5034    0.3329   -0.5000
   -0.5324    0.0695    0.6287    0.2096   -0.5220   -0.0000
   -0.0710   -0.6424    0.1373    0.4501    0.3329    0.5000

Check that Q is orthogonal.

   QQ = Q'*Q

QQ =
    1.0000   -0.0000    0.0000    0.0000   -0.0000   -0.0000
   -0.0000    1.0000   -0.0000   -0.0000   -0.0000    0.0000
    0.0000   -0.0000    1.0000   -0.0000   -0.0000    0.0000
    0.0000   -0.0000   -0.0000    1.0000   -0.0000   -0.0000
   -0.0000   -0.0000   -0.0000   -0.0000    1.0000    0.0000
   -0.0000    0.0000    0.0000   -0.0000    0.0000    1.0000

And that Q*R regenerates our magic square.

   QR = Q*R

QR =
   35.0000    1.0000    6.0000   26.0000   19.0000   24.0000
    3.0000   32.0000    7.0000   21.0000   23.0000   25.0000
   31.0000    9.0000    2.0000   22.0000   27.0000   20.0000
    8.0000   28.0000   33.0000   17.0000   10.0000   15.0000
   30.0000    5.0000   34.0000   12.0000   14.0000   16.0000
    4.0000   36.0000   29.0000   13.0000   18.0000   11.0000

References

Alston S. Householder, "Unitary Triangularization of a Nonsymmetric Matrix", Journal of the ACM 5, 339-242, 1958. <http://dl.acm.org/citation.cfm?doid=320941.320947>

G. W. Stewart, Matrix Algorithms: Volume 1: Basic Decompositions, SIAM, xix+458, 1998. <http://epubs.siam.org/doi/book/10.1137/1.9781611971408>

Get the MATLAB code

Published with MATLAB® R2016a