A numerical method for polynomial eigenvalue
problems using contour integral
Junko Asakuraa Tetsuya Sakuraib Hiroto Tadanob
Tsutomu Ikegamic Kinji KimuradaGraduate School of Systems and Information Engineering, University of TsukubabDepartment of Computer Science, University of TsukubacdGrid Technology Research Center, AIST

Department of Applied Mathematics and Physics, Kyoto UniversityaE-mail: asakura@mma.cs.tsukuba.ac.jp

Abstract
We propose a numerical method using contour integral to solve polynomial eigenvalue problems (PEPs). The method ﬁnds eigenvalues contained in a certain domain
which is deﬁned by a surrounding integral path. By evaluating the contour integral
numerically along the path, the method reduces the original PEP into a small generalized eigenvalue problem, which has the identical eigenvalues in the domain. Error
analysis indicates that the error of the eigenvalues is not uniform: inner eigenvalues
are less erroneous. Four numerical examples are presented, which conﬁrm the theoretical predictions.
Keywords: polynomial eigenvalue problem, matrix polynomial, contour integral,projection method

1

Introduction

The present paper is concerned with a numerical method to solve polynomial eigenvalue problems. The polynomial eigenvalue problem (PEP) [1,4,7] involves ﬁnding aneigenvalue∑λl andi corresponding nonzero eigenvector x that satisfy F (λ)x = 0, whereF (λ) = i=0 λ Ai with real or complex coeﬃcient matrices. The eigenvalue problem of the matrix polynomial is considered to be a generalization of linear eigenvalue
problems.
Polynomial eigenvalue problems can be used in a variety of problems in scienceand engineering. For example, quadratic eigenvalue problems arise in oscillation
1

analysis with damping [9, 18] and stability problems in ﬂuid dynamics [8], and thethree-dimensional (3D) Schrödinger equation can result in a cubic eigenvalue problem[10]. Similarly, the study of higher-order systems of diﬀerential equations leads to amatrix polynomial of degree greater than one [5]. However, its applications are more
complicated than standard and generalized eigenvalue problems.
One reason is in the diﬃculty in solving the PEPs. Polynomial eigenvalue problems are typically solved by linearization [6, 13, 15], which promotes the k-th ordern × n matrix polynomial into the larger kn × kn linear eigenvalue problem. Other
methods, such as Arnoldi shift and invert strategy [2], can be used when several eigenvalues are desired. A disadvantage of the shift-and-invert Arnoldi methods is that achange of the shift parameter requires a new Krylov subspace to be built. Anotherapproach is a direct solution obtained by means of the Jacobi-Davidson method [17],although this method has been investigated far less extensively.
We herein propose an eigensolver, which prevents the inﬂation of the matrix dimension. Indeed, we can even reduce the dimension of the problem by focusing on
only the eigenvalues of physical interest. For linear problems, the Sakurai-Sugiura(SS) method [16] can ﬁnd certain eigenvalues in a given domain, is capable to solve
non-Hermitian systems, and is suitable to modern distributed parallel computers. In
the present paper, we propose a solution of polynomial eigenvalue problems by meansof the SS method. The proposed method enables us to obtain the eigenvalues of thematrix polynomial by solving the generalized eigenvalue problem, which is derived bysolving systems of linear equations. Since these linear systems are independent foreach equation, they can be solved in parallel. The proposed method is discussed fromthe theoretical point of view. In order to give the theoretical proof of the present SS
method, we have redeﬁned PEP in terms of the Smith form which is sometimes called
the canonical form for polynomial matrix. In addition, some results of numerical examples are reported to conﬁrm our theoretical observations.
The remainder of the paper is organized as follows. In the next section, we introduce the Smith form for matrix polynomials, which is sometimes called a canonicalform, and its related results. In Section 3, we present the numerical method forsolving the PEP by means of the SS method and present theoretical results for the
proposed method. In Section 4, we present the algorithm of the SS method for thecase in which the domain is given by a circle and the numerical integration is evaluated via the trapezoidal rule. We also provide the preliminary error analysis for theapproximated eigenvalues of matrix polynomials evaluated by the SS method. Somenumerical examples are shown in Section 5. Finally, conclusions and suggestions for
future study are presented in Section 6.

2

Matrix polynomials

In this section, we ﬁrst deﬁne the term matrix polynomial. Then, we introduce theSmith form for matrix polynomials, which is sometimes called a canonical form, andits related results.

2

Deﬁnition 2.1 (matrix polynomial) If A0 , A1 , . . . , Al ∈ Cn×n and Al 6= O, thenthe matrix-valued function deﬁned on C by
F (z) =l∑z i Ai ,z∈C

i=0

