<< . .

. 30
( : 95)



. . >>

Program 26 - powerm : Power method
function [nu1,x1,niter,err]=powerm(A,z0,toll,nmax)
198 5. Approximation of Eigenvalues and Eigenvectors

q=z0/norm(z0); q2=q; err=[]; nu1=[]; res=toll+1; niter=0; z=A*q;
while (res >= toll & niter <= nmax)
q=z/norm(z); z=A*q; lam=q™*z; x1=q;
z2=q2™*A; q2=z2/norm(z2); q2=q2™;
y1=q2; costheta=abs(y1™*x1);
if (costheta >= 5e-2),
niter=niter+1; res=norm(z-lam*q)/costheta;
err=[err; res]; nu1=[nu1; lam];
else
disp(™ Multiple eigenvalue ™); break;
end
end

A coding of the inverse power method is provided in Program 27. The
input parameter mu is the initial approximation of the eigenvalue. In output,
the vectors sigma and err contain the sequences {σ (k) } and r(k) 2 /|(w(k) )H
q(k) | (see (5.29)). The LU factorization (with partial pivoting) of the ma-
trix Mµ is carried out using the MATLAB intrinsic function lu.

Program 27 - invpower : Inverse power method
function [sigma,x,niter,err]=invpower(A,z0,mu,toll,nmax)
n=max(size(A)); M=A-mu*eye(n); [L,U,P]=lu(M);
q=z0/norm(z0); q2=q™; err=[]; sigma=[]; res=toll+1; niter=0;
while (res >= toll & niter <= nmax)
niter=niter+1; b=P*q; y=L\b; z=U\y;
q=z/norm(z); z=A*q; lam=q™*z;
b=q2™; y=U™\b; w=L™\y;
q2=(P™*w)™; q2=q2/norm(q2); costheta=abs(q2*q);
if (costheta >= 5e-2),
res=norm(z-lam*q)/costheta; err=[err; res]; sigma=[sigma; lam];
else,
disp(™ Multiple eigenvalue ™); break;
end
x=q;
end



Example 5.3 The matrix A in (5.30)
®  ® 
’2 ’0.944 ’0.088
15 2 0.393
   
A= ’3  , V =  ’0.312 0.309 
1 10 0.919 (5.30)
° » ° »
’2 1 0 0.112 0.013 0.947

has the following eigenvalues (to ¬ve signi¬cant ¬gures): »1 = 14.103, »2 = 10.385
and »3 = 0.512, while the corresponding eigenvectors are the vector columns of
matrix V.
5.3 The Power Method 199

0
0
10
10

’2
10
’2
10

’4
(NS)
10
’4
10
’6
10
’6
10
’8
10
’8
10 (S)
’10
10

’10
10 ’12
10

’14
’12
10
10
0 6 12 18
0 10 20 30 40 50 60 70 80



FIGURE 5.2. Comparison between the a posteriori error estimate and the actual
absolute error for matrix A in (5.30) (left); convergence curves for the power
method applied to matrix A in (5.31) in its symmetric (S) and nonsymmetric
(NS) forms (right)

To approximate the pair (»1 , x1 ), we have run the Program 26 with initial datum
z(0) = (1, 1, 1)T . After 71 iterations of the power method the absolute errors are
|»1 ’ ν (71) | = 7.91 · 10’11 and x1 ’ x1 ’11
(71)
∞ = 1.42 · 10 .
(0)
In a second run, we have used z = x2 + x3 (notice that with this choice
±1 = 0). After 215 iterations the absolute errors are |»1 ’ ν (215) | = 4.26 · 10’14
’14
(215)
and x1 ’ x1 ∞ = 1.38 · 10 .
Figure 5.2 (left) shows the reliability of the a posteriori estimate (5.26). The
sequences |»1 ’ ν (k) | (solid line) and the corresponding a posteriori estimates
(5.26) (dashed line) are plotted as a function of the number of iterations (in
abscissae). Notice the excellent agreement between the two curves.

