1.1146 1.2648

has eigenvalues {±i, 1, 2}, i being the imaginary unit. Running Program 36 with

toll equal to the roundo¬ unit yields after 100 iterations

®

2 1.1999 0.5148 4.9004

0 0.7182

’0.0001 ’0.8575

T(101) = .

0 ’0.8186

1.1662 0.0001

° »

0 0 0 1

The obtained matrix is the real Schur form of A, where the 2—2 block T(101) (2:3,

2:3) has complex conjugate eigenvalues ±i. These eigenvalues cannot be computed

•

by the algorithm (5.52)-(5.53) since µ is real.

The problem with this example is that working with real matrices neces-

sarily yields a real shift, whereas a complex one would be needed. The QR

iteration with double shift is set up to account for complex eigenvalues and

allows for removing the 2—2 diagonal blocks of the real Schur form of A.

Precisely, suppose that the QR iteration with single shift (5.52) detects

(k)

at some step k a 2—2 diagonal block Rkk that cannot be reduced into

upper triangular form. Since the iteration is converging to the real Schur

(k)

form of the matrix A the two eigenvalues of Rkk are complex conjugate

¯

and will be denoted by »(k) and »(k) . The double shift strategy consists of

5.7 The QR Iteration with Shifting Techniques 219

the following steps:

determine Q(k) , R(k) such that

Q(k) R(k) = T(k’1) ’ »(k) I (¬rst QR factorization);

then, let

T(k) = R(k) Q(k) + »(k) I;

(5.55)

determine Q(k+1) , R(k+1) such that

¯

Q(k+1) R(k+1) = T(k) ’ »(k) I (second QR factorization);

then, let

¯

T(k+1) = R(k+1) Q(k+1) + »(k) I.

Once the double shift has been carried out the QR iteration with single shift

is continued until a situation analogous to the one above is encountered.

The QR iteration incorporating the double shift strategy is the most e¬ec-

tive algorithm for computing eigenvalues and yields the approximate Schur

form of a given matrix A. Its actual implementation is far more sophisti-

cated than the outline above and is called QR iteration with Francis shift

(see [Fra61], and, also, [GL89], Section 7.5 and [Dem97], Section 4.4.5). As

for the case of the QR iteration with single shift, quadratic convergence can

also be proven for the QR method with Francis shift. However, special ma-

trices have recently been found for which the method fails to converge (see

for an example Exercise 14 and Remark 5.13). We refer for some analysis

and remedies to [Bat90], [Day96], although the ¬nding of a shift strategy

that guarantees convergence of the QR iteration for all matrices is still an

open problem.

Example 5.12 Let us apply the QR iteration with double shift to the matrix

A in Example 5.11. After 97 iterations of Program 37, with toll equal to the

roundo¬ unit, we get the following (approximate) Schur form of A, which displays

on its diagonal the four eigenvalues of A

®

’2.33 + 0.86i

2 1 + 2i 4.90

0 5.02 · 10’14 + i ’2.02 + 6.91 · 10’14 i

0.72

T(97) = (97)

t ’1.78 · 10’14 ’ i ’0.82

0

° 31 »

(97) (97)

t41 t42 0 1

where t31 = 2.06 · 10’17 + 7.15 · 10’49 i, t41 = ’5.59 · 10’17 and t42

(97) (97) (97)

=

’4.26 · 10’18 , respectively. •

Example 5.13 Consider the pseudo-spectral di¬erentiation matrix (10.73) of

order 5. This matrix is singular, with a unique eigenvalue » = 0 of algebraic

220 5. Approximation of Eigenvalues and Eigenvectors

multiplicity equal to 5 (see [CHQZ88], p. 44). In this case the QR method with

double shift provides an inaccurate approximation of the spectrum of the ma-

trix. Indeed, using Program 37, with toll=eps, the method converges after 59

iterations to an upper triangular matrix with diagonal entries given by 0.0020,

0.0006 ± 0.0019i and ’0.0017 ± 0.0012i, respectively. Using the MATLAB intrin-

sic function eig yields instead the eigenvalues ’0.0024, ’0.0007 ± 0.0023i and

0.0019 ± 0.0014i. •

A basic implementation of the QR iteration with double shift is provided

in Program 37. The input/output parameters are the same as those of

Program 36. The output matrix T is the approximate Schur form of matrix

A.

Program 37 - qr2shift : QR iteration with double shift

