Eﬃcient parameter estimation and implementationof a contour integral-based eigensolver
Tetsuya Sakurai†,‡ , Yasunori Futamura† and Hiroto Tadano††
Department of Computer Science, University of Tsukuba‡
CREST, Japan Science and Technology Agency

Abstract: We consider an eigensolver for computing eigenvalues in agiven domain and the corresponding eigenvectors of large-scale matrix
pencils. The Sakurai-Sugiura (SS) method is an eigensolver based oncomplex moments given by contour integrals of matrix inverses withseveral shift points. This method has good parallel scalability, and issuitable for massively parallel computing environments. The SS method
has several parameters, and the choice of these parameters is crucial
for achieving high accuracy and good parallel performance. We discuss
some numerical properties of the method, and present eﬃcient parameter
estimation techniques. We demonstrate the eﬃciency of our method withnumerical experiments.

1

INTRODUCTION

Large-scale eigenvalue problems arise in many scientiﬁc and engineering areas such asnano simulation, vibration analysis, data analysis, etc. Massively parallel computers
are used to solve such large-scale problems, and they require eﬃcient algorithms forparallel computing.
In this paper, we consider the problem of ﬁnding several eigenvalues in a given
domain and their corresponding eigenvectors of the generalized eigenvalue problem
Ax = λBx,
where A, B ∈ Cn×n . In [14], an eigensolver for generalized eigenvalue problems using complex moments is proposed. This method is called the Sakurai-Sugiura (SS)method. In the SS method, contour integrals with a source vector v are used togenerate a subspace spanned by a set of eigenvectors with respect to the eigenvalues
in a target domain. A large-scale eigenvalue problem is reduced to a small eigenvalue problem with Hankel matrices constructed from complex moments. In [16], aninterpretation for ﬁltering of spectrum is used to discuss numerical properties of acontour integral approximated by numerical quadrature. An inﬂuence of approximation by numerical quadrature is considered as a contamination of eigencomponents,
and the choice of an appropriate subspace size provides accurate eigenpairs in atarget domain.
1

A variant of the SS method that improves numerical accuracy by using theRayleigh-Ritz procedure is presented in [15]. Ikegami, et al. [6, 7] presented a block
version of the SS method that uses multiple source vectors instead of the single
source vector for the contour integrals. The block SS method improves numericalstability when the target domain contains many eigenvalues. Moreover, this method
can treat multiple eigenvalues. In [1, 2], the SS method is extended to nonlinear
eigenvalue problems. As related works of eigensolvers using contour integrals, Polizzi[13] proposed an iterative reﬁnement of a contour integral method for symmetric or
Hermitian positive deﬁnite eigenvalue problems. Beyn [3] proposed a method for
nonlinear eigenvalue problems using contour integrals with a singular value decomposition of a matrix with a Hankel type structure. Yokota, et al. [21] proposed aRayleigh-Ritz type method using contour integrals for nonlinear eigenvalue problems. In this method, a subspace that includes target eigenvectors are generated bycontour integrals, and a large-scale nonlinear eigenvalue problem is projected to asmall nonlinear eigenvalue problem, and the projected problem is solved by Hankel
type nonlinear eigensolver using contour integrals.
The SS method computes a set of eigenvalues by computing the solutions tosystems of linear equations(zj B − A)Yj = BV,

j = 1, . . . , N,(1)where V is a matrix with L column vectors and zj is a shift point on the complex
plane. The method computes the desired eigenvalues inside of a border deﬁned
by the set of shifts {zj }. The ﬁrst step of the SS method is the construction of
a subspace that includes the eigenvectors corresponding to the eigenvalues locatedinside the given domain. In this step, solutions of linear systems at several shift
points are used. The second step is to solve the projected problem in the subspace
and to extract the approximate eigenvalues and the corresponding eigenvectors forthe original problem. Since the size of the projected subspace is assumed to be small
compared with the original matrix size, the computational costs of the ﬁrst step isdominant.
Krylov subspace methods for multiple right-hand sides are eﬃcient for solving
the linear systems (1). In [11, 18], methods to improve numerical stability and convergence for block Krylov subspace methods are presented. In the case of standard
eigenvalue problems, the linear systems (1) are shifted linear systems, and a shift
invariance of the Krylov subspace reduces computational costs to obtain solutionsof linear systems at several shift points [12]. The application of the SS method withthe shifted CG method for shell model calculations is reported in [10]. Futamura,et al. [5] presented an eﬃcient implementation techniques for solving such shifted
linear systems in parallel computing, and a report of application for band structure
calculations with the real space density functional theory on the K computer is presented. Yamazaki, et al. [20] implemented a nonlinear version of the SS method,and evaluated parallel performances of the method.
Each of the linear systems is independent with respect to the other shifts, so eachcan be solved without any consideration of the nodes assigned to diﬀerent shifts indistributed computing. Therefore, the method provides coarse-grained parallelismof computation. By employing a parallel linear solver for each shift point, the total
2