The symmetric matrix A in (5.31)
®  ® 
1 3 4 8 1 6
   
A= 3 2 , T= 3 7
1 5 (5.31)
° » ° »
4 2 1 4 9 2

has the following spectrum: »1 = 7.047, »2 = ’3.1879 and »3 = ’0.8868 (to ¬ve
signi¬cant ¬gures).
It is interesting to compare the behaviour of the power method when computing
»1 for the symmetric matrix A and for its similar matrix M = T’1 AT, where T
is the nonsingular (and nonorthogonal) matrix in (5.31).
Running Program 26 with z(0) = (1, 1, 1)T , the power method converges to the
eigenvalue »1 in 18 and 30 iterations, for matrices A and M, respectively. The
sequence of absolute errors |»1 ’ ν (k) | is plotted in Figure 5.2 (right) where (S)
and (NS) refer to the computations on A and M, respectively. Notice the rapid
error reduction in the symmetric case, according to the quadratic convergence
properties of the power method (see Section 5.3.1).
We ¬nally employ the inverse power method (5.28) to compute the eigenvalue
of smallest module »3 = 0.512 of matrix A in (5.30). Running Program 27 with

q(0) = (1, 1, 1)T / 3, the method converges in 9 iterations, with absolute errors
|»3 ’ σ (9) | = 1.194 · 10’12 and x3 ’ x3 ∞ = 4.59 · 10’13 .
(9)

200 5. Approximation of Eigenvalues and Eigenvectors

5.4 The QR Iteration
In this section we present some iterative techniques for simultaneously ap-
proximating all the eigenvalues of a given matrix A. The basic idea consists
of reducing A, by means of suitable similarity transformations, into a form
for which the calculation of the eigenvalues is easier than on the starting
matrix.
The problem would be satisfactorily solved if the unitary matrix U of the
Schur decomposition theorem 1.5, such that T = UH AU, T being upper
triangular and with tii = »i (A) for i = 1, . . . , n, could be determined in a
direct way, that is, with a ¬nite number of operations. Unfortunately, it is
a consequence of Abel™s theorem that, for n ≥ 5, the matrix U cannot be
computed in an elementary way (see Exercise 8). Thus, our problem can
be solved only resorting to iterative techniques.
The reference algorithm in this context is the QR iteration method, that is
here examined only in the case of real matrices. (For some remarks on the
extension of the algorithms to the complex case, see [GL89], Section 5.2.10
and [Dem97], Section 4.2.1).

Let A ∈ Rn—n ; given an orthogonal matrix Q(0) ∈ Rn—n and letting
T(0) = (Q(0) )T AQ(0) , for k = 1, 2, . . . , until convergence, the QR iteration
consists of:
determine Q(k) , R(k) such that
Q(k) R(k) = T(k’1) (QR factorization);
(5.32)
then, let
T(k) = R(k) Q(k) .

At each step k ≥ 1, the ¬rst phase of the iteration is the factorization of
the matrix T(k’1) into the product of an orthogonal matrix Q(k) with an
upper triangular matrix R(k) (see Section 5.6.3). The second phase is a
simple matrix product. Notice that

T(k) = R(k) Q(k) = (Q(k) )T (Q(k) R(k) )Q(k) = (Q(k) )T T(k’1) Q(k)
(5.33)
k ≥ 0,
= (Q(0) Q(1) . . . Q(k) )T A(Q(0) Q(1) . . . Q(k) ),

i.e., every matrix T(k) is orthogonally similar to A. This is particularly
relevant for the stability of the method, since, as shown in Section 5.2, the
conditioning of the matrix eigenvalue problem for T(k) is not worse than it
is for A (see also [GL89], p. 360).
A basic implementation of the QR iteration (5.32), assuming Q(0) = In ,
is examined in Section 5.5, while a more computationally e¬cient version,
starting from T(0) in upper Hessenberg form, is described in detail in Sec-
tion 5.6.
5.5 The Basic QR Iteration 201