function [T,iter]=qr2shift(A,toll,itmax)

n=max(size(A)); iter=0; [T,Q]=houshess(A);

for k=n:-1:2

I=eye(k);

while abs(T(k,k-1)) > toll*(abs(T(k,k))+abs(T(k-1,k-1)))

iter=iter+1; if (iter > itmax), return, end

mu=T(k,k); [Q,R,c,s]=qrgivens(T(1:k,1:k)-mu*I);

T(1:k,1:k)=R*Q+mu*I;

if (k > 2),

Tdiag2=abs(T(k-1,k-1))+abs(T(k-2,k-2));

if abs(T(k-1,k-2)) ¡= toll*Tdiag2;

[lambda]=eig(T(k-1:k,k-1:k));

[Q,R,c,s]=qrgivens(T(1:k,1:k)-lambda(1)*I);

T(1:k,1:k)=R*Q+lambda(1)*I;

[Q,R,c,s]=qrgivens(T(1:k,1:k)-lambda(2)*I);

T(1:k,1:k)=R*Q+lambda(2)*I;

end

end

end, T(k,k-1)=0;

end

I=eye(2);

while (abs(T(2,1)) > toll*(abs(T(2,2))+abs(T(1,1)))) & (iter <= itmax)

iter=iter+1; mu=T(2,2);

[Q,R,c,s]=qrgivens(T(1:2,1:2)-mu*I); T(1:2,1:2)=R*Q+mu*I;

end

5.8 Computing the Eigenvectors and the SVD of a Matrix 221

5.8 Computing the Eigenvectors and the SVD of a

Matrix

The power and inverse iterations described in Section 5.3.2 can be used

to compute a selected number of eigenvalue/eigenvector pairs. If all the

eigenvalues and eigenvectors of a matrix are needed, the QR iteration can be

pro¬tably employed to compute the eigenvectors as shown in Sections 5.8.1

and 5.8.2. In Section 5.8.3 we deal with the computation of the singular

value decomposition (SVD) of a given matrix.

5.8.1 The Hessenberg Inverse Iteration

For any approximate eigenvalue » computed by the QR iteration as de-

scribed in Section 5.7.2, the inverse iteration (5.28) can be applied to the

matrix H = QT AQ in Hessenberg form, yielding an approximate eigenvec-

tor q. Then, the eigenvector x associated with » is computed as x = Qq.

Clearly, one can take advantage of the structure of the Hessenberg matrix

for an e¬cient solution of the linear system at each step of (5.28). Typi-

cally, only one iteration is required to produce an adequate approximation

of the desired eigenvector x (see [GL89], Section 7.6.1 and [PW79] for more

details).

5.8.2 Computing the Eigenvectors from the Schur Form of a

Matrix

Suppose that the (approximate) Schur form QH AQ=T of a given matrix

A∈ Rn—n has been computed by the QR iteration with double shift, Q

being a unitary matrix and T being upper triangular.

Then, if Ax=»x, we have QH AQQH x= QH »x, i.e., letting y=QH x, T

y=»y holds. Therefore y is an eigenvector of T, so that to compute the

eigenvectors of A we can work directly on the Schur form T.

Assume for simplicity that » = tkk ∈ C is a simple eigenvalue of A. Then

the upper triangular matrix T can be decomposed as

®

T11 v T13

T= wT ,

0 »

° »

0 0 T33

where T11 ∈ C(k’1)—(k’1) and T33 ∈ C(n’k)—(n’k) are upper triangular

matrices, v∈ Ck’1 , w∈ Cn’k and » ∈ σ(T11 ) ∪ σ(T33 ).

222 5. Approximation of Eigenvalues and Eigenvectors

Thus, letting y = yk’1 , y, yn’k , with yk’1 ∈ Ck’1 , y ∈ C and yn’k ∈

T T

Cn’k , the matrix eigenvector problem (T - »I) y=0 can be written as

±

(T11 ’ »Ik’1 )yk’1 + vy+ T13 yn’k =0

wT yn’k =0 (5.56)

(T33 ’ »In’k )yn’k = 0.

Since » is simple, both matrices T11 ’ »Ik’1 and T33 ’ »In’k are nonsin-

gular, so that the third equation in (5.56) yields yn’k = 0 and the ¬rst

equation becomes

(T11 ’ »Ik’1 )yk’1 = ’vy.

Setting arbitrarily y = 1 and solving the triangular system above for yk’1