number of nodes is the product of the number of nodes assigned for each linear
system and the number of shift points.
The SS method has several parameters, and the choice of these parameters is
crucial for achieving high accuracy and good parallel performance. In this paper, weshow some numerical properties of the method. The contour integral for a matrix
inverse is regarded as a ﬁlter for an eigensubspace. When the contour integral isapproximated by numerical quadrature, the quadrature error causes contaminationof the eigencomponents corresponding to the eigenvalues located outside of the contour path. Based on these properties, we propose eﬃcient parameter estimation
techniques for the SS method.
In [4], a method for stochastic estimation of number of eigenvalues in a given
domain is proposed. This estimation can be used for predicting appropriate parameters. Maeda, et al. [9] extended this eigenvalue count method to nonlinear
eigenvalue problems.
The rest of this paper is organized as follows. In Section 2, we brieﬂy introduce
the SS method. In Section 3, the properties of numerical quadrature applied for amatrix inverse are discussed. In Section 4, eﬃcient parameter estimation methods
are presented. Some numerical experiments are shown in Section 5. The last section
concludes the paper.

2

A CONTOUR INTEGRAL-BASED EIGENSOLVER

In this section, we brieﬂy introduce the SS method. For matrices A, B ∈ Cn×n ,let λ1 , . . . , λn be eigenvalues of the matrix pencil A − λB, and let x1 , . . . , xn becorresponding eigenvectors. Let Γ be a positively oriented closed Jordan curve inthe complex plane, and let G be a domain for which the border is given by Γ. We willﬁnd the eigenvalues inside Γ and the corresponding eigenvectors by using contour
integrals.

2.1

Eigensubspace obtained by contour integrals

Suppose that m eigenvalues λ1 , . . . , λm are located inside Γ, and other eigenvalues
are located outside Γ. Deﬁne a sequence of matrices F0 , F1 , . . . as∫1Fk =z k (zB − A)−1 B dz, k = 0, 1, . . . .
(2)
2πi Γ
For a matrix V ∈ Rn×L with a positive integer L, let∫1z k (zB − A)−1 B dz V, k = 0, . . . , M − 1,
Sk = Fk V =2πi Γwhere M is chosen such that LM ≥ m, and set
F = [F0 , F1 , . . . , FM −1 ]andS = [S0 , S1 , . . . , SM −1 ].
3(3)

According to [14], the column vectors of S are given by linear combinations of
the eigenvectors with respect to the eigenvalues located inside Γ, and thusspan(S) = span(x1 , . . . , xm ),if the column space of V includes x1 , . . . , xm . V is called a source matrix for the contour integral. In practice, the elements of V are set by a random number generator.
The eigenvectors x1 , . . . , xm are obtained from S when the maximum multiplicityof the eigenvalues in Γ is less than or equal to L.
Using the Rayleigh-Ritz procedure with S, we can extract the eigenpairs. Let
the singular value decomposition of S be
S = U ΣW H ,where Σ = diag(σ1 , . . . , σLM ), U ∈ Cn×LM and W ∈ CLM ×LM . Since the rank of S
is m, σm 6= 0 and σm+1 = · · · = σLM = 0. Setting Um = U (:, 1 : m), we calculate
the projected matrices as
H
Am = Um
AUm ,H
Bm = Um
BUm .

(4)

Let ω1 , . . . , ωm be the eigenvalues of the matrix pencil Am − λBm , and let r1 , . . . , rm
be the corresponding eigenvectors. Then the eigenvalues inside Γ of the matrix pencil
A − λB are given by
λi = ωi , i = 1, . . . , m,and the corresponding eigenvectors are given by
xj = Um rj ,j = 1, . . . , m.

(5)

When the matrices are large, storage of S and computation of the singular value
decomposition restrict the application size of the method. The use of Hankel matrices reduces the memory requirement and computational costs. Let Mk ∈ CL×Lbe∫1
z k V T (zB − A)−1 BV dz.
(6)
Mk =
2πi Γ
<
Let the Hankel matrices HLM , HLM∈ CLM ×LM beM0
M1 · · · MM −1
 M1
M2 · · ·MM 
HLM = 

..
..
..


.
.
.

MM −1 MM
and



M1
 M2

<
HLM
= .
 ..
MM

M2
M3
..
.

···

M2M −2

···
···


MM
MM +1 

.
..

.

MM +1 · · ·4M2M −1

Let the singular value decomposition of H̃m beH̃m = Ũ Σ̃W̃ H ,(7)< = H < (1 : m, 1 : m). Let ω̃ , . . . , ω̃ andwhere H̃m = HLM (1 : m, 1 : m) and H̃m1m
LMq̃1 , . . . , q̃m be the eigenvalues and the corresponding eigenvectors such that<(Σ̃−1 Ũ H H̃mW̃ )q̃i = ω̃i q̃i ,i = 1, . . . , m.

