Introduction

In this section, we extend from PQEP(1.7) to (⋆, ε)-palindromic matrix polyno-mials of even degree 2d

P(λ) ≡

d−1

∑

k=0

λ^2d−kA^⋆_d−k + λ^dA₀+ ε

∑d k=1

λ^d−kA_k, (2.1)

where d ≥ 2, ε = ±1 and ⋆ = H (Hermitian) or T (transpose), Ak ∈ C^n×n (k = 0, 1, . . . , d) and A^⋆₀ = εA₀. The corresponding polynomial eigenvalue prob-lemP(λ)x = 0, with λ ∈ C and x ∈ Cⁿ\{0} being the eigenvalue and the associated eigenvector respectively, is called a (⋆, ε)-palindromic polynomial eigenvalue problem ((⋆, ε)-PPEP). The equation (2.1) is also called a ⋆-PPEP if ε = 1 or a ⋆-anti-PPEP if ε =−1.

The underlying matrix polynomialP(λ) in (2.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(λ) = ελ^2dP(1/λ)^⋆ (2.2)

and explains the world “(anti-)palindromic” [40]. Consequently, taking the (conju-gate) transpose of (2.1), we easily see that the eigenvalues of P(λ) satisfy a “recip-rocal” property, that is, they appear in the pairs of the form (λ, 1/λ^⋆).

The (⋆, ε)-PPEPs arise in solving higher order systems of ordinary or partial dif-ferential equations. In the beginning, a T-palindromic quadratic eigenvalue problem (T-PQEP) is raised in the vibration analysis for fast trains in Germany [20, 21] and then in the study of surface acoustic wave filters [65]. 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 [23]. Furthermore, a ⋆-PPEP of even degree is obtained by solving the lin-ear quadratic discrete-time optimal control problem for higher order systems [2, 63].

A standard approach for computing the eigenpairs ofP(λ) in (2.1) is to linearize it to a 2dn× 2dn linear matrix pencil by the companion linearization and compute its generalized Schur form [55]. Nevertheless, the reciprocal property of eigenvalues of P(λ) is not preserved and this has resulted in large numerical errors [6, 24, 29].

Recently, in order to preserve the reciprocity of eigenvalues, some forerunner’s works [40, 41] propose some good linearizations that linearize P(λ) to palindromic linear pencils of the form λZ^⋆ + Z . This does lead to a great improvement over pre-vious unstructured approaches, keeping the palindromic structure in the original polynomial and enabling structure-preserving numerical methods to be designed.

Later, a QR-like algorithm [53] and a hybrid method [42] 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 QR-like algorithm typically requires O(n⁴) flops and the hybrid method re-quires O(n³log(n)) flops. Recently, a URV-decomposition based structured method of cubic complexity is developed in [54] to solve T-palindromic linear eigenvalue problems, producing eigenvalues which are reciprocally paired to working precision.

A new structure-preserving doubling algorithm with cubic complexity for solving

⋆-palindromic linear eigenvalue problems is developed in [7]. On the other hand, for solving a (⋆, ε)-PQEP, a structure-preserving doubling algorithm is developed in [6, 8] 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⁻¹)-transform [33] and Patel’s algorithm [45], is proposed in [24]

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.

The purpose of this section 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 [24] 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⁻¹)-transform and enlarge the H-skew-Hamiltonian 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 [54, 58] to the linearized T-palindromic linear pencil.

This section is organized as follows. In Section 2.2, we propose a palindromic quadratization for a (⋆, ε)-palindromic matrix polynomial of even degree. In Sec-tion 2.3, we develop a structure-preserving algorithm for solving the H-PQEPs. We consider the structured backward stability in Section 2.5. After that We develop bal-ancing techniques for PPEPs and PQEPs in Section 2.5. Comparisons of numerical results computed by the palindromic quadratization, the palindromic linearization and the standard companion linearization are presented in Section 2.6. Conclusions are given in Section 2.7.

2.2 P-Quadratization of (⋆, ε)-PPEP

In [24], 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 Section 2.3. As the 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 [24] is applicable.

Next we present definitions of quadratization and P-quadratization of a general matrix polynomial and a palindromic matrix polynomial, respectively.

Definition 2.1. (Quadratization/P-Quadratization)