is called a matrix polynomial of degree l.
The matrix polynomial F (z) is called regular when the determinant of F (z) is not
identically zero for all values of z.
Deﬁnition 2.2 (vector polynomial) If u0 , u1 , . . . , ul ∈ Cn and ul 6= 0, then thevector-valued function deﬁned on C byv(z) =l∑z i ui ,z∈C

i=0

is called a vector polynomial of degree l.
Theorem 2.3 (the Smith form [7]) Let F (z) be an n × n regular matrix polynomial. Then, F (z) admits the representation
P (z)F (z)Q(z) = D(z),
where




D(z) = 

d1 (z)

0

(1)

···
...

0
..
.

0

0
dn (z)0..
.

d2 (z)... ...

0

···







is a diagonal polynomial matrix with monic scalar polynomials dj (z) such that dj (z)are divisible by dj−1 (z) for j = 2, ..., n. In addition, P (z) and Q(z) are n × n matrix
polynomials with constant nonzero determinants.
Representation (1), as well as the diagonal matrix D(z), is called the Smith form of
matrix polynomial F (z). D(z) is deﬁned uniquely.
The eigenpairs of the PEP F (λ)x = 0 is formally derived from the Smith form.
Let qj (z) be the column vectors of Q(z),Q(z) = (q1 (z) . . . qn (z)),(2)

P (z)H = (p1 (z) . . . pn (z)).

(3)

and pj (z) be
Let λi be a set of zero points of dn (z), where i = 1, . . . , ˜l ≤ ln. Because dj (z) is
divisible by dj−1 (z), the invariant polynomial dj (z) is represented in terms of λi asl̃∏
dj (z) =(z − λi )αji ,j = 1, . . . , n,i=1where αji ∈ Z+ (non-negative integer) and αji ≤ αj 0 i for j < j 0 (see [7, p. 320]). The
eigenpairs of the PEP is deﬁned by the the zero-points of dj (z) and qj (z).
3Lemma 2.4 Let qj (z) be a vector polynomial of (2) , and λi be a zero point of dj (z).
Then,F (λi )qj (λi ) = 0,
that is, qj (λi ) is the eigenvector of F (z) corresponding to eigenvalue λi .
Proof. Because P (z) and Q(z) are invertible,
F (λi )qj (λi ) = P (λi )−1 D(λi )Q(λi )−1 (Q(λi )ej )= P (λi )−1 D(λi )ej= dj (λi )P (λi )−1 ej .
2

Since dj (λi ) = 0, we have the result of the lemma.

Note that if the eigenvalue λi is simple and not degenerated, i.e., λi is the simple
zero of det F (z), the invariant polynomials satisfy αji = 0 for j 6= n and αni = 1.

3

An eigensolver using contour integral

In this section, we propose an eigensolver using contour integral. This method enables us to obtain the eigenvalues of the matrix polynomial by solving a generalized
eigenvalue problem derived in this section.
Let F (z) be an n × n regular matrix polynomial of degree l, and let λ1 , . . . , λl̃ beﬁnite eigenvalues of the matrix polynomial F (z). For nonzero vector u, v, we deﬁnef (z) := uH F (z)−1 v.
The existence of the Smith form allows us to prove the following.
Theorem 3.1 Let D(z) = diag(d1 (λ), . . . , dn (λ)) be the Smith form for F (z), and
let P (z) and Q(z) be deﬁned by (1). Let χj (z) = uH qj (z)pj (z)H v, 1 ≤ j ≤ n. Then,f (z) admits the representationf (z) =n∑χj (z)

d (z)
j=1 j

.

Proof. It follows from (1) thatf (z) = uH F (z)−1 v
 −1d1 (z)
0
···
0
..
...


0
d−1
.


H
2 (z)
= u Q(z)  P (z)v
..
...
...


0
.
0
···
0 d−1
n (z)
n
∑ uH qj (z)pj (z)H v=
dj (z)j=1=n∑
χj (z)j=1

dj (z)

.
4

(4)

2

This proves the theorem.

Let Γ be a positively oriented closed Jordan curve in the complex plane. Let
λ1 , . . . , λm be distinct eigenvalues that lie in the interior of Γ. Assume that these
eigenvalues are simple and not degenerated. Then we can suppose that αjl = 0 and
αnl = 1 for j = 1, . . . , n − 1 and l = 1, . . . m.
Deﬁnition 3.2 For a non-negative integer k, the moment µk is deﬁned as∫
1
z k f (z)dz, k = 0, 1, ....
µk :=
2πi Γ

(5)

<
and Hm are deﬁned as
Deﬁnition 3.3 Two m × m Hankel matrices Hmµ0µ1 · · · µm−1
 µ1
µ2 · · ·µm Hm :=  ..

..
..
 .

.
.
µm−1 µm · · · µ2m−2



and

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

µ2
µ3
..
.

···
···

µm µm+1 · · ·µmµm+1 

.
..

.
µ2m−1