Then the eigenvalues of the matrix pencil A − λB in Γ are given by
λi = ω̃i .
The eigenvectors are given by
xi = S(:, 1 : m)W̃ qi ,i = 1, . . . , m.

In this computation, the singular value decomposition of S does not required. A
disadvantage using the Hankel matrices with the moment matrices Mk is numericalinstability comparing with the Rayleigh-Ritz procedure in the case of numericalcomputation with large m.
In the case of the nonlinear eigenvalue problem T (λ)x = 0 with a matrix valued
function T (λ), the integrand V T (zB − A)−1 BV in (6) is replaced by V T T (z)−1 V[1, 2]. Note that the derived eigenvalue problem with Hankel matrices are linear even
if the original problem is nonlinear. In [3], the integrand in the contour integral (6)is replaced by T (z)−1 V instead of V T T (z)−1 V .

2.2

Approximation by a numerical quadrature

The contour integral in (2) is approximated by an N -point numerical quadrature.
Suppose that a Jordan curve Γ is represented by scaling and shifting from a Jordan
curve Γ0 with a scaling factor ρ and sa hift γ. Without any loss of generality, we
assume that Γ0 encloses the origin. Let ζ(θ) be a point on Γ0 with a parameter θ,
0 ≤ θ ≤ 2π, and let z on Γ be given byz(θ) = γ + ρζ(θ).
Then the contour integral of a function f (z) is given by
1
2πi

∫

1
f (z)dz =
2π
Γ

∫ 2π
0

1
f (z)(−iρζ (θ))dθ =
2π
0

∫ 2π
ρf (z)w(θ)dθ,(8)0where w(θ) = −iζ 0 (θ). The integral (8) is approximated by the N -point quadrature
rule
∫
N
∑
1
f (z)dz ≈ρwj f (zj ),
(9)
2πi Γ
j=1

where wj = w(θj )∆j /(2π), ζj = ζ(θj ) and zj = γ + ρζj with appropriate θj and ∆j ,j = 1, . . . , N .
5Since
1
2πi

