A filter diagonalization for generalized
eigenvalue problems based on the
Sakurai-Sugiura projection method
Tsutomu Ikegami a,∗ Tetsuya Sakurai b Umpei Nagashima ca Grid Technology Research Center, AIST, 1-1-1 Umezono, Tsukuba 305-8568,

JAPANb Department of Computer Science, University of Tsukuba, Tsukuba 305-8573,

JAPAN
c Research Institute for Computational Science, AIST, 1-1-1 Umezono, Tsukuba

305-8568, JAPAN

Abstract
The Sakurai-Sugiura projection method, which solves a generalized eigenvalue problem to find certain eigenvalues in a given domain, was reformulated by using theresolvent theory. A new interpretation based on the filter diagonalization was given,and the corresponding filter function was derived explicitly. The block version of
the method was also proposed, which enabled to resolve degenerated eigenvalues.
Two numerical examples were provided to illustrate the method.
Key words: Generalized eigenvalue problem, Spectral projection, Filter
diagonalization

1

Introduction

There is a wide area of scientific and engineering applications that are reduced
to the eigenvalue problem [3,11,2,6]. Generally speaking, the dimension of thesystem increases with the refinement of the corresponding physical model,
which can reach few millions. Fortunately, it is often the case that not all of
the eigenvalues are necessary, and those in the physically meaningful range
may suffice. Iterative methods are typically used for those cases [1,10,5,14],∗ Corresponding author.
Email address: t-ikegami@aist.go.jp (Tsutomu Ikegami).

Preprint submitted to Elsevier

3 December 2008

which are suitable to extract eigenvalues with large absolute value. In some
applications, however, it is required to extract eigenvalues embedded in the
middle of the distribution.
In the last decade, a filter diagonalization [4,7] is developed to extract eigenvalues contained in a specified region G. It was first formulated based on time
propagation of quantum wave packets [9], and soon was reformulated as afilter operator. Let H be a Hermitian matrix and v be an arbitrary vector.
P
A filter function f (x) = i ci xi is designed, such that f (x) ' 0 on outsideof G. The filter operator is defined by replacing x with H as f (H). Whenthe operator is applied on v, eigen-elements of H whose eigenvalues are contained in G dominate in f (H)v. Starting from a set of initial vectors vi , we
can thus span an eigen-subspace localized on G by using a set of f (H)vi . In
typical applications, a gaussian function is chosen for f (x), and Chebyshev
series expansion is used instead of the Taylor series. Because the power series
is employed in the expansion, however, it is inevitable for f (x) to diverge as
x → ±∞. Therefore, we have to shift and scale H a priori such that all the
eigenvalues of H are confined in the domain of f (x).
The Sakurai-Sugiura (SS) method [12] also extract eigenvalues in the specified
domain, though it takes a different approach. In the filter diagonalization, thedomain G is determined numerically by the asymptotic behavior of f (x). On
the other hand, the SS method defines G mathematically by a Cauchy integral
path that surrounds G. As a result, the SS method can deal with the nonHermitian systems, where the eigenvalues may be located on the complex
plane. The method is also applicable to the generalized eigenvalue problems,and is suitable to the modern distributed parallel computer.
In the present paper, we will reformulate the SS method in the context of theresolvent theory. The block version of the SS method is also proposed, in which
the relation to the filter diagonalization is discussed. Next section is devoted
to the derivation based on the resolvent theory. In Sec. 3, the SS method isinvestigated in the light of the filter diagonalization, and the algorithm of theblock SS method is given. Two numerical examples are shown in Sec. 4, andSec. 5 conclude.

2

Block Sakurai-Sugiura method

In this section, we reformulate the Sakurai-Sugiura (SS) projection method
based on the resolvent of a matrix pencil, and propose a block version of the
SS method. Let A, B ∈ CN ×N form a regular matrix pencil zB − A. The
regular pencil can be transformed into Weierstrass canonical form [15].

2Theorem 1 (Weierstrass canonical form). Let zB − A be a regular pencil oforder N . Then there exist nonsingular matrices P̃, Q ∈ CN ×N such thatzI − J1
 k1








P̃(zB − A)Q = 
















 , (1)








