A projection method for generalized
eigenvalue problems
Tetsuya Sakurai†

Hiroshi Sugiura‡

July 24, 2002
ISE-TR-02-189†Institute of Information Sciences and Electronics
University of Tsukuba
Tsukuba, Ibaraki 305-8573, Japan‡

Graduate School of Engineering
Nagoya University
Furocho, Chikusaku, Nagoya 464-5501, JapanA projection method for generalized
eigenvalue problems
Tetsuya Sakuraia and Hiroshi Sugiuraba Institute of Information Sciences and Electronics,

University of Tsukuba, Tsukuba 305-8573 (Japan)
sakurai@is.tsukuba.ac.jp
b Graduate School of Engineering,

Nagoya University, Nagoya 464-8603 (Japan)sugiura@nuie.nagoya-u.ac.jp

Abstract.
In this paper, we propose a method for finding certain eigenvalues of ageneralized eigenvalue problem that lie in a given domain of the complex plane.
The proposed method projects the matrix pencil onto a subspace associatedwith the eigenvalues that are located in the domain via numerical integration.
The projection produces a small pencil with Hankel matrices.

1Introduction

The generalized eigenvalue problem
Ax = λBx,
where A and B are n × n real or complex matrices, arises in many applicationsof scientific computations. Often A and B are large and sparse, and only a few
of the eigenvalues are desired. Due to the size of the problem, several methods
for such eigenvalue problems are building sequences of subspaces that containthe desired eigenvectors, for example, see [14, 15, 16, 18, 19]. A comprehensive
report on computational methods for eigenvalue problems was given by Goluband van der Vorst [6].
In this paper, we propose a method for determining certain eigenvalues thatlie in a given domain of the complex plane. The proposed method projects
the matrix pencil A − λB onto a subspace associated with the eigenvalues
that are located in the domain. The proposed approach is based on the root
finding method for an analytic function proposed by Kravanja, Sakurai and
Van Barel [11]. This method finds all of the zeros that lie in a circle by using
numerical integration. Methods for the determination of zeros of analytic
functions that are based on the numerical evaluation of integrals are called
quadrature methods [5, 10, 13].
1

In Section 2, we introduce a function that is derived via the resolvent of
the matrix pencil. We show a projection method using numerical integrationof this function. In Section 3, we present an algorithm for the case in which
the domain is given by a circle and the numerical integration is evaluated via
the trapezoidal rule. We present an error analysis of the method using theresults described in [12, 17]. In Section 4, a number of numerical examples areused to illustrate the properties of the proposed method.

2A projection method

Let A, B ∈ Cn×n , and let λ1 , . . . , λd be finite eigenvalues of the matrix pencil
A − λB. The pencil A − λB is called regular if det(A − λB) is not identicallyzero for λ ∈ C.
For nonzero vectors u, v ∈ Cn , we define
f (z) := uH (zB − A)−1 v.
The function f (z) is analytic when zB − A is nonsingular.
We use the following result about Weierstrass’s canonical form (see e.g.,
Gantmacher [7]).
Theorem 1 Let A − zB be a regular pencil of order n. Then there exist
nonsingular matrices P, Q ∈ Cn×n such thatµ¶zId − JdO
P (zB − A)Q =,(1)OzJn−d − In−dwhere Jd and Jn−d are in Jordan canonical form, Jn−d is nilpotent, and Id
denotes the identity matrix of order d.
We assume that Jd is the diagonal matrix Jd := diag(λ1 , . . . , λd ). Let
p1 , . . . , pn ∈ Cn be vectors such thatP T = (p1 , . . . , pn ),and let q1 , . . . , qn ∈ Cn be
Q = (q1 , . . . , qn ).
Letνj := uH qj pHj v,1 ≤ j ≤ d.

Theorem 2 Let K be the maximum size of Jordan blocks of Jn−d . If A − λB
is regular and A is diagonalizable thenf (z) =dX

νj+ g(z),z − λjj=1where g(z) is a polynomial of degree K − 1.
2(2)Proof. Let P be partitioned in the formµ¶
P1,P =P2
where P1 ∈ Cd×n and P2 ∈ C(n−d)×n , and let
Q = (Q1 , Q2 ),where Q1 ∈ Cn×d and Q2 ∈ Cn×(n−d) .
According to Theorem 1, we havef (z) = uH (zB − A)−1 vµ¶−1zId − JdOHPv= u QOzJn−d − In−d= uH Q1 (zId − Jd )−1 P1 v + uH Q2 (In−d − zJn−d )−1 P2 v.
Since

