Backward error bounds for polynomial eigenvalue problem solved by a Rayleigh–Ritz type contour integral-based eigensolver

Journal Pre-proof Backward error bounds for polynomial eigenvalue problem solved by a Rayleigh-Ritz type contour integral-based eigensolver

Hongjia Chen, Lei Du

PII: DOI: Reference:

S0893-9659(19)30446-X https://doi.org/10.1016/j.aml.2019.106122 AML 106122

To appear in:

Applied Mathematics Letters

Received date : 30 September 2019 Revised date : 1 November 2019 Accepted date : 1 November 2019 Please cite this article as: H. Chen and L. Du, Backward error bounds for polynomial eigenvalue problem solved by a Rayleigh-Ritz type contour integral-based eigensolver, Applied Mathematics Letters (2019), doi: https://doi.org/10.1016/j.aml.2019.106122. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Elsevier Ltd. All rights reserved.

Journal Pre-proof *Manuscript Click here to view linked References

Backward error bounds for polynomial eigenvalue problem solved by a Rayleigh-Ritz type contour integral-based eigensolver Hongjia Chena , Lei Dub1 a

Department of Mathematics, School of Science, Nanchang University, Nanchang, 330031, China b

School of Mathematical Sciences, Dalian University of Technology, Dalian 116024, China

na lP repr oo f

Abstract The contour integral-based eigensolvers have attracted much attention in recent years. In this paper, we consider solving a polynomial eigenvalue problem (PEP) by a contour integral-based eigensolver named the Sakurai-Sugiura method with Rayleigh-Ritz projection (SS-RR method). We derive a backward error bound of PEP solved by the SS-RR method. This bound can be used to show the accuracy of the computed approximate eigenpairs of PEP. The accuracy of the derived bounds is demonstrated by several examples.

Keywords: contour integral, polynomial eigenvalue problem, SS-RR method, backward error 2010 MSC: 65F15, 15A18

1. Introduction

Given a matrix polynomial P (λ) =

Pm

i=0

λi Ai , where Ai ∈ Cn×n , the polynomial eigenvalue

problem (PEP) is to find eigenvalues λ ∈ C and their corresponding eigenvectors x ∈ Cn \{0} such that

P (λ)x =

m X i=0

i

λ Ai

!

x = 0.

(1)

ur

PEP lies at the heart of computation science and engineering. In some applications, such as structural dynamics and structural-acoustic interactions, it is unnecessary to compute all eigen-

Jo

pairs (λ, x); we only focus on the eigenpairs of physical interest. The SS-RR method can extract all eigenvalues inside a Jordan curve Γ using a subspace constructed by a contour integral approach [1]. For solving the PEP, the SS-RR method firstly converts the original PEP (1) into a smaller projected PEP as ! m X i λ Ri y = V H P (λ)V y = 0, R(λ)y :=

(2)

i=0

5

where y ∈ Cℓ and the matrix V ∈ Cn×ℓ , ℓ ≪ n, has orthonormal columns consisting of basis vectors for the subspace constructed by the SS-RR method. 1 Corresponding

author. E-mail: [email protected]

A

Journal Pre-proof

Then the PEP (2) of small or medium-sized can be to solved by linearization as follows L(λ)z = (X − λY )z = 0,

(3)

where X, Y ∈ Cmℓ×mℓ and z ∈ Cmℓ . L(λ) and R(λ) have the same spectrum. The resulting eigenproblem (3) is usually solved by the state-of-the-art stable eigensolvers, such as the QZ method. Infinitely many different linearizations can be constructed such that they share the same eigenstructure with R(λ). One of the most commonly used linearizations is the companion linearization [2] as follows

10

· · · Rm−1





O

O .. . .. .

· · · −Rm

na lP repr oo f

 R0    L(λ) =    

R1 Iℓ

..

.

Iℓ

     I  − λ ℓ      

O

Iℓ

O

Iℓ ∈ Cℓ×ℓ is the identity matrix and z is the eigenvector of L(λ).



   ,   



y



     λy   z=  ..  ,  .    m−1 λ y

(4)

Let (λj , zj ) be approximate eigenpairs of L(λ), the approximate eigenpairs (λj , yj ) of R(λ) are obtained from

(λj , yj ) = (λj , zj (1 : ℓ)),