(i) LetP(λ) 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 I_q−ν



 . (2.3)

(ii) Let P(λ) be an arbitrary ν × ν (⋆, ε)-palindromic matrix polynomial of degree p ≥ 2 with pν = 2q (i.e., P(λ) = ελ^pP(1/λ)^⋆ as in (2.2)). A quadratization Q(λ) of P(λ) having the (⋆, ε)-palindromic structure is called a P-quadratization of P(λ).

Theorem 2.2. 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(

[λ²Q(_λ¹)]

λ=0

) = 0) if and only if ∞ is

an eigenvalue of P(λ) (i.e., det(

[λ^pP(_λ¹)]

λ=0

) = 0).

Proof. (i) The factorization (2.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 λ₀.

(ii) Since

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 (2.3) implies the existence of a more wide class of one-sided factorization

where F (λ) and G(λ) are matrix rational functions of size q× ν. From the factor-ization (2.4) a close connection between eigenpairs of P(λ) and eigenpairs of Q(λ) has been shown in [17].

Theorem 2.3. [17] Assume that (2.4) holds at λ₀ ∈ C with F (λ0) and G(λ₀) being of full column rank. Then F (λ₀)z₁ is an eigenvector of Q(λ) if and only if z1 is an eigenvector of P(λ), both corresponding to eigenvalue λ0.

In Definition 2.1(i), we give a new definition of quadratization for a general matrix polynomial. In Theorems 2.2 and 2.3, 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.

Theorem 2.4. Let P(λ) be an n × n (⋆, ε)-palindromic matrix polynomial of degree 2d as in (2.1) with (⋆, ε) ̸= (T, −1). Then P(λ) can be P-quadratized into a (⋆,

ε)-palindromic quadratic pencil of the form

(ii) (For 2d = 4m + 2) A1 and A0 are given by

and for j = 1, . . . , m− 1:

and let

and

for j = 1, . . . , m− 1. Finally, letting

for j = 1, . . . , m. Via careful calculation we get

Note that the P-quadratization of a (⋆, ε)-palindromic matrix polynomial of odd degree with even matrix dimension can be defined as in Definition 2.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 (2.5) and P(λ) in (2.1).

Theorem 2.5. LetQ(λ) in (2.5) be a P-quadratization of P(λ) in (2.1) with (⋆, ε) ̸=

(T,−1). Denote z = [

z^T₁ · · · zd^T

]T

with z_j ∈ Cⁿ (j = 1, . . . , d). Then

(i) For λ₀ ̸= 0, (λ0, z₁) is an eigenpair ofP(λ) if and only if (λ0, z) is an eigenpair of Q(λ) with zj = F_j,1(λ₀)z₁ (j = 2, . . . , d), where {Fj,1(λ)}^dj=2 are given in (2.9) (for d = 2m) or (2.14) (for d = 2m + 1).

(ii) (0, z₁) 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)







z_2m=−(d^⋆m)⁻¹A₁z₁, z_2j+1 = 0, z_2j =−(d^⋆j)⁻¹A_2m−2j+1z₁;

(2.15a)

(or, for d = 2m + 1)







z_2m+1=−(d^⋆m)⁻¹A₁z₁, z_2m = 0, z2j = 0, z2j+1 =−(d^⋆j)⁻¹A2m−2j+1z1.

(2.15b)

(iii) (∞, z2) is an eigenpair of P(λ) if and only if (∞, z) is an eigenpair of Q(λ), with z₁ = z₃ =· · · = zd= 0.

Proof. (i) From Theorem 2.4, 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 (λ₀) and G(λ₀) defined in (2.4) are of full rank. The assertion in (i) follows

imme-diately from Theorem 2.3 and the relation

and for d = 2m + 1:







q1 = A^⋆₁z1+∑m

k=1A^⋆_2m−2k+3z2k, q2 =−εd1z3, q3 = d^⋆₁z4,

· · · , q2m =−εdmz_2m+1, q_2m+1= d^⋆_mz₁.

(2.17b)

From (2.17), it follows that (∞, z2) is an eigenpair of P(λ); i.e., A^⋆dz₂ = 0 if and only if (∞, z) is an eigenpair of Q(λ), or [

λ²Q(₁

λ

)]

