Novel interpretation ofcontour integral spectral projection methods forsolving generalized eigenvalue problems⋆Akira Imakura⋆,⋆ , Lei Du⋆,† , Tetsuya Sakurai⋆,†⋆University of Tsukuba, Tennodai 1-1-1, Tsukuba-city, Ibaraki, Japan⋆⋆†

JST/CREST

E-mail: imakura@cs.tsukuba.ac.jp

Abstract: For generalized eigenvalue problems, we consider computing all
eigenvalues located in a certain region and their corresponding eigenvectors.
Recently, contour integral spectral projection methods have been proposed forsuch problems. In this study, from an analysis of the relationship between
the contour integral spectral projection and the Krylov subspace, we provide anovel interpretation of these methods. We also propose a new algorithm basedon this interpretation.

1

Introduction

We consider computing all eigenvalues located in a certain region and their corresponding
eigenvectors (λi , xi ), i = 1, 2, . . . , m for the generalized eigenvalue problem in the form
Axi = λi Bxi ,xi ∈ Cn \ {0},λi ∈ Ω ⊂ C,(1)where A, B ∈ Cn×n and the matrix pencil (A, B) is assumed to be regular, i.e., A −zB is nonsingular for some z ∈ C. Such eigenvalue problems arise in many areas ofcomputational science and engineering.
A powerful algorithm for solving the generalized eigenvalue problem (1) is the contour
integral spectral projection method proposed by Sakurai and Sugiura in 2003 [5]. Since
then, several authors have actively studied and proposed improved techniques based on
the concepts of Sakurai and Sugiura [1, 2, 4, 6].
In this study, we analyze the relationship between the contour integral spectral projection and the Krylov subspace, and provide a novel interpretation of contour integral
spectral projection methods. We also propose a new algorithm based on this interpretation.
The remainder of this paper is organized as follows. In Section 2, we brieﬂy describe
the basic concepts of contour integral spectral projection methods. In Section 3, we
provide a novel interpretation and propose a new algorithm based on this interpretation.
In Section 4, we show results from numerical experiments using our proposed algorithm,and ﬁnally we draw some conclusions in Section 5.
Throughout this paper, the following notations are used. Let A = [a1 , a2 , . . . , am ] ∈Cn×m , then span{A} := span{a1 , a2 , . . . , am }. Also, let V, W be subspaces, then V + W
denotes the sum of the subspaces of V, W, i.e., V + W := {v + w|v ∈ V, w ∈ W}.

1

Algorithm 1 The Rayleigh-Ritz procedure
Input: m dimensional subspace Sm
Output: Eigenpairs (λi , xi ) for i = 1, 2, . . . , m1: Compute the orthonormal basis W = [w1 , w2 , . . . , wm ] of the subspace Sm
2: Compute eigenpairs (λi , ui ) of the matrix pencil (W H AW, W H BW ),and xi = W ui for i = 1, 2, . . . , m2

Contour integral spectral projection methods

For solving the generalized eigenvalue problem (1), a contour integral spectral projection
method was ﬁrst introduced by Sakurai and Sugiura in 2003 [5]. For solving (1), they
introduced the rational function
f (z) := uH (zB − A)−1 v,u, v ∈ Cn \ {0},(2)whose poles are eigenvalues λ of the generalized eigenvalue problem; they then considered
computing all poles located in Ω.
All poles located in a certain region of an analytic function can be computed by thealgorithm in [3], which is based on the Cauchy’s integral formula. Applying the algorithmin [3] to the rational function (2), the eigenpairs (λi , xi ) of the generalized eigenvalue
problem (1) can be obtained from the generalized eigenvalue problem with small Hankel
matricesµ0µ1 · · · µM −1µ1µ2···µM µ1 µ2µ2 · · ·µM µ3· · · µM +1 
<HM :=  ..
 , HM :=  ..
,..
..
..
..
.
.
.
.
 .
 .