...
zIkd − Jd
zNd+1 − Ikd+1..

.
zNr − Ikr

where Ji , Ni ∈ Cki ×ki are Jordan blocks, Ni is nilpotent, and Ik denotes theidentity matrix of order k.
−1

Because P̃ and Q are the regular matrices, we can define P = P̃ and Q̃ =Q −1 . According to the Jordan block structures in Eq. (1), we will partitionrow vectors in P̃ and Q̃ into P̃i , Q̃i ∈ Cki ×N , and column vectors in P and Q
into Pi , Qi ∈ CN ×ki , respectively, for i = 1, 2, . . . , r.
Theorem 2. Resolvent of the regular pencil (zB − A)−1 is decomposed into(zB − A)−1 =

d
X
i=1





kX
r
m
i −1
i −1

kX

X
(Ji − αi Iki )
m
m
P̃
−
Q
Qi
z
N
P̃ , (2)
i
i
i


 i
(z − αi )m+1 
m=0

i=d+1

m=0

where αi is an eigenvalue of Jordan block Ji .
Proof. Let W = P̃(zB − A)Q. According to Theorem 1, we have(zB − A)−1 = QW −1 P̃
=

d
X

r
X

Qi (zIki − Ji )−1 P̃i +i=1

Qi (zNi − Iki )−1 P̃i .

(3)

i=d+1

Using the resolvent of the Jordan block,R(z, Ji ) ≡ (zIki − Ji )−1 =

kX
i −1

(Ji − αi Iki )m,m+1
m=0 (z − αi )(4)and (zNi − Iki )−1 = −z −1 R(z −1 , Ni ), we get the result.
Definition 3. Let Γ be a positively oriented closed Jordan curve and G be
the inside of Γ. For a non-negative integer n, the n-th order moment matrix
3

of the pencil zB − A, localized on G, is defined by
1 I n
Mn =
z (zB − A)−1 dzB.
2πi Γ

(5)

Note that the moment matrix can be defined more generally, where a polynomial of z appears in place of z n in Eq. (5). Such a matrix, however, can be
represented by a linear combination of Mn .
Theorem 4. The localized moment matrix is written asX

Mn =Qi Ji n Q̃i .

(6)i;αi ∈G

Proof. Let Eq. (2) to Eq. (5),

Mn =

d
X

½

Qi

i=1

ji (z) = znni (z) = z

n

¾

½

kX
i −1

(Ji − αi Iki )m,m+1
m=0 (z − αi )

kX
i −1

¾

r
X
1 I
1 I
ji (z)dz P̃i B −ni (z)dz P̃i B, (7)
Qi
2πi Γ
2πi Γ
i=d+1

(8)

z m Ni m .

(9)

m=0

Using the relation
n

z =

n
X

Ã !

j=0

we have

Ã

Res
z=α

n

z
(z − α)m+1

³ ´

!

n
(z − α)n−j αj ,
j

(10)



 0,
m > n,= ³´n
αn−m , m ≤ n,
n−m

(11)

where nj is a binomial coefficient and Res(g(z)) gives a residue of g(z) at az=αmpole z = α. Because (Ji − αi Iki ) = 0 for m ≥ ki , the residue of ji (z) is givenbymin(n,ki −1) ÃRes (ji (z)) =

z=αi

X

m=0

!

n
αin−m (Ji − αi Iki )mn−m= Ji n .

(12)

Because ji (z) is regular for z 6= αi and ni (z) is regular for z ∈ C, we have,using the Cauchy’s integral theorem,
Mn =Xi,αi ∈G4

Qi Ji n P̃i B.

(13)

Meanwhile, by definition, the matrix B is expanded as
B=

d
X

r
X

Pi Q̃i +i=1

Pi Ni Q̃i .

(14)

i=d+1

Combining Eqs. (13) and (14), we get the result.
As shown below, we can use the localized moment matrix to extract Jordan
blocks whose eigenvalues are contained in G. Here we employ the following
collective notations: Those Jordan blocks Ji ; αi ∈ G are collected to fromPkΓ × kΓ Jordan matrix JΓ , where kΓ = i;αi ∈G ki . Similarly, the corresponding
eigenvectors Qi and Q̃i are collected to form QΓ ∈ CN ×kΓ and Q̃Γ ∈ CkΓ ×N ,respectively. Let C and D be arbitrary N × m matrices, where N > m ≥ kΓ .
A size-reduced moment matrix is defined as
Mn = C H Mn D

