QR Factorization Chapter4
QR Factorization Chapter4
4. QR-DECOMPOSITION
4.1 Introduction
In this chapter we show that given any high matrix (m n ) A, that has rank n,
it can be written as a product of an mxm unitary matrices Q and an mxn upper
triangular matrix R:
A = QR.
QR decomposition has various applications. It is used, for example, to find the
least squares solution for a system of linear equations Ax=b (multiplying by a
unitary matrix Q does not change the norm) and surprisingly the best known
method for finding the eigenvalues of a general matrix is based on QR
decomposition. QR decomposition is also used to construct an orthonormal basis
for a subspace S of F
m
. We take the matrix A consisting of the original basis
vectors of the subspace, the first n column vectors of Q now form the basis
required and the rest of the columns form an orthonormal basis for S
.
4.2 Householder matrix
In chapter 2 we noticed that left multiplication of a matrix A by suitable
Gaussian matrices transforms A to an upper triangular matrix. The same kind of
transformation can be done by using unitary Householder matrices.
Multiplication by a unitary matrix does not change the length of a vector, only
the direction, so they are safe to use in numerical computation.
Let v 0 be a vector in F
m
. A Householder matrix is a Hermitian, unitary,
mxm matrix, of the form
P = I 2
vv
v
2
*
.
Exercise 4.2.1 Prove, that P is Hermitian and unitary.
Matrix Algebra 1 QR-decomposition 87
v
<v,x>
v
2
<v,x>
v
v
2
x
v
v
2
<v,x>
{
}
v
P
x
x
Picture 4.2.1 The transformation carried out by a Householder matrix.
The vector x is transformed into a mirror image of itself with respect to the
orthogonal complement {v}
of the vector v.
Let x 0 be a vector in F
m
. We try to choose a vector v such that its associated
Householder matrix gives us,
Px=e
1
.
We have shown in exercise 4.2.1 that P is unitary, so || Px ||=|| x ||, which means
that || =||x ||. We write the solution in the following form. The motivation for
this will become evident later.
Px= e
i
||x ||e
1
.
When we write P in terms of its components, we get the equation
Px = x 2
vv*x
v
2
= x 2
< v,x >
v
2
v = -e
i
x e
1
and further
2
< v,x >
v
2
v = x + e
i
x e
1
.
The left and right sides of the equation must be parallel for a solution to exist.
This means that
v = x + e
i
x e
1
( )
Matrix Algebra 1 QR-decomposition 88
where is an unknown scalar for the time being. Substituting v into the previous
equation gives
1 = 2
< x + e
i
x e
1
( )
,x >
x + e
i
x e
1
( )
2
= 2
< x + e
i
x e
1
( )
,x >
< x + e
i
x e
1
( )
, x + e
i
x e
1
( )
>
= 2
< x + e
i
x e
1
( )
,x >
< x + e
i
x e
1
( )
, x + e
i
x e
1
( )
>
=
2 < x + e
i
x e
1
( )
,x >
< x + e
i
x e
1
( )
, x + e
i
x e
1
( )
>
because does not occur in the equation anymore we must find a coefficient so
that the equation holds. We compute all the inner products and get the equation
1 =
2 x
2
+ e
i
x x
1
x
2
+ x
2
+ e
i
x x
1
+ e
i
x x
1
where x
1
is the first element of the vector x and <e
1
, x>= x
1
. Dividing the
equation by ||x|| gives
1 =
2 x + e
i
x
1
( )
2 x + e
i
x
1
+ e
i
x
1
2 x + e
i
x
1
+ e
i
x
1
= 2 x + e
i
x
1
( )
e
i
x
1
+ e
i
x
1
= 2e
i
x
1
e
i
x
1
= e
i
x
1
In general the element x
1
is a complex number and it can be written in the form
x
1
= e
i
|x
1
|. (x
1
0) Substituting this into the equation above reveals that the
solutions are of the form:
Matrix Algebra 1 QR-decomposition 89
= + k, k= ..-n,..-1,0,1,...,n,....
Only the solutions = and = + give different values for
e
i
. Lets compare the
solutions. The formula for the Householder matrix is easy to compute
numerically if the norm of the vector v does not get too small. This is not a big
problem because
v
2
=
2
2 x
2
+ e
io
x x
1
+ e
io
x x
1
1
]
1
If =, then
v
2
=
2
2 x
2
+ 2 x x
1
If = + , then
v
2
=
2
2 x
2
2 x x
1
.
The first solution of these two is numerically well behaved. If we choose the latter
solution, the vector v could be zero or the vector could become very small. So we
choose = and the vector required is
v = x + e
i
x e
1
, x
1
= e
i
x
1
.
We have also shown that the selection of v does not depend on the scalar , so we
can choose it to be one. If |x
1
|=0, we choose =0.
NB. If the vector x is real and
x
1
is its first element, then
a) if
x
1
0, then
x
1
= e
i0
x
1
, so
= 0.
b) if
x
1
< 0, then
x
1
= e
i
x
1
, so = .
Example 4.2.1 Let x = [3,1,5,1]
T
. We form the Householder matrix P such that
Px span e
1
{ }
.
because 3>0, we choose
= 0 and
v = x + x e
1
=
3
1
5
1
+
6
0
0
0
=
9
1
5
1
.
Matrix Algebra 1 QR-decomposition 90
Substitution into the formula gives us
P = I 2
vv
T
v
2
=
1
54
27 9 45 9
9 53 5 1
45 5 29 5
9 1 5 53
.
and a little computation shows us that
Px =
1
54
27 9 45 9
9 53 5 1
45 5 29 5
9 1 5 53
3
1
5
1
=
-6
0
0
0
= -6e
1
.
4.3 Constructing the QR-decomposition by using Householder
matrices
4.3.1 Introduction
In this section we transform a matrix A to upper triangular form as we did when
constructing the LU-decomposition. The main difference being that we use
Householder matrices instead of Gaussian matrices. Householder matrices are
unitary, so they have excellent numerical properties.
4.3.2 QR decomposition
Let
A = a
1
, a
2
,...., a
n
be an mxn matrix, that has rank rank(A) = n and m n .
1
st
step. rank(A) = n, a
1
0, and so we can choose an mxm Householder matrix
P
1
to satisfy the equation
1
a
1
= e
i
a
1
e
1
= r
11
e
1
.
So
Matrix Algebra 1 QR-decomposition 91
P
1
A = r
11
e
1
, P
1
a
2
,...., P
1
a
n
=
r
11
0
|
|
|
|
|
|
r
12
A
2
( )
r
1n
.
because a
1
0,
r
11
0 .
2
nd
step. We write the (m-1)x(n-1) matrix A
(2)
in terms of its columns,
A
2
( )
= a
2
2
( )
, a
3
2
( )
,...., a
n
2
( )
.
Our previous assumption about the rank says that a
2
(2)
0. (because if a
2
(2)
=0,
then rank(A) = rank(P
1
A) < n , a contradiction.) We can thus form the (m-1)x(m-1)
Householder matrix B
2
such that
B
2
a
2
2
( )
= e
i
a
2
2
( )
1, 0,..., 0
T
= r
22
e
1
.
because a
2
0,
r
22
0 .
Let
P
2
=
1 0
T
0 B
2
.
and left multiply P
1
A with it. Now we have
P
2
P
1
A =
r
11
r
12
r
1n
0 r
22
r
2n
0 0 A
3
( )
.
It is easy to see that P
2
is unitary.
3
rd
step. Again we can form the Householder matrix B
3
such that
B
3
a
3
(3)
= r
33
e
1
,
since a
3
(3)
0 where a
3
(3)
is the first column of the matrix A
(3)
, so we have
a
3
(3)
0 and r
33
0. We next construct
Matrix Algebra 1 QR-decomposition 92
P
3
=
I
2
0
0 B
3
.
Now the first three columns of P
3
P
2
P
1
A are in upper triangular form.
n
th
step. The matrix A
(n)
is a m-n+1-vector. According to the rank assumption
A
(n)
0. If m = n, we can stop. If m > n, we can find a Householder matrix B
n
such that
B
n
A
n
( )
= r
nn
e
1
where r
nn
0. We let
P
n
=
I
n1
0
0 B
n
.
We denote Q* = P
n
.....P
2
P
1
which gives us the solution we wanted
Q* A = R =
R
n
0
,
where R is mxn upper triangular matrix ( R
n
is nxn upper triangular matrix) ,
and Q is unitary as it is a product of unitary mxm matrices.
Example 4.3.1 The QR decomposition of the matrix
A =
1 1 1
1 2 1
1 2 1
is
1
3
+
2
3
0
1
3
1
6
1
2
1
3
1
6
+
1
2
3 3
5
3
+
3
3
0
2
3
2
3
0 0 2
.
Exercise 4.3.2. Assume that A is a square matrix and rank(A) = n. Show that Q
can be chosen so that the diagonal elements of R are positive.
Matrix Algebra 1 QR-decomposition 93
(Hint: If the vectors in a orthogonal set {q
i
} are multiplied by complex numbers
i
, that have absolute value |
i
|=1, then the new set of vectors {
i
q
i
} is
orthogonal as well.)
Use this result to transform the diagonal elements of the matrix R to positive
numbers in the QR decomposition of the example 4.2.1.
Answer is
QR =
1
3
2
3
0
1
3
1
6
1
2
1
3
1
6
1
2
3 3
5
3
3
3
0
2
3
2
3
0 0 2
.
Example 4.3.2 Computing the QR decomposition with a pen and paper is a
laborious task. Programming it is easy. Let
A =
0
1
0
0
3
2
.
Lets construct its QR decomposition.
1st step
v =
1
1
0
, P
1
=
0 1 0
1 0 0
0 0 1
P
1
A =
1
0
0
3
0
2
.
2nd step
A
2
( )
=
0
2
, v =
2
2
, B
2
=
0 1
1 0
,
P
2
=
1 0 0
0 0 1
0 1 0
.
P
2
P
1
A =
1
0
0
3
2
0
=
R
2
0
Q= P
2
P
1
( )
1
= P
1
*
P
2
*
=
0 0 1
1 0 0
0 1 0
.
Matrix Algebra 1 QR-decomposition 94
Check, that
A = QR.
If the diagonal elements of R are chosen to be positive then
R =
1
0
0
3
2
0
, Q=
0 0 1
1 0 0
0 1 0
.
Computation times. Constructing the QR decomposition of an mxn matrix
takes n
2
(m-n/3) flops.
Exercise 4.3.3 Let
A =
0 0
1 2
0 1
.
Construct the QR decomposition by using Householder matrices and give the
solution where the diagonal elements of R are positive.
Exercise 4.3.4 Let A be an nxn matrix and rank(A)=n. Show that the QR
decomposition is unique if the diagonal elements of R are chosen to be positive.
Exercise 4.3.5 a) Use QR decomposition to solve the least square problem. Let
mn and the columns of A be linearly independent. Find the vector x such that
the error
|| Ax-b ||
is as small as possible. How big is the error? Hint: Multiplication by a unitary
matrix does not change the norm of a vector.
b) Find the least squares solution, when b=[1, 1, 1]
T
and
A =
0 0
1 3
0 2
Hint: Example 4.3.2.
4.4 Orthonormal basis and decomposition of space
Matrix Algebra 1 QR-decomposition 95
QR decomposition can be used to transform a given linearly independent set of
vectors to an orthonormal set. The orthonormal set spans the same subspace as
the original vector set. An orthonormal basis for the orthogonal complement is
constructed at the same time.
Let {x
1
,x
2
,...,x
n
} be a linearly independent set of vectors in F
m
(m n) and
R = span x
1
, x
2
,..., x
n
{ }
.
Form the mxn matrix X = [x
1
,x
2
,...,x
n
]. It has rank rank(X) = n (why). Assume
that X has the following QR decomposition
n
X = x
1
,..., x
n
= Q
R
n
0
= q
1
, q
2
,..., q
m
r
11
r
12
r
1n
0 r
22
r
23
..
r
2n
.. .. ... ...
0 .. 0 r
nn
0 0 0
.. .. .. ..
0 0... 0 0
m
Direct computation shows, that
X = x
1
,..., x
n
= q
1
, q
2
,..., q
n
R
n
,
where R
n
is a square upper triangular matrix. It has nonzero diagonal elements,
so by theorem 3.4.1
span x
1
,..., x
n
{ }
= span q
1
, q
2
,..., q
n
{ }
.
The set {q
1
,q
2
,...,q
n
} is orthonormal because it consists of the column vectors of a
unitary matrix (exercise 1.20.3). We note also that
x
1
= r
11
q
1
,
so x
1
and q
1
are parallel.
Further we notice that the column vectors of Q
q
1
, q
2
,..., q
n
, q
n+1
, ..., q
m
{ }
,
form an orthonormal basis for F
m
(as any m linearly independent vectors form a
basis). In addition
Matrix Algebra 1 QR-decomposition 96
R = span q
1
, q
2
, ..., q
n
{ }
.
A method to find coordinates in an orthonormal basis can be derived from the
general formula for coordinates, which was presented in chapter 3. Because this
result is frequently used we give it here.
Theorem 4.4.1. Let {q
1
,q
2
,...,q
m
} be an orthonormal basis of F
m
. Every vector in
the space F
m
has the following representation as a linear combination of basis
vectors
x = < q
i
i=1
m
, x > q
i
.
Also
x
2
= < q
i
, x >
2
i=1
m
.
Proof. Every vector is a linear combination of basis vectors
x = c
i
i=1
m
q
i
= q
1
, q
2
,..., q
m
[c
1
,...., c
m
]
T
= Qc.
We can solve the coefficient vector c.
c =
c
1
c
m
= Q* x=
< q
1
, x >
< q
m
, x >
.
Furthermore, because Q is unitary,
x
2
= Qc
2
= c
2
= < q
i
, x >
2
i=1
m
Theorem 4.4.2. Let R F
m
be a subspace. Now
F
m
= R R
.
Proof. According to theorem 4.5.1 a basis {x
1
,...,x
n
} of R can be extended to an
orthonormal basis {q
1
,q
2
,...,q
n
,q
n+1
,...,q
m
} of F
m
so that R = span {q
1
,q
2
,...,q
n
}.
We show that
R
= span q
n+1
, q
n+2
,..., q
m
{ }
.
Matrix Algebra 1 QR-decomposition 97
We can easily see that span{q
n+1
,q
n+2
,...,q
m
} R
. Let x be a vector in R
,
because {q
1
,q
2
,...,q
n
,q
n+1
,...,q
m
} is a basis, x has the following representation:
x = < q
i
i=1
m
, x > q
i
= < q
i
i=1
n
, x > q
i
+ < q
i
, x > q
i
i=n+1
m
.
and since x R
span q
n+1
, q
n+2
,..., q
m
{ }
.
{q
1
,q
2
,...,,q
m
} is orthonormal basis of F
m
and so we know that
F
m
= R R
.
Exercise 4.4.1 Let S be a subspace. Show, that (S
= S.
Hint: It is easy to prove, that
S (S
.
Exercise 4.4.2 Let A be an mxn matrix. Show, that
F
n
= R (A*) N A
( )
, F
m
= R (A) N A*
( )
.
This result verifies the separation of the domain and codomain spaces into
orthogonal sums of subspaces as represented in picture 3.6.1.
Exercise 4.4.3 A well known method for transforming a given set of vectors to
an orthonormal set is the Gram-Schmidt algorithm. Deduce it from the QR
decomposition where the diagonal elements of R are positive by researching the
interaction between matrices X and Q. In numerical computation the Gram-
Schmidt algorithm is not as good as QR decomposition. The Gram-Schmidt
algorithm might not work properly in a situation where the basis vectors are
almost linearly dependent.