.
.
.
.
.
.
µM −1 µM · · · µ2M −2µM µM +1 · · · µ2M −1
whose entries consist of the momentsI1(z − γ)k f (z)dz,µk :=2πi Γk = 0, 1, . . . ,where Γ is the positively oriented Jordan curve (the boundary of Ω), and γ is located inΩ. For details, we refer to [5].
Since then, an improvement on the Sakurai and Sugiura’s algorithm has been proposed for obtaining more accurate eigenpairs; this is called the CIRR method [6] andis based on using the Rayleigh-Ritz procedure (Algorithm 1). Block variants of the algorithms in [5] and in [6] have also been proposed [1, 2] for the higher stability of thealgorithms, especially when many eigenvalues exist in Ω.
Let V ∈ Cn×L \{O}, e.g., a random matrix. Also, let S = [S0 , S1 , . . . , SM −1 ] ∈ Cn×M L ,whereI
1(1)(L)(z − γ)k (zB − A)−1 BV dz, k = 0, 1, . . . , M − 1. (3)Sk = [sk , . . . , sk ] :=2πi ΓThen, the block CIRR method [1] constructs the M L-dimensional subspace
(1)

(L)

(1)

(L)

(1)

(L)

SM := span{S} = span{s0 , . . . , s0 , s1 , . . . , s1 , . . . , sM −1 , . . . , sM −1 },2(4)

Algorithm 2 The block CIRR method
Input: N, M, L ∈ N, V ∈ Cn×L , γ, (zj , ωj ) for j = 1, 2, . . . , NOutput: Eigenpairs (λi , xi ) for i = 1, 2, . . . , M L
1: Solve Yj = (zj B − A)−1 V for j = 1, 2, . . . , N∑N
2: Compute Sk = (1/N ) j=1 ωj (zj − γ)k Yj for k = 0, 1, . . . , M − 1
3: Compute the QR decomposition of S = [S0 , S1 , . . . , SM −1 ] : S = W R
4: Compute eigenpairs (θi , ui ) of the matrix pencil (W H AW, W H BW ),and (λi , xi ) = (θi , W ui ) for i = 1, 2, . . . , M Lfrom the Rayleigh-Ritz procedure. The subspace SM contains all eigenvectors of (1):SM = span{x1 , x2 , . . . , xm } for m ≤ M L. The block CIRR method with L = 1 reduces
to the CIRR method.
In practice, the computation of the contour integral in (3) is approximated by somenumerical integration rule such as the N -point trapezoidal ruleN1 ∑Sk ≈ Sek :=ωj (zj − γ)k (zj B − A)−1 BV,
N j=1

where zj are the quadrature points and ωj are the corresponding weights. The simple
formulation of the block CIRR method is shown in Algorithm 2; for details, we refer
to [2, 6].

3A novel interpretation and a new algorithm

For a generalization of the Jordan canonical form to the matrix pencil, we have thefollowing theorem.
Theorem 1 (Weierstrass canonical form). Let zB − A be regular. Then, there exist
nonsingular matrices Pe, Q ∈ Cn×n such thatPe(zB − A)Q = diag(zIn1 − J1 , . . . , zInd − Jd , zNd+1 − Ind+1 , . . . , zNr − Inr ),(5)where Ji , Ni are Jordan blocks, Ni is nilpotent, and Ini denotes the identity matrix oforder ni .
Based on Theorem 1, in Section 3.1, we provide a novel interpretation of the contour
integral spectral projection. In Section 3.2, we propose a new algorithm based on theblock Arnoldi method.

3.1

A novel interpretation of the contour integral spectral projectione = Q−1 . According toSince the matrices Pe, Q are nonsingular, we deﬁne P = Pe−1 , Q
e into Pei , Q
ei ∈ Cni ×nthe Jordan block structure of (5), we partition row vectors in Pe, Qni ×nand column vectors in P, Q into Pi , Qi ∈ C, respectively, for i = 1, 2, . . . , r. Then,
we can derive the following lemma and theorem.
3Lemma 2. Let ψk (z) be a k-degree polynomial. Then, we haveI1ψk (z − γ)(zB − A)−1 dz = ψk (SΓ )PΓ ,2πi ΓwhereSΓ :=∑ei ,Qi (Ji − γI)QPΓ :=i:λi ∈Ω∑

