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