j = 1, . . . , n(Γ),

(5)

where zj (1 : ℓ) is a sub-vector of zj consisting of 1st to the ℓth elements. n(Γ) ≤ ℓ is the number of approximate eigenvalues in the target region Ω.

Finally, the approximate eigenpairs (λj , xj ) of P (λ) are computed by (λj , xj ) = (λj , V yj ),

j = 1, . . . , n(Γ).

(6)

For solving the PEP via the SS-RR method, it is important to estimate how accurate the approximate eigenpairs of the PEP (1) may be when the companion linearization (4) derived from the projected PEP (2). By ”accurate”, we mean the backward errors of eigenpairs.

ur

15

The purpose of this work is to investigate the backward errors of approximate eigenpairs of P (λ) obtained via the SS-RR method. Firstly, we give a few key definitions for the backward

Jo

errors in Section 2. Then we establish the bounds for backward errors of the SS-RR method in Section 3. These bounds give useful prediction about when the approximate eigenpairs of P (λ) 20

can be computed with both small and large backward errors. Next we present a few numerical experiments that confirm our predictions in Section 4. Finally, we give a short conclusion in Section 5. 2. Backward errors for matrix polynomial and linearization form In this section, we firstly introduce the definition of the normwise backward error of an ap-

25

proximate eigenpair of R(λ) and L(λ). Then we introduce the relation between (λ, y) of R(λ) and

Journal Pre-proof

(λ, z) of L(λ). Definition 1 ([3]). The normwise backward error of a finite approximate eigenpair (λ, x) of R(λ) is defined by ηR (λ, y) := min{ǫ : (R(λ) + ∆R(λ))y = 0, k∆Ri k2 ≤ ǫkRi k2 , i = 0 : m}, where ∆R(λ) =

m P

λi ∆Ri with ∆Ri being the perturbation to Ri .

i=0

An analogous definition holds for the backward error of ηL (λ, z) of an approximate eigenpair

na lP repr oo f

(λ, z). Definition 2 ([3]). The normwise backward error of a finite approximate eigenpair (λ, z) of L(λ) is defined by

ηL (λ, z) := min{ǫ : (L(λ) + ∆L(λ))z = 0, k∆Xk2 ≤ ǫkXk2, k∆Y k2 ≤ ǫkY k2 }, 30

where ∆X and ∆Y being the perturbation to X and Y .

The explicit and computable expression for the backward error ηR (λ, y) of an approximate eigenpair of R(λ) is also given in [3] by

kR(λ)yk2 ηR (λ, y) = Pm , ( i=0 |λ|i kRi k2 ) kyk2

(7)

and that of an approximate eigenpair of L(λ) is given by ηL (λ, z) =

kL(λ)zk2 . (kXk2 + |λ|kY k2 ) kzk2

(8)

where X, Y are from coefficient matrices of linearization (4).

To bound the backward error ratio ηR (λ, y)/ηL (λ, z), an appropriate two-sided transformation

ur

for L(λ) is given in [4] by

G(λ)L(λ)z = (eT1 ⊗ R(λ))z = R(λ)(eT1 ⊗ Iℓ )z = R(λ)y,

(9)

Jo

where G(λ) ∈ Cℓ×mℓ , e1 ∈ Cm×1 is the first unit vector and y = (eT1 ⊗ Iℓ )z. For the companion linearization (4), we give an explicit expression for G(λ) such that  m m P 1 P j 1 λ Rj , − λ12 λj Rj , . . . , − λm−1 G(λ) = Iℓ , − λ j=1

j=2

m P

j=m−1

 λj Rj .

(10)

3. Backward error analysis for P (λ) solved by the SS-RR method Now we have all the ingredients for bounding the backward errors of the eigenpairs of P (λ). The following lemma gives an upper bound for the backward error ratio ηP (λ, x)/ηR (λ, y) based 35

on (7).

Journal Pre-proof

Lemma 1 ([5]). Let (λ, y) be the approximate eigenpair of R(λ), with y normalized so that kyk2 = 1, (λ, x) be the approximate eigenpair of P (λ), where x = V y and kV yk2 = 1. Then, the bound for the ratio of ηP (λ, x) to ηR (λ, y) is given by ηP (λ, x) ≤ UP R , ηR (λ, y)