∈ Cm×m .

(15)

Theorem 5. If ranks of both C H QΓ and Q̃Γ D are kΓ , non-singular part of amatrix pencil zM0 − M1 is equivalent to zIk − JΓ .
Γ

Proof. Using the collective notations,
Mn = C H QΓ JΓ n Q̃Γ D.

(16)

Because ranks of C H QΓ and Q̃Γ D are kΓ , there exist m × m regular matricesP̃m and Qm , such that




Ik 
P̃m C H QΓ =  Γ ,
0
³

(17)

´

Q̃Γ D Qm = Ik , 0 .

(18)

Γ

Then, we have

 zIkΓ − JΓ,P̃m (zM0 − M1 ) Qm = 

0

(19)

which is the result of the theorem.
From Theorem 5, the interior eigenvalues of the origial matrix pencil zB − A
are obtained by solving the size-reduced eigenvalue problem zM0 − M1 . The
corresponding eigenvectors are also obtained from the eigenvectors P̃m and Qmof the size-reduced system. Let Qm be split into (QΓ , Q∅ ), where QΓ ∈ Cm×kΓ
5defines the non-singular components of the right-eigenvectors of zM0 − M1 .
Similarly, P̃m is split into (P̃Γ , P̃∅ )T , P̃Γ ∈ CkΓ ×m . Note that P̃∅ and Q∅
define right and left common null spaces of Mn , respectively.
Theorem 6. The right-eigenvectors of the original matrix pencil zB − A isgiven by QΓ = M0 DQΓ , and its adjoint is given by Q̃Γ = P̃Γ C H M0 .
Proof. Because M0 = QΓ Q̃Γ , we have
M0 DQΓ = QΓ Q̃Γ DQΓ = QΓ Ik = QΓ .
Γ

(20)

Similarly, P̃Γ C H M0 = Q̃Γ .
The left-eigenvectors are obtained by solving on the transposed system.
We can now derive the original SS method. Let v ∈ CN be an arbitrary vector
and µn = v H Mn v. Two m × m Hankel matrices Hm and H<
m are defined as




 µ0


 µ1
Hm = 
 .
 ..



µ1 · · · µm−1 µ2 · · ·
..
.










µm
..
.

(21)

µm−1 µm · · · µ2m−2
and





 µ1



 µ2
<
Hm = 
 .
 ..



µ2 · · ·
µ3
..
.

µm 

· · · µm+1 

.
.. 
. 

(22)




µm µm+1 · · · µ2m−1
Theorem 7. If all elements of Q̃Γ v and v H QΓ are non-zero, and there isno degeneracy in JΓ , then non-singular part of a matrix pencil zHm − H<m isequivalent to zIkΓ − JΓ .
Proof. By choosing row vectors of C H and column vectors of D to be(C H )i,∗ = v H QΓ JΓ i−1 Q̃Γ(23)

D∗,i = QΓ JΓ i−1 Q̃Γ v

(24)

and
for i = 1, 2, . . . , m, respectively, we have Hm = M0 and H<m= M1 . As forthe rank of Q̃Γ D, we consider that column vectors of Q̃Γ D form the Krylov
6series of JΓ starting from Q̃Γ v. Because JΓ is not degenerated, and elementsof Q̃Γ v are non-zero, these column vectors are linearly independent, and thus
the rank of Q̃Γ D is kΓ . Similarly, the rank of C H QΓ is kΓ . From Theorem 5,
we have the result of the theorem.
Rigorously, not all of the elements of Q̃Γ v need to be non-zero; the prerequisite
non-zero elements are those corresponding to the last Jordan canonical vector
for each Jordan blocks in JΓ .
A block version of the SS method is located between Theorems 5 and 7. Let V
be an N ×l matrix, whose column vectors are used as initial vectors. The block
SS method is defined by replacing µn in Theorem 7 with matrices V H Mn V .
Apparently from the proof of the theorem, up to the l-th order degeneracy inJΓ can be separated in the block SS method.