(zId − Jd )−1 = (z − λ1 )−1(z − λ2 )−1..


,
.
(z − λd )−1it follows thatuH Q1 (zId − Jd )−1 P1 v




p1 v


= (uH q1 , . . . , uH qd )(zId − Jd )−1  ... pd v=dXuH qj pHj v

z − λjj=1=LetdXνj.
z − λjj=1
D=


0 10

 ∈ Rk×k .
...
1 0..

.

Then(I − zD)−1 = (I + zD + (zD)2 + · · · + (zD)k−1 ).
3(3)Therefore if the maximum size of Jordan blocks of (In−d − zJn−d ) is K, theng(z) := uH Q2 (In−d − zJn−d )−1 P2 v
is a polynomial of degree K − 1.

¤The resolvent form of a pencil has important roles in eigenvalue problems
and matrix analysis (see e.g., [2] and [8]).
Let Γ be a positively oriented closed Jordan curve in the complex plane.
Let λ1 , . . . , λm be distinct eigenvalues that lie in the interior of Γ. We consider
the problem to find all the poles of f (z) inside Γ. LetZ
1µk :=z k f (z)dz, k = 0, 1, . . . .
(4)2πi Γ<
Let the m × m Hankel matrices Hm and Hm
beHm := [µi+j−2 ]mi,j=1<and Hm:= [µi+j−1 ]mi,j=1 .

We have the following theorem:
Theorem 3 If νj 6= 0 for 1 ≤ j ≤ m then the eigenvalues of the pencil<Hm− λHm are given by λ1 , . . . , λm .
Proof. Let Vm be the Vandermonde matrix

11···
 λ1
λ2 · · ·
Vm :=  .
..
.
 .
.

1λm..
.



.


λm−1λm−1· · · λm−1
12m
Let Dm := diag(ν1 , . . . , νm ) and Λm := diag(λ1 , . . . , λm ). The residue theorem
implies thatmXµk =νj λkj , k = 0, 1, . . . .
(5)
j=1

Then, we can verify that
Hm = Vm Dm VmTand<= Vm Dm Λm VmT .
HmTherefore,<− λHm = Vm Dm ΛVmT − λVm Dm VmTHm= Vm Dm (Λ − λI)VmT .
4

Since λ1 , . . . , λm are distinct, Vm is nonsingular, and since Dm is also nonsingular, we have the result of the theorem.
¤Therefore, we can obtain the eigenvalues λ1 , . . . , λm by solving the gener<
alized eigenvalue problem Hmx = λHm x. If we find an appropriate Γ thatincludes a small number of eigenvalues, then the size of the derived eigen<problem is small. Note that the elements of Hm and Hm
are defined via theintegration (4).

3The case in which Γ is given by a circle

In this section, we consider the case in which Γ is given by a circle and theintegration is evaluated via a trapezoidal rule on the circle. Let γ and ρ be
the center and the radius, respectively, of the given circle. Let N be a positive
integer, and let
2π iωj := e N j , j = 0, 1, . . . , N − 1.
If Γ is the unit circle, we obtain the following approximations for µk
N −11 X k+1
ω f (ωj ),µk ≈ µ̂k :=N j=0 jk = 0, 1, . . . ,(6)by approximating the integral of equation (4) via the N -point trapezoidal rule.
Note thatf (ωj ) = v H (ωj B − A)−1 u,and the equation (6) can be evaluated via the first Fourier transform.
When Γ is not the unit circle, we apply the shift and the scaling for f (z),i.e., we use f (γ + ρωj ) instead of f (ωj ).
<Let the m × m Hankel matrices Ĥm and ĤmbeĤm := [µ̂i+j−2 ]mi,j=1<and Ĥm:= [µ̂i+j−1 ]mi,j=1 .

<− λĤm . Then we regard
Let ζ1 , . . . , ζm be the eigenvalues of the pencil Ĥmλ̂j = γ + ρζj , 1 ≤ j ≤ m as the approximations for λ1 , . . . , λm .
Thus, we obtain the following:

