Palindromic quadratization and structure-preserving algorithm for palindromic matrix polynomials of even degree

全文

(1)Numerische Mathematik. Numer. Math. (2011) 118:713–735 DOI 10.1007/s00211-011-0370-7. Palindromic quadratization and structure-preserving algorithm for palindromic matrix polynomials of even degree Tsung-Ming Huang · Wen-Wei Lin · Wei-Shuo Su. Received: 30 December 2009 / Revised: 25 November 2010 / Published online: 8 March 2011 © Springer-Verlag 2011. Abstract In this paper, we propose a palindromic quadratization approach, transforming a palindromic matrix polynomial of even degree to a palindromic quadratic pencil. Based on the (S + S −1 )-transform and Patel’s algorithm, the structurepreserving algorithm can then be applied to solve the corresponding palindromic quadratic eigenvalue problem. Numerical experiments show that the relative residuals for eigenpairs of palindromic polynomial eigenvalue problems computed by palindromic quadratized eigenvalue problems are better than those via palindromic linearized eigenvalue problems or polyeig in MATLAB. Mathematics Subject Classification (2000). 65F15 · 15A18 · 15A57. 1 Introduction In this paper, we consider (, ε)-palindromic matrix polynomials of even degree 2d P(λ) ≡. d−1 k=0. λ2d−k Ad−k + λd A0 + ε. d . λd−k Ak ,. (1). k=1. T.-M. Huang Department of Mathematics, National Taiwan Normal University, Taipei 116, Taiwan e-mail: [email protected] W.-W. Lin (B) Department of Mathematics, National Taiwan University, Taipei 106, Taiwan e-mail: [email protected] W.-S. Su Department of Applied Mathematics, National Chiao Tung University, Hsinchu 300, Taiwan e-mail: [email protected]. 123.

(2) 714. T.-M. Huang et al.. where d ≥ 2, ε = ±1 and = H (Hermitian) or T (transpose), Ak ∈ Cn×n (k = 0, 1, . . . , d) and A0 = ε A0 . The corresponding polynomial eigenvalue problem P(λ)x = 0, with λ ∈ C and x ∈ Cn \{0} being the eigenvalue and the associated eigenvector respectively, is called a (, ε)-palindromic polynomial eigenvalue problem ((, ε)-PPEP). It is also called a -PPEP if ε = 1 or a -anti-PPEP if ε = −1. The underlying matrix polynomial P(λ) in (1) has the property that reversing the order of coefficients, followed by taking the (conjugate) transpose, leads to the original matrix polynomial (anti-)invariant, which satisfies P(λ) = ελ2d P(1/λ). (2). and explains the world “(anti-)palindromic” [20]. Consequently, taking the (conjugate) transpose of (1), we easily see that the eigenvalues of P(λ) satisfy a “reciprocal” property, that is, they appear in the pairs of the form (λ, 1/λ ). The (, ε)-PPEPs arise in solving higher order systems of ordinary or partial differential equations. A T-palindromic quadratic eigenvalue problem (T-PQEP) is first raised in the vibration analysis for fast trains in Germany [11,12] and then in the study of surface acoustic wave filters [28]. An H-palindromic quadratic eigenvalue problem (H-PQEP) arises in the computation of the Crawford number, for detecting definite Hermitian pairs or hyperbolic or elliptic quadratic eigenvalue problems [9]. A -PPEP of even degree is obtained when solving the linear quadratic discrete-time optimal control problem for higher order systems [2,27]. A standard approach for computing the eigenpairs of P(λ) in (1) is to linearize it to a 2dn × 2dn linear matrix pencil by the companion linearization and compute its generalized Schur form [26]. However, the reciprocal property of eigenvalues of P(λ) is not preserved, generally, thus producing large numerical errors [6,13,15]. Recently, some pioneering works [20,21] propose some good linearizations that linearize P(λ) to palindromic linear pencils of the form λZ + Z which preserves the reciprocity of eigenvalues. This does lead to a vast improvement over previous unstructured approaches, reflecting the palindromic structure in the original polynomial and enabling structure-preserving numerical methods to be designed. Later, a Q R-like algorithm [24] and a hybrid method [19] which combines Jacobi-type method with the Laub’s trick, a postprocessing step of the generalized Schur form, are proposed for solving T-palindromic linear eigenvalue problems efficiently. The Q R-like algorithm typically requires O(n 4 ) flops and the hybrid method requires O(n 3 log(n)) flops. Recently, a URV-decomposition based structured method of cubic complexity is developed in [25] to solve T-palindromic linear eigenvalue problems, producing eigenvalues which are reciprocally paired to working precision. A new structurepreserving doubling algorithm with cubic complexity for solving -palindromic linear eigenvalue problems is developed in [4]. On the other hand, for solving a (, ε)-PQEP, a structure-preserving doubling algorithm is developed in [5,6] via the computation of a solvent of a nonlinear matrix equation associated with the (, ε)-PQEP. Lately, a numerically stable structure-preserving algorithm (SPA), based on the (S + S −1 )transform [18] and Patel’s algorithm [22], is proposed in [13] to solve the T-PQEP directly. The numerical results obtained by the SPA algorithm show much promise and the computational cost of SPA is about a half of that of the URV-based method.. 123.

(3) Palindromic quadratization and structure-preserving algorithm. 715. The purpose of this paper is to develop a palindromic quadratization which transforms a (, ε)-palindromic matrix polynomial of even degree with (, ε) = (T, −1) into a (, ε)-palindromic quadratic pencil. If = T and ε = 1, then we can apply the SPA algorithm in [13] to solve the associated quadratized T-PQEP directly. If = H and ε = ±1, we first transform the associated quadratized (H, ε)-PQEP to an H-skew-Hamiltonian pencil by the (S + S −1 )-transform and enlarge the H-skewHamiltonian pencil to a real skew-Hamiltonian pencil, to which the SPA algorithm is applicable. Note that for the case (, ε) = (T, −1), the T-anti-PPEP can then be solved by applying the URV-based method [23,25] to the linearized T-palindromic linear pencil. Throughout this paper, Cn×m is the set of all n × m complex matrices, Cn = Cn×1 and C = C1 ; In and 0n×m (or simply I and 0 if their dimensions are clear from the context) denote the n × n identity matrix and n × m zero matrix, respectively. The jth column of an identity matrix is denoted by e j . The direct sum of two matrices is denoted by “⊕”. This paper is organized as follows. In Sect. 2, we propose a palindromic quadratization for a (, ε)-palindromic matrix polynomial of even degree. In Sect. 3, we develop a structure-preserving algorithm for solving the H-PQEPs. We develop balancing techniques for PPEPs and PQEPs in Sect. 4. Comparisons of numerical results computed by the palindromic quadratization, the palindromic linearization and the standard companion linearization are presented in Sect. 5. Conclusions are given in Sect. 6. 2 P-Quadratization of (, ε)-PPEP In [13], a structure-preserving algorithm is well-developed for solving the T-PQEP. A similar structure-preserving algorithm for solving the H-PQEP will be introduced in Sect. 3. As the H-anti-PQEP can be easily transformed to the H-PQEP, all (, ε)-PQEPs with (, ε) = (T, −1) can be solved by structure-preserving algorithms. Furthermore when (, ε) = (T, −1), we shall propose a new palindromic quadratization (P-quadratization) which can be utilized to transform a (, ε)-PPEP into a (, ε)-PQEP so that the structure-preserving algorithm in [13] is applicable. Next we present definitions of quadratization and P-quadratization of a general matrix polynomial and a palindromic matrix polynomial, respectively. Definition 1 (Quadratization/P-Quadratization) (i) Let P(λ) be an arbitrary ν × ν matrix polynomial of degree p ≥ 2 with pν = 2q. A q × q quadratic matrix polynomial (quadratic pencil) Q(λ) is a quadratization of P(λ) if there are matrix rational functions E(λ) and F(λ) of size q × q with nonzero and constant determinants satisfying the two-sided factorization E(λ)Q(λ)F(λ) =. P(λ) 0 . 0 Iq−ν. (3). 123.