3

Filter diagonalization

From the viewpoint of the filter diagonalization, the localized moment matrix
Mn is considered to be a filter operator. In the following, we focus on the casewhere the matrix pencil zB − A is diagonalizable. In this case, Ji is reduced toa scalar, αi . Comparing Equations (5) and (6), a filter function corresponding
to Mn is derived as

fn (x) =1 I zn
dz
2πi
Γ z−x


 xn ,x ∈ G,

 0,otherwise.

=

(25)

When Mn is operated on the i-th right-eigenvector, its amplitude is modulated
by fn (αi ),Q̃i Mn Qi = fn (αi )Q̃i Qi .
(26)
In the conventional filter diagonalization, a single filter operator f (H) is applied on a set of random vectors to construct a basis set that spans the filteredeigen-subspace. On the other hand, the SS method applies a set of modulated
filter operators Mn on a single vector to build the basis set. From this point ofview, the block SS method is considered naturally as a hybrid method between
the filter diagonalization and the SS method.
In the actual calculations, the path integral in Eq. (5) is evaluated numericallyby using an appropriate quadrature. Let zj be the quadrature points on Γ and
7

n=0

1
n=2

0.5
n=8

-1

-2

1

2

-0.5

n=1

G
-1

Fig. 1. Plots of filter functions f˜n (x) for M = 16, γ = 0, ρ = 1, n = 0, 1, 2, and 8.

wj be the corresponding weights, where j = 1, 2, . . . , M . The localized momentmatrix is then approximated by
M̄n =

M
X

wj zjn (zj B − A)−1 B.

(27)

j=1

Accordingly, the effective filter function f¯n (x) is derived for the quadrature asf¯n (x) =

M
X
j=1

wj

zjn
.
zj − x(28)

Thanks to the explicit form of the filter function, several properties of the SS
method can be analyzed quantitatively. As an example, let Γ be a circle ofradius ρ centered at a point γ. Dividing the circle into M division and applying
the trapezoidal rule, we have zj = γ + ρ exp( 2πi(j − 21 )) and wj = (zj − γ)/M .
M
From the numerical reason, the momental weight of z n is also replaced by ashift-and-scaled one, ((z − γ)/ρ)n . Then, the filter function becomesf˜n (x) =

M
X

Ã

wj

j=1

=

zj − γρ!n1zj − x(29)x̃n
1 + x̃M(30)

where x̃ = (x − γ)/ρ. In Fig. 1, f˜n (x) is plotted along the real axis for M = 16,
γ = 0, ρ = 1, n = 0, 1, 2, and 8. The shape of the effective filter functions
deviates from the ideal rectangle, which exude outside of G. The value and thederivative at the boundary are f˜n (γ+ρ) = 1/2 and f˜n0 (γ+ρ) = (−M/4+n/2)/ρ
for n = 0, 1, . . . , M −1. From here, we notice that n should be kept well smaller
than M/2.
The filtered subspace is spanned by the set of filtered vectors. If there exist
multiple eigenvalues in G, we need as many vectors to span the subspace.
Higher order moment matrices become necessary then, which may, in turn,8require large M . Instead, we can prepare multiple initial vectors and use theblock SS method. The algorithm of the block SS method is shown below.
Algorithm 1 (Block SS method).
Input: V ∈ CN ×l , {zj , wj } for j = 1, 2, . . . , M
Output: αk for k = 1, 2, . . . , K
1. Solve Ṽj = (zj B − A)−1 BV and calculate Vj = V H Ṽj ∈ Cl×lPn
2. Compute µ̄n = M
j=1 wj zj Vj
ml×ml
3. Construct Hankel matrices Hm and H<m ∈ C
4. Perform singular value decomposition, Hm = WsU
H −1/2
5. Construct H = s−1/2 WH H<∈ CK×K
mU s
6. Compute eigenvalues of H to have αkIf eigenvectors Qk are also wanted, let qk be eigenvectors of H,Pn
7. Compute S̄n = M
j=1 wj zj Ṽj8. Compute (Q1 , . . . , QK ) = (S̄0 , . . . , S̄m−1 )UH s−1/2 (q1 , . . . , qK )In the block SS method, the filtered subspace is implicitly spanned by the ml
column vectors of S̄n . At steps 2 and 7, the momental weight zjn may as well
be replaced by the shifted-and-scaled one. Note that the resulting eigenvalues
are also shifted-and-scaled. At step 4, small singular value components areomitted, so that s ∈ RK×K where K < ml. If no ignorable singular values arefound, there may be more eigenvalues than ml, and we should increase either
m or l to span the subspace in higher dimensions. Because the effective filter
function is not perfect, eigenvalues outside of G may be contaminated to be
K ≥ kΓ .