λ=0z = q = 0, with z₁ = z₃ =

· · · = zd= 0.

Remark 2.6. (i) Theorem 2.4 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 [40] to transform it into a T-palindromic linear pencil λZ^T + Z, and then solve it by the QR-like algorithm [53], the hybrid method [42], the URV-based method [54]

or the doubling algorithm [7].

(ii) On the other hand, if we rewrite λZ^T + Z to a T-palindromic quadratic pencil bQ(bλ) ≡ bλ²Z^T + bλ0 + Z by letting bλ² = λ, then the SPA algorithm [24] can also be used to solve its eigenpairs. It is shown in [24] that applying the SPA to solve Q(bλ)y = 0 is mathematically equivalent to applying the URV-based method to solveb λZ^T + Z.

Applying the P-quadratization in Theorem 2.4, an H-anti-PPEP can be quadra-tized 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: (λ²A^H₁ +λA0−A1)x = 0, with A^H₀ =−A0. Then (ıω, x) is an eigenpair of the H-anti-PQEP if and only if (ω, x) is an eigenpair of the H-PQEP: [

ω²(−A1)^H + ω(ıA₀) + (−A1)]

x = 0.

Proof. The result can be easily obtained by setting λ = ıω and using the fact (ıA₀)^H = ıA₀.

We now consider the ⋆-even and ⋆-odd polynomial eigenvalue problems of even degree. Let C(λ) = ∑2d

k=0λ^kC_k, where C_k ∈ C^n×n (k = 0, 1, . . . , 2d) and C_2d ̸= 0.

The polynomial eigenvalue problem C(λ)x = 0 is called a ⋆-even polynomial eigen-value problem, if C_2k^⋆ = C_2k(k = 0, 1, . . . , d) and C_2k^⋆ ₋₁ =−C2k−1(k = 1, . . . , d); and it is called a ⋆-odd polynomial eigenvalue problem, if C_2k^⋆ = −C2k (k = 0, 1, . . . , d) and C_2k^⋆ ₋₁ = C2k−1 (k = 1, . . . , d). By the Cayley transformation, it was shown in [40] that a ⋆-even/odd polynomial eigenvalue problem can be transformed to a (⋆,±1)-PPEP, respectively.

Figure 2.1: Relations between various structured polynomial eigenvalue problems (PEPs).

We illustrate the relationship among various structured polynomial eigenvalue problems in Figure 2.1. We see that even, odd, anti-palindromic and T-palindromic polynomial eigenvalue problems of even degree can be P-quadratized to PQEPs. Thus, the SPA algorithm in [24] can be applied to solve the associated T-PQEPs. On the other hand, we see that even, odd, 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.

2.3 ⋆-palindromic Quadratic Eigenvalue Problems

Consider the ⋆-PQEP

Q(λ)x ≡ (λ²A^⋆₁+ λA₀ + A₁)x = 0, (2.18)

where A₀, A₁ ∈ C^n×n with A^⋆₀ = A₀. The eigenvalues of Q(λ) clearly appear in the

“reciprocal” pairs of the form (λ, 1/ˆλ) where

λ = ¯ˆ λ if ⋆ = H, λ = λ if ⋆ =ˆ ⊤.

Classical linearizations of (2.18) in a companion form, generally, do not preserve the symplectic structure. Fortunately, the special linearization of (2.18) (see [8] or [24])

where y = ¹_λA₁x and multiplying the second equation of (2.18) satisfies

MJ M^⋆ =LJ L^⋆, J = J2n ≡

so that the matrix pair (M, L) has eigenvalues λ and 1/ˆλ, preserving reciprocity.

The pencil M − λL or the matrix pair (M, L) are called H-symplectic.

For any real symplectic matrix pair (M, L) satisfying (2.20), a structure-preserving (S + S⁻¹)-transform for the computation of all eigenvalues was proposed by [33].

The (S + S⁻¹)-transform (Ms,Ls) of an ⋆-symplectic matrix pair (M, L) is defined

by

Ms ≡ MJ L^⋆+LJ M^⋆ ≡ KJ , Ls ≡ LJ L^⋆ ≡ N J . (2.21)

