Householder Reflections and the QR Decomposition
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>
- Category:
- Algorithms,
- History,
- Matrices,
- People
Comments
To leave a comment, please click here to sign in to your MathWorks Account or create a new one.