Download Numerical Linear Algebra, Lecture Notes - Mathematics - 3 and more Study notes Mathematical Methods for Numerical Analysis and Optimization in PDF only on Docsity!
8 The QR algorithm
We are now going to use the QR factorisation to compute the eigenvalues and eigenvectors
of a square matrix A ∈ R
n×n
. Though the method described here works for arbitrary
square matrices, we will concentrate on symmetric ones.
Recall that a symmetric matrix A is called similar to a matrix B if there is a nonsingular
matrix P such that A = P
− 1 BP. Similar matrices have the same eigenvalues, since if
A = P
− 1
BP ,
0 = det(A − λI) = det(P
− 1
(B − λI)P ) = det(P
− 1
) det(P ) det(B − λI),
and hence the characteristic polynomials associated to A and B have the same zeros,
meaning that A and B have the same eigenvalues.
Definition 8.1 The basic QR algorithm works as follows. First set A 1
= A. Then, for
k = 1, 2 , 3 , · · · do the following two steps:
- form the QR factorisation Ak = QkRk
- set A k+
= R
k
Q
k
Proposition 8.2 The symmetric matrices A 1 , A 2 ,... , Ak,... are all similar and thus have
the same eigenvalues.
Proof: Since
A
k+
= R
k
Q
k
= (Q
T
k
Q
k
)R
k
Q
k
= Q
T
k
(Q
k
R
k
)Q
k
= Q
T
k
A
k
Q
k
= Q
− 1
k
A
k
Q
k
we see that Ak+1 is symmetric if Ak is, and that Ak+1 is similar to Ak.!
This basic QR algorithm works since it is possible to show that for a symmetric matrix A,
the sequence of matrices A k
eventually converges to a diagonal matrix for k → ∞. The
entries on the diagonal therefore have to be the eigenvalues of the matrix A. Moreover, if
we collect the products of matrices Q k
we have
Ak+1 = Q
T
k
AkQk = Q
T
k
Q
T
k− 1
Ak− 1 Qk− 1 Qk = · · · = Q
T
k
Qk− 1 · · · Q
T
1
AQ 1 · · · Qk− 1 Qk
=: Q
T
1 ..k
AQ 1 ..k,
such that the columns of Q 1 ..k
are an approximation to the eigenvectors of A.
However, for being computationally efficient, we need two improvements.
The first one follows from the fact that in each step, we have to compute a QR factorisation,
which will cost us O(n
3 ) time. Then, we have to form the new iteration Ak+1 which is
done by multiplying two full matrices, which also costs O(n
3 ) time. Both is too expensive
if a large number of steps is required.
A possible idea to reduce this cost is to first transform A to a sparse form, this may cost
O(n
3 ) time but is required only once.
The idea to reduce this cost is to apply explicit similarity transformations QAQ
− 1
QAQ
T , with Q orthogonal, so that QAQ
T is tridiagonal. We will use Householder trans-
formations to achieve this. It would, of course, be nice if the resulting matrix QAQ
T would
be a diagonal matrix but the following observation shows that, in general, we cannot expect
this.
Suppose we choose a Householder matrix H(w) such that
H(w)A =
α 1
. A 1
then H(w)AH(w)
T is generally full, i.e. all zeros created by pre-multiplication are de-
stroyed by the post-multiplication. This simply comes from the fact that multiplication
from the right means that the first column is replaced by linear combinations of all columns.
If, however, we write
A =
a 11
a ˜ 1
A
with a 11
∈ R, a˜ 1
∈ R
n− 1 and
A ∈ R
n− 1 ×n− 1 , choose a Householder matrix
H
1
∈ R
n− 1 ×n− 1
with
H 1 ˜a 1 = α 1 e 1 , and define
H 1 =
T
H
1
then we have H
T
1
= H
1
and
H
1
AH
1
T
H
1
a 11
a ˜ 1
A
T
H
1
a 11 ∗ · · · ∗
α 1 e 1
H 1
A
T
H 1
a 11
α 1
e 1
H
1
A
H
1
0 A
2
(a 11 , a 21 )
T into a multiple of (1, 0). This matrix
H(1, 2) can be “filled-up” to an n × n
orthogonal matrix H(1, 2) in the usual way. Important here is that the computation of
H(1, 2)A costs constant time since it involves only changing the first two rows and, because
of the tridiagonal structure, only the first three columns. Continuing like this, we can use
the index pairs (2, 3), (3, 4), · · · , (n − 1 , n) to create the QR factorisation
R = H(n − 1 , n)H(n − 2 , n) · · · H(1, 2)A,
such that(because of the symmetry of the Householder matrices)
Q = H(1, 2)H(2, 3) · · · H(n − 1 , n).
In each step, at most six entries are changed. Thus the computational cost is O(n).
This can be visualised as follows:
(1,2)
(2,3)
(3,4)
(4,5)
Next, to compute the next iteration A k+
= R
k
Q
k
= R
k
H(1, 2)H(2, 3) · · · H(n − 1 , n) we
see that, the first matrix H(1, 2) acts only on the first two columns, the next matrix H(2, 3)
on columns 2 and 3 and so on:
(1,2)
(2,3)
(3,4)
(4,5)
The new matrix is indeed a Hessenberg matrix again. But since it is also symmetric, all
entries above the super-diagonal have to be zero. This ensures that the new matrix is a
tridiagonal matrix and that the computation can again be done in linear time.!
A further improvement can be achieved by introducing shifts into the algorithm. The
iteration then takes the form that for k = 1, 2 ,...
- form the QR factorization of Ak − μkI = QkRk
- set A k+
= R
k
Q
k
I
It is easy to see that for any chosen sequence of values of μ k
∈ R the matrices A k
are
symmetric and tridiagonal if A 1 has these properties, and they are all similar to A 1.
The simplest shift to use is a n,n
, which leads rapidly in almost all cases to
Ak =
Tk 0
T λ
where T k
is n − 1 by n − 1 and tridiagonal, and λ is an eigenvalue of A 1
. Inductively, once
this form has been found, the QR algorithm with shift an− 1 ,n− 1 can be concentrated only
on the n − 1 by n − 1 leading submatrix T k
. This process is called deflation.
The overall algorithm for calculating the eigenvalues of an n by n symmetric matrix is
hence:
Algorithm 1: QR algorithm with shifts
Input : Symmetric matrix A ∈ R
n×n
Output : Approximate eigenvalues and eigenvectors
reduce A to tridiagonal form
for m = n, n − 1 ,... 2 do
while a m− 1 ,m
> tol do
[Q, R] = qr(A − a m,m
∗ I)
A = R ∗ Q + a m,m
∗ I
record eigenvalue λ m
= a m,m
A ← leading m − 1 by m − 1 submatrix of A
record eigenvalue λ 1 = a 1 , 1