It is easily seen that K and N are both ⋆-skew-Hamiltonian, i.e., KJ = J K^⋆ and N J = J N^⋆. Hence, if µ is an eigenvalue of (K, N ), so is ˆµ.

The relationship between eigenpairs of an ⋆-symplectic matrix pair and its (S + S⁻¹)-transform has been given in [33].

Theorem 2.7. [33] Let (M, L) be ⋆-symplectic and (Ms,Ls) be its (S + S⁻¹ )-transform as in (2.21). Suppose z_s is an eigenvector of (Ms,Ls) corresponding to µ = ν +¹_ν. If (L^⋆−_ν¹M^⋆)z_s̸= 0 or (L^⋆− νM^⋆)z_s̸= 0, then (ν, J (L^⋆−¹_νM^⋆)z_s) or (1/ν,J (L^⋆− νM^⋆)z_s) is an eigenpair of (M, L), respectively.

We have some different derivative theorems between ⋆ =⊤ [24] and ⋆ = H here.

Then the relationship between eigenpairs of (Ms,Ls) and Q(λ), ⋆ = H is given as follows.

Theorem 2.8. Let (M, L) be H-symplectic of the form in (2.19) and (Ms,Ls) be its (S +S⁻¹)-transform. Suppose that z_s is an eigenvector of (Ms,Ls) corresponding to the eigenvalue µ̸= ±2, and denote zs ≡ [zs1^T, z^T_s2]^T with z_s1, z_s2 ∈ Cⁿ. Let ν be a root of the quadratic equation λ + ¹_λ = µ. Then

(i) at least one of vectors z_s1+_ν¹z_s2 and z_s1+ νz_s2 is nonzero;

(ii) if z_s1+_ν¹z_s2 ̸= 0, then zs1+ ¹_νz_s2 is an eigenvector of Q(λ) corresponding to ν;

(iii) if z_s1+¹_νz_s2 = 0, then z_s2 is an eigenvector of Q(λ) corresponding to ¹_ν.

Proof. (i) Suppose that zs1+¹_νzs2 = 0 and zs1+νzs2 = 0. It implies that (ν−¹_ν)zs2 = 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, zs1+ ¹_νz_s2 ̸= 0 or zs1+ νz_s2 ̸= 0.

(ii) Since MJ M^H =LJ L^H, by (2.21) it holds that

0 = (Ms− µLs) z_s= (M − νL)J (L^H − 1

νM^H)z_s. (2.22)

From (2.22), we obtain

(M − νL)



 z_s1+¹_νz_s2 x_ν



 = 0, (2.23)

where

x_ν ≡ 1

νA^H₁ z_s1− 1

νA₀z_s2− A1z_s2. (2.24)

Substituting (M, L) of (2.19) into (2.23), we have

xν = 1

νA1(zs1+ 1

νzs2) (2.25)

and

A₀(z_s1+ 1

νz_s2) + x_ν + νA^H₁ (z_s1+ 1

νz_s2) = 0. (2.26)

Substituting x_ν of (2.25) into (2.26) and multiplying (2.26) by ν, we get Q(ν)(zs1+

1

νz_s2) = 0.

(iii) Since z_s1+ ¹_νz_s2 = 0, it follows that z_s1 = −_ν¹z_s2 ̸= 0 and xν = 0 in (2.25).

Substituting these results into (2.24), it holds that

0 = xν =− (1

ν )2

A^H₁ zs2− 1

νA0zs2− A1zs2.

Therefore, z_s2 is an eigenvector of Q(λ) corresponding to the eigenvalue _ν¹.

In [24], a structure-preserving algorithm (SPA) based on Patel’s algorithm [45]

has been developed for solving T-PQEPs. In order to solve an H-PQEP, we apply the (S + S⁻¹)-transform to the H-symplectic pair (M, L) of the form (2.19) and get the generalized eigenvalue problemKzs= µN zs, where K and N are defined in (2.21). Substituting (M, L) in (2.19) into (2.21), the H-skew-Hamiltonian K and N can be represented as

K =

