# Householder Reflections and the QR Decomposition 1

Posted by **Cleve Moler**,

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

**503**views (last 30 days) | |

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

## Recent Comments