(4) 716. T.-M. Huang et al.. (ii) Let P(λ) be an arbitrary ν × ν (, ε)-palindromic matrix polynomial of degree p ≥ 2 with pν = 2q (i.e., P(λ) = ελ p P(1/λ) as in (2)). A quadratization Q(λ) of P(λ) having the (, ε)-palindromic structure is called a P-quadratization of P(λ). Theorem 1 Let Q(λ) be a q × q quadratization of a ν × ν matrix polynomial P(λ) of degree p with pν = 2q. Then (i) λ0 ∈ C is a finite eigenvalue of Q(λ) (i.e., det(Q(λ0 )) = 0) if and only if λ0 is a finite eigenvalue of P(λ) (i.e., det(P(λ 0 )) = 0). (ii) ∞ is an eigenvalue of Q(λ) (i.e., det [λ2 Q( λ1 )]λ=0 = 0) if and only if ∞ is an eigenvalue of P(λ) (i.e., det λ p P λ1 λ=0 = 0). Proof. (ii). (i) The factorization (3) implies that det(Q(λ)) = c det(P(λ)) for some nonzero constant c, so that Q(λ) and P(λ) are singular or nonsingular for precisely the same values of λ0 . Since. 1 1 1 2q 2q = λ det Q = cλ det P det λ Q λ λ λ . 1 , = c det λ p P λ . 2. both Q(λ) and P(λ) have or have no infinite eigenvalues.. . Since both det(E(λ)) and det(F(λ)) are nonzero and constant, it is easily seen that the two-sided factorization (3) implies the existence of a more wide class of one-sided factorization Q(λ)F(λ) ≡ Q(λ)F(λ). In 0. . = E(λ)−1. . In P(λ) ≡ G(λ)P(λ), 0. (4). where F(λ) and G(λ) are matrix rational functions of size q × ν. From the factorization (4) a close connection between eigenpairs of P(λ) and eigenpairs of Q(λ) has been shown in [8]. Theorem 2 [8] Assume that (4) holds at λ0 ∈ C with F(λ0 ) and G(λ0 ) being of full column rank. Then F(λ0 )z 1 is an eigenvector of Q(λ) if and only if z 1 is an eigenvector of P(λ), both corresponding to eigenvalue λ0 . In Definition 1(i), we give a new definition of quadratization for a general matrix polynomial. In Theorems 1 and 2, we show the connection between eigenpairs of a general matrix polynomial and its quadratization. We next present a P-quadratization for a palindromic matrix polynomial of even degree explicitly.. 123.

(5) Palindromic quadratization and structure-preserving algorithm. 717. Theorem 3 Let P(λ) be an n × n (, ε)-palindromic matrix polynomial of degree 2d as in (1) with (, ε) = (T, −1). Then P(λ) can be P-quadratized into a (, ε)palindromic quadratic pencil of the form Q(λ) ≡ λ2 A1 + λA0 + εA1. (5). with A0 = εA0 , where (i) (For 2d = 4m) A1 and A0 are given by

(6) A1 =. d1 I. A1. √ ε εd1 A2. 0.

(7). , A0 =. A0 −. √. εI −. √ ε A2 A2. 0. √ − εd1 d1 I. 0. (6). if m = 1; otherwise, ⎡. A1. 0. ···. ⎢ √ ⎢ ε εd1 A2m 0 ··· ⎢ ⎢ ⎢ ⎢ A2m−1 d I . . . 1 ⎢ ⎢ ⎢ 0 −εd2 I A1 = ⎢ ⎢ ⎢ ⎢ .. ⎢ . ⎢ ⎢ ⎢ ⎢ A3 ⎣ 0 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ A0 = ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣. A0 −. ···. ···. 0. ···. ···. ···. ... .. ... .. √ √ ε I − ε A2m A2m 0. ... dm I. .. dm−1 I. 0 −εdm I. 0. 0. 0 .. . ε A2 0. ···. ⎥ 0 ⎥ ⎥ ⎥ .. ⎥ . ⎥ ⎥ ⎥ .. ⎥ . ⎥ ⎥, ⎥ .. ⎥ . ⎥ ⎥ ⎥ .. ⎥ . ⎥ ⎦ 0. (7a). A2m−2 0 · · · A2 0. √ − εd1 d1 I. ε A2m−2. ⎤. ···. ⎤. ⎥ .. ⎥ .⎥ ⎥ ⎥ .. ⎥ .⎥ ⎥ ⎥ .. ⎥ ; 0 .⎥ ⎥ ⎥ .. ⎥ .. . .⎥ ⎥ ⎥ 0 0⎥ ⎦ ··· ··· ··· 0 (7b). 123.

(8) 718. T.-M. Huang et al.. (ii) (For 2d = 4m + 2) A1 and A0 are given by ⎡. A1. ⎢ A2m+1 ⎢ ⎢ ⎢ 0 ⎢ ⎢ ⎢ A2m−1 A1 = ⎢ ⎢ ⎢ 0 ⎢ ⎢ . ⎢ .. ⎢ ⎣ A3 ⎡. ··· ··· . −εd1 I . .. ··· ···. d1 I. .. 0 0. ... ··· ···. . −εd2 I . . .. .. ··· ···. ... .. dm−1 I. 0 A2m. A0 0 ⎢ ⎢ ε A2m 0 ⎢ ⎢ ⎢ 0 0 ⎢ ⎢ ⎢ A0 = ⎢ ε A2m−2 ⎢ ⎢ 0 ⎢ ⎢ .. ⎢ . ⎢ ⎣ εA 2 0 ··· ···. A2m−2. 0. ···. 0 ···. 0 ···. A2. 0 −εdm I ⎤. 0 .. ⎥ .⎥ ⎥ .. ⎥ .⎥ ⎥ .. ⎥ .⎥ ⎥; .. ⎥ 0 .⎥ ⎥ .. ⎥ .. . .⎥ ⎥ 0 0⎦ ··· ··· ··· 0. ⎤ dm I 0 ⎥ ⎥ .. ⎥ . ⎥ ⎥ .. ⎥ . ⎥ ⎥ .. ⎥ , . ⎥ ⎥ .. ⎥ . ⎥ ⎥ 0 ⎦. (8a). 0. (8b). in which d1 , . . . , dm are arbitrary nonzero constants. Proof. (i) (For 2d = 4m) We define the n × n matrix rational functions for j = 1, . . . , m: ⎛ ⎞ 2 j−2 −1 √ ⎝ε ελ2m In + F2 j,1 (λ) = λ2 j−1 d j λk A2m−k ⎠ , (9a) k=0. and for j = 1, . . . , m − 1: F2 j+1,1 (λ) = λ2m−2 j In .. (9b). Let F1 (λ) =. In F2,1 (λ)T F3,1 (λ)T · · · F2m,1 (λ)T 0 I(2m−1)n. T .. (10). Routine but tedious calculation in terms of F1 (λ) in (10), Q(λ) in (5) and A0 , A1 in (6) or (7) leads to Q 1 (λ) In = , (11) Q(λ)F1 (λ) 0(d−1)n,n 0(d−1)n,n. 123.