Algorithm:Input: u, v ∈ Cn , N , m, γ, ρ
Output: λ̂1 , . . . , λ̂m
1. Set ωj ← exp(2πij/N ), j = 0, . . . , N − 1
2. Form yj = ((γ + ρωj )B − A)−1 u, j = 0, . . . , N − 1
3. Set fj ← v H yj , j = 0, . . . , N − 1
4. Compute µ̂k , k = 0, . . . , 2m − 1
5<5. Compute the eigenvalues ζ1 , . . . , ζm of the pencil Ĥm− λĤm .
6. Set λ̂j ← γ + ρζj , j = 1, . . . , m.

Let us now investigate the influence of the quadrature error of µ̂k in the<algorithm. Error analyses for the eigenvalues of the pencil Hm− λHm arepresented in [12] and [17]. Let η be|λj − γ|.
j>mρη := minThen, we have the following error estimation for the computational results:|λ̂j − λj | = O(η 2m−N ),1 ≤ j ≤ m,if N ≥ 2m.

4

Numerical examples

We provide various numerical examples of the proposed method. The algorithm was implemented in Matlab Ver. 6.0 for a Linux operating system.
Computation was done with double precision arithmetic. The vectors u and v
were chosen randomly using the Matlab function ’rand’. In Example 1, 2 and 3,
the linear equations that appear in the calculation for f (z) were solved via theMatlab command ’\’. In Example 4, the linear equations were solved by thepreconditioned Bi-CGSTAB method [20] in conjugation with the incomplete
LU factorization. The stopping criterion for the residual was 10−12 .
Example 1 Let us illustrate the results in the previous sections using a simple
example [18]. Let
 99

10 ··· 0100100


..
98
 0 100

.
0

.
.
... ...
A=
(7)
.. 
 ..
,


11 
 0 · · · 0 100
100
00 ··· 00 100
and let B = I. The eigenvalues of the pencil A − λB are λj = (j − 1)/100 for
j = 1, . . . , 100.
The eigenvalues in the circle with center γ = 0.015 and radius ρ = 0.02
are λ = 0.0, 0.01, 0.02, and 0.03. The computational results obtained using theparameters m = 4 and N = 64 are given byλ̂1 = −0.0000000220414688,
λ̂2 = 0.0099996628360414,
λ̂3 = 0.0199994536835002,
λ̂4 = 0.0299998807651666.
6

The maximum error of {λ̂j } was 5.5 × 10−7 . In this case, η = 1.25, andη 2m−N = η −56 ≈ 3.7 × 10−6 . The computational results for N = 128 are givenby
λ̂1 = −0.000000000000082,
λ̂2 = 0.009999999999377,
λ̂3 = 0.019999999999253,
λ̂4 = 0.029999999999877.
The maximum error of {λ̂j } was 7.5 × 10−13 , and η 2m−N = η −120 ≈ 2.3 × 10−12 .
Example 2 A is given by (7). B is an n × n matrixµ¶
0 0
B=0 I4
where I4 is an identity matrix of size four. n this case, B is nonsingular, andthe the pencil has only four eigenvalues at 0, 0.01, 0.02 and 0.03. All of theeigenvalues are located inside the circle.
The computational results obtained using N = 16 and m = 4 are given byλ̂1 = 0.0000000000000000,
λ̂2 = 0.0099999999999999,
λ̂3 = 0.0200000000000001,
λ̂4 = 0.0300000000000000.
The maximum error of {λ̂j } was 3.2 × 10−16 .
Example 3 A problem that arises in the modal analysis of dissipative magnetohydrodynamics (MHD) [3]. Both A and B are nonsymmetric, and n = 416.
In Figure 1, we show the results obtained using parameters γ = −0.2+0.6i,
ρ = 0.05, m = 5, and N = 64. The + symbol denotes {λ∗j } that are evaluated
via the Matlab command ‘eig’. The circle formed by the dotted line shows Γ,and the ◦ symbol denotes {λ̂j }. The presented method found five eigenvalues
inside the circle. The maximum error of {λ̂j } was 8.5 × 10−7 .
Example 4 A problem that is derived from the finite element method for amolecular electronic state [9]. Both A and B are real symmetric, and n = 9264.
The number of nonzero elements is 240968 in both A and B.
We computed four eigenvalues in the interval [−10, −8]. The results obtained using γ = −9, ρ = 1, N = 32 and m = 4 are given in Table 1. In thetable, λ∗j was evaluated via the Matlab command ’eigs’ with the starting value
λ̂j .
The results obtained using γ = −9, ρ = 0.1, N = 32 and m = 3 are givenin Table 2.
7