where UP R =

kP (λ)xk2 . kV H P (λ)V yk2

Remark 1. From [5], we know that kP (λ)xk2 ≈ kV H P (λ)V yk2 for the SS-RR method. Therefore, UP R ≈ 1. To bound the backward error of R(λ) relative to the backward error of L(λ), we have the

na lP repr oo f

following lemma. Lemma 2. Let (λ, y) and (λ, z) be approximate eigenpairs of R(λ) and companion linearization L(λ), respectively, with y and z normalized so that kyk2 = 1 and kzk2 = 1. The bound for the ratio of ηR (λ, y) to ηL (λ, z) is given by ηR (λ, y) ≤ URL , ηL (λ, z)

where URL = m

3/2

max{ max kRj k2 , 1}2 j=0:m

min{kR0 k2 , kRm k2 }

.

Proof. From (9), we have

kR(λ)yk2 ≤ kG(λ)k2 kL(λ)zk2 ,

which, along with (7) and (8), give

kG(λ)k2 (kXk2 + |λ|kY k2 ) ηR (λ, y) Pm ≤ . j ηL (λ, z) j=0 |λ| kRj k2

(11)

Based on Lemma 3.5 in [4], we have

kXk2 ≤ m max{ max (kRj k2 ), 1}, j=0:m−1

and

ur

kG(λ)k2 ≤ GU (λ) =

kY k2 ≤ m max{kRm k2 , 1},

o n X √ m−1 |λ|j max 1, max kRj k2 . m j=1:m

j=0

Jo

According to (11),

ηR (λ, y) ≤ ηL (λ, z)

mGU (λ)(|λ| + 1) max{ max kRj k2 , 1} j=0:m

(|λ|m

+ 1) min{kR0 k2 , kRm k2 }

.

Based on the upper bound GU , we have ηR (λ, y) ≤ ηL (λ, z)

m−1 P

|λ|j × max{ max kRj k2 , 1}2 j=0:m j=0 . (|λ|m + 1) min{kR0 k2 , kRm k2 }

m(|λ| + 1)

By the Theorem 3.6 in [4], note that (|λ| + 1)

m−1 P j=0

|λ|j

(|λ|m + 1)

≤ m1/2 .

Journal Pre-proof

Therefore, we have ηR (λ, y) ≤ ηL (λ, z)

m3/2 max{ max kRj k2 , 1}2 j=0:m

min{kR0 k2 , kRm k2 }

.

40

Theorem 1. Let (λ, x) and (λ, z) be the approximate eigenpairs of the matrix polynomial P (λ) and a companion linearization L(λ), respectively. Then the upper bound for the backward error ratio ηP (λ, x)/ηL (λ, z) is given by UP L = m

3/2

max{ max kRj k2 , 1}2 j=0:m

min{kR0 k2 , kRm k2 }

na lP repr oo f

ηP (λ, x) . UP L , ηL (λ, z)

.

(12)

Proof. Based on Lemmas 1 and 2, we have

ηP (λ, x) ηR (λ, y) ηP (λ, x) = ≤ UP R URL . URL = UP L . ηL (λ, z) ηR (λ, y) ηL (λ, z)

UP L only depends on the norms of the coefficient matrices of R(λ). From the bound (12), we investigate the condition when the backward error of the eigenpairs of P (λ) solved by the SS-RR method is not much larger than that of the companion linearization L(λ).

QZ method is a backward stable method for solving the GEP L(λ)z = 0, that is, ηL (λ, z) = O(mℓǫ), where ǫ is the machine precision and mℓ is the size of matrix coefficients of the GEP L(λ)z = 0 [3]. Therefore, we have the upper bound for backward error of approximate eigenpairs of the SS-RR method as

ηP (λ, x) . UP L O(mℓǫ).

45

(13)

If UP L = O(1), the backward error ηP (λ, x) is small and if UP L ≫ O(1), the backward error

ur

ηP (λ, x) may be very large.

4. Numerical experiments

Jo

In this section, we illustrate the bounds for the backward error of computed approximate eigenpairs of PEP via the SS-RR method. All test problems (shown in Table 1) are from the 50