{

∫
k

ζ dζ =
Γ0

1, k = −10, otherwise

for integer k, the quadrature points ζ1 , . . . , ζN on Γ0 and the corresponding weights
w1 , . . . , wN are set to satisfyN∑{
wj ζjk =j=1

ν 6= 0, k = −1,0,
k = 0, . . . , N − 2(10)where ν is a nonzero constant.
In particular, when Γ is a circle with center γ and radius ρ, and the quadraturepoints are set aszj = γ + ρ(cos θj + i sin θj ),j = 1, . . . , N,

where θj = (2π/N ) × (j − 1/2), j = 1, . . . , N , then Γ0 is the unit circle and thequadrature weights are given by
wj = cos θj + i sin θj ,j = 1, . . . , N.

In the case that all the eigenvalues are located on the real axis, it might be betterto put the quadrature points closer to the real axis as follows:
zj = γ + ρ(cos θj + iα sin θj ),j = 1, . . . , N(11)with a vertical scaling factor 0 < α < 1. The corresponding quadrature weights are
given bywj = α cos θj + i sin θj , j = 1, . . . , N.
(12)
In [12], quadrature points are set on straight lines to reuse solutions of linear systems.
The Gauss-Legendre quadrature rule on a circle is used for the numerical quadraturein [13].
Using the quadrature rule (9), Fk and Sk are approximated by
Fk ≈ F̂k =N∑ρwj ζjk (zj B − A)−1 B(13)ρwj ζjk (zj B − A)−1 BV.

(14)

j=1

and
Ŝk = F̂k V =

N
∑
j=1

Matrices F and S are approximated by F̂ = [F̂0 , . . . , F̂M −1 ] and Ŝ = [Ŝ0 , . . . , ŜM −1 ].
The Rayleigh-Ritz procedure for Ŝ gives the approximate eigenvalues λ̂i and theeigenvectors x̂i . Let the singular value decomposition of Ŝ be
Ŝ = Û Σ̂Ŵ H ,6

where Σ̂ = diag(σ̂1 , . . . , σ̂LM ). Let K be the number of singular values of Ŝ thatsatisfy σ̂i ≥ δ, 1 ≤ i ≤ K with small δ > 0. We calculate the projected matrices as
Â = Û (:, 1 : K)H (A − γB)Û (:, 1 : K),

B̂ = Û (:, 1 : K)H B Û (:, 1 : K).

(15)

Let ω̂1 , . . . , ω̂K be the eigenvalues of the matrix pencil Â − λB̂, and let r̂1 , . . . , r̂K
be the corresponding eigenvectors. Then the approximate eigenvalues inside Γ aregiven byλ̂i = γ + ω̂i , i = 1, . . . , K,
and the corresponding eigenvectors are given byx̂j = Û (:, 1 : K)r̂j ,3

j = 1, . . . , K.

(16)

FILTERING FOR A SUBSPACEIn this section, we discuss the properties of the subspace obtained by the numericalquadrature (14) from the view-point of a ﬁlter for a subspace.
Here, for simplicity, we consider the case that all the eigenvalues inside Γ aresimple, and the inverse of the matrix zB − A is expanded asn∑xi y H(zB − A)−1 =i

i=1z − λi,(17)where xi and yi are the right and left eigenvectors corresponding to the eigenvalue λi .
This expansion can be generalized to non-Hermitian and nonlinear cases ([2, 3, 14]).
Let Pi = xi yiH B, 1 ≤ i ≤ n. With the expansion (17), from the residue theorem,
we have∫1Fk =z k (zB − A)−1 Bdz
2πi Γ
)
∫
n (
∑
1
z k Pidz=
2πi Γ z − λi
=

i=1
m
∑

λki Pi ,i=1

and
Sk = F k V =m∑

λki Pi V.

i=1

Deﬁne a function Fk (λ) as
1
Fk (λ) =
2πi
{

Then
Fk (λi ) =

∫

zk
dz.
Γ z−λ

λki λi ∈ G,
0, otherwise
7and Sk is represented as
Sk =

n
∑

Fk (λi )Pi V.

i=1

This equation shows that a projected component associated with Pi in V is ﬁlteredwith the factor Fk (λi ). Therefore the function Fk (λ) is regarded to give the factor
of ﬁltering with respect to λ.
For the case that the contour integral is approximated by the numerical quadrature, we deﬁne the corresponding ﬁlter function byF̂k (λ) =N∑ρwj ζjkzj − λ

j=1.

The following result is obtained.
Theorem 1 Let λ be a complex number that is located outside Γ. Then the following
holds:
∞∑
νN +p −p F̂k (λ) = −νN −1 η −N +k 1 + η −1
η
,
(18)
νN −1
p=0

where η = (λ − γ)/ρ and νp =∑N

pj=1 wj ζj .

Proof. Since |η| = |(λ − γ)/ρ| > |ζj | for 1 ≤ j ≤ N , we haveN∑
ρwj ζjk
j=1

zj − λ

=

N
∑
j=1

=

N
∑
wj ζjkwj ζjk=(zj − γ)/ρ − (λ − γ)/ρ
ζj − η

N (
∑
j=1

= −

∞
∑

−1
η


)

j=1

wj ζjk1 − ζj /ηη −p−1

N
∑

p=0



wj ζjp+k  .

j=1

Since the quadrature weights w1 , . . . wN satisfyN∑

k = 0, . . . , N − 2,

wj ζjk = 0, 1,

j=1

we haveF̂k (λ) =N∑ρwj ζjk
j=1

zj − λ∞∑= −η −p−1

p=N −1−k= − νN −1 η −N +k +

N
∑
j=1

∞
∑
p=0

8


wj ζjp+k νN +p η −N +k−1−p  .

Thus we have eqn (18).
If |(λ − γ)/ρ| is suﬃciently large then the ﬁlter F̂k (λ) is approximated byF̂k (λ) =N∑ρwj ζjk
j=1

zj − λ(≈ −νN −1λ−γρ)−N +k.

(19)

This implies that the eigencomponents corresponding to the eigenvalues locatedoutside Γ in each column vector of Ŝk = F̂k V are reduced in proportion to the(−N + k)-th power of magnitude of the scaled distance |(λ − γ)/ρ|.
Suppose that the integer m0 is taken as)(
λi − γ −N +M −1≤ δ, m0 < i ≤ n
(20)
νN −1
ρ
with small δ > 0. Then, from eqn (18), we haveŜk = F̂k V =n∑0F̂k (λi )Pi V =

i=1

m
∑

F̂k (λi )Pi V + O(δ).

i=1

4

EFFICIENT PARAMETER ESTIMATION AND IMPLEMENTATION

4.1

Selection of subspace size

The SS method has some parameters, and the choice of these parameters aﬀects
the accuracy and performance of the method. The number of quadrature points N
determines the number of systems of linear equations to solve, and consequently N
speciﬁes the number of computing nodes to use in parallel computing. Therefore
we assume that N is ﬁxed in advance. In practice, N is chosen as N = 16 or 32
depending on the number of computing nodes or memory requirements, and it is
not necessary to take a large N to reduce the quadrature error as was observed in
the previous section.
The parameter M speciﬁes the upper bound of the degree of moments. Increasing
M gives a larger subspace size LM . However, the decay factor of the ﬁlter depends
on −N + k with 0 ≤ k ≤ M − 1, and a large M diminishes the performance of
the ﬁlter. Considering a performance of the ﬁler and computational costs, we set
M = N/4.
The number of column vectors LM of Ŝ should be taken such that the minimum
singular value of Ŝ becomes suﬃciently small. Since M depends on N , we shallextend the number of column vectors of Ŝ by increasing the number of source vectors
L. Since m0 is larger than or equal to m, an approximation for m can be used as alower bound of m0 . To predict m, we use the stochastic estimation method proposedin [4, 9]. The number of eigenvalues inside Γ is given by∫1
m=
trace((zB − A)−1 B)dz.
2πi Γ9However, the cost of computing the trace of the matrix inverse is large if the matrices
are large. The trace of (zB − A)−1 B is approximated by−1

trace((zB − A)

(
B) ≈

1
L0

)∑
L0

viT (zB − A)−1 Bvi ,i=1

with some integer L0 , where the elements of the sample vectors vi ∈ Rn are taken as−1 or 1 with equal probability. The contour integral is approximated by the N -point
numerical quadrature. Thus, the estimated number of eigenvalues m̃ is given by
m ≈

N
∑

(
)
ρwj trace((zj B − A)−1 B)

j=1

≈

N
∑

((
ρwj

j=1

1
L0

)∑
L0

)
viT (zj B − A)−1 Bvi

= m̃.

i=1

Note that
m̃ =

N
∑

((
ρwj

j=1

1
L0

)∑
L0

)
viT (zj B − A)−1 Bvi= trace(V T Ŝ0 )/L0 ,i=1

and Ŝ0 can be used for the stochastic estimation of the number of eigenvalues by
setting V = [v1 , . . . , vL0 ] in the SS method.
Using m̃, we set the approximation of m0 as κm̃ with a parameter κ ≥ 1, andconsequently we set L = dm0 /M e ≈ dκm̃/M e, where dxe returns the smallest integer
not less than x. When the subspace size LM is not suﬃciently large, the minimum
singular value σmin of Ŝ is not small. In this case, we increment L until σmin satisﬁes
the condition σmin ≤ δ ×σ1 with small δ > 0. The computation of the singular valuesof Ŝ is rather expensive, so we may use the Hankel matrix Ĥ instead of Ŝ.

4.2

Iterative refinement of a subspace

After setting appropriate L, we apply the Rayleigh-Ritz procedure with Ŝ. The
increase of L causes an increase in the size of the projected subspace. It causes
an increase in the cost for computing the singular value decomposition of Ŝ andthe solution of the projected eigenvalue problem with matrices Â and B̂. To avoid
increasing the size of the projected space, we restrict the size of L, and apply therecurrence reﬁnement described below.
(0)
Setting Ŝ0 = Ŝ0 , and recurrently applying F̂0 , we have
(r−1)

Ŝ0
(r−1)

Using Ŝ0

(r−2)

= F̂0 Ŝ0(0)= · · · = (F̂0 )r−1 Ŝ0 .

(21), the output matrix with r reﬁnements is given by
(r)

(r−1)

Ŝk = F̂k Ŝ0,

k = 0, . . . , M − 1,

10

(22)

(r)

(r)

and Ŝ (r) = [Ŝ0 , . . . , ŜM −1 ]. The corresponding ﬁlter is given by (Fk (λ))r and isapproximated by()
λ − γ −r(N −k)(F̂k (λ))r ≈ (−νN −1 )r.
ρTherefore the recurrence application of the ﬁlter process makes the decay factor of
the ﬁlter smaller. The reﬁnement is terminated if the smallest singular value of Ŝ (r)becomes suﬃciently small with a threshold δ > 0.
In the case that some residuals of the obtained approximate eigenpairs are notsmall enough for a given tolerance, we can brush up the resulting approximate
eigenpairs by setting the source matrix of the SS method as
V = [x̂1 , . . . , x̂m̂ ]C,where C ∈ Rm̂×L for which the elements are given by random numbers, and x̂1 , . . . , x̂m̂
are the selected eigenvectors that are regarded as the approximate eigenvectors withrespect to the eigenvalues inside Γ. This reﬁnement technique using approximate
eigenvectors for the source matrix V is used in [13].

4.3Linear solvers for a complex shift

When A and B are real symmetric, the shifted matrix C = zB − A with a complex
shift z is complex symmetric. Therefore, a linear solver for complex symmetricsystems is used to solve the system(zB − A)Y = BV.
For a direct solver, the modiﬁed Cholesky factorization saves computational costs
for factorization. For an iterative solver, Krylov subspace methods for complex
symmetric systems, such as the COCG method, can be used.
When Γ is symmetric with respect to the real axis, the quadrature points are set
as zN −j+1 = z j , j = 1, . . . , N/2. Then, for real matrices A and B, the solutions at
zN −j+1 are obtained by
YN −j+1 = (zN −j+1 B − A)−1 BV = Y j
without any computations on zN −j+1 .
When A and B are Hermitian, we use the property(zj B − A)H = z j B H − AH = z j B − A.
If the LU factorization at zj is calculated as zj B − A = LU then we have
YN −j+1 = (z j B − A)−1 V = (U H LH )−1 V.
Therefore the LU factorization at zj can be used for the calculation at zN −j+1 .

11

5

NUMERICAL EXAMPLESIn this section we show some numerical examples. The computations are performed
in MATLAB 8.0.0. in double precision arithmetic. Random numbers are generated
by the function rand, and the projected small eigenvalue problems are solved by eig.
The systems of linear equations are solved by lu. The factorized matrices are held
during the computation, and only triangular solves are applied in the recurrence
reﬁnements.
In the following examples, the quadrature points are set by eqn (11) and thecorresponding weights are set by eqn (12) with α = 0.1. The relative residual forthe eigenpair (λ̂i , x̂i ) is calculated byresi =

kAx̂i − λ̂B x̂i k2
kAx̂i k2 + |λ̂i |kB x̂i k2.

We removed the eigenvalues with resi ≥ 10−2 inside Γ as spurious eigenvalues.
Example 1. The matrices A and B are taken from BCSSTK11 and BCSSTM11of the BCS Structural Engineering Matrices in Matrix Market [8]. A and B are
real symmetric and B is positive deﬁnite. The matrix dimension is n = 1, 473 with34, 241 nonzero entries. The parameters are set as N = 16 and L = 16. The domain
is set as γ = 103 and ρ = 5 × 102 . In this example, L is ﬁxed, and the iterative
reﬁnement is not applied.
The results are shown in Table 1. The number of singular values that are greaterthan δ = 10−12 is K = 18. Therefore 18 eigenvalues are obtained from the projected
problem, of which 7 eigenvalues are located inside Γ. The residuals of the eigenvalues
located inside Γ are small, however the residuals of the eigenvalues located outside
Γ are related to the scaled distance |ηi | = |(λi − γ)/ρ|.
Example 2. In this example, we apply the iterative reﬁnement deﬁned by (21)and (22). The matrices A and B are same as Example 1. The parameters are set as
N = 16, L = 16, and the domain is set as γ = 2 × 105 and ρ = 2 × 104 .
In Figure 1, The singular values of Ŝ (r) at r-th reﬁnement are shown. We cansee that the ratio of the minimum singular value and the maximum singular value
increases by the iterative reﬁnement. After two reﬁnement, the minimum singularvalue becomes small enough. Table 2 shows the residuals of the calculated eigenvalues located inside Γ. In the table, the notation mean(resi ) is given by the geometricmean of the residuals deﬁned bymean(resi ) =

( m̂
∏

)1/m̂
resi

,

i=1

where m̂ is the number of calculated eigenvalues located inside Γ.
Example 3. In this example, we use the stochastic estimation of the number of
eigenvalues in Γ to set the initial L, and the iterative reﬁnement of Ŝ is also used.
The matrices A and B are same as Example 1. The parameters are set as N = 16
12

Table 1: Results of Example 1.
λ̂i2345.08723030540
2398.81729572773
2628.94468521146
2723.54384863656
3383.97540832681
3501.25383608303
3561.62085364923
3629.33212408543
3796.50112783802
4022.39762561787
4100.71462746484
4175.86741050601
4770.43635520514
5071.04303115872
5185.64239506030
5325.06302301902
5608.24863853754
5874.78406307974

resi
3.0 × 10−01
3.3 × 10−01
4.9 × 10−02
1.4 × 10−02
3.8 × 10−08
9.0 × 10−11
2.7 × 10−11
4.0 × 10−11
4.8 × 10−11
3.1 × 10−11
1.5 × 10−11
3.4 × 10−11
5.1 × 10−06
1.6 × 10−04
3.5 × 10−03
1.6 × 10−02
1.0 × 10−01
6.6 × 10−01

|ηi |
3.3
3.2
2.7
2.6
1.2
–
–
–
–
–
–
–
1.5
2.1
2.4
2.7
3.2
3.8

r=0

log10 (σi /σ1 )

i
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

r=1
r=2

i
Figure 1: Singular values in r-th iterative reﬁnement.

13

Table 2: Results in Example 2.
]reﬁnement012min(resi )1.8 × 10−07
3.8 × 10−12
1.9 × 10−14mean(resi )9.1 × 10−06
1.1 × 10−10
7.2 × 10−13

max(resi )1.7 × 10−04
1.1 × 10−09
1.2 × 10−11

and δ = 10−12 . The domain is set as γ = 2 × 105 and ρ = 2 × 104 . The numberof sample vectors for the stochastic estimation of the number of eigenvalues in Γ isset as L0 = 16. The initial guess of the number of column vectors of V is given by
L = d2m̃/M e, i.e. κ = 2.
In Table 3, we show the residuals of the eigenvalues located inside Γ. The numberof eigenvalues in Γ is m = 30 and the estimated number of eigenvalues is m̃ = 32.7.
The number of iterative reﬁnement is 2.
Example 4. The matrices A and B are taken from BCSSTK13 and BCSSTM13.
A and B are real symmetric and B is positive semi-deﬁnite. The matrix dimension isn = 2, 003 with 83, 883 nonzero entries. The parameters are the same as in Example
3. The domain is set as γ = 106 and ρ = 4 × 105 .
In Table 4, we show the residuals of the eigenvalues located inside Γ. The numberof eigenvalues in Γ is m = 73 and the estimated number of eigenvalues is m̃ = 77.7.
The number of column vectors of V is L = 55 and the number of iterative reﬁnementis 2. The maximum, mean and minimum residuals are 2.1 × 10−10 , 8.6 × 10−12 and2.7 × 10−13 , respectively. We can obtain the eigenpairs in the given domain with thesame initial parameters.
Example 5.
The matrices A and B are derived from molecular orbital calculations for a model DNA [19]. A and B are real symmetric and B is positivedeﬁnite. The matrix dimension is n = 1, 980 with 728, 080 nonzero entries. The
parameters are the same as in Example 3 and 4. The domains are given by the intervals [−0.20, −0.15], [−0.25, −0.15], [−0.30, −0.15], [−0.35, −0.15], [−0.40, −0.15],[−0.45, −0.15] and [−0.50, −0.15].
In Table 5, we show the number of eigenvalues in the given interval (]ev), theestimated number of eigenvalues (Est. ]ev), the number of column vectors of V(L), the number of iterative reﬁnement (]reﬁnement) and the maximum residuals of
eigenvalues in the interval (max(resi )). In the results, the maximum residuals are
suﬃciently small by estimating appropriate L and the number of iterative reﬁnement
for each domain.