4

Numerical examples

Two numerical examples are given in this section. The first one depicts how
the localized moment matrix works as a filter operator. The other one demonstrates the improvement achieved by the block version of the SS method. The
algorithm was implemented on Mathematica 5.2. The linear equation solver
LinearSolve[] was used to evaluate (zj B − A)−1 BV numerically. The calculated eigenvalues are compared against those obtained by Eigensystem[].
Example 8. Two random matrices, A, B ∈ R100×100 , were generated. The
elements of A were taken randomly in [−1, 1], while B was diagonal dominant,where random noises of [−0.1, 0.1] were added to every entry of an identity
matrix. The initial vector was setup artificially as a sum of all the rightP
eigenvectors, v = i Qi . We define the i-th eigen-amplitude in v by Q̃i v, so
9~(a) Qv~–(b) QM0v

1

1
5

0.5
0

0

-5
0
Re

5

0.5
0

Im

0

-5
0
Re

-5
5

Im

-5
5

Fig. 2. Example 8: the eigen-amplitudes of a vector (a) before and (b) after theoperation of the localized moment matrix.
(a)

(b)

2
α8

10

~ –
Amplitude, | Qi Mnv |

α6
1

Im

Γ
0

α4

-1

-2

α3
α1

α5

α2

10

10

10

-1

0
Re

1

2

α4

α2,3

-2

α5,6

-4

-6

α1

α7
10

-2

0

α7,8

-8

0

5

10

15

Order of momentum matrix, nFig. 3. Example 8: (a) distribution of the eigenvalues around the origin, and (b)eigen-amplitudes of M̄n v for α1∼8 .

that all the amplitudes in the initial vector were unity. The integral path Γ
was taken to be a unit circle centered at origin. The M -point trapezoidal rule
was employed for the quadrature, where M = 64, γ = 0.0, and ρ = 1.0.
The eigen-amplitudes of v and M̄0 v are plotted in Fig. 2. Each box is located at
the eigenvalue point on the complex plane, and the real part of the amplitudeis plotted as its height; the imaginary part was negligibly small. As can beseen from the figure, four primary eigen-elements remained in M̄0 v, which arelocated inside of Γ. Detailed location of these eigenvalues, α1∼4 , are plotted inFig. 3(a), along with four more eigenvalues, α5∼8 , located on the periphery of
Γ. The eigen-amplitudes are modulated as we change the order of the momentmatrix n. In Fig. 3(b), absolute values of the amplitudes are plotted as afunction of n. The amplitudes of other eigen-elements than α1∼8 were lower
than the plotted area. For the primary eigen-elements, α1∼4 , the amplitude
became smaller as n increases, which is roughly proportional to |αi |n . On theother hand, the amplitude of the contaminant eigen-elements, α5∼8 , increased
with n. It is because the tail of the filter function f˜n (x) extends more outward
for larger n. In the present example, f˜15 (x) becomes lower than the thresholdof 10−14 at x > 1.9 along the real axis.
1016×16
To obtain the eigenvalues, Hm , H<
were calculated. The singularm ∈ Cvalue decomposition was performed on Hm to extract the non-singular partfrom the size-reduced matrix pencil. The singular values were compared withthe largest one, and those components relatively smaller than 1.1 × 10−14 were
removed. Two component were rejected at this step, and the approximatedeigenvalues were calculated with the rest of 14 components. The accuracy ofthus obtained eigenvalues are not uniform. The maximum error in the primaryeigenvalues α1∼4 was 1.7 × 10−14 . The corresponding eigenvectors are mainly
spanned by the large singular value components, which are calculated in highaccuracy. The smaller components also participate, though the contributions
are so small that their inferior accuracy doesn’t affect. Next two eigenvalues, α5and α6 , are located near the periphery of G, the error of which was 2.4 × 10−12 .
The errors in the following 7 eigenvalues were 3.2 × 10−8 ∼ 1.5 × 10−6 . Theseeigenvectors are spanned by the smaller singular value components. The rest
eigenvalue could not be assigned to any exact correspondents, the eigenvectorof which is spanned by the trace singular value components.

