Contents lists available atScienceDirect
Linear Algebra and its Applications
j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / l a aThe palindromic generalized eigenvalue problem
A
∗
x
= λ
Ax: Numerical solution and applications
Tiexiang Li
a, Chun-Yueh Chiang
b, Eric King-wah Chu
c,∗, Wen-Wei Lin
d aDepartment of Mathematics, Southeast University, Nanjing 211189, People’s Republic of ChinabCenter for General Education, National Formosa University, Huwei 632, Taiwan cSchool of Mathematical Sciences, Building 28, Monash University, VIC 3800, Australia
dCMMSC and NCTS, Department of Applied Mathematics, National Chiao Tung University, Hsinchu 300, Taiwan
A R T I C L E I N F O A B S T R A C T
Article history: Received 29 June 2009 Accepted 19 November 2009 Available online 12 January 2010 Submitted by V. Mehrmann AMS classification: 15A18 15A22 65F15 Keywords:
Palindromic generalized eigenvalue problem
Doubling algorithm Singular descriptor system
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 eigenvalues 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.
© 2009 Elsevier Inc. All rights reserved.
1. Introduction
In this paper, we develop the palindromic doubling algorithm (PDA) for the numerical solution of the palindromic generalized eigenvalue problem (PGEP)
∗Corresponding author. Tel.: +61 3 99054480; fax: +61 3 99054403.
E-mail addresses: [email protected] (T. Li), [email protected] (C.-Y. Chiang), [email protected] (E.K.-w. Chu), [email protected] (W.-W. Lin).
0024-3795/$ - see front matter © 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.laa.2009.12.020
A∗x
= λ
Ax, (1.1) where A is a real or complex N×
N matrix,λ ∈
C
and x∈
C
N\{
0}
are eigenvalue and the correspond-ing eigenvector of (1.1), respectively. Here, the symbol “∗
”=
(Transpose) or H (Hermitian). The pencilA∗
− λ
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 analysis [2,
4,12]. A PGEP with real coefficient matrices can also be shown to be equivalent 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., byqz 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 of cubic complexity, a URV-decomposition based structured method [9] and a structure-preserving algorithm [3] for PGEPs 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 multiplicities two. In appli-cation 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,
C
m×n andR
m×n denote the sets of m×
n complex and real matrices,respectively. For convenience, we denote
C
n=
C
n×1,C = C
1,R
n=
R
n×1andR = R
1. The open left-half plane and the open unit disk are denoted byC
−andD
1, respectively; 0m×n(
0m)
and Imare 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·
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/∈ σ(
A∗, A)
(the assumption can be removed later in Remark 3.1). We then haveA∗
(
A∗+
A)
−1A= ((
A∗+
A) −
A)(
A∗+
A)
−1(
A+
A∗−
A∗)
= (
I−
A(
A∗+
A)
−1)((
A∗+
A) −
A∗)
=
A(
A∗+
A)
−1A∗.
(2.1)From (2.1), it is easily seen that
A
(
A∗+
A)
−1, A∗(
A∗+
A)
−1 A−
∗ A
We now define the doubling transform A
→
A byA
=
A(
A∗+
A)
−1A.
(2.3)Theorem 2.1. The matrix pair
(
A∗, A)
has the doubling property;
i.
e.
, ifA∗U
=
AUS, (2.4)where U
∈
C
N×and S∈
C
×, thenA∗U
=
AUS2.
(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 Theorem2.1, we see that the doubling transform (2.3) preserves the palindromic structure. So, for a palindromic matrix pairA∗0, A0
with A0
∈
C
N×N, we can develop the PDA to generate thesequence A∗k, Ak
if no breakdown occurs in the iterative process.
PDA Algorithm Given A0
∈
C
N×N,τ
(a small tolerance),for k
=
0, 1, 2,. . .
, compute Ak+1=
Ak A∗k+
Ak −1 Ak, (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=
H∗k, Kk=
1 2 Ak−
A∗k= −
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+
KkHk−1Kk+
Kk.
The iteration (2.6) in the PDA can be simplified toHk+1
=
1 2 Hk+
K0Hk−1K0 , Kk+1=
Kk= · · · =
K0.
3. Convergence of PDALet A0
∈
C
N×N. Suppose the eigenvalue “1" of
A∗0, A0
(if exists) has partial multiplicity one or two, and the other unimodular eigenvalues ofA∗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˜ J0∗⊕
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 “⊕
" denotes the direct sum of matrices.Since C0D0
=
D0C0, from (3.1b) we have that A∗0ZD0=
A0ZC0. From Theorem2.1and steps in thePDA, it follows that
A∗kZD20k
=
AkZC2k
0
.
(3.2)Substituting (3.1b) into (3.2), we get
A∗kZ In 0n,n˜ 0n,n˜ J0∗ 2k
⊕
Ir=
AkZ ⎡ ⎣J2 k 0⊕
Ω 2k 0 0,˜⊕
k 0˜n,n I⊕
Ω2 k 0 ⎤ ⎦ , (3.3) whereΓk=
2kΩ2 k−10 . On the other hand, we can interchange the role of
A∗0, A0 by considering the pairA0, A∗0
which has the same Kronecker structure asA∗0, A0
. Therefore, there are nonsingular P and Y such that
PA0Y
=
J∗0⊕
Ω0∗ 0,˜⊕
Ir 0n,n˜ I˜⊕
Ω0∗≡
E0, (3.4a) PA∗0Y=
In 0n,˜n 0n,n˜ J0⊕
Ir≡
F0, (3.4b)Since E0F0
=
F0E0, we deduce that A0YF0=
A∗0YE0. Using the similar arguments as in (3.2) and (3.3),we obtain AkY In 0n 0n J2 k 0
⊕
Ir=
A∗kY ⎡ ⎣ J0∗2k⊕
Ω0∗2k 0,˜⊕
Γk∗ 0n,n˜ I˜⊕
Ω∗ 0 2k ⎤ ⎦.
(3.5)We partition Ak, Hkand K0in (2.7a) into four sub-blocks as in
Ak
=
Ak1 Ak3 Ak2 Ak4 , Hk=
Hk1 H∗k2 Hk2 Hk4 , K0=
K01−
K02∗ K02 K04 , (3.6a) where Ak1, Hk1, K01∈
C
n×n, A∗k2, Ak3, Hk2∗, K02∗∈
C
nטnand Ak4, Hk4, K04
∈
C
˜nטn. From (2.7a) and (3.6a),we also have
Ak1
=
Hk1+
K01, Ak2=
Hk2+
K02, (3.6b)Ak3
=
H∗k2−
K02∗, Ak4=
Hk4+
K04.
(3.6c)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
∈
C
n×n; Z∗2, Z3, Y2∗, Y3∈
C
nטnand Z4, Y4∈
C
n˜×˜n. For convenience, we denoteZi,a
≡
Zi(:
, 1: )
, Yi,a≡
Yi(:
, 1: )
, i=
3, 4, (3.8)Theorem 3.1. Let A0
∈
C
N×N.
Suppose that the eigenvalue “1" of
A∗0, A0
(
if exists)
has partial multiplicity one or two, the other unimodular eigenvalues ofA∗0, A0
(
if exist)
have exactly partial multiplicities two, and(
3.1b)
and(
3.4b)
hold withn˜
2(
i.
e.
, r+
m).
Suppose that Z1and Y1in(
3.7)
are invertible, and W≡ [
Z3,a−
Z4,a|
Y3,a−
Y4,a] ∈
C
n˜×2is of full row rank, whereΦ≡
Z2Z1−1,Ψ≡
Y2Y1−1.
If−
1/∈
rj=1e2kiωj, k 0, then the sequence A∗k, Ak
generated by the PDA is well defined and satisfies A∗k Z1 Z2
→
0, linearly as k→ ∞
, (3.10a) Ak Y1 Y2→
0, linearly as k→ ∞
, (3.10b)with convergence rate 1
/
2, where span Z 1 Z2 and span Y 3 Y4form the weakly stable and the unstable invariant subspaces ofA∗0, A0 corresponding to
(
J0⊕
Ω0, In)
and In, J0∗⊕
Ω0∗ , respectively.
Proof. Since−
1/∈
rj=1 e2kiωj, k 0, from (2.6) we see that
−
1/∈ σ
A∗k, Ak, thus, A∗k
+
Akis in-vertible for all k.From (3.6a), (3.3) and (3.7), we have
A∗k1Z1
+
A∗k2Z2=
Ak1Z1 J02k⊕
Ω02k+
Ak3Z2 J02k⊕
Ω02k , (3.11) A∗k3Z1+
A∗k4Z2=
Ak2Z1 J02k⊕
Ω02k+
Ak4Z2 J02k⊕
Ω02k , (3.12) A∗k1Z3 J∗ 0 2k⊕
Ir+
A∗k2Z4 J∗ 0 2k⊕
Ir=
Ak1 Z1(
0,˜⊕
Γk) +
Z3 I˜⊕
Ω02k+
Ak3 Z2(
0,˜⊕
Γk) +
Z4 I˜⊕
Ω02k , (3.13) A∗k3Z3J0∗ 2k⊕
Ir+
A∗k4Z4J0∗ 2k⊕
Ir=
Ak2 Z1(
0,˜⊕
Γk) +
Z3 I˜⊕
Ω02k+
Ak4 Z2(
0,˜⊕
Γk) +
Z4 I˜⊕
Ω02k.
(3.14) Post-multiplying (3.13) by 0˜,⊕
Γk−1Ω02k, we get A∗k1Z3 0˜,⊕
Γk−1Ω02k+
A∗k2Z4 0˜,⊕
Γk−1Ω02k= (
Ak1Z1+
Ak3Z2)
0⊕
Ω02k+ (
Ak1Z3+
Ak3Z4)
0˜,⊕
Ω02kΓk−1Ω02k.
(3.15)Substituting (3.15) into (3.11) and usingΩ02kΓk−1Ω02k
=
2−kΩ02k+1, we haveA∗k1Z1
+
A∗k2Z2= (
Ak1Z1+
Ak3Z2)
J02k⊕
0r+ (
Ak1Z1+
Ak3Z2)
0⊕
Ω02k= (
Ak1Z1+
Ak3Z2)
J02k⊕
0r+
A∗k1Z3+
A∗k2Z4 0˜,⊕
Γk−1Ω02k− (
Ak1Z3+
Ak3Z4)
0˜,⊕
2−kΩ02k+1.
(3.16)Using (3.6a) and re-arranging (3.16), we get
Hk1Z Z1 In
−
J02k⊕
0r−
Z3 0˜,⊕
2−kΩ0 Ir−
Ω2 k 0+
Hk2∗ Z2 In−
J02k⊕
0r−
Z4 0˜,⊕
2−kΩ0 Ir−
Ω2 k 0=
K01 Z1 In+
J02k⊕
0r−
Z3 0˜,⊕
2−kΩ0 Ir+
Ω2 k 0−
K02∗ Z2 In+
J02k⊕
0r−
Z4 0˜,⊕
2−kΩ0 Ir+
Ω2 k 0.
(3.17) Denotek
≡
maxρ(
J0)
2 k , 2−k→
0, as k→ ∞.
(3.18) SinceΩ2k0 is bounded and Z1is invertible, by lettingΦ
≡
Z2Z1−1, (3.17) can be simplified toHk1
= −
H∗k2(
Φ+
O(
k)) +
K01−
K02∗Φ+
O(
k).
(3.19) Post-multiplying (3.13) by I˜⊕
Γk−1, we have A∗k1Z3J0∗ 2k⊕
Γ−1 k+
A∗k2Z4J0∗ 2k⊕
Γ−1 k=
Ak1 Z1(
0,˜⊕
Ir) +
Z3 I˜⊕
Ω02kΓk−1+
Ak3 Z2(
0,˜⊕
Ir) +
Z4 I˜⊕
Ω02kΓk−1.
(3.20)From (3.6a) 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⊕
2Im⊕
0r) +
Z1(
0,˜⊕
Ir) +
O(
k)]
+
K02∗[
Z4(
I⊕
2Im⊕
0r) +
Z2(
0,˜⊕
Ir) +
O(
k)].
(3.21) Substituting (3.19) into (3.21) we getH∗k2
(
Φ+
O(
k))[
Z3(
I⊕
0m+r) +
Z1(
0,˜⊕
Ir) +
O(
k)] −
Z4(
I⊕
0m+r)
−
Z2(
0,˜⊕
Ir) +
O(
k)
=
O(
1).
SinceΦZ1,b=
Z2,b, it holds thatH∗k2
([
ΦZ3,a−
Z4,a] +
O(
k)) =
O(
1) ∈
C
n˜×.
(3.22)On the other hand, from (3.6a), (3.5) and (3.7), we have
Ak1Y1
+
Ak3Y2=
A∗k1Y1 J∗02k⊕
Ω0∗2k+
A∗k2Y2 J0∗2k⊕ (
Ω0∗)
2k , (3.23) Ak2Y1+
Ak4Y2=
A∗k3Y1 J∗02k⊕
Ω0∗2k+
A∗k4Y2 J∗02k⊕
Ω0∗2k , (3.24) Ak1Y3 J2k 0⊕
Ir+
Ak3Y4 J2k 0⊕
Ir=
A∗k1Y3+
A∗k2Y4 I˜⊕ (
Ω0∗)
2k+
A∗k1Y1+
A∗k2Y2 0,˜⊕
Γk∗ , (3.25) Ak2Y3 J2k 0⊕
Ir+
Ak4Y4 J2k 0⊕
Ir=
A∗k3Y3+
A∗k4Y4 I˜⊕ (
Ω0∗)
2k+
A∗k3Y1+
A∗k4Y2 0,˜⊕
Γk∗.
(3.26)As in (3.15) and (3.17), post-multiplying (3.25) by 0˜,
⊕
Γk∗−1Ω0∗2kand substituting it into (3.23), we have Ak1Y1+
Ak3Y2=
A∗k1Y1+
A∗k2Y2 J0∗2k⊕
0r+
A∗k1Y3+
A∗k3Y4 0˜,⊕
2−kΩ0∗−
A∗k1Y3+
A∗k2Y4 0˜,⊕
2−k(
Ω0∗)
2k+1.
(3.27)From (3.6a) and (3.18), (3.27) becomes
Hk1 Y1 In
−
J0∗2k⊕
0r−
Y3 0˜,⊕
2−kΩ0∗ Ir−
Ω0∗ 2k+
Hk2∗ Y2 In−
J∗02k⊕
0r−
Y4 0˜,⊕
2−kΩ0∗ Ir−
Ω0∗ 2k= −
K01 Y1 In+
J0∗2k⊕
0r−
Y3 0˜,⊕
2−kΩ0∗(
Ir+ (
Ω0∗)
2k)
+
K02∗ Y2 In+
J0∗ 2k⊕
0r−
Y4 0˜,⊕
2−kΩ0∗ Ir+ (
Ω0∗)
2k , (3.28) and then Hk1= −
H∗k2(
Ψ+
O(
k)) −
K01+
K02∗Ψ+
O(
k)
, whereΨ=
Y2Y1−1.Post-multiplying (3.25) by I˜
⊕
Γk∗−1and substituting (3.28) into it, we haveHk2∗ Y4
(
I⊕
0m+r) +
Y2(
0,˜⊕
Ir) +
O(
k) − (
Ψ+
O(
k))[
Y3(
I⊕
0m+r)
+
Y1(
0,˜⊕
Ir) +
O(
k)]
=
O(
1).
SinceΨY1,b=
Y2,b, it holds thatHk2∗
([
Ψ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) ∈
C
n˜×2.
(3.30) By the assumption that W≡ [
ΦZ3,a−
Z4,a|
ΨY3,a−
Y4,a] ∈
C
n˜×2is of full row rank, it follows thatH∗k2is uniformly bounded on k. Consequently, (3.19) implies that
Hk1, and in turnAk1andA∗k2, are uniformly bounded on k. From (3.16), it follows thatA∗k1Z1
+
A∗k2Z2=
O(
k) →
0, as k→ ∞.
(3.31)Applying the similar argument as in (3.15) and (3.17) to (3.24) and (3.26), we deduce that
Hk4
= −
Hk2Y2Y1−1
+
O(
k)
+
K04−
K02Y2Y1−1+
O(
k).
Thus, (3.30) implies that
Hk4, and in turnAk4, are uniformly bounded on k. To show A∗k3Z1+
A∗k4Z2=
O(
k)
, we post-multiply (3.14) by 0˜,⊕
Γk−1Ω2 k 0 and obtain A∗k3Z3 0˜,⊕
Γk−1Ω02k+
A∗k4Z4 0˜,⊕
Γk−1Ω02k= (
Ak2Z1+
Ak4Z2)
0⊕
Ω02k+ (
Ak2Z3+
Ak4Z4)
0˜,⊕
2−kΩ02k+1.
(3.32)Substituting (3.32) into (3.12), as in (3.16) we have A∗k3Z1
+
A∗k4Z2= (
Ak2Z1+
Ak4Z2)
J02k⊕
0r+
A∗k3Z3+
A∗k4Z4 0˜,⊕
2−kΩ0− (
Ak2Z3+
Ak4Z4)
0˜,⊕
2−kΩ02k+1=
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) and (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∗2 k+1=
O(
k).
(3.34)Using the boundedness of
Aki, 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)
2kin (3.18) for sufficiently large values of k.
Remark 3.1. Consider the assumption
−
1/∈
U≡
rj=1e2kiωj, k 0in Theorem 3.1. Since U is a countable set (possibly dense on the unit circle only when r→ ∞
), there exist an−
eiθ0/∈
U. WithA∗new
≡
e−iθ0/2A∗0, we have A∗new
+
Anew=
eiθ0/2A0+
e−iθ0/2A0=
e−iθ0/2
A∗0
+
eiθ0A 0
being invert-ible. It is unclear how the “optimal"
θ
0can be found.Theorem 3.2. SupposeA∗0, A0
has no unimodular eigenvalues
.
The sequence A∗k, Ak generated by the PDA satisfies A∗k Z1 Z2→
0, Ak Y1 Y2→
0,quadratically, as k
→ ∞
, with convergence rateρ(
J0).
Proof. SinceA∗0, A0
has no unimodular eigenvalues, Theorem3.1impliesA∗k, Ak
has no unimodular eigenvalues andA∗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+
J02k I−
J02k −1.
Therefore, Hk1Z1+
H∗k2Z2 O(
1).
This implies(
Hk1+
K01)
Z1+
H∗k2
−
K02∗Z2=
Ak1Z1+
Ak3Z2 O(
1).
(3.37)From (3.35) and (3.37), we have
A∗k1Z1
+
A∗k2Z2=
O
ρ(
J0)
2k
→
0, as k→ ∞.
Similarly, from (3.36), we obtain
Hk2Z1
+
Hk4∗Z2=
K02Z1+
K04∗Z2 I+
J02k I−
J02k −1 ,which is uniformly bounded on k. This implies
A∗k3Z1
+
A∗k4Z2=
O
ρ(
J0)
2k
→
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 Y
1
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 PGEP, 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 eigen-values, 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∈
C
2n×2n. First, weapply the PDA to A0until convergence to Ak. Then we compute the bases Zs,Ys
∈
C
2n×nfor the right and left null spaces of A∗k, respectively, satisfyingA∗kZs
=
0, Ys∗A∗k=
0.
This implies that there are S and T
∈
C
n×nwithρ(
S)
1 andρ(
T)
1 such thatA∗0Zs
=
A0ZsS, A0Ys=
A∗0YsT.
(4.1)From (4.1), S and T can be computed by
S
=
Ys∗A0Zs −1 Ys∗A∗0Zs≡
S−11S2, T=
Zs∗A∗0Ys −1 Zs∗A0Ys≡
S−∗1 S∗2.
Rewrite the second equation of (4.1) as
A0 YsS−∗1
=
A∗0 YsS−∗1 S2∗S1−∗=
A∗0 YsS1−∗ S∗.
Compute Sgj
= λ
jgjand S∗hj= λ
j∗hj, as well as zj=
Zsgjand yj=
YsS−∗1
hj, for j
=
1,. . .
, n. It holds thatTable 4.1
Results for Example 1.
ITs 20 21 22 23 24
Errk 1.26e-6 6.29e-7 3.15e-7 1.58e-7 8.01e-8
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 Theorem3.1shows the PDA converges linearly with rate 1
/
2 when all unimodular eigenvalues ofA∗0, A0
have partial multiplicities two.
Example 4.1. Given
α =
cos(θ)
andβ =
sin(θ)
withθ =
0.
62. LetJ0
=
02 Γ I2 I2 , Js=
03 diag(λ
1,λ
2,λ
3)
I3 03 ,whereΓ
=
α −ββ α , and|λ
i| <
1 for i=
1, 2, 3. We constructA0
=
Q∗(
J0⊕
Js)
Q ,where Q is an unitary matrix. It is easily seen thatA∗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 assumptions in Theorem3.1. The kth absolute error as in (3.10a) computed by the PDA is defined byErrk
≡
A∗kZs,k,where Zs,kis an orthogonal basis corresponding to the five smallest singular values of Ak.
In Table4.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 ofA∗0, A0
have partial multiplies two. Furthermore, the residualA∗0Zs
−
A0Zssis given by 8.
07e−
8, where Zs≡
Zs,24andsis the corresponding approximate stable eigenvalue matrix.4.2. GCARE
In this subsection, we are interested in finding the c-stabilizing solution of the generalized continuous-time algebraic Riccati equation (GCARE)
AcXcEc
+
Ec XcAc−
Nc+
Ec XcBc R−c1 Nc+
EcXcBc+
Mc=
0, (4.2)which solves the continuous-time linear-quadratic control problem min u 1 2 ∞ 0 x u Mc Nc Nc Rc x u dt (4.3a)
subject to the descriptor linear system
Ecx
˙
=
Acx+
Bcu, x(
0) =
x0, (4.3b)where Ec, Ac, Mc
=
Mc, Xc=
Xc∈
R
n×n, Bc, Nc
∈
R
n×mand Rc=
Rc∈
R
m×mwith Ecand Rcbeing nonsingular. Furthermore, the c-stabilizing closed-loop matrix pencil of (4.3b) is given by Ac
+
BcKc−
λ
Ecwith theσ(
Ac+
BcKc, Ec) ⊆
C
−, whereKc
≡ −
R−c1 BcXcEc+
Nc.
Let Mc− λ
Lc≡
⎡ ⎢ ⎣ 0 Ac Bc Ac Mc Nc Bc Nc Rc ⎤ ⎥ ⎦− λ
⎡ ⎣−
0Ec E0c 00 0 0 0 ⎤ ⎦.
(4.4)One common approach to solve (4.2) is to compute the n-dimensional, c-stable invariant subspaceUcof the symmetric/skew-symmetric pencilMc
− λ
Lccorresponding to the eigenvalue matrix pair(
Sc, Ec)
withσ(
Sc, Ec) ⊆
C
−, whereUc is the column space of Uc∈
R
(2n+m)×n which satisfiesMcUcEc=
LcUcSc.We assume that the matrix pencilMc
− λ
Lchas no eigenvalues on the imaginary axis. The gener-alized eigenvalues of(M
c,Lc)
can be arranged byλ
1,. . .
,λ
2n; ¯λ
1,. . .
,¯λ
2n; ∞
,. . .
!,∞
" m,
where
λ
i∈
C
−, for 1 i 2n. The m trivial infinity eigenvalues are from the nonsingularity of Rc. With Uc=
⎡ ⎣XcInEc Kc ⎤ ⎦}
}
nn}
m ,Xcis the c-stabilizing solution of GCARE (4.2) and Kcis the optimal controller for (4.3b) [11]. In order to utilize the PDA to compute an orthogonal basis V
=
V1, V2, V3
forUcwith V1, V2
∈
R
n×n, we consider the Cayley transformation A0− λ
A0= (
Mc+
Lc) − λ(M
c−
Lc)
, where A0=
Mc−
Lc=
⎡ ⎢ ⎣ 0 Ac−
Ec Bc Ac+
Ec Mc Nc Bc Nc Rc ⎤ ⎥ ⎦.
(4.5)Then the c-stabilizing solution Xcfor GCARE (4.2) can be obtained by Xc
=
V1V2−1E−c1.To measure the accuracy of the computed solution Xc for the GCARE, we use the “normalized" residual
(
NRc)
NRc≡
A cXcEc+
EcXcAc−
Nc+
EcXcBc R−c1 Nc+
EcXcBc+
Mc AcXcEc+
EcXcAc+ (
Nc+
EcXcBc)
R−c1 Nc+
EcXcBc+
Mc.
In the following example, we compare the residuals NRcof Xccomputed bycare in MATLAB, Newton’s method (NTM) [1] and the PDA.
Example 4.2 (Example 4.3 of [1]). Let n
=
m=
8. GivenAc
=
diag 0 0 0 0 , 0 1−
1 0 , 0 2−
2 0 ,−
1 1 0−
1 , Rc=
I8, Ec=
I8, Gc≡
BcR−c1Bc=
trid(
1, 2, 1) +
e8e1+
e1e8, 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 eigen-values having linear elementary divisors. We apply the NTM method to GCARE (4.2) with X0=
I8, andapply the PDA to A0
=
Gc Ac−
Ec Ac+
Ec Mc ,which is a degenrate form of (4.5) with Nc
=
0. The toleranceτ
in the NTM and the PDA is chosen to be 10−10. The numerical results are given in Table4.2.Table 4.2
Results for Example 3.
care NTM PDA
NRc ∗ 5.25×10−10 6.61×10−10
Iter. no. - 10 27
From Table4.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 Sylvester equation in each step.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)
AdXdAd
−
EdXdEd−
Nd+
AdXdBd Rd+
BdXdBd −1 Nd+
AdXdBd+
Md=
0, (4.6) which solves the discrete-time linear-quadratic control problemmin uk 1 2 ∞ # k=0 xk uk M d 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∈
R
n×mand Rd=
Rd∈
R
m×mwith E dbeing singu-lar. Furthermore, the d-semi-stabilizing closed-loop matrix pencil of (4.7b) is given by Ad+
BdKd− λ
Ed with theσ(
Ad+
BdKd, Ed) ⊆
D
1∪ {∞}
, whereKd
≡ −
Rd+
BdXdBd −1 BdXdAd+
Nd.
Let Md− λ
Ld≡
⎡ ⎢ ⎢ ⎣ 0 Ad Bd Ed Md Nd 0 Nd Rd ⎤ ⎥ ⎥ ⎦− λ
⎡ ⎢ ⎢ ⎣ 0 Ed 0 Ad 0 0 Bd 0 0 ⎤ ⎥ ⎥ ⎦.
One common approach to solve (4.6) is to compute the n-dimensional, d-semi-stable invariant sub-spaceUd of the matrix pencilMd
− λ
Ldcorresponding to the eigenvalue matrix pair(
Sd, Ed)
withσ (
Sd, Ed) ⊆
D
1∪ {∞}
, whereUdis the column space of Ud∈
R
(2n+m)×nwhich satisfiesMdUdEd=
LdUdSd. With Ud=
⎡ ⎣XdInEd Kd ⎤ ⎦}
}
nn}
m ,Xdis the d-semi-stabilizing solution of GDARE (4.6) and Kdis the optimal controller for (4.7b) [11]. Assume that the matrix pencilMd
− λ
Ldhas no eigenvalues on the unit circle, rd=
nullity(
Ed)
and ind∞(
Ad, Ed)
1. From [11] we see thatσ (M
d,Ld) = σ(
Ad+
BdKd, Ed) ∪ σ
Ed,
(
Ad+
BdKd)
∪ σ (
0m, Im).
(4.8)Table 4.3
Correspondence amongλd,μandλ.
λd 0< |λd| <1 |λd| >1 0 ∞ m trivial∞ μ Re(μ) <0 Re(μ) >0 −1 1 m trivial∞ λ λ = λd λ = λd 0 ∞ m trivial−1
{λ
1,. . .
,λ
n−rd,∞
,. . .
!,∞
" rd} ∪
⎧ ⎪ ⎨ ⎪ ⎩0, ! ". . .
r , 0 d ,λ
−11,. . .
,λ
−n−1rd ⎫ ⎪ ⎬ ⎪ ⎭∪ {∞
,. . .
!m ,∞
"}
, (4.9) whereλ
i∈
D
1(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 subspaceUd.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, RcsatisfyEc 0 Ac Bc
= χ
Ed 0 Ad Bd Wd=
√
1 2 Ad+
Ed Bd Ad−
Ed Bd Wd, (4.10a) Mc Nc Nc Rc=
Wd Md Nd Nd Rd Wd, (4.10b) in whichχ =
√1 2 I n In −In In, and [Ad+Ed Bd]
=
[H 0]Wdis the QR-factorization with Wdbeing orthogonal and H being lower triangular.By the important property of fWin [11], it is assumed that rank
(
Ad+Ed Bd) =
n and(M
d,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
(M
c,Lc)
by(
Ec, Ac, Bc, Mc, Nc, Rc)
as in (4.4) which satisfiesMc
+
Lc=
W−1MdW , (4.11)where W
≡
diag(
√
2In, Wd)
. Let(A
0,A0) = (M
c+
Lc,Mc−
Lc)
(4.12)be the Cayley transformation of
(M
c,Lc)
. From (4.10b) and (4.11), we see that the eigenvaluesλ
d∈
σ (M
d,Ld)
,μ ∈ σ(M
c,Lc)
andλ ∈ σ(A
0,A0)
satisfy the relationship in Table4.3, in whichμ =
(λ −
1)(λ +
1)
−1. From Table4.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
(M
d,Ld)
for the computation of the d-semi-stabilizing solution Xdof GDARE (4.6).Firstly, we apply the PDA to the matrixA0until convergence toAk. Then we compute the orthogonal basesNrandN
∈
R
(2n+m)×nfor the right and left null spaces ofAk; i.e.,AkNr
=
0, AkN=
0, (4.13)which form orthogonal bases for the d-stable invariant subspaces of A0,A0 and A0,A0 , respec-tively.
We then compute the QR-factorizationA0Nr
=
Q1R1, where Q1 is orthogonal and R1 is uppertriangular. Next compute
S
=
QsA0Nr, T=
QsA0Nr, (4.14)where Qs
=
Q1(:
, 1:
n)
. We see that(
S, T)
forms the d-stable eigenvalue matrix pair ofA0,A0
associated withNr, and T is clearly nonsingular.
We would like to separate the invariant subspaces of A0,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Φ∈
R
n×nsuch thatΦGΦ
=
G11 G12 0 G22 , where G11∈
R
s×s withσ(
G11) =
λ ∈ σ
A0,A0|
0< |λ| <
1 and G22∈
R
(n−s)×(n−s) withσ (
G22) = {
0}
. Sinceσ(
G11)
+σ(
G22) = φ
, there is aΨ=
I s Ψ12 0 In−s such that Ψ−1 ΦGΦΨ=
G11 0 0 G22 ,whereΨ12solves the Sylvester equation G11Ψ12
−
Ψ12G22=
G12uniquely. ThenV0
=
NrΦΨ(:
, s+
1:
n)
, Vs=
NrΦΨ(:
, 1:
s)
(4.15)span the invariant subspaces of
(A
0,A0)
corresponding to the zero and nonzero d-stable eigenvalues,respectively.
Let
ˆζ
spans the left null space of Ed. Thenζ
=
ˆζ
, 0∈
R
(2n+m)×rdcontains the r deigenvectors of(M
d,Ld)
corresponding to the trivial zeros. From the transformation (4.11), we see that W−1ζ
contains the rd eigenvectors ofA0,A0
corresponding to trivial zeros. Now, we want to extract
W−1
ζ
from span{
V0}
.Compute the QR-factorization,W −1ζ V0-
=
Q0R0, where Q0is orthogonal and R0is uppertri-angular. Let
V0
=
Q0(:
, rd+
1:
n−
s)
, (4.16)which forms the eigenvectors of
A0,A0
corresponding to zero eigenvalues of
(
Sd, Ed)
. We will next find the invariant space U∞of A0,A0
corresponding to the infinity eigenvalues. Compute the QR-factorizationA0N
=
Q∞R∞, where Q∞is orthogonal and R∞is upper triangular.Let N∞
=
NQ∞(:
, s+
1:
n) ≡
⎡ ⎣N∞N∞,1,2 N∞,3 ⎤ ⎦}
}
nn}
m , Vs=
,Vs V0-≡
⎡ ⎣VVs1s2 Vs3 ⎤ ⎦}
}
nn}
m.
From the Cayley transform, there is a full rank matrix Z∈
R
(n−s)×rdso thatV
=
⎡ ⎣VVs1s2 N∞N∞,1,2ZZ Vs3 N∞,3Z ⎤ ⎦≡
⎡ ⎣VV12 V3 ⎤ ⎦ (4.17)is a basis of an invariant subspace of A0,A0
, satisfying Span
{
V} =
Span . XcEc In Kc / . To determine Z, (4.17) and the fact Xc=
XcimplyVs2 ZN∞,2 Ec ,Vs1 N∞,1Z
-=
Vs1 ZN∞,1 Ec , Vs2 N∞,2Z-.
That is,
Vs2EcVs1
=
Vs1EcVs2, (4.18a)ZN∞,2EcVs1
=
ZN∞,1EcVs2, (4.18b)Vs2EcN∞,1Z
=
Vs1EcN∞,2Z, (4.18c)ZN∞,2EcN∞,1Z
=
ZN∞,1EcN∞,2Z.
(4.18d)Since Ecis nonsingular, from the isotropic property of V s1 Vs2 andN∞1 N∞2
, (4.18a) and (4.18d) hold au-tomatically. Since (4.18b) is the transpose of (4.18c), the matrix Z is solved by finding the basis of Null Vs2EcN∞,1
−
Vs1EcN∞,2.
Finally, we have the d-semi-stabilizing solution Xdfor GDARE (4.6) can be obtained by
Xd
=
Xc=
V1V2−1E− 1c
.
(4.19)We summarize the computational steps (4.12)–(4.19) for Xdin 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 Xdof (4.6).1. ConstructA0via (4.12). 2. Apply PDA to A0,A0
until dist
(
Null(A
k)
, Null(A
k−1)) < τ
.3. Compute basesNr,Nfor the right and left null spaces ofAk as in (4.13).
4. Compute bases V0, Vsfor d-stable invariant subspaces of A0,A0 as in (4.15). 5. Compute eigenvectors V0of A0,A0 corresponding to zeros as in (4.16). 6. Determine Z by (4.18c). 7. Compute Xd
=
V1V2−1Ec−1as in (4.19).In the following, we apply Algorithm 4.1 for a discrete-time descriptor system with Edbeing singular. To measure the accuracy of the computed solution Xdfor the GDARE, we use the “normalized" residual:
NRd≡ A dXdAd−EdXdEd− Nd+AdXdBd Rd+BdXdBd −1 Nd+AdXdBd +Md AdXdAd+EdXdEd+ Nd+AdXdBd Rd+BdXdBd −1 Nd+AdXdBd + Md .
Example 4.3. With n
=
10, m=
6, we let Ad=
rand(
n)
, Bd=
rand(
n, m)
, Nd=
rand(
n, m)
. ConstructEd
=
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(M
d,Ld)
is also 3. Algorithm 4.1 givesNRd
=
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