However, Patel’s algorithm can only be applied to (K, N ) of (2.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 (2.27) into an enlarged real T-skew-Hamiltonian pair so that Patel’s algorithm can be applied. We extend (K, N ) in (2.27) to a real 4n× 4n matrix pair (K2,N2) by

Theorem 2.9. The multiplicities of eigenvalues of (K2,N2) are all even.

Proof. Define eK2 ≡ ΠK2Π and eN2 ≡ ΠN2Π, where Π = I_n⊕

From (2.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).

Theorem 2.10. (i) If (α + ıβ, x + ıy) is an eigenpair of (K, N ), then

Proof. (i) SinceK(x+ıy) = (α+ıβ)N (x+ıy), comparing the real and the imaginary parts of both sides leads to

K2

eigenvalue α + ıβ, it holds that

KRx− KIy = α(NRx− NIy)− β(NIx +NRy), KIx +KRy = β(NRx− NIy) + α(NIx +NRy),

by setting x = x₁ − y2 and y = x₂ + y₁. Thus, (α + ıβ, x + ıy) is an eigenpair of (K, N ).

From Theorem 2.10, the eigenpairs of (K, N ) can be computed from the eigen-pairs of ( eK2, eN2). Since eK2 and eN2 are both real skew-Hamiltonian, based on Patel’s

approach [24, 45], the pair ( eK2, eN2) can be reduced to block upper triangular forms are upper Hessenberg and upper triangular, respectively.

From Theorem 2.9 and (2.29), we see that the pair (K₁₁, N₁₁) has the same and the associated eigenvectors of Q(λ) by Theorem 2.8. We present the structure-preserving algorithm for solving H-PQEP in Algorithm 2.1.

Algorithm 2.1 Structure-Preserving Algorithm (SPA) for H-PQEP

Input: An H-palindromic quadratic pencil Q(λ) ≡ λ²A^H₁ + λA₀+ A₁ with A₀, A₁ ∈ C^n×n and A^H₀ = A₀;

Output: All eigenvalues and eigenvectors of Q(λ).

1: Form the matrix pair ( eK2, eN2) = (ΠK2Π, ΠN2Π) as in (2.28); it is an eigenvector of Q(λ) corresponding to _ν¹_j;

2.4 Structured Backward Perturbation Analysis

Typically, algorithms would approach the right solution in the limit, if there were no round-off or truncation errors. However, depending on the specific computational method, errors can be magnified and causing the error to grow exponentially.Let {µ, z} be a computed eigenpair of

( _d

Theoretically, we would like to have (∑_d

ℓ=0

A_ℓµ^ℓ )

z = 0. But practically we have (∑_d

ℓ=0

Aℓµ^ℓ )

z = −r with residual r ̸= 0 but usually tiny. Backward error analysis asks if the computed eigenpair {µ, z} is an exact eigenpair of a nearby PEP such

that ( _d

where△Aℓare called backward perturbation matrices and ρ_ℓare scaling parameters.

Tisseur [56] developed a backward error perturbation analysis for PEP generally, where△Aℓ is no structure. However, as A are Hermitian, we would like to enforcing that△Aℓ should be Hermitian, too. This is the reason why the structured backward perturbation analysis be developed. That is to say, we consider PPEP with

( _d

We have mentioned the fast train application [29, 43? ] which yield a problem of this form with d = 2,⋆ =⊤, and ε = 1 before. Let ∥·∥ be either the spectral norm

∥·∥₂ or the Frobenius norm ∥·∥_F.Now, we are interested in knowing the structrued

backward error

△ = min

√∑_⌊d/2⌋

ℓ=0 ∥△Aℓ∥²,

subject to (∑d

ℓ=0(A_ℓ+ ρ_ℓ△Aℓ) µ^ℓ )

z = 0, (2.31)

scaling parameters ρ_d−ℓ = ρ_ℓ ≥ 0

A_d−ℓ+ ρ_d−ℓ△Ad−ℓ = ε (A_ℓ+ ρ_ℓ△Aℓ)^⋆ for ℓ = 0, 1, ...⌊d/2⌋

The constraints in (2.31) require ρ_d−ℓ = ρ_ℓ and △Ad−ℓ = ε△A^⋆d−ℓ just like the palindromic form, and we decompose (2.30) which is equivalent to

( _d

∑

ℓ=0

ρ_ℓ△Aℓµ^ℓ )

z = r def

= − ( _d

∑

ℓ=0

A_ℓµ^ℓ )

z, (2.32)

where △Aℓ ∈ C^n×n(ℓ = 0, 1, ...d) . We will seek if (2.32) has a solution {△Aℓ, ℓ = 0, 1, ...⌊d/2⌋} and if it does, we’ll seek △Aℓ such that

min

⌊d/2⌋∑

ℓ=0

∥△Aℓ∥² (2.33)

where ∥·∥ is either ∥·∥₂ or∥·∥_F.

The key to solve (2.32) and (2.33) is a reduction technique that has been used in [4, 25, 35, 57]. The technique allows us to consider (2.32) in the case of 2×2 reduced matrix when we deal with backward errors in a 2-dimensional subspace spanned by {z, r} or by {z, ¯r}.

Let Q∈ C^n×n be a unitary matrix (QQ^H = I_n), such that

Such Q always can be generated, and the equation (2.34) also implies

|α| = ∥z∥₂, γ = z^Hrˆ

We can rewrite (2.32) by multiplying Q , then

Q^⋆

β = β and ˆˆ γ = γ if ⋆ = H, and ˆβ = ¯β and ˆγ = ¯γ if ⋆ = ⊤.

Since Q is unitary, (2.31) and (2.36) have the same solvability property. It means that if the former is solvable, so is the latter, and moreover

⌊d/2⌋∑

ℓ=0

∥△Bℓ∥² =

⌊d/2⌋∑

ℓ=0

∥△Aℓ∥². (2.37)

Thus we can focus on the optimal solution{△Bℓ, ℓ = 0, 1, ...,⌊d/2⌋} to (2.37) which means

min

⌊d/2⌋∑

ℓ=0

∥△Bℓ∥². (2.38)

That also generate one solution {△Aℓ, ℓ = 0, 1, ...,⌊d/2⌋} to (2.32) in the sence of (2.33), and vice versa. Furthemore, we transfrom (2.31) to the reduced structrued backward error

△p = min





 vu ut^⌊d/2⌋∑

ℓ=0

∥△Bℓ∥²_p : (2.36) satisfied





 for p =2,F

It follows from (2.35) that

|δ1| = z^Hrˆ

∥z∥²₂ ,|δ2| =

√∥ˆr∥²₂∥z∥²₂− |z^Hˆr|²

∥z∥²₂ ,

√

|δ1|²+|δ2|² = ∥r∥₂

∥z∥₂ (2.39)

The complex calculations and technical operation in Theorem 2.12 [47] are omit and we present the significant results.We first define a few parameters in term of a

given approximate eigenpair {µ, z} of PPEP (2.30): for even d

Theorem 2.11. [47] Let {µ, z} be a given approximate eigenpair of PPEP (2.30).

Suppose ⋆ = H and ε = ±1 in (2.31), and δ1 and δ₂ are as in (2.39) with ˆr = r which is defined in (2.32). Let

ϕ =







ϕeven for even d ϕodd, for odd d

Φ =







Φeven for even d Φodd, for odd d

Ψ =







Ψeven for even d Ψodd, for odd d

Theorem 2.12. For the structure backward error △F defined in (2.31), we have

1. If |µ| = 1 and z^Hr/(√

εµ^d/2) /∈ R, then △F = +∞

2. If |µ| = 1 and z^Hr/(√

εµ^d/2)∈ R, then

△F =

√

|δ1|²

Φ + |δ2|² Ψ

3. If |µ| ̸= 1, then

△F ≤

√|δ1|²

ϕ + |δ2|² Ψ

Let us look back at our original problems (2.1) and apply the Theorem 2.12 in our PPEP. Taking d = 4, ε = 1 for example, we consider

P(λ) = λ⁴A^H₂ + λ³A^H₁ + λ²A0+ λA1 + A2. (2.40)

From the quadration we mentioned in Theorem 2.4 where A1 and A0 are given by

A1 =



 A₁ I A₂ 0



 , A⁰ =



 A₀− I − A^H2 A₂ 0

0 −I





which satisties

Q(λ) = λ²A^H1 + λA0+A1. (2.41)

In section 2, we propose a method to solve the equation (2.41) with structure-preserving algorithm. In this case, it is quite reasonable to use structure-structure-preserving backward error to estimate the stability. First, we select the equation (2.40) and transform it to (2.41). If (λ₀, z₁) is an eigenpair ofP(λ), then (λ0, z) is an eigenpair of Q(λ) with z = [z₁^⊤, z₂^⊤]^⊤ where

z₂ = 1 λ₀

(λ²₀A^H2 z₁+ z₁) .

By the structure-preserving backward error analysis, we use the equivalent

no-tation on Q(λ) and let {µ, z} = {

λ₀, [z^⊤₁, z₂^⊤]^⊤}

. The theorem 2.12 gives us an estimated backward error upper bound and the bound is △F ≤ √

|δ1|²

ϕ +^|δ_Ψ^2|2. The terms ϕ and Ψ just relate to λ0 which is preserved in our algorithm. In spite of this, the term δ1 contains the eigenvector z which is changed in our process (|δ1| = |^zHrˆ|

∥z∥²2

).

More precisely, the norm ofA2cause considerable impact on our structure-preserving backward error. In general, the norm of Aℓ, ℓ = 0, 1, ...⌊d/2⌋ become perturbative factors that we have to face it. Next section we provide a classical technique that modify our backward error and easily implement in our algorithm.

2.5 Balancing of P(λ) and Q(λ)

Scaling [3, 9, 15, 34] is a commonly used technique for standard eigenvalue prob-lems for the improvement of the sensitivity of eigenvalues. In this section, we first propose a diagonal scaling forP(λ) in (2.1). Then, we determine the free parameters d₁, . . . , d_m in (2.7) and (2.8) to improve the backward errors of eigenpairs for P(λ) as in [26, 27, 47].

In order to balance the entries of coefficient matrices inP(λ), 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−kA^⋆_d−k+ λ^dA₀+ ε

∑d k=1

λ^d−kA_k )

D^⋆

are close to one as much as possible. That is, we determine α₁, . . . , α_nand β₂, . . . , β_n

so that

2^αj+ıβjA_k(j, ℓ)2^{αℓ−ıβℓ} ≈ 1, (2.42)

for j, ℓ = 1, 2, . . . , n and k = 0, 1, . . . , d, where A_k(j, ℓ) is the (j, ℓ)-th entry of A_k. By taking logarithm of (2.42), the parameters, α₁, . . . , α_n and β₂, . . . , β_n can be determined by solving the least square problems

αj + αℓ =−Re(log2(Ak(j, ℓ))), βj− βℓ =−Im(log2(Ak(j, ℓ))),

where Re(c) and Im(c) represent the real and imaginary parts of c, respectively.

Then, the parameters α₁, . . . , α_n and β₂, . . . , β_n are determined by the associated normal equations

B^TB [α₁,· · · , αn]^T = B^Tb, C^TC [β₂,· · · , βn]^T = C^Tc.

We now determine d₁, . . . , d_m in (2.7) or (2.8), to balance the magnitudes of entries of A0 and A1 in Q(λ). For convenience, we define

d_i =







d⁽¹⁾_i , if 2d = 4m, d⁽²⁾_i , if 2d = 4m + 2;

for i = 1, . . . , m.

From the row balancing of A1 in (2.7a) or (2.8a), we first set

η_i^(s)= max{1, max{∥A2m−k+s−1∥1 : k = 0, 1, . . . , 2i− 3 + s}} .

Then we take δ_i^(s)to be the geometric average of η_i^(s)and the average of the absolute

magnitudes of entries of A_2m−2i+1 in A1; i.e.,

δ_i^(s) = vu utη^(s)_i

( _n

∑

j=1

∑n ℓ=1

|A2m−2i+1(j, ℓ)|/n² )

for i = 1, . . . , m and s = 1, 2. Although the value of d^(s)_i can be set to δ^(s)_i to balance the entries of A1, we also need to consider the balance of the entries of both A0

and A1 in (2.7) or (2.8). As a result, we take the values of d^(s)₁ , . . . , d^(s)m to be the geometric average of the maximal values of δ₁^(s), . . . , δm^(s) 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

d^(s)_i :=

vu

utρ^(s) max

0≤k≤d

{ _n

∑

j=1

∑n ℓ=1

|Ak(j, ℓ)|/n² }

with

ρ^(s)= max{δi^(s); i = 1, . . . , m}.

2.6 Numerical Results

In [24], 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 performance and accuracy for solving H-PPEP of even degree by using structure-preserving algorithms and companion linearization.

For solving an n× n H-PPEP of even degree 2d, we apply the P-quadratization in Section 2 to transform it into a dn × dn H-PQEP. We then apply the SPA (Algorithm 2.1) in Section 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 [40, 41] to transform the H-PPEP into a palindromic linear pencil λZ^H + Z, and then utilize SPA to solve the H-PQEP:

(ˆλ²Z^H + ˆλ0 + Z)x = 0 with λ = ˆλ². The combination of the “good” linearization and SPA is called the PL_SPA algorithm. As mentioned in Remark 2.6 (ii), we see that applying the SPA to ˆλ²Z^T + ˆλ0 + Z is mathematically equivalent to applying the URV-based method [54] to λZ^T + Z.

2.6.1 Computational Cost

For making PQ_SPA more efficient, we reorder the submatrices of ( eK2, eN2) in step 1 of Algorithm 2.1 by the permutations

Π₁ =

where D₁ is a (2d− 2)n × (2d − 2)n diagonal matrix and V3, V₄ ∈ R^n×n. Set

Ke2 :=



 Π₁ 0 0 Π2



 eK2



 Π^T₂ 0 0 Π^T₁



 .

Now, we compare the computational costs for PQ_SPA and PL_SPA:

• (QR factorization and updating with real arithmetic operations) Compute Q₁ and Z₁ such that Q^T₁Ne2Z₁ = diag

(

N₁₁⁽¹⁾, (N₁₁⁽¹⁾)^T )

, where N₁₁⁽¹⁾ is upper triangular, and update Q^T₁Ke2Z₁. It requires (80/3n³ + 32dn³) and 341¹₃d³n³ flops for PQ_SPA and PL_SPA, respectively.

• (Given’s rotations and updating with real arithmetic operations) Reducing the new pair ( eK2, eN2) produced by above step to block upper triangular forms of (2.29), it requires 232d³n³− (296d²− 24d)n² and 1856d³n³− 800d²n² flops for PQ_SPA and PL_SPA, respectively.

• (Computing eigenvalues of (K11, N₁₁)) Computing eigenvalues of the real upper Hessenberg and triangular pair (K₁₁, N₁₁) by QZ algorithm, it requires 176d³n³ and 1408d³n³ flops for PQ_SPA and PL_SPA, respectively, to obtain the upper quasi-triangular and triangular pair.

• The eigenvectors of ( eK2, eN2) can be computed by an additional (408d³ + 32d)n³ − (332d² − 16d)n² and 1920d³n³ − 1088d²n² flops for PQ_SPA and PL_SPA, respectively.

We summarize the computational flops of PQ_SPA and PL_SPA in Table 2.1.

Eigenvalues Eigenvectors Total

PQ_SPA (408d³+ 32d + 80/3)n³ (408d³+ 32d)n³ (816d³+ 64d + 80/3)n³ PL_SPA 3605¹₃d³n³ 1920d³n³ 5525¹₃d³n³

Table 2.1: Computational flops of PQ_SPA and PL_SPA

2.6.2 Numerical Experiments

For an approximate eigenpair (λ, x) of the palindromic matrix polynomialP(λ), we define the associated relative residual by

RRes≡ RRes(λ, x) := ∥P(λ)x∥2

[∑d

k=1(|λ|^d+k+|λ|^d−k)∥Ak∥2+|λ|^d∥A0∥2

]∥x∥2

.

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 (2.1)).

As mentioned before, theoretically, the eigenvalues of H-PPEP appear in recipro-cal pairs (λ, 1/¯λ). So, if we sort the eigenvalues in ascending order by modulus, the product 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

r_i ≡ |λiλ¯_2dn+1−i− 1| (i = 1, . . . , dn).

All numerical experiments are carried out using MATLAB 2008b with machine precision eps ≈ 2.22 × 10⁻¹⁶.

LetCn,bdenote the set of n×n complex matrices which real and imaginary parts

LetCn,bdenote the set of n×n complex matrices which real and imaginary parts

在文檔中 保結構迴紋式二次化求解迴紋式方程 (頁 15-71)