Qi Pei .

i:λi ∈Ω

Proof. From [2, Theorem 4], we have the following relationI∑1(z − γ)k (zB − A)−1 dz =Qi (Ji − γI)k Pei ,2πi Γi:λ ∈Ωk = 0, 1, . . . .

iei Qi = I and Q
ei Qj = O for i ̸= j,Here, since Q∑(Qi (Ji − γI)k Pei =i:λi ∈Ω∑)k (eiQi (Ji − γI)Qi:λi ∈Ω∑)
Qi Pei= SΓk PΓ .

i:λi ∈Ω

Therefore, Lemma 2 is proved.
Definition. Let A ∈ Cn×n , B = [b1 , b2 , . . . , bl ] ∈ Cn×l . Also, let Km (A, b1 ) be theKrylov subspace deﬁned byKm (A, b1 ) := span{b1 , Ab1 , . . . , Am−1 b1 }□Then, the subspace Km(A, B) is deﬁned by the sum of the Krylov subspaces, i.e.,□Km(A, B) := Km (A, b1 ) + Km (A, b2 ) + · · · + Km (A, bl ) = span{[B, AB, . . . , Am−1 B]}.

Theorem 3. Let SM be the subspace of the block CIRR method defined by (4). Then,we have□
SM = K M(SΓ , PΓ BV ).
Proof. From the deﬁnition of Sk (3) and Lemma 2, we have Sk = SΓk PΓ V . Therefore,□(SΓ , PΓ BV ).
SM = span{[S0 , S1 , . . . , Sk ]} = span{[PΓ BV, SΓ PΓ BV, . . . , SΓM −1 PΓ BV ]} = KM

Therefore, Theorem 4 is proved.
Remark. According to Theorem 4, the block CIRR method can be regarded as the□Rayleigh-Ritz procedure based on the block Krylov subspace KM(SΓ , PΓ BV ). Here, we□note that in the block CIRR method, the basis vectors of KM (SΓ , PΓ BV ) are explicitly
computed by (3) and the QR decomposition of S (Algorithm 2).

43.2

Proposal of a new algorithm based on the block Arnoldi
methodIn Section 3.1, we analyzed the relationship between the contour integral spectral projection and the block Krylov subspace. In this section, we consider applying the block□Arnoldi procedure to the block Krylov subspace KM(SΓ , PΓ BV ), and we propose a newalgorithm based on the block Arnoldi method.
For the eigenpairs of the generalized eigenvalue problem (1) and matrix multiplications for SΓ , we have the following two theorems.
Theorem 4. Let HM ∈ CM L×M L be the banded Hessenberg matrix obtained from the ma□(SΓ , PΓ BV ):trix formula of the block Arnoldi procedure for the block Krylov subspace KM
W H SΓ W = HM ,□where columns of W are the orthonormal basis of KM(SΓ , PΓ BV ). Also, let (θi , ui ) be
the standard eigenpairs of HM . Then, in the case of m ≤ M L, the eigenpairs (λi , xi ) of
the generalized eigenvalue problem (1) are obtained by (λi , xi ) = (θi + γ, W ui ).

Proof. From the deﬁnitions of SΓ and W , Theorem 5 is proved.
□Theorem 5. Let U = [u1 , u2 , . . . , uL ], ui ∈ KM
(SΓ , PΓ BV ). Then, there exist (M − 1)L×Ldegree polynomials ψM −1 (z) and Ψ ∈ Csuch thatI1
ψM −1 (z − γ)(zB − A)−1 BV Ψdz.
U = ψM −1 (SΓ )PΓ V Ψ =2πi Γ