The original SS method has a difficulty to resolve nearly degenerated eigenvalues congested around the center of Γ. It is because the filter functions of largen are flat at the center, which makes the filtered vectors linearly dependent.
The next example shows that the difficulty can be removed by the block SS
method.
Example 9. A real symmetric matrix A ∈ R400×400 was prepared, which has
five primary eigenvalues, -10.03, -10.02, -10.01, -10.00, and -9.99, in the rangeof [−10.5, −9.5]. Other eigenvalues were taken randomly in the range [−40, 40],
and a random unitary matrix was prepared to construct A. An identity matrix was used for B. The circular integral path and the trapezoidal rule was
employed, with the center and radius as γ = −10 and ρ = 0.5, respectively.
The shifted-and-scaled momental weight was employed as described in Sec. 3.
The block SS method was examined on this system with varying the numberof initial vectors l. Those initial vectors of R400 were generated randomly
and ortho-normalized. The result is summarized in Table 1. To maintain thecomputational amount nearly the same, the number of trapezoidal points M
was reduced as l increased. The order of the Hankel matrix m was taken
to be large enough for the size-reduced matrix pencil z−γ
Hm − H<m to haveρa singular part. This process certifies the over-completeness of the basis set
that spans the filtered subspace. The linear dependence of the basis set wasremoved by the singular value decomposition on Hm , where the threshold of
the singularity was taken to be 1.1 × 10−14 . For l = 1 ∼ 3, a total of 12 vectors
were (implicitly) generated, while 16 vectors were necessary for l = 4, due tothe duller filter function with M = 32.
Starting from a single initial vector, only four dimensions out of five could11Table 1
Example 9: errors of eigenvalues calculated by the block SS method. The columns,
from left to right, are the number of initial vectors, the number of trapezoidal points,
the order of the Hankel matrix, the number of obtained eigenvalues, errors of theeigenvalues in G, and errors of the contaminant eigenvalues, respectively.
l

M

m

Nα

|∆α∈G |

|∆α∈G
/ |

1

128

12

6

1.3 × 10−6 ∼ 4.1 × 10−3#

8.0 × 10−11 ∼ 6.9 × 10−9

2

64

6

8

7.7 × 10−12 ∼ 9.5 × 10−8

3.0 × 10−12 ∼ 6.5 × 10−5

3

42

4

10

4.5 × 10−13 ∼ 1.1 × 10−10

8.1 × 10−13 ∼ 2.3 × 10−2

4

32

4

14

1.8 × 10−15 ∼ 1.5 × 10−12

4.8 × 10−14 ∼ 8.1 × 10−2

# only four eigenvalues out of five were obtained.

be extracted for the primary subspace. Because the primary eigenvalues arelocated near the origin, the corresponding eigen-elements diminished quickly in
Mn v, as n ≥ 1 increased. As a result, peripheral eigen-elements overwhelmedin Mn v, leaving the basis set incomplete in the primary subspace. The situationis the same for l > 1, but in these cases the column vectors of M0 V alreadyspan l-dimensions with moderate linear independency. The eigenvalues became
more accurate as l increases, reflecting the improved linear independency inthe primary subspace. With the large l, the order of the moment m can also
be kept smaller, where the practical filter functions are expected to behave
better (see Fig. 1).
In this example, the number of linear equations to be solved were kept nearly
constant to average the computation cost among different l. It is often thecase, however, that the cost to solve a linear equation A(x1 , x2 ) = (v1 , v2 ) ismuch smaller than twice the cost of Ax = v. This observation makes the block
version of the SS method preferable, especially when the size of the matrix
is large and preconditioning is necessary. The block SS method also inherits
the adaptability to the parallel computing in the original, which makes themethod attractive from a computational point of view.