The following theorem is one of the main result of the present paper.
<− λHm are given by λ1 , . . . , λm , if
Theorem 3.4 The eigenvalues of the pencil HmHHχn (λl ) = u qn (λl )pn (λl ) v 6= 0 for 1 ≤ l ≤ m.

Proof. Let D(z) = diag(d1 (z), . . . , dn (z)) be the Smith form of F (z). Then, from
Theorem 3.1, we have∫1µk =z k f (z)dz
2πi Γ
∫
n
∑
1
χj (z) k
=
z dz
2πi
d
(z)
j
Γ
j=1
(
)
m
∑
χn (z) k=Resz , λl
d
n (z)
l=1
=

m
∑
l=1

=

m
∑

∏

χn (λl )
∏
λk
αni l
(λ
−
λ
)
i6=l (λl − λi )
l
i
i>m

i≤m

νl λkl .

l=15Here, νl are given as follows:νl = ∏χn (λl )∏,αnii6=l (λl − λi )i>m (λl − λi )i≤mfor l = 1, . . . , m.
Let Vm be the Vandermonde matrix

1
1
···
 λ1
λ
···
2

Vm :=  ..
..
 .
.
m−1
m−1
λ1
λ2
···

1
λm
..
.




.


λm−1
m

Let Dm := diag(ν1 , . . . , νm ), Λm := diag(λ1 , . . . , λm ). One can easily verify that
Hm = Vm Dm VmT ,<= Vm Dm Λm VmT ,Hmtherefore,<− λl Hm = Vm Dm (Λm − λl I )VmT .
Hm(6)

If χn (λl ) 6= 0 for 1 ≤ l ≤ m，then νl 6= 0 for 1 ≤ l ≤ m. Since λ1 , . . . , λm are
distinct, Vm is nonsingular. Since Dm is also nonsingular, we have the result of thetheorem.
2
Therefore, we can obtain the eigenvalues λ1 , . . . , λm of the matrix polynomial F (z)<x = λHm x.
by solving the generalized eigenvalue problem HmNext, we consider a method to evaluate eigenvectors. For the nonzero vectorv ∈ Cn , let
∫
1
sk :=
z k F (z)−1 vdz, k = 0, 1, . . . ,2πi Γand let
S := [s0 , . . . , sm−1 ].
Note that µk = uH sk . Then, we have the following theorem.
<− λHm and let
Theorem 3.5 Let (λl , wl ) be the eigenpairs of the pencil Hm

xl := Swl ,l = 1, . . . , m.

Then xl are the eigenvectors of F (z) corresponding to eigenvalue λl for l = 1, . . . , m.
For the proof of Theorem 3.5, we need the following lemma.
Lemma 3.6 Let qn (z) be a vector polynomial of (2) and let (λl , wl ) be the eigenpairs<− λHm . Then,of the pencil Hm
qn (λl ) = cl Swl ,for l = 1, . . . , m.
6