We can also compute the matrix multiplications SΓ U based on the contour integral asfollows:
I1SΓ U =(z − γ)ψM −1 (z − γ)(zB − A)−1 BV Ψdz.
(6)
2πi Γ
□
Proof. From Lemma 2 and the deﬁnition of KM(SΓ , PΓ BV ), Theorem 6 is proved.

Remark. According to Theorems 5 and 6, the generalized eigenvalue problem (1) canbe solved by the block Arnoldi method for computing eigenpairs of SΓ . Here, matrix
multiplications for SΓ can be eﬃciently computed based on the contour integral (6).
Consequently, we propose a new algorithm based on the block Arnoldi method, as
shown in Algorithm 3. Here, the computation of the contour integral is approximatedby the N -point trapezoidal rule, just as in the block CIRR method (Algorithm 2).
When the contour integrals are computed without approximation (without numerical integration), the proposed method is mathematically equivalent to the block CIRR
method and to the block variant of Hankel-based method proposed in [2]. The algorithmin [4] can also be regarded as a special case of the restarted version of the proposed
algorithm with M = 1.
We show the contributions to the computational costs of the block CIRR method(Algorithm 2) and the proposed method (Algorithm 3) in Table 1. Here, we assume thatn ≫ L, M, N and that the QR decomposition of an n × m matrix requires 2nm2 ﬂops.
5

Algorithm 3 Proposed algorithm based on the block Arnoldi method
Input: N, M, L ∈ N, V ∈ Cn×L , γ, (zj , ωj ) for j = 1, 2, . . . , N
Output: Eigenpairs (λi , xi ) for i = 1, 2, . . . , M L
1: Solve Yj = (zj B − A)−1 BV for i = 1, 2, . . . , N∑N
2: W0 = (1/N ) j=1 ωj Yj
3: Compute the QR decomposition of W0 : W0 = W1 R
4: Set Ψ1,j = R−1 for j = 1, 2, . . . , N
5: for k = 1, 2, . . . , M do:e k,j = (zj − γ)Ψk,j for j = 1, 2, . . . , N
6:Ψ
fk = (1/N ) ∑N ωj Yj Ψe k,j
7:
W
j=1
8:
for i = 1, 2, . . . , k do:
fk
9:Hi,k = WiH W
e k,j = Ψe k,j − Ψi,j Hi,k for j = 1, 2, . . . , N
10:Ψ
fk = W
fk − Wi Hi,k11:
W
12:
end forfk : W
fk = Wk+1 Hk+1,k13:
Compute the QR decomposition of We k,j H −1 for j = 1, 2, . . . , N
14:
Ψk+1,j = Ψk+1,k15: end for
16: Set W = [W1 , W2 , . . . , WM ] and HM = {Hi,j }1≤i,j≤M
17: Compute eigenpairs (θi , ui ) of HM , and (λi , xi ) = (θi + γ, W ui ) for i = 1, 2, . . . , M L
Table 1: Main contributions to the computational costs of the block CIRR method
and the proposed algorithm. Here, GEP and SEP denote the generalized and standard
eigenvalue problems, respectively.
Operation type
Solve (zj B − A)Yj = BV
Mat-Vec for A and B
Solve small eigenvalue problem
Construct the orthonormal basis
Set small eigenvalue problem
Compute eigenvectors

Block CIRR
NM L (for A) + M L (for B)
GEP: 12nLM N + 2n(LM )24n(LM )22nM L × m[times][time][time][ﬂops][ﬂops][ﬂops]

Proposed algorithm
N [times]
0 [times]
SEP: 1[time]
222nL M N + 2n(LM )[ﬂops]
0[ﬂops]2nM L × m
[ﬂops]