5

Conclusion

The SS method was reformulated based on the resolvent theory, and its block
version was proposed. The numerical example indicated that the block version
has a potential not only to resolve degeneracies, but also to achieve higher
accuracy under the same computational cost. The method was also interpreted
in the context of the filter diagonalization. The explicit form of the filter
function was derived for the SS method, which helps to investigate several
properties of the method quantitatively.
12

In the Rayleigh-Ritz approach of the SS method [13], we don’t use the Hankel
matrices, but construct a size-reduced matrix pencil directly in the filteredsubspace. In this approach, the quadrature for the path integral needs not
be accurate enough, as far as the corresponding filter function f¯n (x) is well
localized. More specifically, f¯0 (x) need not be flat in x ∈ G, as far as f¯0 (x) '0 for x ∈/ G. Indeed, Murakami proposed a novel filter operator similar toEq. (27) [8], where the momental weight is replaced by constants and {zj , wj }are chosen specifically. The corresponding filter function does oscillate in G,
and falls off quickly outside of G. Liberation from the accuracy allows us tochoose the integral path and quadrature points more freely, and to tweak thequadrature weight to make the filter function tail off faster. For example, if
the quadrature points can be arranged on an even grid, the survey of theeigenvalue distribution becomes much easier, because we can reuse most of
the quadrature points for different integral paths. Designing such an optimized
quadrature is a part of our future work.

Acknowledgements
The author TI thanks Dr. H. Murakami and Dr. R. J. Harrison for the fruitful
discussions and valuable suggestions.

References[1] W. E. Arnoldi, The principle of minimized iteration in the solution of the matrix
eigenproblem, Q. Appl. Math. 9 (1951) 17–29.
[2] Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, H. van der Vorst, Templates for
the Solution of Algebraic Eigenvalue Problems: a Practical Guide, SIAM,
Philadelphia, 2000.
[3] F. Chatelin, Eigenvalues of Matrices, John Wiley & Sons, 1993.
[4] R. Chen, H. Guo, A general and efficient filter-diagonalization method without
time propagation, J. Chem. Phys. 105 (1996) 1311–1317.
[5] E. R. Davidson, The iterative calculation of a few of the lowest eigenvalues andcorresponding eigenvectors of large real-symmetric matrices, J. Comp. Phys. 17(1975) 87–94.
[6] G. H. Golub, H. van der Vorst, Eigenvalue computation in the 20th century, J.
Comp. Appl. Math. 123 (2000) 35–65.
[7] V. A. Mandelshtam, H. S. Taylor, A low-strage filter diagonalization methodfor quantum eigenenergy calculation or for spectral analysis of time signals, J.
Chem. Phys. 106 (1996) 5085–5090.

13

[8] H. Murakami, in: Proceedings of High Performance Computing Symposium,
Tsukuba, 2007.
[9] D. Neuhauser, Bound state eigenfunctions from wave packets: Time→energy
resolution, J. Chem. Phys. 93 (1990) 2611–2616.
[10] C. Paige, Computational variants of the Lanczos method for the eigenproblems,
J. Inst. Math. Appl. 10 (1972) 373–381.
[11] Y. Saad, Numerical Methods for Large Eigenvalue Problems, Halstead Press,
Manchester, 1992.
[12] T. Sakurai, H. Sugiura, A projection method for generalized eigenvalue
problems using numerical integration, J. Comp. Appl. Math. 159 (2003) 119–128.
[13] T. Sakurai, H. Tadano, CIRR: a Rayleigh-Ritz type method with contour
integral for generalized eigenvalue problems, Hokkaido Math. J. 37 (2007), in
print.
[14] G. Sleijpen, H. van der Vorst, A Jacobi-Davidson iteration method for linear
eigenvalue problems, SIAM J. Matrix Anal. Appl. 17 (1996) 401–425.
[15] H. W. Turnbull, A. C. Aitken, An Introduction to the Theory of CanonicalMatrices, Blackie & Sons, London, 1952.

14