(9) Palindromic quadratization and structure-preserving algorithm. 719. where √ √ Q 1 (λ) = λ2 A1 + λ A0 − ε I − ε A2m A2m + ε A1 m−1 √ √ +λ ε A2m λ2m ε ε I + A2m + λ2m−2k λ2 A2m−2k+1 + λA2m−2k +ελ. 1−2m. k=1. √. ε ελ. 2m. I+. 2m−2 . . λ A2m−k. = λ1−2m P(λ).. k. (12). k=0. Next, we define for j = 1, . . . , m: E 1,2 j+1 (λ) = λ2 j−2m In , for j = 1, . . . , m − 1, ⎛ ⎞ j−2 −1 2 √ ⎝ λ2m−k A2m−k + ε In ⎠ , E 1,2 j (λ) = λ2m−2 j+1 d j k=0. and let E 1 (λ) =. In E 1,2 (λ) E 1,3 (λ) · · · E 1,2m (λ) 0 I(2m−1)n. .. (13). From (12) and the definition of E 1 (λ) in (13), we have . In 0n,(d−1)n E 1 (λ)Q(λ)F1 (λ) = λ1−2m P(λ) 0n,(d−1)n .. Then from (10) and (13), it follows for m = 1 that . 0 In. . 0 E 1 (λ)Q(λ)F1 (λ) In. . √ = −λ εd1 d1 In ;. or for m ≥ 2: . . 0. . 0 I(d−1)n E 1 (λ)Q(λ)F1 (λ) I(d−1)n ⎡ √ ⎤ 2 −λ εd1 d1 In λ d1 In 0 ··· ··· 0 .. ⎢ ⎥ .. ⎢ ⎥ . . 0 −λ2 εd2 In εd1 In ⎢ ⎥ ⎢ ⎥ .. .. .. .. ⎢ ⎥ . . . 0 −d2 In . ⎢ ⎥. =⎢ ⎥ . . . . .. .. .. .. 2 ⎢ ⎥ 0 λ dm−1 In ⎢ ⎥ ⎢ ⎥ .. .. 2 ⎣ . 0 −λ εdm In ⎦ . εdm−1 In 0 ··· ··· 0 −dm In 0. 123.

(10) 720. T.-M. Huang et al.. Using the factorizations

(11). I 0 √n ε I I n λd j n. . √

(12) In −λ εd j d j In λ2 d j In εd j In 0 0. √ λ In εd j. In. √ 0 −λ εd j d j In √ = 0 λ ε In . and

(13). In d j+1 √ In ελ. =. 0 In. √ λ ε In 0. . √ √ In λ εd j+1 In λ ε In −λ2 εd j+1 In 0 In −d j+1 In 0 √ 0 , −λ εd j+1 d j+1 In. it holds that E 2m−1 (λ) · · · E 2 (λ)E 1 (λ)Q(λ)F1 (λ)F2 (λ) · · · F2m−1 (λ) √ √ √ √ = diag λ1−2m P(λ), −λ εd1 d1 In , λ ε In , . . . , λ ε In , −λ εdm dm In , where

(14) E 2 j (λ) = I(2 j−1)n ⊕

(15) E 2 j+1 (λ) = I2 jn ⊕. I √n ε λd j In. In d j+1 √ In ελ. 0 In. 0 In. ⊕ I(2m−2 j−1)n ,. ⊕ I(2m−2 j−2)n ,. and

(16) F2 j (λ) = I(2 j−1)n ⊕ F2 j+1 (λ) = I2 jn. In. √ λ In εd j. ⊕ I(2m−2 j−1)n , 0 In √ In λ εd j+1 In ⊕ I(2m−2 j−2)n , ⊕ 0 In. for j = 1, . . . , m − 1. Finally, letting √ √ √ E 2m (λ) := diag λ2m−1 In , −(λ εd1 d1 )−1 In , (λ ε)−1 In , −(λ εd2 d2 )−1 In , √ √ . . . , (λ ε)−1 In , −(λ εdm dm )−1 In , E(λ) := E 2m (λ) · · · E 1 (λ) and F(λ) := F1 (λ) · · · F2m−1 (λ), one can easily verify that E(λ)Q(λ)F(λ) = diag(P(λ), I(2m−1)n ). Furthermore, it holds that det(E(λ)) = √ (−ε)m /( ε mj=1 d j d j ) and det(F(λ)) = 1.. 123.

(17) Palindromic quadratization and structure-preserving algorithm. 721. (For 2d = 4m + 2) Let. (ii). . 2 j = I(2 j−1)n. 0 In ⊕ In 0. ⊕ I(2m−2 j)n , for j = 1, . . . , m.. Similar to part (i), we define n × n matrix rational functions E 1 (λ) and F1 (λ) by E 1 (λ) = 2 F1 (λ) =. In E 1,2 (λ) E 1,3 (λ) · · · E 1,2m+1 (λ) 0 I2mn. In F2,1 (λ)T F3,1 (λ)T · · · F2m+1,1 (λ)T 0 I2mn. , T. with E 1,2 j (λ) = λ2 j−2m−2 In ,. E 1,2 j+1 (λ) = d −1 j. 2 j−1. λk+1 A2m+k−2 j+2 ,. k=0. F2 j,1 (λ) = λ. 2m−2 j+2. In ,. j−1 −1 2 F2 j+1,1 (λ) = d j λ2 j λk A2m−k+1. (14). k=0. for j = 1, . . . , m. Via careful calculation we get Q(λ)F1 (λ). In 0. . =. λ−2m P(λ) 0. . and . . In 0 E 1 (λ)Q(λ)F1 (λ) = λ−2m P(λ) 0 .. Letting. E j (λ) = 2 j. I2( j−1)n ⊕ ⎡ In F j (λ) = I(2 j−3)n ⊕ ⎣ 0 0. . In. λ−2 In. 0 In ⎤. . ⊕ I(2m−2 j+1)n ,. 0 λ2 I n In In ⎦ ⊕ I(2m−2 j+1)n , 0 In. for j = 2, . . . , m, and E m+1 (λ). −1 −1 −1 −1 −1 −1 = diag λ2m In , In , 2 In , In , 2 In , . . . , In , 2 In , d1 λ εd1 d2 λ εd2 dm λ εdm. 123.