6

CONCLUSIONSIn this paper, we considered an eigensolver for computing the eigenvalues in agiven domain and the corresponding eigenvectors of large-scale matrix pencils. The

14

Table 3: Results in Example 3.
i
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

λ̂i
181301.355856
181353.297523
185810.063953
185856.309721
189076.069885
190580.274469
191916.768828
192249.997887
192450.352262
195110.875562
195362.147280
195522.864186
196453.465229
196779.318796
203358.448118

resi
3.0 × 10−12
8.2 × 10−13
3.1 × 10−12
2.2 × 10−12
1.3 × 10−12
1.7 × 10−12
4.6 × 10−12
6.2 × 10−12
8.8 × 10−12
8.9 × 10−13
1.6 × 10−12
2.1 × 10−12
9.5 × 10−13
1.1 × 10−12
5.6 × 10−13

i
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

λ̂i
206423.180896
207887.176182
209720.799807
211359.608331
211525.005509
211778.728062
211798.736010
214623.208612
215071.649241
216638.323804
216782.856683
216875.914785
217120.082795
217475.120411
217803.381541

resi
2.2 × 10−12
4.5 × 10−12
1.2 × 10−12
1.6 × 10−12
1.2 × 10−12
1.0 × 10−12
1.4 × 10−12
1.7 × 10−12
1.2 × 10−12
1.1 × 10−12
5.0 × 10−13
4.9 × 10−13
1.4 × 10−12
1.3 × 10−13
5.8 × 10−13