cl ∈ C\{0}Proof. From equation(6),<− λj Hm )wl
0 = (Hm= Vm Dm (Λm − λl I )VmT wl .

Since Vm and Dm are nonsingular, and Λm el = λl el , VmT wl admits the following
representation:(7)
VmT wl = αl el , αj ∈ C\{0}.
Here, el is the l-th unit vector. Let pn (z) be a vector polynomial of (3) and let
σl := ∏pn (λl )H v∏,
αni
i6=l (λl − λi )i>m (λj − λi )l = 1, . . . , m.

l≤mNote that σl 6= 0 if χn (λl ) 6= 0. As in the proof of Theorem 3.4, we can derive thefollowing equation.
S = [s0 , . . . , sm−1 ] = [σ1 qn (λ1 ), . . . , σm qn (λm )]VmT .
Therefore,qn (λl ) =1 1
1
SVm−T el = S wl = cl Swl .
σl
σl αl

Here, cl are
cl =

1
σl αl
2

for l = 1, . . . , m.
From Lemma 2.4 and Lemma 3.6, we obtain
0 = F (λl )qn (λl ) = F (λl )(cl Swl )for l = 1, . . . , m. Hence xl = cl Sxl is the eigenvector of F (z) corresponding to λl .
This proves Theorem 3.5.

4Case in which Γ is given by a circle

In actual applications, the contour integrals of equation (5) is evaluated numerically.
In this section, we will introduce a detailed algorithm where a trapezoidal rule isemployed on a circular integral path.
Let γ and ρ be the center and the radius, respectively, of the given circle. Let N
be a positive integer, and let
ωj = γ + ρ exp(2πi(j + 1/2)/N ),7

j = 0, 1, . . . , N − 1.

The integral of (5) is evaluated via the N -point trapezoidal rule, and the momental
weight is shifted and scaled based on a numerical consideration. Therefore, we obtainthe following approximations for µk .
N −1 (

1 ∑µk ≈ µ̂k :=N j=0

ωj − γρ)k+1
f (ωj ),

k = 0, 1, . . . .

Note that the resulting eigenvalues are also shifted and scaled.
To calculate the eigenvectors, letyj := F (ωj )−1 v,

j = 0, 1, . . . , N − 1,and approximate sk as
N −1 (

1 ∑sk ≈ ŝk :=N j=0

ωj − γ
ρ

)k+1
yj ,

k = 0, 1, . . . .

(8)Note that µ̂k = uH ŝk , becausef (ωj ) = uH F (ωj )−1 v = uH yj .
<andBecause of the numerical approximation, the derived Hankel matrices Hm
Hm contain numerical errors. Let us start from a generalized eigenvalue problem
Aw = λBw, where A and B are perturbed by εF and εG, respectively, wherekF k2 ≤ 1, kGk2 ≤ 1 and |ε| ¿ 1. The perturbed problem becomes(A + εF )w(ε) = λ(ε)(B + εG)w(ε).
Again, we assume λ to be a simple eigenvalue.
Let vector u be the corresponding left eigenvectors. Then, the following estimation
is given in [12].
uH (F − λG)w
λ(ε) − λ = ε+ O(ε2 ).
(9)uH Bw
The ﬁrst term of the right-hand side of equation (9) is evaluated asεkuk2 kwk2|uH (F − λG)w|≤ ε(1 + |λ|) H.
H|u Bw||u Bw|(10)

From equation (9) and equation (10), we obtain|λ(ε) − λ| ≤ ε(1 + |λ|)kuk2 kwk2+ |O(ε2 )|.
H|u Bw|(11)

From equation (11), we can deﬁne the following sensitivity factor for each eigenvalue.

8

Deﬁnition 4.1 Let λi be an eigenvalue of the matrix pencil A − λB, and let xi andwi be the corresponding right and left eigenvectors. Then we deﬁne the sensitivity
factor of the eigenvalue λi as
τi := (1 + |λi |)kui k2 kwi k2,|uH
i Bwi |

i = 1, 2, . . . .

(12)

<
and
Now, we will investigate the sensitivity factor in the SS method. Let A = Hm<
B = Hm . From equation (7), the right eigenvector of the matrix pencil Hm − λHm
becomes
wi = αi Vm−T ei

for the eigenvalue λi . Since the pencil is symmetric, the corresponding left eigenvector
is ui = w̄i Then we can derive the following equation
kui k2 kwi k2
kVm−T ei k22
=
−T
T
−1
|uH
|eT
i Hm w i |i Vm (Vm Dm Vm )Vm ei |kVm−T ei k22
1=
kV −T ei k22 ,=T|νi | m|ei Dm ei |where Dm = diag(ν1 , . . . , νm ).
<− λHm asTherefore, we obtain the sensitivity factor of the eigenvalues of Hmτi =1
kV −T ei k22 (1 + |λi |).
|νi | m(13)

When µk is approximated numerically, eigenvalues on the periphery of Γ start tocontaminate. Therefore, we usually take the dimension of the Hankel matrix m̂ larger
than m. Assuming that those peripheral eigenvalues are also simple, Then, we can
write
m̂
∑
µ̂k =
ν̂i λki + O(ε).
i=1Here ν̂i for i ≤ m, ν̂i for i > m, λ̂i and ε approach to νi , 0, λi and 0 respectively, as theintegrator is improved. Indeed, in the case of the N -point trapezoidal rule applied on
the circular path treated here, we can show that
ν̂i =

ci
1
·
,
ρ 1 + ηiN

(14)

where ηi = (λi − γ)/ρ and ci is a constant determined by F (z). From equation(13), we can see that the inner eigenvalues λi , i = 1, . . . , m are less sensitive to theperturbation, while the peripheral eigenvalues may ﬂuctuate largely.
Note that Hankel matrices tend to be very ill-conditioned [19]. Also in the present
case, the condition number of the Hankel matrix can be enormous because of the
trace |ν̂i | for i > m. As shown in the sensitivity factor analysis, however, the innereigenvalues can be calculated accurately.
The algorithm of the SS method is shown below.
9

Algorithm 1 (SS method)
Input: u, v ∈ Cn , N, K, δ, γ, ρ
Output: λ̂1 , . . . , λ̂m̂ , x̂1 , . . . , x̂m̂
1. Set ωj ← γ + ρ exp(2πi(j + 1/2)/N ),−1

2. Compute yj = F (ωj ) v,

j = 0, . . . , N − 1

j = 0, . . . , N − 1

3. Compute ŝk , k = 0, . . . , 2K − 1 by (8)4. Form µ̂k = uH ŝk , k = 0, . . . , 2K − 1<
K
5. Construct ĤK = [µ̂i+j−2 ]Ki,j=1 and ĤK = [µ̂i+j−1 ]i,j=1

6. Perform singular value decomposition of ĤK
7. Omit small singular value components s < δ so that<<Ĥm̂ = ĤK (1 : m̂, 1 : m̂), Ĥm̂= ĤK
(1 : m̂, 1 : m̂),where m̂ ≤ K<
8. Compute the eigenpairs (ζ1 , w1 ), . . . , (ζm̂ , wm̂ ) of the pencil Ĥm̂
− λĤm̂9. Construct S = [ŝ0 , . . . , ŝm̂−1 ]10. Form x̂j = Swj , j = 1, . . . , m̂
11. Set λ̂j ← γ + ρ ζj ,j = 1, . . . , m̂In practice, we choose vectors u and v as random vectors. The dimension of the
Hankel matrices m̂ is determined from the singular value decomposition of HK withthe singular value threshold of δ. Therefore, K must be taken large enough thatK ≥ m̂, where m̂ is the number of inner and peripheral eigenvalues.
The block version of the SS method for generalized eigenvalue problems was proposed in [11]. The numerical examples in [11] indicate that the block SS method has
the potential to achieve higher accuracy.
Let U and V be n × L matrices, the column vectors of which are linearly independent and used as initial vectors. The block SS method is deﬁned by replacing f (z)in equation (5) with matrices U H F (z)−1 V . For a non-negative integer k, the matrix
Mk is deﬁned as∫1
z k U H F (z)−1 V dz, k = 0, 1, . . . .
(15)
Mk :=
2πi Γ
<
Then, m̃L × m̃L block Hankel matrices Hm̃Land Hm̃L are deﬁned asM0M1 · · · Mm̃−1
 M1
M2 · · ·Mm̃ Hm̃L := 

..
..
..


.
.
.
Mm̃−1 MK · · · M2m̃−2

10



and

M1
 M2

<
:=  ..
Hm̃L
 .

···
···

M2
M3
..
.

Mm̃ Mm̃+1 · · ·
Mm̃
Mm̃+1 .
..
.
M2m̃−1

Here, the order of the block Hankel matrices m̃ can be smaller than m, where m̃L ≥ m.
Let
∫
1
Sk :=
z k F (z)−1 V dz, k = 0, 1, . . . ,2πi Γand let S := [S0 , ..., Sm̃−1 ].
By approximating the integral of equation (15) via the N -point trapezoidal rule,
we obtain the following approximations for Mk .
N −1 (

1 ∑Mk ≈ M̂k :=N j=0

ωj − γρ)k+1
U H F (ωj )−1 V,

j = 0, 1, . . . .

Similarly, we obtain the follow approximations for Sk .
N −1 (

1 ∑Sk ≈ Ŝk :=N j=0

ωj − γρ)k+1
F (ωj )−1 V,

j = 0, 1, . . . .

(16)

Note that Mk = U H Sk . Therefore, the algorithm of the block SS method is shown
below.

Algorithm 2 (Block SS method)
Input: U, V ∈ Cn×L , N, K, L, δ, γ, ρ
Output: λ̂1 , . . . , λ̂m̂ , x̂1 , . . . , x̂m̂
1. Set ωj ← γ + ρ exp(2πi(j + 1/2)/N ),−1

2. Compute F (ωj ) V ,j = 0, . . . , N − 1

j = 0, . . . , N − 1

3. Compute Ŝk , k = 0, . . . , 2K − 1 by (16)
4. Form M̂k = U H Ŝk , k = 0, . . . , 2K − 1
<
KL×KL
K
5. Construct ĤKL = [M̂i+j−2 ]Ki,j=1 and ĤKL = [M̂i+j−1 ]i,j=1 ∈ C

6. Perform singular value decomposition of ĤKL
7. Omit small singular value components s < δ so that<<
= ĤKL(1 : m̂, 1 : m̂),Ĥm̂ = ĤKL (1 : m̂, 1 : m̂), Ĥm̂
where m̂ ≤ KL<
− λĤm̂8. Compute the eigenpairs (ζ1 , w1 ), . . . , (ζm̂ , wm̂ ) of the pencil Ĥm̂

9. Construct S = [Ŝ0 , . . . , Ŝm̂−1 ]10. Compute xj = Swj , j = 1, . . . , m̂
11. Set λ̂j ← γ + ρ ζj ,j = 1, . . . , m̂11

In practice, we choose matrices U and V as random matrices.
Note that Algorithm 2 with L = 1 is equivalent to Algorithm 1. At Step 2, we needto solve the systems of linear equations with L right hand side vectors. Practically,
however, we can take smaller N and larger ρ in the block version, so that the total
computational cost is similar to, or even smaller than the single case. The rigorous
proof of the block SS method, which includes the PEP with degeneracy, is underway.

5

Numerical examples

In this section, we investigate the theory using several examples of polynomial eigenvalue problems. The algorithm was implemented in MATLAB 7.4. We generated amatrix V := [v1 , . . . , vL ] using the MATLAB function rand and set U = V . The
MATLAB command mldivide\ was used to evaluate F (z)−1 V numerically. We compared the eigenvalues {λ̂j } that were obtained by the Block SS method and theeigenvalues {λ∗j } that were computed by calling the MATLAB function polyeig. We
deﬁne the residual of λj as kF (λ̂j )x̂j k2 , where x̂j is the corresponding eigenvector
and kx̂j k2 = 1.

Example 1
In Example 1, we conﬁrm the error estimation of the approximated eigenvalues in
Section 4. We consider the quadratic eigenvalue problem
F (λ) = λ2 A2 + λ(A1 + βIn ) + A0 + βA1
given in [3, p.10]. Here, Ai are matrices of order n = 10 and A0 is positive deﬁnite.
We set β = 0.8.
We took the integral path Γ as
Γ = γ + ρ eiθ

(γ = −0.8, ρ = 7).

so that 4 eigenvalues lie in Γ. The distribution of eigenvalues and the integral path
are shown in Figure 1. We took N = 32, K = 16, L = 1 and δ = 10−12 .
We obtained 14 eigenvalues of F (z) and computed approximation of their sensitivity factors τ̂k . Note that from equation (13) and equation (14), we can approximate
sensitivity factors τk byτ̂k = ρ(1 + |ηk |N )(1 + |λk |)kVm−T ek k22 ,where we set λk ← λˆk and m ← m̂.
The numerical results are shown in Table 1. In Table 1, underlines denote disagreement of digits between λ̂k and λ∗k . As shown in Table 1, the sensitivity factors
τ̂k are exponentially ampliﬁed if λk is outside Γ. We can also see that the relative
errors for the outer eigenvalues actually increase exponentially. Meanwhile, the inner eigenvalues were calculated accurately, even though the condition number of the
Hankel matrix was very large such that cond(Hm̂ ) = 8.25 × 1012 .
12

30
20
10
0
−10
−20
−30
−40

−30

−20

−10

0

10

20

30

40

Figure 1: Eigenvalues in Example 1.

Table 1: Relative errors and sensitivity factors in Example 1.
k
1
2
3
4
5
6
7
8
9
10
11
12
13
14

λ̂k
|λ̂k − λ∗k |/|λ̂k |−0.799999999999723 + 3.037712186344210 i
1.09×10−12− 0.800000000000111 − 3.037712186344154 i
1.07×10−12− 0.799999999999710 + 6.231420704050912 i
7.44×10−13− 0.799999999999999 − 6.231420704050882 i
7.38×10−13− 0.799999988132865 + 9.389825194259766 i
2.98×10−8− 0.799999997431621 − 9.389825196314405 i
3.00×10−8− 0.799979432271814 + 12.540406718026578 i
6.12×10−5
6.18×10−5− 0.799992884762234 −12.540413682454103 i− 0.799811460749113 + 15.699437889062162 i
8.54×10−4− 0.799917135958713 −15.699515527393906 i
8.59×10−4
9.98×10−3− 0.771833428737543 + 21.756907623843826 i− 0.784076104391534 −21.770287407471734 i
9.32×10−3
2.45×10−2− 0.788213898295595 + 28.952844300096686 i− 0.792216122017860 −28.958853070430731 i
2.47×10−213τ̂k
1.54×103
1.54×103
3.03×103
3.03×103
1.46×107
1.46×107
2.01×1010
2.00×1010
8.82×1011
8.78×1011
2.49×1013
2.50×1013
7.76×1013
7.79×1013

Table 2: Relative errors and residuals in Example 2.
|λ̂k − λ∗k |/|λ̂k |k
λ̂k
1
0.333333333333372
1.05 ×10−1320.499999999999953
8.24 ×10−1431.000000000000012
9.10 ×10−15
41.000000000000001 i 1.02 ×10−15
5 − 1.000000000000001 i 1.02 ×10−15

kF (λ̂k )x̂k k2
1.78 ×10−14
1.41 ×10−14
1.53 ×10−14
1.94 ×10−14
1.49 ×10−14

Example 2
In Example 2, we show that the SS method can be used to solve polynomial eigenvalue
problems. We consider the quadratic eigenvalue problem with

λ+1
6λ2 − 6λ0,0
F (λ) =  2λ 6λ2 − 7λ + 1
00
λ2 + 1
given in [18, p.250]. There are six eigenvalues λ1 = 1/3, λ2 = 1/2, λ3 = 1, λ4 = i, λ5 =−i, λ6 = ∞.
The integral path Γ was taken as follows:Γ = γ + ρ eiθ

(γ = 0, ρ = 1.5).

We took N = 10, K = 5, L = 1 and δ = 10−12 .
The numerical results are shown in Table 2. The SS method found ﬁve eigenvalues
inside the circle. Therefore, the SS method can be used for polynomial eigenvalue
problems.

Example 3
Next, we consider the connected damped mass-spring system [18, p.258]
F (λ) = λ2 A2 + λA1 + A0 ,where
A2 = In ,A1 = κ1 tridiag(−1, 3, −1),

A0 = κ0 tridiag(−1, 3, −1).

We took n = 50 degrees of freedom and chose κ0 = 5, κ1 = 3.
The integral path Γ was taken as follows:Γ = γ + ρ eiθ

(γ = −2 + 1.5i, ρ = 1).

All of the eigenvalues lie in the left half-plane shown in Figure 2. The circle showsΓ. A total of 16 eigenvalues lie in Γ.
14

2.5
2
1.5
1
0.5
0
−0.5
−1
−1.5
−2
−14

−12

−10

−8

−6

−4

−2

0

Figure 2: Eigenvalues in Example 3.
We set N = 32, K = 8, L = 24 and δ = 10−12 .
The numerical results are shown in Table 3. As shown in Table 3, we obtained all of
the eigenvalues in Γ. The largest residual of the computed eigenvalues is 2.72 × 10−12 .

Example 4
In Example 4, we show that the SS method can be applied to higher-order matrix
polynomials. We consider the quartic eigenvalue problem
F (λ) = λ4 A4 + λ3 A3 + λ2 A2 + λA1 + A0 ,given in [14]. Here, Ai are matrices of order n = 64, obtained by a tensor productconstruction. We constructed Ai by calling the MATLAB function kron.
The integral path Γ was taken as follows:Γ = γ + ρ eiθ

(γ = 1 − i, ρ = 0.5).

The distribution of eigenvalues and the integral path are shown in Figure 3. A
total of 13 eigenvalues lies in Γ. We set N = 64, K = 8, L = 24 and δ = 10−12 .
The numerical results are shown in Table 4. As shown in Table 4, we obtained all of
the eigenvalues in Γ. The largest residual of the computed eigenvalues is 2.32 × 10−12 .

15Table 3: Relative errors and residuals in Example 3.
k
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

λ̂k
|λ̂k − λ∗k |/|λ̂k | kF (λ̂k )x̂k k2−1.505690013788871 + 1.658874406037429 i
3.86×10−15
4.45×10−13−1.522738470984195 + 1.660440439240554 i
5.31×10−16
1.94×10−13−15−1.551080700948297 + 1.662653801098332 i
2.42×102.30×10−13−1.590609191894973 + 1.664930340377379 i
5.05×10−15
4.41×10−13−1.641173998718531 + 1.666471692426448 i
3.08×10−15
5.42×10−13−15−1.702583311786932 + 1.666279620106145 i
5.20×10
2.22×10−13−1.774604184541429 + 1.663167841804553 i
1.13×10−15
3.80×10−13−15−1.856963417142646 + 1.655767170992366 i
1.33×106.00×10−13−1.949348592811159 + 1.642519012495447 i
8.58×10−16
6.12×10−13−2.051409262931337 + 1.621650675206375 i
4.42×10−15
4.51×10−13−15−2.162758276422983 + 1.591122525925358 i
6.99×10
7.57×10−13−2.282973248338001 + 1.548529617142127 i
9.00×10−15
2.72×10−12−2.411598162111222 + 1.490924157988098 i
2.07×10−15
1.42×10−12−15−2.548145099387273 + 1.414487027309264 i
1.07×10
9.79×10−14−2.692096090862231 + 1.313876810729122 i
1.06×10−15
9.84×10−14−15−2.842905081118472 + 1.180779814421682 i
8.54×10
3.94×10−13

2
1.5
1
0.5
0
−0.5
−1
−1.5
−2
−1.5

−1

−0.5

0

0.5

1

Figure 3: Eigenvalues in Example 4.

16

1.5

Table 4: Relative errors and residuals in Example 4.
k
1
2
3
4
5
6
7
8
9
10
11
12
13

6

|λ̂k − λ∗k |/|λ̂k |λ̂k0.744283786331637 − 0.646541111813926 i 3.76 ×10−140.738844823886407 − 0.731658534005984 i 3.87 ×10−140.864617980453669 − 0.651815654480539 i 1.23 ×10−14
0.994127888031147 − 0.535135868221425 i 8.85 ×10−15
0.953854040217050 − 0.611439886362259 i 6.30 ×10−15
0.863349700394677 − 0.797929809342592 i 1.68 ×10−14
1.026189973208218 − 0.685703044215538 i 5.96 ×10−15
0.971854722649304 − 0.783539836463599 i 1.18 ×10−140.848570953056571 − 0.925677807336431 i 2.05 ×10−14
1.056265535074977 − 0.904134007343116 i 7.36 ×10−15
0.970370449857825 − 1.001776965449533 i 4.73 ×10−15
0.930660687304600 − 1.240183199928941 i 8.54 ×10−15
1.054414864515329 − 1.244513158205443 i 1.53 ×10−14

kF (λ̂k )x̂k k2
9.57 ×10−13
1.20 ×10−12
1.43 ×10−12
9.69 ×10−13
5.88 ×10−13
8.96 ×10−13
9.48 ×10−13
2.32 ×10−12
1.33 ×10−12
1.45 ×10−12
6.59 ×10−13
1.08 ×10−12
1.86 ×10−12

Conclusions

In the present paper, we have proposed a novel method to solve polynomial eigenvalue
problems. By using contour integral, the method converts the original problem into asmall linear eigenvalue problem that has identical eigenvalues in the domain deﬁned by
the integral path. The numerical evaluation of the contour integral can be performed
highly in parallel, which makes the method suitable to the modern cluster computerarchitecture. The present method is closely related to the contour integral method
proposed for generalized eigenvalue problems [16], though the theories behind them
are not identical. Relevance and irrelevance between them are still under study.
From an application point of view, the potential applicability of the polynomial
eigenvalue problem has not been scrutinized yet, especially for the higher order polynomials. It is probably due to the lack of eﬀective numerical methods to solve highdimensional and high-order polynomial eigenvalue problems. The proposed method
can deal with high-order problems without pain, which may derive quite a new algorithm in the ﬁeld of science and engineering. We are now under collaboration withphysicists and chemists in this direction.

References[1] Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, and H. Van Der Vorst, Templates forthe solution of algebraic eigenvalue problems: A practical guide, SIAM, Philadelphia (2000).
[2] Z. Bai and Y. Su, SOAR: A second order Arnoldi method for the solution of thequadratic eigenvalue problem, SIAM J. Matrix Anal. Appl., 26 (2005) 640–659.
17

[3] T. Betcke, N. J. Higham, V. Mehrmann, C. Schröder, and F. Tisseur, NLEVP:
A collection of nonlinear eigenvalue problem, MIMS EPrint 2008.40 (2008).
[4] F. R. Gantmacher, The theory of matrices, vol. 1-2, Chelsea Publishing Company(1977).
[5] I. M. Gel’fand, Lecture on linear algebra, Dover Publications, INC (1961).
[6] I. Gohberg, M. A. Kaashoek and P. Lancaster, General theory of regular matrix
polynomials and band toeplitz operators, Integr. Eq. Oper. Theory, 11 (1988)
776–882.
[7] I. Gohberg, P. Lancaster, F. Rodman, Matrix polynomials, Academic Press, NewYork (1982).
[8] R. S. Heeg, Stability and transition of attachment-line ﬂow, Ph.D.thesis, Universiteit Twente, Enschede, the Netherlands (1998).
[9] N. J. Higham and F. Tisseur, Bounds for eigenvalues of matrix polynomials, Lin.
Alg. Appl., 358 (2003) 5–22.
[10] T. Hwang, W. Lin, J. Liu and W. Wang, Jacobi-Davidson methods for cubic
eigenvalue problems, Numer. Lin. Alg. Appl., 12 (2005) 605–624.
[11] T. Ikegami, T. Sakurai, U. Nagashima, A ﬁlter diagonalization for generalized
eigenvalue problems based on the Sakurai-Sugiura projection method, Technical
Report CS-TR-08-13, Tsukuba (2008).
[12] P. Kravanja, T. Sakurai, H. Sugiura, M. Van Barel, A perturbation result for generalized eigenvalue problems and its application to error estimation in a quadrature method for computing zeros of analytic functions, J. Comput. Appl. Math.,
161 (2003) 339–347.
[13] D. S. Mackey, N. Mackey, C. Mehl and V. Mehrmann, Vector spaces of linearizations for matrix polynomials, SIAM J. Matrix Anal. Appl., 28 (2006) 971–1004.
[14] V. Mehrmann and D. Watkins, Polynomial eigenvalue problems with Hamiltonian structure, Electron. Trans. Numer. Anal., 13 (2002) 106–118.
[15] E. Pereira, On solvents of matrix polynomials, Appl. Numer. Math., 47 (2003)
197–208.
[16] T. Sakurai, H. Sugiura, A projection method for generalized eigenvalue problems,
J. Comput. Appl. Math., 159 (2003) 119–128.
[17] G. L. G. Sleijpen, A. G. L. Booten, D. R. Fokkema and H. A. Van Der Vorst,
Jacobi-Davidson type methods for generalized eigenproblems and polynomial
eigenproblems, BIT, 36 (1996) 595–633.

18

[18] F. Tisseur, K. Meerbergen, The quadratic eigenvalue problem, SIAM Rev., 43(2001) 235–286.
[19] E. E. Tyrtyshnikov, How bad are Hankel matrices?, Numer. Math., 67 (1994)
261–269.

19