(18) 722. T.-M. Huang et al.. one can also verify that E(λ)Q(λ)F(λ) = diag(P(λ), I2mn ), where E(λ) = E m+1 (λ) · · · E 1 (λ) and F(λ) = F1 (λ) · · · Fm (λ). Furthermore, it holds that det(E(λ)) =. ε−m /( mj=1 (d j d j )) and det(F(λ)) = 1. Note that the P-quadratization of a (, ε)-palindromic matrix polynomial of odd degree with even matrix dimension can be defined as in Definition 1. However, to the best of our knowledge, a P-quadratization of this type has not been found. We now show the relationship between eigenpairs of Q(λ) in (5) and P(λ) in (1). Theorem 4 Let Q(λ) in (5) be a P-quadratization of P(λ) in (1) with (, ε) = T with z j ∈ Cn ( j = 1, . . . , d). Then (T, −1). Denote z = z 1T · · · z dT (i) For λ0 = 0, (λ0 , z 1 ) is an eigenpair of P(λ) if and only if (λ0 , z) is an eigenpair of Q(λ) with z j = F j,1 (λ0 )z 1 ( j = 2, . . . , d), where {F j,1 (λ)}dj=2 are given in (9) (for d = 2m) or (14) (for d = 2m + 1). (ii) (0, z 1 ) is an eigenpair of P(λ) if and only if (0, z) is an eigenpair of Q(λ), where, for j = 1, . . . , m − 1: (for d = 2m) (or, for d = 2m + 1). z 2m = −(dm )−1 A1 z 1 , z 2 j+1 = 0, z 2 j = −(d j )−1 A2m−2 j+1 z 1 ;. (15a). z 2m+1 = −(dm )−1 A1 z 1 , z 2m = 0, z 2 j = 0, z 2 j+1 = −(d j )−1 A2m−2 j+1 z 1 .. (15b). (iii) (∞, z 2 ) is an eigenpair of P(λ) if and only if (∞, z) is an eigenpair of Q(λ), with z 1 = z 3 = · · · = z d = 0. Proof (i) From Theorem 3, there are dn × dn matrix rational functions E(λ) and F(λ) with nonzero and constant determinants such that E(λ)Q(λ)F(λ) = diag(P(λ), I(d−1)n ). Since λ = 0 is the only pole of E(λ) and F(λ), the matrices F(λ0 ) and G(λ0 ) defined in (4) are of full rank. The assertion in (i) follows immediately from Theorem 2 and the relation ⎤ ⎡ ⎤ ⎡ z1 z1 ⎢ F (λ )z ⎥ ⎢ z ⎥ 2,1 0 1 z1 z1 ⎥ ⎢ 2⎥ ⎢ = F1 (λ0 ) =⎢ F(λ0 )z 1 = F(λ0 ) ⎥ = ⎢ .. ⎥ . .. 0 0 ⎦ ⎣ . ⎦ ⎣ . Fd,1 (λ0 )z 1 zd T (ii) By the definition of Q(λ) in (5), we have Q(0)z = εq ≡ ε q1T · · · qdT , where for d = 2m: √ q1 = A1 z 1 + dm z 2m , q2 = ε εd1 A2m z 1 , q3 = A2m−1 z 1 + d1 z 2 , (16a) z 2m−2 , q2m = −εdm z 2m−1 ; q4 = −εd2 z 3 , . . . , q2m−1 = A3 z 1 + dm−1 and for d = 2m + 1: q1 = A1 z 1 + dm z 2m+1 , q2 = A2m+1 z 1 , q3 = −εd1 z 2 , q4 = A2m−1 z 1 + d1 z 3 , z 2m−1 , q2m+1 = −εdm z 2m . q5 = −εd2 z 4 , . . . , q2m = A3 z 1 + dm−1 (16b). 123.

(19) Palindromic quadratization and structure-preserving algorithm. 723. From (16), we see that (0, z 1 ) is an eigenpair of P(λ); i.e., Ad z 1 = 0 if and only if (0, z) is an eigenpair of Q(λ); i.e., Q(0)z = εq = 0, where {z j }dj=1 satisfy (15). (iii) By the definition of Q(λ) in (5) again, we have λ2 Q λ1 λ=0 z = q ≡ T T q1 · · · qdT , where for d = 2m: . !m−1 √ q1 = A1 z 1 + εd1 A2m z 2 + k=1 A2m−2k+1 z 2k+1 , q2 = d1 z 3 , q3 = −εd2 z 4 , . . . , q2m−2 = dm−1 z 2m−1 , q2m−1 = −εdm z 2m , q2m = dm z 1 ; (17a). and for d = 2m + 1: ! q1 = A1 z 1 + m k=1 A2m−2k+3 z 2k , q2 = −εd1 z 3 , q3 = d1 z 4 , . . . , q2m = −εdm z 2m+1 , q2m+1 = dm z 1 .. (17b). From (17), it follows that (∞, z 2 ) is an eigenpair of P(λ); i.e., Ad z 2 = 0 if and only if (∞, z) is an eigenpair of Q(λ), or λ2 Q λ1 λ=0 z = q = 0, with z 1 = z 3 = · · · =. z d = 0. Remark 1 (i) Theorem 3 cannot be applied to the case (, ε) = (T, −1). In fact, up to now, there is no structure-preserving algorithm to solve the T-anti-PQEP directly. So it is pointless to transform a T-anti-PPEP to a T-anti-PQEP. Thus, for a (T, −1)palindromic matrix polynomial, we can apply the palindromic linearization [20] to transform it into a T-palindromic linear pencil λZ T + Z , and then solve it by the Q R-like algorithm [24], the hybrid method [19], the URV-based method [25] or the doubling algorithm [4]. (ii) On the other hand, if we rewrite λZ T + Z to a T-palindromic quadratic pencil "" λ0 + Z by letting " λ2 = λ, then the SPA algorithm [13] can also Q( λ) ≡ " λ2 Z T + " be used to solve its eigenpairs. It is shown in [13] that applying the SPA to solve "" Q( λ)y = 0 is mathematically equivalent to applying the URV-based method to solve λZ T + Z . Applying the P-quadratization in Theorem 3, an H-anti-PPEP can be quadratized into an H-anti-PQEP whose eigenpairs can then be computed from an H-PQEP by the following relationship. Proposition 1 Given an H-anti-PQEP: (λ2 A1H + λA0 − A1 )x = 0, with A0H = −A0 . Then (ıω, x) is an eigenpair of the H-anti-PQEP if and only if (ω, x) is an eigenpair of the H-PQEP: ω2 (−A1 ) H + ω(ı A0 ) + (−A1 ) x = 0. Proof The result can be easily obtained by setting λ = ıω and using the fact (ı A0 ) H = ı A0 .. We now consider and -odd polynomial eigenvalue problems of even de! the -even k n×n (k = 0, 1, . . . , 2d) and C = 0. The gree. Let C(λ) = 2d 2d k=0 λ C k , where C k ∈ C polynomial eigenvalue problem C(λ)x = 0 is called a -even polynomial eigenvalue =C problem, if C2k 2k (k = 0, 1, . . . , d) and C 2k−1 = −C 2k−1 (k = 1, . . . , d); and = −C (k = 0, 1, . . . , d) it is called a -odd polynomial eigenvalue problem, if C2k 2k. 123.

(20) 724. T.-M. Huang et al.. H−even PEP. [20]. H−palindromic PEP. Theorem 3. H−palindromic QEP. T−even PEP. [20]. T−palindromic PEP. Theorem 3. T−palindromic QEP. H−odd PEP. [20]. H−anti−palindromic. Theorem 3. H−anti−palindromic QEP. [20, 21]. T−palindromic GEP. T−odd PEP. PEP [20]. T−anti−palindromic PEP. Proposition 1. H−palindromic QEP. Remark 1. T−palindromic QEP. Fig. 1 Relations between various structured polynomial eigenvalue problems (PEPs) and C2k−1 = C2k−1 (k = 1, . . . , d). By the Cayley transformation, it was shown in [20] that a -even/odd polynomial eigenvalue problem can be transformed to a (, ±1)-PPEP, respectively. We illustrate the relationship among various structured polynomial eigenvalue problems in Fig. 1. We see that T-even, T-odd, T-anti-palindromic and T-palindromic polynomial eigenvalue problems of even degree can be P-quadratized to T-PQEPs. Thus, the SPA algorithm in [13] can be applied to solve the associated T-PQEPs. On the other hand, we see that H-even, H-odd, H-anti-palindromic and H-palindromic polynomial eigenvalue problems of even degree can be P-quadratized to H-PQEPs. Therefore, we are motivated to develop a structure-preserving algorithm in the next section to solve the H-PQEP.. 3 H-palindromic quadratic eigenvalue problems Consider the H-PQEP Q(λ)x ≡ λ2 A1H + λA0 + A1 x = 0,. (18). where A0 , A1 ∈ Cn×n with A0H = A0 . The eigenvalues of Q(λ) clearly appear in the ¯ “reciprocal” pairs of the form (λ, 1/λ). Classical linearizations of (18) in a companion form, generally, do not preserve the symplectic structure. Fortunately, the special linearization. (M − λL)z ≡. 0 A1 −A0 −I. . . 0 I −λ A1H 0. x =0 y. (19). of (18) (see [5] or [13]) satisfies 0 In , ≡ −In 0 . MJ M. H. = LJ L , J = J2n H. (20). ¯ preserving reciprocity. The so that the matrix pair (M, L) has eigenvalues λ and 1/λ, pencil M − λL or the matrix pair (M, L) are called H-symplectic.. 123.