Table 4: Results in Example 4.
i
1
2
3
4
5
..
.

λ̂i
602514.527692
605178.148251
616657.672408
623758.141144
641859.031825
..
.

33
34
35
36
37

924036.280859
927854.750782
941218.254886
942132.221466
960716.560772

resi
1.2 × 10−12
2.1 × 10−11
5.2 × 10−12
2.2 × 10−11
1.3 × 10−12..
.

7.0 × 10−11
4.3 × 10−12
9.8 × 10−12
1.5 × 10−12
1.0 × 10−11

i
38
39
40
41
42
..
.

λ̂i
964884.799128
971058.404128
973436.179279
981630.285398
985027.771304
..
.

70
71
72
73

1332026.80482
1348423.99041
1372139.51897
1379152.51378

resi
2.7 × 10−11
3.0 × 10−11
9.5 × 10−12
3.0 × 10−11
6.4 × 10−11..
.
7.5 × 10−12
4.7 × 10−13
8.4 × 10−12
1.1 × 10−12

Table 5: Results in Example 5.
Interval[−0.20, −0.15][−0.25, −0.15][−0.35, −0.15][−0.40, −0.15][−0.45, −0.15][−0.50, −0.15]

]ev
22
78
198
262
333
406

Est. ]ev
23.9
80.0
196.3
270.1
327.9
410.5