If A has real eigenvalues, distinct in module, it will be seen in Section 5.5
that the limit of T(k) is an upper triangular matrix (with the eigenvalues of
A on the main diagonal). However, if A has complex eigenvalues the limit
of T(k) cannot be an upper triangular matrix T. Indeed if it were T would
necessarily have real eigenvalues, although it is similar to A.
Failure to converge to a triangular matrix may also happen in more
general situations, as addressed in Example 5.9.
For this, it is necessary to introduce variants of the QR iteration (5.32),
based on de¬‚ation and shift techniques (see Section 5.7 and, for a more
detailed discussion of the subject, [GL89], Chapter 7, [Dat95], Chapter 8
and [Dem97], Chapter 4).
These techniques allow for T(k) to converge to an upper quasi-triangular
matrix, known as the real Schur decomposition of A, for which the following
result holds (for the proof we refer to [GL89], pp. 341-342).

Property 5.8 Given a matrix A ∈ Rn—n , there exists an orthogonal ma-
trix Q ∈ Rn—n such that
® 
R11 R12 . . . R1m
 
0 R22 . . . R2m 
 
QT AQ =  . , (5.34)
. 
. .
..
. 
. .
..
.
° »
0 0 . . . Rmm

where each block Rii is either a real number or a matrix of order 2 having
complex conjugate eigenvalues, and

Q = lim Q(0) Q(1) · · · Q(k) (5.35)
k’∞

Q(k) being the orthogonal matrix generated by the k-th factorization step of
the QR iteration (5.32).


The QR iteration can be also employed to compute all the eigenvectors
of a given matrix. For this purpose, we describe in Section 5.8 two possi-
ble approaches, one based on the coupling between (5.32) and the inverse
iteration (5.28), the other working on the real Schur form (5.34).



5.5 The Basic QR Iteration
In the basic version of the QR method, one sets Q(0) = In in such a way
that T(0) = A. At each step k ≥ 1 the QR factorization of the matrix T(k’1)
202 5. Approximation of Eigenvalues and Eigenvectors

can be carried out using the modi¬ed Gram-Schmidt procedure introduced
in Section 3.4.3, with a cost of the order of 2n3 ¬‚ops (for a full matrix A).
The following convergence result holds (for the proof, see [GL89], Theorem
7.3.1, or [Wil65], pp. 517-519).

Property 5.9 (Convergence of QR method) Let A ∈ Rn—n be a ma-
trix such that
|»1 | > |»2 | > . . . > |»n |.
Then
® 
»1 t12 ... t1n
 
 ... 
0 »2 t23
 
= 
lim T(k) . . (5.36)
 . . ..
 .
. .
k’+∞
.
. . .»
°
0 0 ... »n

As for the convergence rate, we have
k
»i
(k)
|ti,i’1 | =O for k ’ +∞.
, i = 2, . . . , n, (5.37)
»i’1

Under the additional assumption that A is symmetric, the sequence {T(k) }
tends to a diagonal matrix.
If the eigenvalues of A, although being distinct, are not well-separated, it
follows from (5.37) that the convergence of T(k) towards a triangular matrix
can be quite slow. With the aim of accelerating it, one can resort to the
so-called shift technique, which will be addressed in Section 5.7.

Remark 5.2 It is always possible to reduce the matrix A into a triangular
form by means of an iterative algorithm employing nonorthogonal similarity
transformations. In such a case, the so-called LR iteration (known also as
Rutishauser method, [Rut58]) can be used, from which the QR method has
actually been derived (see also [Fra61], [Wil65]). The LR iteration is based
on the factorization of the matrix A into the product of two matrices L
and R, respectively unit lower triangular and upper triangular, and on the
(nonorthogonal) similarity transformation

L’1 AL = L’1 (LR)L = RL.

The rare use of the LR method in practical computations is due to the loss
of accuracy that can arise in the LR factorization because of the increase
in module of the upper diagonal entries of R. This aspect, together with
the details of the implementation of the algorithm and some comparisons
with the QR method, is examined in [Wil65], Chapter 8.
5.6 The QR Method for Matrices in Hessenberg Form 203