yields (formally)

«

’(T11 ’ »Ik’1 )’1 v

¬ ·

y=¬ ·.

1

0

The desired eigenvector x can then be computed as x=Qy.

An e¬cient implementation of the above procedure is carried out in the

intrinsic MATLAB function eig. Invoking this function with the format [V,

D]= eig(A) yields the matrix V whose columns are the right eigenvectors

of A and the diagonal matrix D contains its eigenvalues. Further details can

be found in the strvec subroutine in the LAPACK library, while for the

computation of eigenvectors in the case where A is symmetric, we refer to

[GL89], Chapter 8 and [Dem97], Section 5.3.

5.8.3 Approximate Computation of the SVD of a Matrix

In this section we describe the Golub-Kahan-Reinsch algorithm for the

computation of the SVD of a matrix A ∈ Rm—n with m ≥ n (see [GL89],

Section 5.4). The method consists of two phases, a direct one and an iter-

ative one.

In the ¬rst phase A is transformed into an upper trapezoidal matrix of

the form

B

U T AV = , (5.57)

0

where U and V are two orthogonal matrices and B ∈ Rn—n is upper bidi-

agonal. The matrices U and V are generated using n + m ’ 3 Householder

matrices U1 , . . . , Um’1 , V1 , . . . , Vn’2 as follows.

5.8 Computing the Eigenvectors and the SVD of a Matrix 223

The algorithm initially generates U1 in such a way that the matrix A(1) =

(1)

U1 A has ai1 = 0 if i > 1. Then, V1 is determined so that A(2) = A(1) V1

(2)

has a1j = 0 for j > 2, preserving at the same time the null entries of the

previous step. The procedure is repeated starting from A(2) , and taking U2

(3)

such that A(3) = U2 A(2) has ai2 = 0 for i > 2 and V2 in such a way that

(4)

A(4) = A(3) V2 has a2j = 0 for j > 3, yet preserving the null entries already

generated. For example, in the case m = 5, n = 4 the ¬rst two steps of the

reduction process yield

® ®

•••• ••00

0 • • • 0 • • •

A = U1 A = 0 • • • ’’ A = A V1 = 0 • • • ,

(1) (2) (1)

°0 • • •» °0 • • •»

• • • • • •

0 0

having denoted by • the entries of the matrices that in principle are di¬erent

than zero. After at most m ’ 1 steps, we ¬nd (5.57) with

U = U1 U2 . . . Um’1 , V = V1 V2 . . . Vn’2 .

In the second phase, the obtained matrix B is reduced into a diagonal

matrix Σ using the QR iteration. Precisely, a sequence of upper bidiagonal

matrices B(k) are constructed such that, as k ’ ∞, their o¬-diagonal

entries tend to zero quadratically and the diagonal entries tend to the

singular values σi of A. In the limit, the process generates two orthogonal

matrices W and Z such that

W T BZ = Σ = diag(σ1 , . . . , σn ).

The SVD of A is then given by

Σ

UT AV = ,

0

with U = Udiag(W, Im’n ) and V = VZ.

The computational cost of this procedure is 2m2 n + 4mn2 + 9 n3 ¬‚ops,

2

which reduces to 2mn2 ’ 2 n3 ¬‚ops if only the singular values are computed.

3

In this case, recalling what was stated in Section 3.13 about AT A, the

method described in the present section is preferable to computing directly

the eigenvalues of AT A and then taking their square roots.

As for the stability of this procedure, it can be shown that the computed

σi turn out to be the singular values of the matrix A + δA with

¤ Cmn u A 2 ,

δA 2

224 5. Approximation of Eigenvalues and Eigenvectors

Cmn being a constant dependent on n, m and the roundo¬ unit u. For other

approaches to the computation of the SVD of a matrix, see [Dat95] and

[GL89].

5.9 The Generalized Eigenvalue Problem

Let A, B ∈ Cn—n be two given matrices; for any z ∈ C, we call A ’ zB a

matrix pencil and denote it by (A,B). The set σ(A,B) of the eigenvalues of

(A,B) is de¬ned as

σ(A, B) = {µ ∈ C : det(A ’ µB) = 0} .

The generalized matrix eigenvalue problem can be formulated as: ¬nd » ∈

σ(A,B) and a nonnull vector x ∈ Cn such that

Ax = »Bx. (5.58)