L
16
40
99
136
164
206

15

]reﬁnement
1
2
2
2
2
2

max(resi )2.8 × 10−13
2.1 × 10−12
8.5 × 10−12
1.7 × 10−12
9.0 × 10−12
9.4 × 10−12

Sakurai-Sugiura (SS) method is an eigensolver based on complex moments given by
the contour integrals of the matrix inverses with several shift points.
Some numerical properties of the method were presented from the view-point ofa ﬁlter for a subspace. According to the results, eﬃcient parameter estimation techniques were shown. The contour integral for a matrix inverse is regarded as a ﬁlter
for an eigensubspace. When the contour integral is approximated by a numericalquadrature, the quadrature error causes contamination of the eigencomponents corresponding to the eigenvalues located outside of the contour path. We demonstrated
the eﬃciency of our method with numerical experiments.
In the numerical experiments, we used a sparse direct solver. The use of iterative
linear solvers for multiple right-hand sides such as block Krylov subspace solvers areuseful because our eigensolver requires very small number of iterative reﬁnement.
The combination of block iterative linear solvers and performance evaluation withsuch solvers in practical problems remaining as future work. We are implementing asoftware package of the SS method for large-scale parallel computing. Applicationsof practical problems and performance evaluations of the presented method in a
massively parallel computing environment also remain as our future work.