Example 5.4 We apply the QR method to the symmetric matrix A∈ R4—4 such
that aii = 4, for i = 1, . . . , 4, and aij = 4 + i ’ j for i < j ¤ 4, whose eigenvalues
are (to three signi¬cant ¬gures) »1 = 11.09, »2 = 3.41, »3 = 0.90 and »4 = 0.59.
After 20 iterations, we get
® 
6.44 · 10’10 ’3.62 · 10’15 9.49 · 10’15
11.09
 
 
’10 ’11 ’16
 6.47 · 10 1.43 · 10 4.60 · 10 
3.41
 .
T(20) = 
 1.74 · 10’21 
1.43 · 10’11 1.16 · 10’4
0.90
° »
2.32 · 10’25 2.68 · 10’15 1.16 · 10’4 0.58

Notice the “almost-diagonal” structure of the matrix T(20) and, at the same
time, the e¬ect of rounding errors which slightly alter its expected symmetry.
Good agreement can also be found between the under-diagonal entries and the

estimate (5.37).

A computer implementation of the basic QR iteration is given in Program
28. The QR factorization is executed using the modi¬ed Gram-Schmidt
method (Program 8). The input parameter niter denotes the maximum
admissible number of iterations, while the output parameters T, Q and R
are the matrices T, Q and R in (5.32) after niter iterations of the QR
procedure.

Program 28 - basicqr : Basic QR iteration

function [T,Q,R]=basicqr(A,niter)
T=A;
for i=1:niter,
[Q,R]=mod grams(T);
T=R*Q;
end




5.6 The QR Method for Matrices in Hessenberg
Form
The naive implementation of the QR method discussed in the previous
section requires (for a full matrix) a computational e¬ort of the order of
n3 ¬‚ops per iteration. In this section we illustrate a variant for the QR
iteration, known as Hessenberg-QR iteration, with a greatly reduced com-
putational cost. The idea consists of starting the iteration from a matrix
(0)
T(0) in upper Hessenberg form, that is, tij = 0 for i > j + 1. Indeed, it can
be checked that with this choice the computation of T(k) in (5.32) requires
only an order of n2 ¬‚ops per iteration.
204 5. Approximation of Eigenvalues and Eigenvectors

To achieve maximum e¬ciency and stability of the algorithm, suitable
transformation matrices are employed. Precisely, the preliminary reduc-
tion of matrix A into upper Hessenberg form is realized with Householder
matrices, whilst the QR factorization of T(k) is carried out using Givens
matrices, instead of the modi¬ed Gram-Schmidt procedure introduced in
Section 3.4.3.
We brie¬‚y describe Householder and Givens matrices in the next section,
referring to Section 5.6.5 for their implementation. The algorithm and ex-
amples of computations of the real Schur form of A starting from its upper
Hessenberg form are then discussed in Section 5.6.4.


5.6.1 Householder and Givens Transformation Matrices
For any vector v ∈ Rn , let us introduce the orthogonal and symmetric
matrix

P = I ’ 2vvT / v 2 . (5.38)
2


Given a vector x ∈ Rn , the vector y = Px is the re¬‚ection of x with respect
to the hyperplane π = span{v}⊥ formed by the set of the vectors that are
orthogonal to v (see Figure 5.3, left). Matrix P and the vector v are called
the Householder re¬‚ection matrix and the Householder vector, respectively.


xk
v
y
x θ
x
y
π
xi

FIGURE 5.3. Re¬‚ection across the hyperplane orthogonal to v (left); rotation by
an angle θ in the plane (xi , xk ) (right)


Householder matrices can be used to set to zero a block of components of
a given vector x ∈ Rn . If, in particular, one would like to set to zero all the
components of x, except the m-th one, the Householder vector ought to be
chosen as

<< . .

. 30
( : 95)



. . >>