<< . .

. 33
( : 95)

. . >>

’0.1814 ’3.2330
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
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
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;
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 
° 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

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
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);
if (k > 2),
if abs(T(k-1,k-2)) ¡= toll*Tdiag2;
end, T(k,k-1)=0;
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;
5.8 Computing the Eigenvectors and the SVD of a Matrix 221

5.8 Computing the Eigenvectors and the SVD of a
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

5.8.2 Computing the Eigenvectors from the Schur Form of a
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 ∈

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=¬ ·.
 

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
U T AV = , (5.57)

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) =
U1 A has ai1 = 0 if i > 1. Then, V1 is determined so that A(2) = A(1) V1
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
such that A(3) = U2 A(2) has ai2 = 0 for i > 2 and V2 in such a way that
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 = ,

with U = Udiag(W, Im’n ) and V = VZ.
The computational cost of this procedure is 2m2 n + 4mn2 + 9 n3 ¬‚ops,
which reduces to 2mn2 ’ 2 n3 ¬‚ops if only the singular values are computed.
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

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
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
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-
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

<< . .

. 33
( : 95)

. . >>