ACKNOWLEDGMENT
This research was partially supported by the Core Research for Evolutional Scienceand Technology (CREST) Project of JST and the Grant-in-Aid for Scientiﬁc Research of Ministry of Education, Culture, Sports, Science and Technology, Japan,
Grant number: 23105702.

References[1] J. Asakura, T. Sakurai, H. Tadano, T. Ikegami and K. Kimura, A numerical method for nonlinear eigenvalue problems using contour integrals, JSIAM
Letters 1 (2009) 52–55.
[2] J. Asakura, T. Sakurai, H. Tadano, T. Ikegami and K. Kimura, A numericalmethod for polynomial eigenvalue problems using contour integral, Japan J.
Indust. Appl. Math. 27 (2010) 73–90.
[3] W. J. Beyn, An integral method for solving nonlinear eigenvalue problems, Lin.
Alg. Appl. 436 (2012) 3839–3863.
[4] Y. Futamura, H. Tadano and T. Sakurai, Parallel stochastic estimation methodof eigenvalue distribution, JSIAM Letters 2 (2011) 27–30.
[5] Y. Futamura, T. Sakurai, S. Furuya and J.-I. Iwata, Eﬃcient algorithm forlinear systems arising in solutions of eigenproblems and its application toelectronic-structure calculations, Proc. 10th International Meeting on HighPerformance Computing for Computational Science (VECPAR 2012) (accepted).

16

[6] 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. 233 (2010) 1927–1936.
[7] T. Ikegami and T. Sakurai, Contour integral eigensolver for non-Hermitiansystems: a Rayleigh-Ritz-type approach, Taiwanese J. Math. 14 (2010) 825–
837.
[8] Matrix Market, http://math.nist.gov/MatrixMarket/.
[9] Y. Maeda, Y. Futamura, T. Sakurai, Stochastic estimation method of eigenvalue
density for nonlinear eigenvalue problem on the complex plane, JSIAM Letters
3 (2011) 61–64.
[10] T. Mizusaki, K. Kaneko, M. Honma and T. Sakurai, Filter diagonalization ofshell-model calculations, Phys. Rev. C 82 (2010) 024310.
[11] M. Naito, H. Tadano, T. Sakurai, A modiﬁed block IDR(s) method for computing high accuracy solutions, JSIAM Letters 4 (2012) 25–28.
[12] H. Ohno, Y. Kuramashi, H. Tadano and T. Sakurai, A quadrature-based eigensolver with a Krylov subspace method for shifted linear systems for Hermitian
eigenproblems in lattice QCD, JSIAM Letters 2 (2010) 115–118.
[13] E. Polizzi, Density-matrix-based algorithms for solving eigenvalue problems,
Phys. Rev. B 79 (2009) 115112.
[14] T. Sakurai and H. Sugiura, A projection method for generalized eigenvalue
problems, J. Comput. Appl. Math. 159 (2003) 119–128.
[15] T. Sakurai and H. Tadano, CIRR: a Rayleigh-Ritz type method with contour
integral for generalized eigenvalue problems, Hokkaido Math. J. 36 (2007) 745–757.
[16] T. Sakurai, J. Asakura, H. Tadano and T. Ikegami, Error analysis for a matrix
pencil of Hankel matrices with perturbed complex moments, JSIAM Letters 1(2009) 76–79.
[17] T. Sakurai, H. Tadano and Y. Kuramashi, Application of block Krylov subspace algorithms to the Wilson-Dirac equation with multiple right-hand sides
in lattice QCD, Comput. Phys. Commun. 181 (2010) 113–117.
[18] H. Tadano, T. Sakurai and Y. Kuramashi, Block BiCGGR: A new block Krylov
subspace method for computing high accuracy solutions, JSIAM Letters 1(2009) 44–47.
[19] T. Watanabe, Y. Inadomi, H. Umeda, K. Fukuzawa, S. Tanaka, T. Nakano and
U. Nagashima, Fragment molecular orbital (FMO) and FMO-MO calculationsof DNA: accuracy validation of energy and interfragment interaction energy, J.
Comput. Theo. Nanoscience 6 (2009) 1328–1337.

17

[20] I. Yamazaki, H. Tadano, T. Sakurai and T. Ikegami, Performance comparison ofparallel eigensolvers based on a contour integral method and a Lanczos method,
Parallel Computing (in press).
[21] S. Yokota and T. Sakurai, A projection method for nonlinear eigenvalue problems using contour integrals, JSIAM Letters (accepted).

18