collection of nonlinear eigenvalue problems NLEVP [6]. The degree of problem plasma−drift is 3, and the degree of other test problems is 2. The dimensions of all problems are n = 400. For each problem, the Jordan curve Γ is a circle with center γ and radius ρ whose values are given in Table 1. In our experiments, the eigenpairs of companion linearizations L(λ) are computed using MATLAB’s eig function and the eigenpairs of R(λ) and P (λ) are recovered by (5) and (6). The

55

backward errors of L(λ), R(λ) and P (λ) are obtained using (7) and (8), while the upper bounds for the backward errors of P (λ) are calculated using (13).

Journal Pre-proof

Table 1: The polynomial eigenvalue problem and parameters of the SS-RR method [5, 6]. Problem

n

center γ −2 + 2.6 ×

106 i

radius ρ 3×

#eigs

105

22

400

shaft

400

2 × 105 i

9 × 104

18

wiresaw1

400

−180i

40

26

wiresaw2

400

140i

40

26

sleeper

400

−16

0.2

29

spring

400

−12

1

26

dirac

400

−5

0.7

24

plasma− drift

400

10

1

10

na lP repr oo f

damped− beam

Table 2: Maximum backward error ηP (λ, x) and upper bound (13) .

Problems

max ηP (λ, x)

UP L

Upper bound (13)

damped−beam

3.5 × 10−6

7.9 × 1022

1.6 × 107

1.3 × 10−12

8.0 × 107

1.8 × 10−5

2.9 × 102

6.4 × 10−14

1.0 × 106

2.9 × 10−10

shaft

wiresaw1 wiresaw2 sleeper spring dirac

plasma−drift

1.6 × 10−10

7.6 × 1020

1.2 × 10−12

8.6 × 1010

4.0 × 10−15

2.0 × 103

6.7 × 10−15 2.2 × 10−15 5.6 × 10−13

4.2 × 104

2.1 × 105

1.9 × 10−5

4.4 × 10−13 1.2 × 10−12

The largest backward error ηP (λ, x) among all the computed eigenpairs and corresponding upper bound are shown in Table 2. Note that the backward error ηP (λ, x) is small when the upper bound (13) is not much far from O(ǫ). This is confirmed by Table 2, such as problems 60

sleeper, spring, plasma−drift and dirac. From (12), we know that the value of UP L depends

ur

on the norms of coefficient matrices of R(λ). Here, note that UP L only governs the backward errors and has nothing to do with the conditioning of the problem. Hence, if the coefficient matrices Ri vary widely in norms (for example, when kR0 k2 ≫ kR2 k2 and kR2 k2 ≫ kR1 k2 in QEP), the 65

Jo

upper bound (13) is too far from O(ǫ). In this case, the backward error ηP (λ, x) may be greatly larger than O(ǫ) and the SS-RR method is not stable for solving PEP. This is reflected in Table 2, such as problems damped−beam, shaft, wiresaw1 and wiresaw2. 5. Conclusion In this article, we have established an upper bound for the backward error of P (λ) solved by the SS-RR method. This bound can give useful information to indicate whether the backward 70

error of P (λ) is small or large.

Journal Pre-proof

References References [1] S. Yokota, T. Sakurai, A projection method for nonlinear eigenvalue problems using contour integrals, JSIAM Letters 5 (2013) 41–44. 75

[2] I. Gohberg, P. Lancaster, L. Rodman, Matrix polynomials, volume 58 of classics in applied mathematics, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA. [3] F. Tisseur, Backward error and condition of polynomial eigenvalue problems, Linear Algebra

na lP repr oo f

and its Applications 309 (1-3) (2000) 339–361.

[4] N. J. Higham, R.-C. Li, F. Tisseur, Backward error of polynomial eigenproblems solved by 80

linearization, SIAM journal on matrix analysis and applications 29 (4) (2007) 1218–1241. [5] H. Chen, A. Imakura, T. Sakurai, Improving backward stability of sakurai-sugiura method with balancing technique in polynomial eigenvalue problem, Applications of Mathematics 62 (4) (2017) 357–375.

[6] T. Betcke, N. J. Higham, V. Mehrmann, C. Schr¨oder, F. Tisseur, Nlevp: A collection of nonlinear eigenvalue problems, ACM Transactions on Mathematical Software (TOMS) 39 (2)

ur

(2013) 7.

Jo

85