(21) Palindromic quadratization and structure-preserving algorithm. 725. For any real symplectic matrix pair (M, L) satisfying (20), a structure-preserving (S +S −1 )-transform for the computation of all eigenvalues was proposed by [18]. The (S + S −1 )-transform (Ms , Ls ) of an H-symplectic matrix pair (M, L) is defined by Ms ≡ MJ L H + LJ M H ≡ KJ ,. Ls ≡ LJ L H ≡ N J .. (21). It is easily seen that K and N are both H-skew-Hamiltonian, i.e., KJ = J K H and ¯ N J = J N H . Hence, if μ is an eigenvalue of (K, N ), so is μ. The relationship between eigenpairs of an H-symplectic matrix pair and its (S + S −1 )-transform has been given in [18]. Theorem 5 [18] Let (M, L) be H-symplectic and (Ms , Ls ) be its (S + S −1 )to μ = transform z s is an eigenvector of (Ms , Ls ) corresponding as in (21). Suppose ν + ν1 . If L H − ν1 M H z s = 0 or (L H −νM H )z s = 0, then ν, J L H − ν1 M H z s or (1/ν, J (L H − νM H )z s ) is an eigenpair of (M, L), respectively. Then the relationship between eigenpairs of (Ms , Ls ) and Q(λ) is given as follows. Theorem 6 Let (M, L) be H-symplectic of the form in (19) and (Ms , Ls ) be its (S + S −1 )-transform. Suppose that z s is an eigenvector of (Ms , Ls ) corresponding T , z T ]T with z , z ∈ Cn . Let ν be to the eigenvalue μ = ±2, and denote z s ≡ [z s1 s1 s2 s2 1 a root of the quadratic equation λ + λ = μ. Then (i) at least one of vectors z s1 + ν1 z s2 and z s1 + νz s2 is nonzero; (ii) if z s1 + ν1 z s2 = 0, then z s1 + ν1 z s2 is an eigenvector of Q(λ) corresponding to ν; (iii) if z s1 + ν1 z s2 = 0, then z s2 is an eigenvector of Q(λ) corresponding to ν1 . Proof (i) Suppose that z s1 + ν1 z s2 = 0 and z s1 +νz s2 = 0. It implies that ν − ν1 z s2 = 0. If z s2 = 0, then z s1 = 0 and z s = 0 which contradicts the fact that z s is an eigenvector. Hence, ν = ±1 and then μ = ±2 which contradicts the assumption that μ = ±2. Therefore, z s1 + ν1 z s2 = 0 or z s1 + νz s2 = 0. (ii) Since MJ M H = LJ L H , by (21) it holds that. 1 H H zs . 0 = (Ms − μLs ) z s = (M − νL)J L − M ν. (22). From (22), we obtain (M − νL). z s1 + ν1 z s2 xν. = 0,. (23). where xν ≡. 1 H 1 A1 z s1 − A0 z s2 − A1 z s2 . ν ν. (24). 123.

(22) 726. T.-M. Huang et al.. Substituting (M, L) of (19) into (23), we have xν =. 1 1 A1 z s1 + z s2 ν ν. (25). and. 1 1 A0 z s1 + z s2 + xν + ν A1H z s1 + z s2 = 0. ν ν. (26). Substituting xν of (25) into (26) and multiplying (26) by ν, we get Q(ν) z s1 + ν1 z s2 =0. (iii) Since z s1 + ν1 z s2 = 0, it follows that z s1 = − ν1 z s2 = 0 and xν = 0 in (25). Substituting these results into (24), it holds that. 2 1 1 0 = xν = − A1H z s2 − A0 z s2 − A1 z s2 . ν ν Therefore, z s2 is an eigenvector of Q(λ) corresponding to the eigenvalue ν1 .. . In [13], a structure-preserving algorithm (SPA) based on Patel’s algorithm [22] has been developed for solving T-PQEPs. In order to solve an H-PQEP, we apply the (S + S −1 )-transform to the H-symplectic pair (M, L) of the form (19) and get the generalized eigenvalue problem Kz s = μN z s , where K and N are defined in (21). Substituting (M, L) in (19) into (21), the H-skew-Hamiltonian K and N can be represented as −A1 0 A1H − A1 A0 , N = . (27) K= 0 −A1H A1 − A1H A0 However, Patel’s algorithm can only be applied to (K, N ) of (27) in the real case, but cannot be directly applied to (K, N ) in the complex conjugate case. In the following, we convert (K, N ) of (27) into an enlarged real T-skew-Hamiltonian pair so that Patel’s algorithm can be applied. We extend (K, N ) in (27) to a real 4n × 4n matrix pair (K2 , N2 ) by K2 =. K R −K I KI KR. . , N2 =. N R −N I NI NR. ∈ R4n×4n ,. (28). where K = K R + ıK I and N = N R + ıN I . From (28), it is easily seen that if μ is an eigenvalue of (K, N ), then μ and μ¯ are eigenvalues of (K2 , N2 ). Theorem 7 The multiplicities of eigenvalues of (K2 , N2 ) are all even. 0 In # # ⊕ In . It is Proof Define K2 ≡ K2 and N2 ≡ N2 , where = In ⊕ In 0 #2 and N #2 are real skew-Hamiltonian; i.e., (K #2 J4n )T = −K #2 J4n easy to check that K T # # and (N2 J4n ) = −N2 J4n . Therefore, from the result of [18], it follows that the #2 ) are all even. #2 , N. multiplicities of eigenvalues of (K. 123.

(23) Palindromic quadratization and structure-preserving algorithm. 727. From (21), we see that μ is an eigenvalue of (K, N ) if and only if μ¯ is also an eigenvalue. We now give the relationship between eigenpairs of (K, N ) and (K2 , N2 ). x y Theorem 8 (i) If (α + ıβ, x + ı y) is an eigenpair of (K, N ), then ±ı y −x are eigenvectors (K2 , N2 ) corresponding to the eigenvalues α ± ıβ. of y x (ii) If 1 +ı 1 is an eigenvector of (K2 , N2 ) corresponding to the eigenvalue x2 y2 α + ıβ, then (x1 − y2 ) + ı(x2 + y1 ) is an eigenvector of (K, N ) corresponding to α + ıβ. Proof (i) Since K(x +ı y) = (α +ıβ)N (x +ı y), comparing the real and the imaginary parts of both sides leads to K2. . x x y y ±ı ±ı = (α ± ıβ)N2 . y y −x −x. y x1 +ı 1 is an eigenvector of (K2 , N2 ) corresponding to the eigenx2 y2 value α + ıβ, it holds that . (ii) Since. K R x − K I y = α(N R x − N I y) − β(N I x + N R y), K I x + K R y = β(N R x − N I y) + α(N I x + N R y), by setting x = x1 − y2 and y = x2 + y1 . Thus, (α + ıβ, x + ı y) is an eigenpair of (K, N ).. From Theorem 8, the eigenpairs of (K, N ) can be computed from the eigenpairs of #2 ). Since K #2 and N #2 are both real skew-Hamiltonian, based on Patel’s approach #2 , N (K # # [13,22], the pair (K2 , N2 ) can be reduced to block upper triangular forms #2 := Q T K #2 Z = K. . K 11 K 12 , T 0 K 11. #2 := Q T N #2 Z = N. . N11 N12 T , 0 N11. (29). T Z J , and K , N where Q, Z ∈ R4n×4n are orthogonal satisfying Q = J4n 4n 11 11 ∈ 2n×2n are upper Hessenberg and upper triangular, respectively. R From Theorem 7 and (29), we see that the pair (K 11 , N11 ) has the same spec#2 ). We then apply the QZ algorithm to (K 11 , N11 ) to compute all #2 , N trum as (K . 2n z˜ j 2n μ j , Z are the 2n eiits eigenpairs {(μ j , z˜ j )} j=1 . Consequently, 0 j=1 %T %T $ $ genpairs of (K2 , N2 ). Let μ j = α j + ıβ j and Z z˜ Tj , 0T ≡ x1Tj , x2Tj + %T $ ı y1Tj , y2Tj with α j , β j ∈ R and x1 j , x2 j , y1 j , y2 j ∈ R2n . From Theorem 8, '2n & α j + ıβ j , J T (x1 j − y2 j +ı(x2 j + y1 j )) j=1 are eigenpairs of (Ms , Ls ). Finally, we compute all eigenvalues and the associated eigenvectors of Q(λ) by Theorem 6. We present the structure-preserving algorithm for solving H-PQEP in Algorithm 1.. 123.

(24) 728. T.-M. Huang et al.. Algorithm 1 Structure-Preserving Algorithm (SPA) for H-PQEP Input: An H-palindromic quadratic pencil Q(λ) ≡ λ2 A1H +λA0 + A1 with A0 , A1 ∈ Cn×n and A0H = A0 ; Output: All eigenvalues and eigenvectors of Q(λ). #2 , N #2 ) = (K2 , N2 ) as in (28); 1: Form the matrix pair (K #2 , N #2 ) to block upper triangular forms as in (29); 2: Reduce (K 3: Compute eigenpairs {(μ j , z˜ j )}2n of (K 11 , N11 ) defined in (29) by the QZ algorithm; j=1 y x1 j z˜ j + ı 1 j , j = 1, 2, . . . , 2n; ≡ 4: Compute Z x2 j y2 j 0 5: Compute the eigenpair (μ j , z j ), for j = 1, 2, . . . , 2n, of (Ms , Ls ) by $ %T z j = J T (x1 j − y2 j + ı(x2 j + y1 j )) ≡ z Tj1 , z Tj2 ; 6: Compute ν j and ν1 by solving ν 2 − μ j ν + 1 = 0; Compute x j1 ≡ z j1 + ν1 z j2 and x j2 ≡ z j1 + ν j z j2 j j for j = 1, 2, . . . , 2n; 7: If x j1 = 0, then it is an eigenvector of Q(λ) corresponding to ν j ; If x j2 = 0, then it is an eigenvector of Q(λ) corresponding to ν1 ; j. 4 Balancing of P(λ) and Q(λ) Scaling [1,3,7,16] is a commonly used technique for standard eigenvalue problems for the improvement of the sensitivity of eigenvalues. In this section, we first propose a diagonal scaling for P(λ) in (1). Then, we determine the free parameters d1 , . . . , dm in (7) and (8) to improve the backward errors of eigenpairs for P(λ) as in [10,14,17]. In order to balance the entries of coefficient matrices in P(λ), we define a complex diagonal matrix D ≡ diag 2α1 , 2α2 +ıβ2 , . . . , 2αn +ıβn with α j , β j ∈ R so that the magnitudes of entries of coefficient matrices in the new (, ε)-palindromic matrix polynomial D. d−1 k=0. λ. 2d−k. Ad−k. + λ A0 + ε d. d . λ. d−k. Ak. D. k=1. are close to one as much as possible. That is, we determine α1 , . . . , αn and β2 , . . . , βn so that 2α j +ıβ j Ak ( j, )2α −ıβ ≈ 1,. (30). for j, = 1, 2, . . . , n and k = 0, 1, . . . , d, where Ak ( j, ) is the ( j, )-th entry of Ak . By taking logarithm of (30), the parameters, α1 , . . . , αn and β2 , . . . , βn can be determined by solving the least square problems α j + α = −Re(log2 (Ak ( j, ))), β j − β = −Im(log2 (Ak ( j, ))),. 123.

(25) Palindromic quadratization and structure-preserving algorithm. 729. where Re(c) and Im(c) represent the real and imaginary parts of c, respectively. Then, the parameters α1 , . . . , αn and β2 , . . . , βn are determined by the associated normal equations B T B [α1 , . . . , αn ]T = B T b, C T C [β2 , . . . , βn ]T = C T c. We now determine d1 , . . . , dm in (7) or (8), to balance the magnitudes of entries of A0 and A1 in Q(λ). For convenience, we define ( di =. di(1) , if 2d = 4m, for i = 1, . . . , m. (2) di , if 2d = 4m + 2;. From the row balancing of A1 in (7a) or (8a), we first set (s). ηi. = max {1, max{

(26) A2m−k+s−1

(27) 1 : k = 0, 1, . . . , 2i − 3 + s}}. (s). (s). Then we take δi to be the geometric average of ηi magnitudes of entries of A2m−2i+1 in A1 ; i.e., (s). δi. and the average of the absolute. ) ⎛ ⎞ * n n * * (s) ⎝ = +ηi |A2m−2i+1 ( j, )|/n 2 ⎠ j=1 =1. for i = 1, . . . , m and s = 1, 2. Although the value of di(s) can be set to δi(s) to balance the entries of A1 , we also need to consider the balance of the entries of both A0 and (s) (s) A1 in (7) or (8). As a result, we take the values of d1 , . . . , dm to be the geometric (s) (s) average of the maximal values of δ1 , . . . , δm and the maximal average for the absolute magnitudes of entries of Ak (k = 0, . . . , d); i.e., for i = 1, . . . , m and s = 1, 2, we set ) ⎧ ⎫ * n n * ⎨ ⎬ * |Ak ( j, )|/n 2 di(s) := +ρ (s) max ⎭ 0≤k≤d ⎩ j=1 =1. with (s). ρ (s) = max{δi ; i = 1, . . . , m}. 5 Numerical results In [13], an SPA is proposed for solving T-PQEPs. Numerical experiments show that SPA performs well on the T-PQEP arising from a finite element model of high-speed trains and rails. In this section, we shall focus on the numerical comparison of the. 123.