The pair (», x) satisfying (5.58) is an eigenvalue/eigenvector pair of the

pencil (A,B). Note that by setting B=In in (5.58) we recover the standard

matrix eigenvalue problem considered thus far.

Problems like (5.58) arise frequently in engineering applications, e.g.,

in the study of vibrations of structures (buildings, aircrafts and bridges)

or in the mode analysis for waveguides (see [Inm94] and [Bos93]). Another

example is the computation of the extremal eigenvalues of a preconditioned

matrix P’1 A (in which case B = P in (5.58)) when solving a linear system

with an iterative method (see Remark 4.2).

Let us introduce some de¬nitions. We say that the pencil (A,B) is regular

if det(A-zB) is not identically zero, otherwise the pencil is singular. When

(A,B) is regular, p(z) = det(A ’ zB) is the characteristic polynomial of the

pencil; denoting by k the degree of p, the eigenvalues of (A,B) are de¬ned

as:

1. the roots of p(z) = 0, if k = n;

2. ∞ if k < n (with multiplicity equal to n ’ k).

Example 5.14 (Taken from [Par80], [Saa92] and [GL89])

’1 0 0 1

σ(A, B) = ±i

p(z) = z 2 + 1

A= , B= =’

0 1 1 0

’1 0 0 0

σ(A, B) = {0, ∞}

A= , B= p(z) = z =’

0 0 0 1

1 2 1 0

σ(A, B) = C.

A= , B= p(z) = 0 =’

0 0 0 0

The ¬rst pair of matrices shows that symmetric pencils, unlike symmetric ma-

trices, may exhibit complex conjugate eigenvalues. The second pair is a regular

pencil displaying an eigenvalue equal to in¬nity, while the third pair is an example

•

of singular pencil.

5.9 The Generalized Eigenvalue Problem 225

5.9.1 Computing the Generalized Real Schur Form

The de¬nitions and examples above imply that the pencil (A,B) has n ¬nite

eigenvalues i¬ B is nonsingular.

In such a case, a possible approach to the solution of problem (5.58) is

to transform it into the equivalent eigenvalue problem Cx = »x, where the

matrix C is the solution of the system BC = A, then apply the QR iteration

to C. For actually computing the matrix C, one can use Gauss elimination

with pivoting or the techniques shown in Section 3.6. This procedure can

yield inaccurate results if B is ill-conditioned, since computing C is a¬ected

by rounding errors of the order of u A 2 B’1 2 (see [GL89], p. 376).

A more attractive approach is based on the following result, which gen-

eralizes the Schur decomposition theorem 1.5 to the case of regular pencils

to (for a proof, see [Dat95], p. 497).

Property 5.10 (Generalized Schur decomposition) Let (A,B) be a

regular pencil. Then, there exist two unitary matrices U and Z such that

UH AZ = T, UH BZ = S, where T and S are upper triangular. For i =

1, . . . , n the eigenvalues of (A,B) are given by

»i = tii /sii , if sii = 0,

»i = ∞, if tii = 0, sii = 0.

Exactly as in the matrix eigenvalue problem, the generalized Schur form

cannot be explicitly computed, so the counterpart of the real Schur form

(5.34) has to be computed. Assuming that the matrices A and B are real, it

˜ ˜

can be shown that there exist two orthogonal matrices U and Z such that

˜ ˜ ˜ ˜ ˜ ˜

T = UT AZ is upper quasi-triangular and S = UT BZ is upper triangular.

This decomposition is known as the generalized real Schur decomposition

of a pair (A,B) and can be computed by a suitably modi¬ed version of the

QR algorithm, known as QZ iteration, which consists of the following steps

(for a more detailed description, see [GL89], Section 7.7, [Dat95], Section

9.3):

1. reduce A and B into upper Hessenberg form and upper triangular

form, respectively, i.e., ¬nd two orthogonal matrices Q and Z such

that A = QT AZ is upper Hessenberg and B = QT BZ is upper trian-

gular;

2. the QR iteration is applied to the matrix AB ’1 to reduce it to real

Schur form.

To save computational resources, the QZ algorithm overwrites the matrices

A and B on their upper Hessenberg and triangular forms and requires 30n3

¬‚ops; an additional cost of 36n3 operations is required if Q and Z are

226 5. Approximation of Eigenvalues and Eigenvectors

also needed. The method is implemented in the LAPACK library in the

subroutine sgges and can be invoked in the MATLAB environment with