The Palindromic Generalized Eigenvalue Problem
A
∗
x = λAx: Numerical Solution and Applications
Tiexiang Li
∗Chun-Yueh Chiang
†Eric King-wah Chu
‡Wen-Wei Lin
§June 30, 2009
Abstract
In this paper, we propose the palindromic doubling algorithm (PDA) for the palindromic generalized eigenvalue problem (PGEP) A∗x = λAx. We establish a complete convergence theory of the PDA for PGEPs without unimodular eigenvalues, or with unimodular eigenval-ues of partial multiplicities two (one or two for eigenvalue 1). Some important applications from the vibration analysis and the optimal control for singular descriptor linear systems will be presented to illustrate the feasibility and efficiency of the PDA.
Keywords. palindromic generalized eigenvalue problem, doubling algorithm, singular de-scriptor system
AMS subject classifications. 15A18, 15A22, 65F15
1
Introduction
In this paper, we develop the palindromic doubling algorithm (PDA) for the numerical solution of the palindromic generalized eigenvalue problem (PGEP)
A∗x = λAx, (1.1)
where A is a real or complex N × N matrix, λ ∈ C and x ∈ CN\{0} are eigenvalue and the
corresponding eigenvector of (1.1), respectively. Here, the symbol “∗00 = > (Transpose) or H
(Hermitian). The pencil A∗− λA and the pair (A∗, A) are usually called a palindromic linear pencil and a palindromic matrix pair, respectively. It is easily seen that the eigenvalues of (1.1) satisfy the reciprocal property, i.e., they appear in pairs as in {λ, 1/λ∗}.
The PGEPs with complex coefficient matrices were firstly suggested as “good” linearizations [5, 6] of palindromic polynomial/quadratic matrix pencils, arising from the study of vibration
∗Department of Mathematics, Southeast University, Nanjing, 211189, People’s Republic of China;
†Center for general education, National Formosa University, Huwei 632, Taiwan; [email protected] ‡School of Mathematical Sciences, Building 28, Monash University, VIC 3800, Australia; Phone:
+61-3-99054480, FAX: +61-3-99054403; [email protected]
§Department of Applied Mathematics, National Chiao Tung University, Hsinchu 300, Taiwan;
[email protected]. This research work is partially supported by the National Science Council and the National Center for Theoretical Sciences in Taiwan.
analysis [2, 4, 12]. A PGEP with real coefficient matrices can also be shown to be equiva-lent to the generalized continuous/discrete-time algebraic Riccati equations, associated with the continuous/discrete-time, linear-quadratic optimal control problems (see [11] for details).
The standard approach for solving the PGEP is to compute its generalized Schur form (e.g., by qz in MATLAB), ignoring its symmetric or palindromic structure in (A∗, A). However, the reciprocal property of eigenvalues of (1.1) is not preserved by computation generally, producing large numerical errors [7]. Recently, a QR-like algorithm [8] and a hybrid method [7] (which combines Jacobi-type method with the Laub’s trick) were proposed for the PGEP. The QR-like algorithm generally requires O(N4) flops and the hybrid method requires O(N3log(N )) flops. Alternatively, for methods cubic complexity, a URV-decomposition based structured method [9] and a structure-preserving algorithm [3] of are developed to solve PGEP were proposed, producing eigenvalues which are paired to working precision. Unfortunately for PGEPs, these more efficient (and equivalent) methods require the transformation of the PGEP to the quadratic form (µ2A∗+
µ · 0 + A)x = 0, leading to operations in larger 2N × 2N matrices. The PDA is a unique and more direct, thus more efficient, algorithm for the PGEP.
The purpose of this paper is to develop the PDA for solving the PGEP structurally. We establish quadratic convergence and linear convergence with rate 1/2 of the PDA, respectively, when (A∗, A) has no unimodular eigenvalues and has unimodular eigenvalues with partial mul-tiplicities two. In application to discrete-time optimal control problems, we especially develop a new algorithm combined with the PDA (as in Algorithm 4.1) for solving the optimal control of singular descriptor linear systems. To our knowledge, the associated generalized discrete-time algebraic Riccati equation (GDARE) has not been solved successfully in a structure-preserving manner.
This paper is organized as follows. In Section 2 we develop the palindromic doubling algorithm (PDA) for solving PGEPs. In Section 3 we establish the convergence theory for the PDA. In Section 4 we use the PDA to compute numerical solutions structurally in different applications in PGEPs, GCAREs and GDAREs. Concluding remarks are given in Section 5.
Throughout this paper, Cm×n
and Rm×n denote the sets of m × n complex and real matrices,
respectively. For convenience, we denote Cn
= Cn×1
, C = C1
, Rn
= Rn×1
and R = R1. The open
left-half plane and the open unit disk are denoted by C− and D1, respectively; 0m×n(0m) and Im
are the m × n(m × m) zero matrix and the m × m identity matrix, respectively. We use σ(A, B) to denote the spectrum of the matrix pair (A, B), and k · k denotes the 2-norm of a matrix.
2
Palindromic Doubling Algorithm
For a given palindromic matrix pair (A∗, A), we shall develop a doubling algorithm for solving the associated PGEP which preserves the palindromic structure at each iterative step.
Suppose −1 6∈ σ(A∗, A) (the assumption can be removed later in Remark 3.1). We then have
A∗(A∗+ A)−1A = ((A∗+ A) − A)(A∗+ A)−1(A + A∗− A∗) = (I − A(A∗+ A)−1)((A∗+ A) − A∗)
From (2.1), it is easily seen that
A(A∗+ A)−1, A∗(A∗+ A)−1 A∗ −A
= 0. (2.2)
We now define the doubling transform A → ˆA by ˆ
A = A(A∗+ A)−1A. (2.3) Theorem 2.1 The matrix pair ( ˆA∗, ˆA) has the doubling property; i.e., if
A∗U V = AU V S, (2.4) where U, V ∈ CN ×` and S ∈ C`×`, then ˆ A∗U V = ˆAU V S2. (2.5)
Proof. Multiplying the both sides of (2.4) by A∗(A∗+ A)−1, and (2.1) and (2.4) imply (2.5). From Theorem 2.1, we see that the doubling transform (2.3) preserves the palindromic struc-ture. So, for a palindromic matrix pair (A∗0, A0) with A0∈ CN ×N, we can develop the PDA to
generate the sequence {(A∗k, Ak)} if no breakdown occurs in the iterative process.
PDA Algorithm Given A0∈ CN ×N, τ (a small tolerance),
for k = 0, 1, 2, . . . , compute
Ak+1= Ak(A∗k+ Ak)−1Ak, (2.6)
if dist(Null(Ak+1), Null(Ak)) < τ , then stop,
end for
Here, “Null(·)” denotes the null space of the given matrix and “dist(·, ·)”denotes the distance between two subspaces.
To develop the PDA further, denote
Ak = Hk+ Kk, (2.7a) where Hk= 1 2(A ∗ k+ Ak) = Hk∗, Kk= 1 2(A ∗ k− Ak) = −Kk∗ (2.7b)
are the ∗-symmetric and ∗-anti-symmetric parts of Ak, respectively. Then the iteration (2.6) can
be rewritten as Ak+1 = Hk+1+ Kk+1= 1 2(Hk+ Kk)H −1 k (Hk+ Kk) = 1 2(Hk+ KkH −1 k Kk) + Kk.
The iteration (2.6) in the PDA can be simplified to
Hk+1= 1 2(Hk+ K0H −1 k K0), Kk+1= Kk= · · · = K0.
3
Convergence of PDA
Let A0∈ CN ×N. Suppose the eigenvalue “100of (A∗0, A0) (if exists) has partial multiplicity one or
two, and the other unimodular eigenvalues of (A∗
0, A0) (if exist) have exactly partial multiplicities
two. By the theorem of Kronecker canonical form there are nonsingular matrices Q and Z such that QA∗0Z =J0⊕ Ω0 0`,˜`⊕ Ir 0n,n˜ I`˜⊕ Ω0 ≡ C0, (3.1a) QA0Z = In 0n,˜n 0n,n˜ J˜0∗⊕ Ir ≡ D0, (3.1b)
where Ω0=diag(eiω1, · · · , eiωr), J0∈ C`×`consists of stable Jordan blocks (i.e., ρ(J0) < 1, where
ρ(·) is the radius of the spectrum) and ˜J0∗= J0⊕Imwith n = `+r, ˜` = `+m, ˜n = n+m = `+r+m
and N = 2n + m. Here “⊕00 denotes the direct sum of matrices.
Since C0D0= D0C0, from (3.1) we have that A∗0ZD0= A0ZC0. From Theorem 2.1 and steps
in the PDA, it follows that
A∗kZD20k = AkZC2
k
0 . (3.2)
Substituting (3.1) into (3.2), we get
A∗kZ In 0n,˜n 0n,n˜ ( ˜J0∗)2 k ⊕ Ir = AkZ " J02k⊕ Ω0 0`,˜`⊕ Γk 0˜n,n I`⊕ Ω2 k 0 # , (3.3) where Γk= 2kΩ2 k−1
0 . On the other hand, we can interchange the role of (A∗0, A0) by considering
the pair (A0, A∗0) which has the same Kronecker structure of (A∗0, A0). Therefore, there are
nonsingular P and Y such that
P A0Y = J∗ 0 ⊕ Ω∗0 0`,˜`⊕ Ir 0n,n˜ I`⊕ Ω∗0 ≡ E0, (3.4a) P A∗0Y = In 0n,˜n 0n,n˜ J˜0⊕ Ir ≡ F0, (3.4b)
Since E0F0 = F0E0, we deduce that A0Y F0 = A∗0Y E0. Using the similar arguments as in
(3.2)–(3.3), we obtain AkY J2k 0 ⊕ Ir 0n 0n In = A∗kY " (J0∗)2k⊕ (Ω∗ 0)2 k 0`,˜`⊕ Γ∗ k 0n,n˜ I`˜⊕ (Ω∗0)2 k # . (3.5)
We partition Ak, Hk and K0in (3.3) into four sub-blocks as in
Ak = Ak1 Ak3 Ak2 Ak4 , Hk = Hk1 Hk2∗ Hk2 Hk4 , K0= K01 −K02∗ K02 K04 , (3.6a) where Ak1, Hk1, K01 ∈ Cn×n, A∗k2, Ak3, Hk2∗ , K ∗ 02 ∈ Cnטn and Ak4, Hk4, K04 ∈ Cnט˜ n. From
(2.7) and (3.6a), we also have
Ak1= Hk1+ K01, Ak2= Hk2+ K02, (3.6b)
Furthermore, we partition Z in (3.3) and Y in (3.5) as in Z =Z1 Z3 Z2 Z4 , Y =Y1 Y3 Y2 Y4 , (3.7)
where Z1, Y1∈ Cn×n; Z2∗, Z3, Y2∗, Y3∈ Cnטn and Z4, Y4∈ Cnט˜ n. For convenience, we denote
Zi,a≡ Zi(:, 1 : `), Yi,a≡ Yi(:, 1 : `), i = 3, 4, (3.8)
Zi,b≡ Zi(:, ` + 1 : n), Yi,b≡ Yi(:, ` + 1 : n), i = 1, 2. (3.9)
Theorem 3.1 Let A0∈ CN ×N. Suppose that the eigenvalue “100of (A∗0, A0) (if exists) has partial
multiplicity one or two, the other unimodular eigenvalues of (A∗0, A0) (if exist) have exactly partial
multiplicities two, and (3.1) and (3.4) hold with ˜n ≤ 2` (i.e., n + m ≤ `). Suppose that Z1 and
Y1 in (3.7) are invertible, and W ≡ [ΦZ3,a− Z4,a|ΨY3,a− Y4,a] ∈ Cn×2`˜ is of full row rank, where
Φ ≡ Z2Z1−1, Ψ ≡ Y2Y1−1. If −1 6∈
Sr
j=1{e 2kiω
j, k ≥ 0}, then the sequence {(A∗
k, Ak)} generated
by the PDA is well defined and satisfies
A∗kZ1 Z2 → 0, linearly as k → ∞, (3.10a) Ak Y1 Y2 → 0, linearly as k → ∞, (3.10b)
with convergence rate 1/2, where spanZ1 Z2
and spanY3 Y4
form the weakly stable and the unstable invariant subspaces of (A∗0, A0) corresponding to (J0⊕ Ω0, In) and (In, J0∗⊕ Ω∗0),
respectively.
Proof. Since −1 6∈Sr
j=1{e
2kiωj, k ≥ 0}, from (2.6) we see that −1 6∈ σ(A∗
k, Ak), thus, A∗k+ Ak
is invertible for all k.
From (3.6a), (3.3) and (3.7), we have
A∗k1Z1+ A∗k2Z2=Ak1Z1(J2 k 0 ⊕ Ω 2k 0 ) + Ak3Z2(J2 k 0 ⊕ Ω 2k 0 ), (3.11) A∗k3Z1+ A∗k4Z2=Ak2Z1(J2 k 0 ⊕ Ω 2k 0 ) + Ak4Z2(J2 k 0 ⊕ Ω 2k 0 ), (3.12) A∗k1Z3((J0∗) 2k ⊕ Ir) + A∗k2Z4((J0∗) 2k ⊕ Ir) =Ak1[Z1(0`⊕ Γk) + Z3(I`⊕ Ω2 k 0 )] + Ak3[Z2(0`⊕ Γk) + Z4(I`⊕ Ω2 k 0 )], (3.13) A∗k3Z3(( ˜J0∗)2 k ⊕ Ir) + A∗k4Z4(( ˜J0∗)2 k ⊕ Ir) =Ak2[Z1(0`,˜`⊕ Γk) + Z3(I`˜⊕ Ω 2k 0 )] + Ak4[Z2(0`,˜`⊕ Γk) + Z4(I`˜⊕ Ω2 k 0 )]. (3.14) Post-multiplying (3.13) by 0˜`,`⊕ Γ−1k Ω2 k 0 , we get A∗k1Z3(0`,`˜ ⊕ Γ −1 k Ω 2k 0 ) + A∗k2Z4(0`,`˜ ⊕ Γ −1 k Ω 2k 0 ) = (Ak1Z1+ Ak3Z2)(0`⊕ Ω2 k 0 ) + (Ak1Z3+ Ak3Z4)(0`,`˜ ⊕ Ω2 k 0 Γ −1 k Ω 2k 0 ). (3.15)
Substituting (3.15) into (3.11) and using Ω2k 0 Γ −1 k Ω 2k 0 = 2−kΩ 2k+1 0 , we have A∗k1Z1+ A∗k2Z2=(Ak1Z1+ Ak3Z2)(J2 k 0 ⊕ 0r) + (Ak1Z1+ Ak3Z2)(0`⊕ Ω2 k 0 ) =(Ak1Z1+ Ak3Z2)(J2 k 0 ⊕ 0r) + (A∗k1Z3+ Ak2∗ Z4)(0`,`˜ ⊕ Γ−1k Ω2 k 0 ) − (A∗k1Z3+ A∗k3Z4)(0`,`˜ ⊕ 2−kΩ2 k+1 0 ). (3.16)
Using (3.6) and re-arranging (3.16), we get
Hk1 n Z1[In− (J2 k 0 ⊕ 0r)] − Z3[0`,`˜ ⊕ 2−kΩ0(Ir− Ω2 k 0 )] o + Hk2∗ nZ2[In− (J2 k 0 ⊕ 0r)] − Z4[0`,`˜ ⊕ 2−kΩ0(Ir− Ω2 k 0 )] o =K01 n Z1[In+ (J2 k 0 ⊕ 0r)] − Z3[0`,`˜ ⊕ 2−kΩ0(Ir+ Ω2 k 0 )] o − K02∗ nZ2[In+ (J2 k 0 ⊕ 0r)] − Z4[0`,`˜ ⊕ 2−kΩ0(Ir+ Ω2 k 0 )] o . (3.17) Denote k ≡ max{ρ(J0)2 k , 2−k} → 0, as k → ∞. (3.18) Since kΩ2k
0 k is bounded and Z1is invertible, by letting Φ ≡ Z2Z1−1, (3.17) can be simplified to
Hk1= −Hk2∗ (Φ + O(k)) + K01− K02∗ Φ + O(k). (3.19) Post-multiplying (3.13) by I`˜⊕ Γ−1k , we have A∗k1Z3(( ˜J0∗) 2k ⊕ Γ−1k ) + A∗k2Z4(( ˜J0∗) 2k ⊕ Γ−1k ) = Ak1[Z1(0˜`,`⊕ Ir) + Z3(I`˜⊕ Ω2 k 0 Γ−1k )] + Ak3[Z2(0`,`˜ ⊕ Ir) + Z4(I`⊕ Ω2 k 0 Γ−1k )]. (3.20)
From (3.6) and (3.18), (3.20) becomes
Hk1[Z3(I`⊕ 0m+r) + Z1(0`,˜`⊕ Ir) + O(k)] + Hk2∗ [Z4(I`⊕ 0m+r) + Z2(0`,˜`⊕ Ir) + O(k)]
= −K01[Z3(I`⊕ 0m+r) + Z1(0`,˜`⊕ Ir) + O(k)] + K02∗ [Z4(I`⊕ 0m+r) + Z2(0`,˜`⊕ Ir) + O(k)].
(3.21)
Substituting (3.19) into (3.21) we get
Hk2∗ n(Φ + O(k))[Z3(I`⊕ 0m+r) + Z1(0`,˜`⊕ Ir) + O(k)] − Z4(I`⊕ 0m+r)
−Z2(0`,˜`⊕ Ir) + O(k))
o
= O(1).
Since ΦZ1,b= Z2,b, it holds that
On the other hand, from (3.6a), (3.5) and (3.7), we have Ak1Y1+ Ak3Y2=A∗k1Y1((J0∗)2 k ⊕ (Ω∗0)2 k ) + A∗k2Y2((J0∗)2 k ⊕ (Ω∗0)2 k ), (3.23) Ak2Y1+ Ak4Y2=A∗k3Y1((J0∗) 2k⊕ (Ω∗ 0) 2k) + A∗ k4Y2((J0∗) 2k⊕ (Ω∗ 0) 2k), (3.24) Ak1Y3( ˜J2 k 0 ⊕ Ir) + Ak3Y4( ˜J2 k 0 ⊕ Ir) =(A∗k1Y3+ Ak2∗ Y4)(I`˜⊕ (Ω∗0) 2k ) + (A∗k1Y1+ A∗k2Y2) (0`,˜`⊕ Γ∗k), (3.25) Ak2Y3( ˜J2 k 0 ⊕ Ir) + Ak4Y4( ˜J2 k 0 ⊕ Ir) =(A∗k3Y3+ Ak4∗ Y4)(I`˜⊕ (Ω∗0) 2k) + (A∗ k3Y1+ A∗k4Y2) (0`,˜`⊕ Γ∗k). (3.26) As in (3.15)–(3.17), post-multiplying (3.25) by 0˜`,`⊕ (Γ∗k)−1(Ω∗0)2 k
and substituting it into (3.23), we have Ak1Y1+ Ak3Y2=(A∗k1Y1+ A∗k2Y2)((J0∗) 2k⊕ 0 r) + (A∗k1Y3+ A∗k3Y4)(0`,`˜ ⊕ 2−kΩ∗0) − (A∗k1Y3+ A∗k2Y4)(0`,`˜ ⊕ 2 −k(Ω∗ 0) 2k+1 ). (3.27)
From (3.6) and (3.18), (3.27) becomes
Hk1 n Y1[In− ((J0∗) 2k ⊕ 0r)] − Y3[0`,`˜ ⊕ 2−kΩ∗0(Ir− (Ω∗0) 2k )]o + Hk2∗ nY2[In− ((J0∗) 2k ⊕ 0r)] − Y4[0`,`˜ ⊕ 2−kΩ∗0(Ir− (Ω∗0) 2k )]o = − K01 n Y1[In+ ((J0∗) 2k ⊕ 0r)] − Y3[0`,`˜ ⊕ 2−kΩ∗0(Ir+ (Ω∗0) 2k )]o + K02∗ nY2[In+ ((J0∗) 2k⊕ 0 r)] − Y4[0`,`˜ ⊕ 2−kΩ∗0(Ir+ (Ω∗0) 2k)]o, (3.28) and then Hk1= −Hk2∗ (Ψ + O(k)) − K01+ K02∗ Ψ + O(k), where Ψ = Y2Y1−1.
Post-multiplying (3.25) by I`˜⊕ (Γ∗k)−1 and substituting (3.28) into it, we have
Hk2∗ nY4(0`⊕ 0m+r) + Y2(0`,˜`⊕ Ir) + O(k) − (Ψ + O(k))[Y3(I`⊕ 0m+r)
+Y1(0`,˜`⊕ Ir) + O(k)]
o
= O(1).
Since ΨY1,b= Y2,b, it holds that
Hk2∗ ([ΨY3,a− Y4,a] + O(k)) = O(1) ∈ C˜n×`. (3.29)
Combining (3.22) and (3.29) we get
Hk2∗ ([ΦZ3,a− Z4,a|ΨY3,a− Y4,a] + O(k)) = O(1) ∈ Cn×2`˜ . (3.30)
By the assumption that W ≡ [ΦZ3,a− Z4,a|ΨY3,a− Y4,a] ∈ C˜n×2` is of full row rank, it follows
that kHk2∗ k is uniformly bounded on k. Consequently, (3.19) implies that kHk1k, and in turn
kAk1k and kA∗k2k, are uniformly bounded on k. From (3.16), it follows that
Applying the similar argument as in (3.15)–(3.17) to (3.24) and (3.26), we deduce that
Hk4= −Hk2(Y2Y1−1+ O(k)) + K04− K02Y2Y1−1+ O(k).
Thus, (3.30) implies that kHk4k, and in turn kAk4k, are uniformly bounded on k.
To show A∗k3Z1+ A∗k4Z2= O(k), we post-multiply (3.14) by 0`,`˜ ⊕ Γ −1 k Ω 2k 0 and obtain A∗k3Z3(0`,`˜ ⊕ Γ −1 k Ω 2k 0 ) + A∗k4Z4(0˜`,`⊕ Γ −1 k Ω 2k 0 ) = (Ak2Z1+ Ak4Z2)(0`⊕ Ω2 k 0 ) + (Ak2Z3+ Ak4Z4)(0˜`,`⊕ 2−kΩ2 k+1 0 ) (3.32)
Substituting (3.32) into (3.12), as in (3.16) we have
A∗k3Z1+ A∗k4Z2=(Ak2Z1+ Ak4Z2)(J2 k 0 ⊕ 0r) + (A∗k3Z3+ Ak4∗ Z4)(0˜`,`⊕ 2−kΩ0) − (Ak2Z3+ Ak4Z4)(0`,`˜ ⊕ 2−kΩ2 k+1 0 ) =O(k) → 0, as → ∞. (3.33)
Combining (3.31) and (3.33), we have shown that
A∗ k1 A∗k2 A∗k3 A∗k4 Z1 Z2 = O(k) → 0, as k → ∞.
Similarly, as in (3.15)-(3.16), from (3.24) and (3.26) we have
Ak2Y1+ Ak4Y2=(A∗k3Y1+ A∗k4Y2)((J0∗) 2k ⊕ 0r) + (Ak2Y3+ Ak4Y4)(0`,`˜ ⊕ 2−kΩ∗0) − (A∗k3Y3+ A∗k4Y4)(0`,`˜ ⊕ 2−k(Ω∗0) 2k+1) = O( k). (3.34)
Using the boundedness of kAkik, i = 1, . . . , 4, and combining (3.27) and (3.34), we have shown
that Ak1 Ak3 Ak2 Ak4 Y1 Y2 = O(k) → 0, as k → ∞, because 1 2k dominates ρ(J0)2 k
in (3.18) for sufficiently large values of k. Remark 3.1 The assumption −1 6∈ Sr
j=1{e 2kiω
j, k ≥ 0} in Theorem 3.1 can be removed. In
fact, since U ≡Sr
j=1{e 2kiω
j, k ≥ 0} is a countable set (possibly dense on the unit circle only when
r → ∞), there exist an −eiθ06∈ U . Let A∗
new ≡ e
−iθ0
2 A∗0. Then A∗new+Anew = eiθ02 A0+e −iθ0
2 A0=
e−iθ02 (A∗
0+ eiθ0A0) is invertible. Consequently, the PDA is well-defined for (A∗new, Anew).
Theorem 3.2 Suppose (A∗0, A0) has no unimodular eigenvalues. The sequence {(A∗k, Ak)}
gen-erated by the PDA satisfies
A∗kZ1 Z2 → 0, Ak Y1 Y2 → 0, quadratically, as k → ∞, with convergence rate ρ(J0).
Proof. Since (A∗0, A0) has no unimodular eigenvalues, Theorem 3.1 implies (A∗k, Ak) has no
unimodular eigenvalues and (A∗k+ Ak) is invertible. So, the PDA is well-defined.
From (3.6a), (3.3) and (3.7), we have
A∗k1Z1+ A∗k2Z2= Ak1Z1J2 k 0 + Ak3Z2J2 k 0 , (3.35) A∗k3Z1+ A∗k4Z2= Ak2Z1J2 k 0 + Ak4Z2J2 k 0 , (3.36)
From (3.6a), it holds that
Hk1Z1+ Hk2∗ Z2= (K01Z1− K01∗ Z2)(I + J2 k 0 )(I − J 2k 0 )−1. Therefore, kHk1Z1+ Hk2∗ Z2k ≤ O(1). This implies k(Hk1+ K01)Z1+ (Hk2∗ − K02∗ )Z2k = kAk1Z1+ Ak3Z2k ≤ O(1). (3.37)
From (3.35) and (3.37), we have
A∗k1Z1+ A∗k2Z2= O(ρ(J0)2
k
) → 0, as k → ∞.
Similarly, from (3.36), we obtain
Hk2Z1+ Hk4∗ Z2= (K02Z1+ K04∗ Z2)(I + J2
k
0 )(I − J2
k
0 )−1,
which is uniformly bounded on k. This implies
A∗k3Z1+ A∗k4Z2= O(ρ(J0)2
k
) → 0, as k → ∞.
This shows that A∗kZ1 Z2
→ 0, quadratically, with convergence rate ρ(J0). Similarly, from (3.6a),
(3.5) and (3.7), we can also show that Ak
Y1
Y2
→ 0 quadratically, with rate ρ(J0).
4
Numerical Solution and Applications
In this section, we want to apply the PDA to find all the eigenpairs of a general GPEP, and solve the c-/d-stabilizing solutions of generalized continuous/discrete-time algebraic Riccati equations (GCARE/GDARE). We especially develop Algorithm 4.1 in subsection 4.3 for the computation of the d-semi-stabilizing solution of GDAREs arising in the optimal control of singular descriptor linear systems. To our knowledge, Algorithm 4.1 is the first structure-preserving algorithm for solving GDAREs associated with singular descriptor systems.
For operation counts or complexity, it depends on the details in the individual applications and whether efficiency can be squeezed from these fine structures. From the PDA, it is suffice to say that the algorithm is of O(N3) complexity per iteration. In addition, for problems without
unimodular eigenvalues, the convergence is quadratic and typically less than ten iterations are required for convergence to machine accuracy.
4.1
PGEP
In this subsection, we apply the PDA to solve the PGEP A∗0x = λA0x, where A0∈ C2n×2n. First,
we apply the PDA to A0 until convergence to Ak. Then we compute the bases Zs,Ys ∈ C2n×n
for the right and left null spaces of A∗k, respectively, satisfying A∗kZs= 0, Ys∗A∗k= 0.
This implies that there are S and T ∈ Cn×n with ρ(S) ≤ 1 and ρ(T ) ≤ 1 such that
A∗0Zs= A0ZsS, A0Ys= A∗0YsT. (4.1)
From (4.1), S and T can be computed by
S = (Ys∗A0Zs)−1(Ys∗A ∗
0Zs) ≡ S1−1S2,
T = (Zs∗A∗0Ys)−1(Zs∗A0Ys) ≡ S1−∗S2∗.
Rewrite the second equation of (4.1) as
A0(YsS1−∗) = A ∗ 0(YsS1−∗)S ∗ 2S −∗ 1 = A ∗ 0(YsS1−∗)S ∗.
Compute Sgj = λjgj and S∗hj = λ∗jhj, as well as zj = Zsgjand yj = (YsS1−∗)hj, for j = 1, . . . , n.
It holds that
A∗0zj = λjA0zj, λ∗jA ∗
0yj = A0yj, j = 1, . . . , n.
In the following example, we report the numerical results of the PDA to illustrate the linear convergence in the critical case. Recall that Theorem 3.1 shows the PDA converges linearly with rate 1/2 when all unimodular eigenvalues of (A∗
0, A0) have partial multiplicities two.
Example 4.1. Given α = cos(θ) and β = sin(θ) with θ = 0.62. Let
J0= 02 Γ I2 I2 , Js= 03 diag(λ1, λ2, λ3) I3 03 , where Γ =α −β β α
, and |λi| < 1 for i = 1, 2, 3. We construct
A0= Q∗(J0⊕ Js)Q,
where Q is an unitary matrix. It is easily seen that (A∗0, A0) has eigenvalues {α + ıβ, α −
ıβ, λ1, λ2, λ3, 1/λ∗1, 1/λ∗2, 1/λ∗3} with partial multiplicities {2, 2, 1, 1, 1, 1, 1, 1} which satisfy the
as-sumptions in Theorem 3.1. The kth absolute error as in (3.10a) computed by the PDA is defined by
Errk ≡ kA∗kZs,kk,
where Zs,kis an orthogonal basis corresponding to the five smallest singular values of Ak.
In Table 4.1, we list the absolute errors from the 20th to 24th iterations computed by the PDA which is observed to be linearly convergent with rate 1/2. Here, the tolerance τ in the PDA is chosen to be the optimal√1e − 16 = 1e − 8, because the unimodular eigenvalues of (A∗0, A0)
have partial multiplies two. Furthermore, the residual kA∗0Zs− A0ZsΛsk is given by 8.07e − 8,
ITs 20 21 22 23 24 Errk 1.26e-6 6.29e-7 3.15e-7 1.58e-7 8.01e-8
Table 4.1. Results for Example 1
4.2
GCARE
In this subsection, we are interested in finding the c-stabilizing solution of the generalized continuous-time algebraic Riccati equation (GCARE)
A>cXcEc+ Ec>XcAc− (Nc+ Ec>XcBc)R−1c (Nc+ E>cXcBc)>+ Mc= 0, (4.2)
which solves the continuous-time linear-quadratic control problem
min u 1 2 Z ∞ 0 x u > Mc Nc Nc> Rc x u dt (4.3a)
subject to the descriptor linear system
Ecx = A˙ cx + Bcu, x(0) = x0, (4.3b)
where Ec, Ac, Mc= Mc>, Xc = Xc> ∈ Rn×n, Bc, Nc∈ Rn×m and Rc= R>c ∈ Rm×m with Ec and
Rc being nonsingular. Furthermore, the c-stabilizing closed-loop matrix pencil of (4.3) is given
by Ac+ BcKc− λEc with the σ(Ac+ BcKc, Ec) ⊆ C− , where
Kc≡ −R−1c (B>cXcEc+ Nc>).
One common approach to solve (4.2) is to compute the n-dimensional, c-stable invariant subspace Sc of the symmetric/skew-symmetric pencil Mc− λLc corresponding to the eigenvalue
matrix pair (Sc, Ec) with σ(Sc, Ec) ⊆ C−, where
Mc− λLc≡ 0 Ac Bc A>c Mc Nc Bc> Nc> Rc − λ 0 Ec 0 −E> c 0 0 0 0 0 . (4.4)
We assume that the matrix pencil Mc − λLc has no eigenvalues on the imaginary axis. The
generalized eigenvalues of (Mc, Lc) can be arranged by
λ1, . . . , λ2n; ¯λ1, . . . , ¯λ2n; ∞, . . . , ∞
| {z }
m
,
where λi ∈ C−, for 1 ≤ i ≤ 2n. The m trivial infinity eigenvalues are from the nonsingularity of
Rc.
If Sc can be spanned by Uc, where
Uc= XcEc In Kc } n } n } m ,
then Xcis the c-stabilizing solution of GCARE (4.2) and Kcis the optimal controller for (4.3)[11].
In order to utilize the PDA to compute an orthogonal basis V =V1>, V2>, V3>> for Sc with
V1, V2∈ Rn×n, we consider the Cayley transformation
A> 0 − λA0= (Mc+ Lc) − λ(Mc− Lc), where A0= Mc− Lc= 0 Ac− Ec Bc A>c + Ec> Mc Nc B> c Nc> Rc . (4.5)
Then the c-stabilizing solution Xc for GCARE (4.2) can be obtained by Xc= V1V2−1E−1c .
To measure the accuracy of the computed solution ˜Xcfor the GCARE, we use the “normalized”
residual (NRc) NRc≡ kA> cX˜cEc+ Ec>X˜cAc− (Nc+ Ec>X˜cBc)R−1c (Nc+ Ec>X˜cBc)>+ Mck kA> cX˜cEck + kEc>X˜cAck + k(Nc+ Ec>X˜cBc)Rc−1(Nc+ Ec>X˜cBc)>k + kMck .
In the following example, we compare the residuals NRc of ˜Xc computed by care in MATLAB,
Newton’s method (NTM)[1] and the PDA.
Example 4.2 (Example 4.3 of [1]) Let n = m = 8. Given
Ac= diag 0 0 0 0 , 0 1 −1 0 , 0 2 −2 0 ,−1 1 0 −1 , Rc= I8, Ec= I8, Gc≡ BcR−1c B > c = trid(1, 2, 1) + e8e>1 + e1e>8, Mc= 08, Nc = 08,
where trid(1, 2, 1) is a 8 × 8 tridiagonal matrix with the sub-, main- and super-diagonal elements being 1, 2 and 1, respectively.
It is readily seen that Xc = 0 and σ(Ac− GcXc) = {−1, 0, ±i, ±2i} with purely imaginary
eigenvalues having linear elementary divisors. We apply the NTM method to GCARE (4.2) with X0= I8, and apply the PDA to
A0=
Gc Ac− Ec
A>c + E>c Mc
which is a degenrate form of (4.5) with Nc = 0. The tolerance τ in the NTM and the PDA is
chosen by 10−10. The numerical results are given in Table 4.2.
From Table 4.2, care in MATLAB dose not work because of the existence of the purely imaginary eigenvalues. We see that the NTM and the PDA almost have the same accuracy. Both methods have linear convergence rate 1/2, but the PDA requires much more iterative steps. However, the PDA only needs to compute a LU-factorization in each step, and NTM is accelerated by some modified technique [1] which needs to solve a more expensive Sylverster equation in each step.
care NTM PDA NRc ∗ 5.25 × 10−10 6.61 × 10−10
Iter. no. - 10 27
Table 4.2. Results for Example 3
4.3
GDARE
In this subsection, we are interested in finding the d-semi-stabilizing solution of the generalized discrete-time algebraic Riccati equation (GDARE)
A>dXdAd− Ed>XdEd− (Nd+ A>dXdBd)(Rd+ Bd>XdBd)−1(Nd+ A>dXdBd)>+ Md= 0 (4.6)
which solves the discrete-time linear-quadratic control problem
min uk 1 2 ∞ X k=0 xk uk > Md Nd Nd> Rd xk uk (4.7a)
subject to the singular descriptor linear system
Edxk+1= Adxk+ Bduk, x0= x0, (4.7b) where Ed, Ad, Md = Md>, Xd = Xd> ∈ R n×n, B d, Nd ∈ Rn×m and Rd = Rd> ∈ R m×m with E d
being singular. Furthermore, the d-semi-stabilizing closed-loop matrix pencil of (4.7) is given by Ad+ BdKd− λEd with the σ(Ad+ BdKd, Ed) ⊆ D1∪ {∞}, where
Kd≡ −(Rd+ Bd>XdBd)−1(Bd>XdAd+ Nd>).
One common approach to solve (4.6) is to compute the n-dimensional, d-semi-stable invariant subspace Sdof the matrix pencil Md− λLd corresponding to the eigenvalue matrix pair (Sd, Ed)
with σ(Sd, Ed) ⊆ D1∪ {∞}, where Md− λLd≡ 0 Ad Bd Ed> Md Nd 0 Nd> Rd − λ 0 Ed 0 A>d 0 0 Bd> 0 0 .
If Sd can be spanned by Ud, where
Ud= XdEd In Kd } n } n } m ,
then Xd is the d-semi-stabilizing solution of GDARE (4.6) and Kd is the optimal controller for
(4.7)[11].
We assume that the matrix pencil Md−λLdhas no eigenvalues on the unit circle, rd=nullity(Ed)
and ind∞(Ad, Ed) ≤ 1. From [11] we see that
So, the generalized eigenvalues of (Md, Ld) corresponding to (4.8) can be arranged by {λ1, . . . , λn−rd, ∞, . . . , ∞ | {z } rd } ∪ {0, . . . , 0 | {z } rd , λ−11 , . . . , λ−1n−r d} ∪ {∞, . . . , ∞ | {z } m }, (4.9)
where λi ∈ D1(can possibly be zero) , i = 1, . . . , n − rd. For convenience, we apply the convention
that 0 and ∞ are mutually reciprocal. The rdinfinity and rdzero eigenvalues in (4.9) are from the
assumption rd= nullity(Ed). The last trivial m infinity eigenvalues are from the last m columns
of Ld. In fact, (Ad+ BdKd, Ed) is an eigenvalue matrix pair associated with the d-semi-stable
invariant subspace Sd.
We now introduce an elegant transformation between the coefficient matrices of the GDARE (4.6) and GCARE (4.2) proposed by [11]. We define
fW : (Ed, Ad, Bd, Md, Nd, Rd) → (Ec, Ac, Bc, Mc, Nc, Rc),
where the matrices Ec, Ac, Bc, Mc, Nc, Rc satisfy
Ec 0 Ac Bc = χEd 0 Ad Bd Wd>= √ 2 2 Ad+ Ed Bd Ad− Ed Bd , (4.10a) Mc Nc Nc> Rc = Wd Md Nd Nd> Rd Wd>, (4.10b) in which χ = √ 2 2 In In −In In
, andAd+ Ed Bd = H 0 Wd is the QR-factorization with Wd
being orthogonal and H being lower triangular.
By the important property of fW in [11], it is assumed that rank Ad+ Ed Bd = n and
(Md, Ld) has no eigenvalue “−1”. Thus, the coefficient matrix tuple (Ec, Ac, Bc, Mc, Nc, Rc)
corresponds to a GCARE (4.2) with Ecand Rcbeing nonsingular. Furthermore, the GDARE (4.6)
and the GCARE (4.2) share the same stabilizing solutions, i.e., Xd= Xc.
We construct (Mc, Lc) by (Ec, Ac, Bc, Mc, Nc, Rc) as in (4.4) which satisfies
Mc+ Lc= W−1MdW, (4.11)
where W ≡ diag(√2In, Wd>). Let
(A>0, A0) = (Mc+ Lc, Mc− Lc) (4.12)
be the Cayley transformation of (Mc, Lc). From (4.10) and (4.11), we see that the eigenvalues
λd∈ σ(Md, Ld), µ ∈ σ(Mc, Lc) and λ ∈ σ(A>0, A0) satisfy the relationship in Table 4.3, in which
µ = (λ − 1)(λ + 1)−1. From Table 4.3, we see that the key property of the transformation λd→ λ
is to transform m trivial infinity eigenvalues to m trivial −1 while preserving other eigenvalues (including nontrivial ∞) unchanged.
In the following, we use the PDA and the special structure of (Md, Ld) for the computation
of the d-semi-stabilizing solution Xd of GDARE (4.6).
Firstly, we apply the PDA to the matrix A0 until convergence to Ak. Then we compute the
orthogonal bases Nrand N`∈ R(2n+m)×n for the right and left null spaces of A>k; i.e.,
λd 0 < |λd| < 1 |λd| > 1 0 ∞ m trivial ∞
µ Re(µ) < 0 Re(µ) > 0 −1 1 m trivial ∞ λ λ = λd λ = λd 0 ∞ m trivial −1
Table 4.3. Correspondence among λd, µ and λ
which form orthogonal bases for the d-stable invariant subspaces of (A>0, A0) and (A0, A>0),
respectively.
We then compute the QR-factorization A0Nr = Q1R1, where Q1 is orthogonal and R1 is
upper triangular. Next compute
S = Q>sA>0Nr, T = Q>sA0Nr, (4.14)
where Qs= Q1(:, 1 : n). We see that (S, T ) forms the d-stable eigenvalue matrix pair of (A>0, A0)
associated with Nr, and T is clearly nonsingular.
We would like to separate the invariant subspaces of (A>0, A0) corresponding to the zero and
nonzero d-stable eigenvalues. Let G = T−1S. By Van Dooren’s algorithm [10], there is an orthogonal matrix Φ ∈ Rn×nsuch that
ΦGΦ> =G11 G12 0 G22
,
where G11 ∈ Rs×s with σ(G11) = {λ ∈ σ(A>0, A0)| 0 < |λ| < 1} and G22 ∈ R(n−s)×(n−s) with
σ(G22) = {0}. Since σ(G11)T σ(G22) = φ, there is a Ψ = Is Ψ12 0 In−s such that Ψ−1Φ>GΦΨ =G11 0 0 G22 ,
where Ψ12 solves the Sylvester equation G11Ψ12− Ψ12G22= G12 uniquely. Then
V0= NrΦΨ(:, s + 1 : n), Vˆs= NrΦΨ(:, 1 : s) (4.15)
span the invariant subspaces of (A>0, A0) corresponding to the zero and nonzero d-stable
eigen-values, respectively.
Let ˆζ` spans the left null space of Ed. Then ζ` =h ˆζ`>, 0
i>
∈ R(2n+m)×rd contains the r
d
eigenvectors of (Md, Ld) corresponding to the trivial zeros. From the transformation (4.11), we
see that W−1ζ
`contains the rd eigenvectors of (A>0, A0) corresponding to trivial zeros. Now, we
want to extract W−1ζ
` from span{V0}.
Compute the QR-factorizationW−1ζ
` V0 = Q0R0, where Q0is orthogonal and R0is upper
triangular. Let
ˆ
V0= Q0(:, rd+ 1 : n − s) (4.16)
We will next find the invariant space U∞ of (A>0, A0) corresponding to the infinity
eigenval-ues. Compute the QR-factorization A0N`= Q∞R∞, where Q∞ is orthogonal and R∞ is upper
triangular. Let N∞= N`Q∞(:, s + 1 : n) ≡ N∞,1 N∞,2 N∞,3 } n } n } m , Vs=Vˆs Vˆ0 ≡ Vs1 Vs2 Vs3 } n } n } m .
From the Cayley transform, there is a full rank matrix Z ∈ R(n−s)×rd so that
V = Vs1 N∞,1Z Vs2 N∞,2Z Vs3 N∞,3Z ≡ V1 V2 V3 (4.17)
is a basis of an invariant subspace of (A>0, A0), satisfying Span{V } = Span
XcEc In Kc .
To determine Z, (4.17) and the fact Xc= Xc> imply
Vs2> Z>N> ∞,2 Ec>Vs1 N∞,1Z = Vs1> Z>N> ∞,1 EcVs2 N∞,2Z . That is, Vs2>Ec>Vs1= Vs1>EcVs2, (4.18a) Z>N∞,2> Ec>Vs1= Z>N∞,1> EcVs2, (4.18b) Vs2>Ec>N∞,1Z = Vs1>EcN∞,2Z, (4.18c) Z>N∞,2> Ec>N∞,1Z = Z>N∞,1> EcN∞,2Z. (4.18d)
Since Ec is nonsingular, from the isotropic property of
Vs1 Vs2 and N∞1 N∞2 , (4.18a) and (4.18d) hold automatically. Since (4.18b) is the transpose of (4.18c), the matrix Z is solved by finding the basis of Null(Vs2>E>c N∞,1− Vs1>EcN∞,2).
Finally, we have the d-semi-stabilizing solution Xd for GDARE (4.6) can be obtained by
Xd= Xc= V1V2−1E −1
c . (4.19)
We summarize the computational steps (4.12)–(4.19) for Xd in Algorithm 4.1.
Algorithm 4.1 (for GDARE (4.6) )
Input: Ed, Ad, Bd, Md, Nd, Rd; τ (a small tolerance);
Output: The d-semi-stabilizing solution Xd of (4.6).
1. Construct A0 via (4.12).
2. Apply PDA to (A>0, A0) until dist(Null(Ak), Null(Ak−1)) < τ .
3. Compute bases Nr, N` for the right and left null spaces of A>k as in (4.13).
4. Compute bases V0, ˆVs for d-stable invariant subspaces of (A>0, A0) as in (4.15).
5. Compute eigenvectors ˆV0 of (A>0, A0) corresponding to zeros as in (4.16) .
6. Determine Z by (4.18c).
In the following, we apply Algorithm 4.1 for a discrete-time descriptor system with Ed being
singular. To measure the accuracy of the computed solution ˜Xd for the GDARE, we use the
“normalized” residual: NRd≡ kA> dX˜dAd− Ed>X˜dEd− (Nd+ A>dX˜dBd)(Rd+ Bd>X˜dBd)−1(Nd+ A>dX˜dBd)>+ Mdk kA> dX˜dAdk + kEd>X˜dEdk + k(Nd+ A>dX˜dBd)(Rd+ Bd>X˜dBd)−1(Nd+ A>dX˜dBd>)>k + kMdk .
Example 4.3. With n = 10, m = 6, we let Ad = rand(n), Bd = rand(n, m), Nd = rand(n, m).
Construct
Ed= Ed1Ed2>, Rd= Rd1+ Rd1>, Md = Md1+ Md1>,
where Ed1= rand(n, n − 3), Ed2= rand(n, n − 3), Rd1= rand(m) and Md1= rand(n). We check
that nullity(Ed) = 3 and the algebraic multiplicity of the zero eigenvalue of (Md, Ld) is also 3.
Algorithm 4.1 gives
NRd= 1.472e − 015.
In this example, the PDA in Step 2 converges quadratically. In addition, Steps 4 and 5 are not required, and Z in Step 6 is chosen to be the obvious I3.
Example 4.4. With n = 20, m = 15, we let
˜
Ad=0n,4 rand(n, n − 4) , Bd= rand(n, m), ˜Nd>=0m,4 rand(m, n − 4) ,
˜
Ed = Ed1Ed2>, ˜Md= diag(04, Md1+ Md1>), Rd= Rd1+ R>d1,
where Ed1 = rand(n, n − 2), Ed2 = rand(n, n − 2), Rd1 = rand(m) and Md1 = rand(n − 4).
Construct
(Ad, Nd, Ed, Md) = ( ˜AdQ, Q>N˜d, Q>E˜d, Q>M˜dQ),
where Q is an arbitrarily orthogonal matrix. We check that nullity(Ed) = 2 and the algebraic
multiplicity of the zero eigenvalue of (Md, Ld) is 6. Then Algorithm 4.1 gives
NRd= 7.501e − 015.
In this example, the PDA in Step 2 converges quadratically. The eigenvectors ˆV0 ∈ R55×4
computed by (4.16) in Step 5 actually correspond to the 4 zero eigenvalues of (Ad+ BdKd, Ed),
and Z in Step 6 is a nontrivial 6 × 2-matrix of full rank as computed by (4.18c).
5
Concluding Remarks
In this paper, we have developed the palindromic doubling algorithm (PDA) for solving the palin-dromic generalized eigenvalue problem (PGEP) A∗x = λAx structurally. We prove quadratic convergence and linear convergence with rate 1/2 of the PDA, when (A∗, A) has no unimodu-lar eigenvalues and has unimoduunimodu-lar eigenvalues with partial multiplicities two (one or two for eigenvalue 1), respectively. Algorithm 4.1 is specially developed for the computation of the d-semi-stabilizing solution of the generalized discrete-time algebraic Riccati equation (GDARE) for the singular descriptor linear system. It is the first structure-preserving algorithm for singular descriptor systems.
Our numerical experience indicates that the PDA is not necessarily better than other spe-cialist algorithms (if exist) for solving the original problem, without linearizing the associated
palindromic matrix polynomials. Such specialist algorithms may be able to better utilize the finer structures of the original problems. Our numerical examples showed selected applications for which the PDA was better or when no specialist structure preserving algorithms exist. For future work, research will be conducted on how the finer structures can be fully utilized for indi-vidual applications. For a general PGEP without finer structures, the PDA is the only structure preserving algorithm which performs reasonably efficiently. Consequently, the “good” vibrations from “good” linearizations [5, 6] can always be computed using the PDA, in the absence of better methods. Of course, numerical solutions from the PDA or other methods may be refined using the finer structures in the original problems, if feasible.
References
[1] C.-H. Guo, and P. Lancaster. Analysis ans Modification of Newton’s Method For Algebraic Riccati Equations. Math. Comp., 67:1089–1105, 2006.
[2] A Hilliges, C. Mehl, and V. Mehrmann. On the solution of palindramic eigenvalue problems. In Proceedings 4th European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS), Jyv¨askyl¨a, Finland, 2004.
[3] T.-M. Huang, W.-W. Lin, and J. Qian. Structured algorithms for palindromic quadratic eigenvalue problems arising from vibration of fast trains. SIAM J. Matrix Anal. Appl., 30:1566–1592, 2009.
[4] C. F. Ipsen. Accurate eigenvalues for fast trains. SIAM News, 37, 2004.
[5] D. S. Mackey, N. Mackey, C. Mehl, and V. Mehrmann. Vector spaces of linearizations for matrix polynomials. SIAM J. Matrix Anal. Appl., 28:971–1004, 2006.
[6] D. S. Mackey, N. Mackey, C. Mehl, and V. Mehrmann. Structured polynomial eigenvalue problems: Good vibrations from good linearizations. SIAM J. Matrix Anal. Appl., 28:1029– 1051, 2006.
[7] D. S. Mackey, N. Mackey, C. Mehl, and V. Mehrmann. Numerical methods for palindromic eigenvalue problems: Computing the anti-triangular Schur form.. Numer. Linear Algebra Appl., 16:63–186, 2009.
[8] C. Schr¨oder. A QR-like algorithm for the palindromic eigenvalue problem. Technical Report 388, TU Berlin, MATHEON, Germany, 2007.
[9] C. Schr¨oder. URV decomposition based structured methods for palindromic and even eigen-value problems. Technical Report 375, TU Berlin, MATHEON, Germany, 2007.
[10] P. Van Dooren. The computation of Kronecker’s canonical form of a singular pencil. Linear Alg. Appl., 27:103–141, 1979.
[11] H. Xu. Transformations between discrete-time and continuous-time algebraic Riccati equa-tions Linear Alg. Appl., 425:77–101, 2007.
[12] S.Zaglmayr. Eigenvalue problems in saw-filter simulations. Diplomarbeit, Institute of Com-putation Mathematics, Johannes Kepler Univerisity Linz, Linz, Austria, 2002.