(28) 730. T.-M. Huang et al.. performance and accuracy for solving H-PPEP of even degree by using structurepreserving algorithms and companion linearization. For solving an n × n H-PPEP of even degree 2d, we apply the P-quadratization in Sect. 2 to transform it into a dn × dn H-PQEP. We then apply the SPA (Algorithm 1) in Sect. 3 to solve the H-PQEP. The combination of the P-quadratization and SPA is called the PQ_SPA algorithm. On the other hand, we can also use the “good” linearization [20,21] to transform the H-PPEP into a palindromic linear pencil λZ H + Z , ˆ + Z )x = 0 with λ = λˆ 2 . and then utilize SPA to solve the H-PQEP: (λˆ 2 Z H + λ0 The combination of the “good” linearization and SPA is called the PL_SPA algorithm. As mentioned in Remark 1 (ii), we see that applying the SPA to λˆ 2 Z T + λˆ 0 + Z is mathematically equivalent to applying the URV-based method [25] to λZ T + Z .. 5.1 Computational cost #2 ) in step 1 #2 , N For making PQ_SPA more efficient, we reorder the submatrices of (K of Algorithm 1 by the permutations ⎡ ⎤ In 0 0 0 0 ⎢ ⎢ 0 0 I(d−1)n 0 ⎢ 0 ⎥ ⎢ ⎢ ⎥ ⎢ ⎥ 0 0 I(d−2)n ⎥ , 2 = ⎢ 1 = ⎢ 0 0 ⎢ ⎢ ⎣0 0 0 In 0 ⎦ ⎣ 0 In 0 0 0 ⎡. ⎤ 0 0 In 0 0 0 0 0⎥ 0 I(d−2)n 0 0 ⎥ 0 0 0 0 0 In ⎥ ⎥. 0 0 0 0 I(d−2)n 0 ⎥ ⎥ 0 0⎦ 0 0 0 In 0 0 0 0 0 In. #2 and get We substitute A1 in Theorem 3 into N 1 0 # 2T 0 # N2 N2 := 0 2 0 1T ⎤ ⎡ ⎡ D1 0 0 0 D1 0 −V1 V2 ⎥ ⎢ ⎢ 0 D1 V2 0 0 0 D V 1 1 ⎥ ⊕⎢ = ⎢ ⎣ 0 0 V3 −V4 ⎦ ⎣ −V1T V2T V3T V4T 0 0 V4 V3 V2T V1T −V4T V3T. ⎤ ⎥ ⎥, ⎦. where D1 is a (2d − 2)n × (2d − 2)n diagonal matrix and V3 , V4 ∈ Rn×n . Set T #2 2 0T . #2 := 1 0 K K 0 2 0 1. 123.

(29) Palindromic quadratization and structure-preserving algorithm. 731. Table 1 Computational flops of PQ_SPA and PL_SPA Eigenvalues. Eigenvectors. Total. PQ_SPA. (408d 3 + 32d + 80/3)n 3. (408d 3 + 32d)n 3. (816d 3 + 64d + 80/3)n 3. PL_SPA. 3605 31 d 3 n 3. 1920d 3 n 3. 5525 31 d 3 n 3. Now, we compare the computational costs for PQ_SPA and PL_SPA: –. –. –. –. (Q R factorization and updating. with real arithmetic operations) Compute Q 1 and T (1) (1) (1) T #2 Z 1 = diag N , N , where N11 is upper trianguZ 1 such that Q 1 N 11 11 #2 Z 1 . It requires (80/3n 3 + 32dn 3 ) and 341 1 d 3 n 3 flops for lar, and update Q 1T K 3 PQ_SPA and PL_SPA, respectively. (Given’s rotations and updating with real arithmetic operations) Reducing the new #2 ) produced by above step to block upper triangular forms of (29), it #2 , N pair (K requires 232d 3 n 3 − (296d 2 − 24d)n 2 and 1856d 3 n 3 − 800d 2 n 2 flops for PQ_SPA and PL_SPA, respectively. (Computing eigenvalues of (K 11 , N11 )) Computing eigenvalues of the real upper Hessenberg and triangular pair (K 11 , N11 ) by QZ algorithm, it requires 176d 3 n 3 and 1408d 3 n 3 flops for PQ_SPA and PL_SPA, respectively, to obtain the upper quasi-triangular and triangular pair. #2 ) can be computed by an additional (408d 3 +32d)n 3 − #2 , N The eigenvectors of (K 2 2 (332d − 16d)n and 1920d 3 n 3 − 1088d 2 n 2 flops for PQ_SPA and PL_SPA, respectively. We summarize the computational flops of PQ_SPA and PL_SPA in Table 1.. 5.2 Numerical experiments For an approximate eigenpair (λ, x) of the palindromic matrix polynomial P(λ), we define the associated relative residual by

(30) P(λ)x

(31) 2 % . d+k + |λ|d−k

(32) A

(33) + |λ|d

(34) A

(35)

(36) x

(37) k 2 0 2 2 k=1 |λ|. RRes ≡ RRes(λ, x) := $! d. We will show numerical results of RRes and the reciprocal property of eigenpair (λ, x) for the H-PPEPs, computed by PQ_SPA, PL_SPA and polyeig in MATLAB (applied directly to (1)). As mentioned before, theoretically, the eigenvalues of H-PPEP appear in reciprocal ¯ So, if we sort the eigenvalues in ascending order by modulus, the prodpairs (λ, 1/λ). uct of the i-th eigenvalue and the conjugate of the (2dn + 1 − i)-th eigenvalue should be one. Therefore, we define the reciprocities of the computed eigenvalues by ri ≡ |λi λ¯ 2dn+1−i − 1| (i = 1, . . . , dn).. 123.

(38) 732. T.-M. Huang et al. −12. (b) 10. polyeig PL_SPA PQ_SPA −13. 10. Reciprocity r i. RRes of eigenpairs. (a) 10−12. −14. 10. −15. −13. 10. −14. 10. 10. −15. −16. 10. 10. 0. 0. 10. 10. | λi |. | λi |. Fig. 2 Relative residuals of eigenpairs and the associated reciprocity for Example 1. All numerical experiments are carried out using MATLAB 2008b with machine precision eps ≈ 2.22 × 10−16 . Let Cn,b denote the set of n × n complex matrices which real and imaginary parts are randomly generated by the normal distribution with zero mean and standard deviation b. Example 1 Consider the H-PPEP with d = 5 and Ak ∈ Cn,100 (k = 0, . . . , 5). Example 2 Consider the H-PPEP with d = 4 and A1 , A3 , A4 ∈ Cn,100 , and A0 and A2 being defined as Ak−1 = B1k · diag ϕ1(k) , . . . , ϕn(k) · B2k ∈ Cn×n (k = 1, 3) where B1k , B2k ∈ Cn,1 , and (. (k). (k). ϕi = 4i+k− , ϕ +i = 4i−k (i = 1, . . . , ), (k) if n is odd; ϕn = 4n/2−k. (31). with = n/2 (if n is even) or = (n − 1)/2 (otherwise). Example 3 Consider the H-PPEP with d = 4 and A0 , A1 , A3 , A4 ∈ Cn,100 , and A2 being defined as (3) A2 = B13 · diag ϕ1 , . . . , ϕn(3) · B23 ∈ Cn×n , (3). where B13 , B23 ∈ Cn,1 , and ϕi. is defined in (31) with k = 3.. We present the relative residuals (RRes) and the reciprocities of eigenpairs computed by the polyeig, PL_SPA and PQ_SPA for Examples 1–3, using the balancing technique in Sect. 4 with n = 30. Numerical results are shown in Figs. 2, 3 and 4. We indicate the results computed by polyeig, PL_SPA and PQ_SPA by “×”, “” and “o”, respectively. For the PL_SPA and PQ_SPA, all reciprocities of eigenvalues. 123.

(39) Palindromic quadratization and structure-preserving algorithm. (a) 10. 0. i. 10. Reciprocity r. RRes of eigenpairs. (b) 102. polyeig PL_SPA PQ_SPA. −5. 733. −10. 10. −2. 10. −4. 10. −6. 10. −15. 10. −8. −4. 10. −2. 10. 0. 10. 2. 10. 10. 4. −4. 10. 10. −2. 10. | λi |. 0. 10. 2. 10. 4. 10. |λ | i. Fig. 3 Relative residuals of eigenpairs and the associated reciprocity for Example 2 with larger

(40) A0

(41) 2 and

(42) A2

(43) 2. (a). −5. −8. polyeig PL_SPA PQ_SPA. −10. −6. 10 i. 10. Reciprocity r. RRes of eigenpairs. (b) 10. 10. −12. 10. −14. 10. −8. 10. −9. −16. 10. 10. −10. −18. 10. −7. 10. −4. 10. −2. 10. 0. 10. |λ | i. 2. 10. 4. 10. 10. −4. 10. −2. 10. 0. 10. 2. 10. 4. 10. |λ | i. Fig. 4 Relative residuals of eigenpairs and the associated reciprocity for Example 3 with larger

(44) A2

(45) 2. are preserved to machine accuracy, which are ignored in Figs. 2b–4b. From Figs. 2–4, we see that most of relative residuals of eigenpairs computed by the PQ_SPA are better than that computed by the PL_SPA, only a few exceptions. Overall, we conclude that applying P-quadratization and SPA (Algorithm 1) to solve PPEPs not only preserves the reciprocal property but also provides higher accuracy than that by PL_SPA and polyeig in MATLAB. 6 Conclusions In this paper, we mainly propose a palindromic quadratization to transform the (, ε)palindromic matrix polynomial of even degree with (, ε) = (T, −1) to a (, ε)-palindromic quadratic pencil, instead of the orthodox palindromic linearization approach. The structure-preserving algorithm for solving palindromic quadratic eigenvalue problem based on (S+S −1 )-transform and Patel’s algorithm can then be applied. Numerical experiments show that relative residuals of approximate eigenpairs for the palindromic polynomial eigenvalue problem computed by the PQ_SPA are better than those by the PL_SPA and polyeig in MATLAB. Moreover, the computational cost for PQ_SPA is much cheaper than that for PL_SPA.. 123.

(46) 734. T.-M. Huang et al.. Acknowledgments The authors are grateful to Professor Tisseur and anonymous referees for many valuable comments and helpful suggestions. This work is partially supported by the National Science Council and the National Center for Theoretical Sciences in Taiwan.. References 1. Betcke, T.: Optimal scaling of generalized and polynomial eigenvalue problems. SIAM J. Matrix Anal. Appl. 30, 1320–1338 (2008) 2. Byers, R., Mackey, D.S., Mehrmann, V., Xu, H.: Symplectic, BVD, and palindromic approaches to discrete-time control problems. Technical report, Preprint 14-2008, Institute of Mathematics, Technische Universität Berlin (2008) 3. Chen, T.-Y., Demmel, J.W.: Balancing sparse matrices for computing eigenvalues. Linear Algebra Appl. 309, 261–287 (2000) 4. Chiang, C.-Y., Chu, E.K.-W., Li, T., Lin, W.-W.: The palindromic generalized eigenvalue problem A∗ x = λAx: numerical solution and applications. Linear Algebra Appl. (2010, to appear) 5. Chu, E.K.-W., Huang, T.-M., Lin, W.-W.: Structured doubling algorithms for solving g-palindromic quadratic eigenvalue problems. Technical report, NCTS Preprints in Mathematics 2008-4-003, National Tsing Hua University, Hsinchu, Taiwan (2008) 6. Chu, E.K.-W., Hwang, T.-M., Lin, W.-W., Wu, C.-T.: Vibration of fast trains, palindromic eigenvalue problems and structure-preserving doubling algorithms. J. Comput. Appl. Math. 219, 237–252 (2008) 7. Fan, H.-Y., Lin, W.-W., Van Dooren, P.: Normwise scaling of second order polynomial matrices. SIAM J. Matrix Anal. Appl. 26, 252–256 (2004) 8. Grammont, L., Higham, N.J., Tisseur, F.: A framework for analyzing nonlinear eigenproblems and parametrized linear systems. Technical report, The MIMS Secretary, School of Mathematics, the University of Manchester. MIMS EPrint: 2009.51 (2009) 9. Higham, N.J., Tisseur, F., Van Dooren, P.M.: Detecting a definite hermitian pair and a hyperbolic or elliptic quadratic eigenvalue problem, and associated nearness problems. Linear Algebra Appl. 351–352, 455–474 (2002) 10. Highman, N.J., Li, R.-C., Tisseur, F.: Backward error of polynomial eigenproblems solved by linearization. SIAM J. Matrix Anal. Appl. 29, 1218–1241 (2007) 11. Hilliges, A.: Numerische Lösung von quadratischen eigenwertproblemen mit Anwendungen in der Schiendynamik. Master’s thesis, Technical University Berlin, Germany, July 2004 12. Hilliges, A., Mehl, C., Mehrmann, V.: On the solution of palindramic eigenvalue problems. In: Proceedings 4th European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS). Jyväskylä, Finland (2004) 13. Huang, T.-M., Lin, W.-W., Qian, J.: Structure-preserving algorithms for palindromic quadratic eigenvalue problems arising from vibration on fast trains. SIAM J. Matrix Anal. Appl. 30, 1566–1592 (2009) 14. Huang, T.-M., Lin, W.-W., Su, W.-S.: Palindromic quadratization and structure-preserving algorithm for palindromic matrix polynomials of even degree. Technical report, NCTS Preprints in Mathematics, National Tsing Hua University, Hsinchu, Taiwan, 2009-6-002 (2009) 15. Ipsen, C.F.: Accurate eigenvalues for fast trains. SIAM News, 37 (2004) 16. Lemonnier, D., Van Dooren, P.: Balancing regular matrix pencils. SIAM J. Matrix Anal. Appl. 28, 253– 263 (2006) 17. Li, R.-L., Lin, W.-W., Wang, C.-S.: Structured backward error for palindromic polynomial eigenvalue problems. Technical report, NCTS Preprints in Mathematics, National Tsing Hua University, Hsinchu, Taiwan, 2008-7-002 (2008) 18. Lin, W.-W.: A new method for computing the closed-loop eigenvalues of a discrete-time algebraic Riccatic equation. Linear Algebra Appl. 96, 157–180 (1987) 19. Mackey, D.S., Mackey, N., Mehl, C., Mehrmann, V.: Numerical methods for palindromic eigenvalue problems: computing the anti-triangular Schur form. Numer. Linear Algebra Appl. 16, 63–86 (2009) 20. Mackey, D.S., Mackey, N., Mehl, C., Mehrmann, V.: Structured polynomial eigenvalue problems: good vibrations from good linearizations. SIAM J. Matrix Anal. Appl. 28, 1029–1051 (2006) 21. Mackey, D.S., Mackey, N., Mehl, C., Mehrmann, V.: Vector spaces of linearizations for matrix polynomials. SIAM J. Matrix Anal. Appl. 28, 971–1004 (2006) 22. Patel, R.V.: On computing the eigenvalues of a symplectic pencil. Linear Algebra Appl. 188, 591– 611 (1993). 123.

(47) Palindromic quadratization and structure-preserving algorithm. 735. 23. Schröder, C.: SKURV: a Matlab toolbox for the skew URV decomposition of a matrix triple. http:// www.math.tu-berlin.de/~schroed/Software/skurv 24. Schröder, C.: A QR-like algorithm for the palindromic eigenvalue problem. Technical report, Preprint 388, TU Berlin, Matheon, Germany (2007) 25. Schröder, C.: URV decomposition based structured methods for palindromic and even eigenvalue problems. Technical report, Preprint 375, TU Berlin, MATHEON, Germany (2007) 26. Tisseur, F., Meerbergen, K.: A survey of the quadratic eigenvalue problem. SIAM Rev. 43, 234–286 (2001) 27. Xu, H.: On equivalence of pencils from discrete-time and continuous-time control. Linear Algebra Appl. 414, 97–124 (2006) 28. Zaglmayr, S.: Eigenvalue problems in saw-filter simulations. Diplomarbeit, Institute of Computational Mathematics, Johannes Kepler University Linz, Linz, Austria (2002). 123.

(48)