Figure 1: Approximate eigenvalues for Example 3j1
23
4

Table 1: Example 4: γ = −9, ρ = 1, N = 32, m = 4λ̂j|λ̂j − λ∗j |−9.01565065746990 − 0.00000000921016i 1.6 × 10−6−8.98265731343177 − 0.00000003453753i 6.6 × 10−6−8.91765400027911 − 0.00000000888371i 2.0 × 10−6−8.57926086161660 − 0.00000000000034i 1.4 × 10−11

Table 2: Example 4: γ = −9, ρ = 0.1, N = 32, m = 3jλ̂j|λ̂j − λ∗j |1 −9.01565230691512 − 0.00000000000001i 2.7 × 10−13
2 −8.98266396293680 − 0.00000000000003i 3.0 × 10−13
3 −8.91765596597219 − 0.00000000000000i 2.3 × 10−138References[1] W. E. Arnoldi, The principle of minimized iteration in the solution of thematrix eigenproblem, Quarterly of Appl. Math. 9 (1951) 17–29.
[2] F. Chatelin, Eigenvalues of Matrices, John Wiley & Sons, Chichester,
1993.
[3] J. Cullum, W. Kerner and R. Willoughby, A generalized nonsymmetricLanczos procedure, Computer Phys. Communications 59 (1989) 19–48.
[4] E. R. Davidson, The iterative calculation of a few of the lowest eigenvalues and corresponding eigenvectors of large real symmetric matrices, J.
Comput. Phys. 17 (1975) 87–94.
[5] L. M. Delves and J. N. Lyness, A numerical method for locating the zerosof an analytic functions, Math. Comput. 21 (1967) 543–560.
[6] G. H. Golub and H. A. van der Vorst, Eigenvalue computation in the 20thcentury, J. Comput. Appl. Math. 123 (2000) 35–65.
[7] F. R. Gantmacher, The theory of Matrices, Chelsea, New York, 1959.
[8] I. Gohberg, P. Lancaster and I. Rodman, Matrix Polynomials, AcademicPress, New York, 1982.
[9] S. Hyodo, Meso-scale fusion: A method for molecular electronic statecalculation in inhomogeneous materials, Proc. the 21st Century, Mikkabi,
2001, in T. Mitsui (Ed.), Special issue of J. Comput. Appl. Math. (to bepublished).
[10] N. I. Ioakimidis, Quadrature methods for the determination of zeros oftranscendental functions–a review, in P. Keast and G. Fairweather (Eds.),
Numerical Integration: Recent Developments, Software and Applications,
Reidel, Dordrecht, 1987, 61–82.
[11] P. Kravanja, T. Sakurai and M. Van Barel, On locating clusters of zerosof analytic functions, BIT 39 (1999) 646–682.
[12] P. Kravanja, T. Sakurai, H. Sugiura and M. Van Barel, An error analysis of
two related quadrature methods for computing zeros of analytic functions,Part II, Tech. Report Leuven TW-338 (2002).
[13] P. Kravanja and M. Van Barel, Computing the Zeros of Analytic Functions, Lecture Notes in Mathematics 1727, Springer, Berlin, 2000.
[14] C. C. Paige, Computational variants of the Lanczos method for the eigenproblems, J. Inst. Math. Appl. 10 (1972), 373–381.
9[15] B. N. Parlett, The Symmetric Eigenvalue Problem, Prentice-Hall, Englewood Cliffs, NJ, 1980.
[16] A. Ruhe, Rational Krylov sequence methods for eigenvalue computation,
Linear Algebra Appl. 58 (1984) 391–405.
[17] T. Sakurai, P. Kravanja, H. Sugiura and M. Van Barel, An error analysis of
two related quadrature methods for computing zeros of analytic functions,J. Comput. Appl. Math. (to appear).
[18] G. De Samblanx and A. Bultheel, Using implicitly filtered RKS for generalized eigenvalue problems, J. Comput. Appl. Math. 107 (1999) 195–218.
[19] G. L. G. Sleipen and H. A. Van Der Vorst, A Jacobi-Davidson iteration
method for linear eigenvalue problem, SIAM J. Matrix Anal. Appl. 17(1996) 401–425.
[20] H. A. van der Vorst, Bi-CGSTAB: A fast and smoothly converging variantof Bi-CG for the solution of nonsymmetric linear systems, SIAM J. Sci.
Statist. Comput. 13 (1992) 631–644.

10