Table 1 shows that the proposed algorithm is more eﬃcient than the block CIRR method.
Note that the proposed algorithm has approximately the same storage requirements as
the block CIRR method.
As another advantage of the proposed algorithm, which is based on the block Arnoldi
method, it is expected that some improvement techniques used for the (block) Arnoldi
method can be applied to the proposed algorithm (Algorithm 3), such as some type ofrestart techniques, the (block) Lanczos method.

6

−8

proposal
Block CIRR

−5

log10 of relative residual

log10 of relative residual−4−6−7−8−9−10proposal
Block CIRR

−9−10−11−12−13−14123451 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

index of eigenpairs

index of eigenpairs(a) The first example: MHD(b) The second example: ESC

Figure 1: Relative residuals of eigenpairs: ∥Axi − λi Bxi ∥2 /(∥Axi ∥2 + |λi |∥Bxi ∥2 ) fromcomputational results for the two test problems.

4Numerical examples and results

To conﬁrm that the proposed algorithm (Algorithm 3) can solve generalized eigenvalueproblems (1) just as well as the block CIRR method, we provide two numerical examples.
The algorithms were implemented in MATLAB R2013a, and the computation of thecontour integral was approximated by the N -point trapezoidal rule. The input matrix
V was set as a random matrix generated by the MATLAB command “rand”, and eachlinear system was solved by the MATLAB command “\”.
The ﬁrst test problem was taken from a model analysis of dissipative magnetohydrodynamics (MHD). The obtained matrices A and B were real nonsymmetric and n = 416.
We set Γ to be the circle with center γ = −0.2 + 0.6i and radius ρ = 0.05, and consider
ﬁnding ﬁve eigenvalues located within the circle. We also set L = 1, M = 5, and N = 64.
The second test problem was taken from an electronic structure calculation (ESC) fora system of silicon. The obtained matrices A and B were Hermitian, and B was positivedeﬁnite, where n = 1726. The eigenvalues are located on the real axis. To ﬁnd sixteen
eigenvalues in the interval [−0.17, 0.17], we set Γ to be the circle with center γ = 0 andradius ρ = 0.17. We also set L = 4, M = 4, and N = 64.
The relative residuals of eigenpairs for the ﬁrst and second test problems are shownin Figures 1(a) and 1(b), respectively. Figure 1 shows that the proposed method (Algorithm 3) can compute eigenpairs with approximately the same accuracy as the block
CIRR method (Algorithm 2).

5Conclusions

In this paper, we have provided a novel interpretation of contour integral spectral projection methods for solving the generalized eigenvalue problem (1). We have also proposeda new algorithm based on the block Arnoldi method.
In future, we will propose and analyze eﬃcient improvements in the proposed algorithm based on some techniques that are applied to the block Arnoldi method.

7

Acknowledgements
This work is supported partly by Strategic Programs for Innovative Research (SPIRE)
Field 5 “The origin of matter and the universe”, JST/CREST and KAKENHI (GrantNo. 25870099).

References[1] T. Ikegami and T. Sakurai, Contour integral eigensolver for non-Hermitian systems:
a Rayleigh-Ritz-type approach, Taiwan. J. Math. 2010; 14(3):825–837.
[2] T. Ikegami, T. Sakurai and U. Nagashima, A ﬁlter diagonalization for generalized
eigenvalue problems based on the Sakurai-Sugiura projection method, J. Comput.
Appl. Math. 2010; 233:1927–1936.
[3] P. Kravanja, T. Sakurai and van M. Barel, On locating clusters of zeros of analyticfunctions, BIT 1999; 39(4):646–682.
[4] E. Polizzi, A density matrix-based algorithm for solving eigenvalue problems, Phys.
Rev. B 2009; 79:115112.
[5] T. Sakurai and H. Sugiura, A projection method for generalized eigenvalue problems
using numerical integration, J. Comput. Appl. Math. 2003; 159:119–128.
[6] T. Sakurai and H. Tadano, CIRR: a Rayleigh-Ritz type method with counter integral
for generalized eigenvalue problems, Hokkaido Math. J. 2007; 36:745–757.

8