Structured doubling algorithms for weakly stabilizing Hermitian solutions of algebraic Riccati equations

(1)

Structured doubling algorithms for weakly stabilizing Hermitian solutions of algebraic

Riccati equations

Tsung-Ming Huang

^a,1

, Wen-Wei Lin

^b,

∗

aDepartment of Mathematics, National Taiwan Normal University, Taipei 11677, Taiwan.

bDepartment of Mathematics, National Tsing Hua University, Hsinchu 30043, Taiwan.

Abstract

In this paper, we propose structured doubling algorithms for the computation of the weakly stabilizing Hermitian solutions of the continuous- and discrete-time algebraic Riccati equations, respectively. Assume that the partial multiplicities of purely imaginary and unimodular eigenvalues (if any) of the associated Hamiltonian and symplectic pencil, respectively, are all even and the C/DARE and the dual C/DARE have weakly stabilizing Hermitian solutions with property (P). Under these assump- tions, we prove that if these structured doubling algorithms do not break down, then they converge to the desired Hermitian solutions globally and linearly. Numerical experiments show that the structured doubling algorithms perform efficiently and reliably.

Key words: algebraic Riccati equation, Hermitian solution, structured doubling algorithm, purely imaginary eigenvalue, unimodular eigenvalue, global and linear convergence

∗ Corresponding author.

Email addresses: min@math.ntnu.edu.tw (Tsung-Ming Huang), wwlin@math.nthu.edu.tw (Wen-Wei Lin).

1 This work is partially supported by the National Science Council and the National Center for Theoretical Sciences of Taiwan.

(2)

1 Introduction

In this paper, we investigate the structured doubling algorithms for the com- putation of the weakly stabilizing Hermitian solution X to

(I) the continuous-time algebraic Riccati equation (CARE):

−XGX + A^HX + XA + H = 0, (1.1)

or

(II) the discrete-time algebraic Riccati equation (DARE):

X = A^HX(I + GX)⁻¹A + H, (1.2)

where A, G, H ∈ C^n×n with G = G^H, H = H^H and I ≡ In is the identity matrix of order n.

Equations (1.1) and (1.2) arise frequently in the pursuit of the “weakly” stabi- lizing controllers of continuous- and discrete-time H_∞-optimal control systems, respectively (13; 15; 17; 25). In addition, several applications in Wiener fil- tering theory (47), network synthesis (3) and Moser-Veselov equations (9; 40) also involve the Hermitian solution of CAREs.

We consider the 2n × 2n Hamiltonian matrix H associated with the CARE:

H =



 A −G

−H −A^H



 (1.3)

which satisfies

HJ = −JH^H, J =



 0 I

−I 0



; (1.4)

and consider the 2n×2n symplectic pair (M, L) (or symplectic pencil M−λL) associated with the DARE:

M =



 A 0

−H I



, L =



I G 0 A^H



 (1.5)

which satisfies

(3)

MJM^H = LJL^H. (1.6) The special symplectic pair (M, L) of the form in (1.5) is referred to as a standard symplectic form (SSF).

Note that (30; 42) λ ∈ σ(H) if and only if −¯λ ∈ σ(H), and λ ∈ σ(M, L) if and only if 1/¯λ ∈ σ(M, L). Here σ(H) and σ(M, L) denote the spectrums of H and (M, L), respectively. It is well-known that (e.g., (30; 32; 42)) the CARE (1.1) has a weakly stabilizing Hermitian solution X if and only if



 A −G

−H −A^H







 I X



=



 I X



Φ, Φ ∈ C^n×n, (1.7)

where σ(Φ) ⊆ R₋ (the closed left half plane); the DARE (1.2) has a weakly stabilizing Hermitian solution X if and only if



 A 0

−H I







 I X



=



I G 0 A^H







 I X



Φ, Φ ∈ C^n×n, (1.8)

where σ(Φ) ⊆ °₁ (the closed unit disk) and (I + GX) is invertible. The particular invariant subspaces spanned by ^hI, X^Tⁱ^T in (1.7) and (1.8) are usually referred to as stable Lagrangian subspaces.

Definition 1.1 A subspace U ⊆ C²ⁿ with dimension n is called an H-stable Lagrangian subspace, if U satisfies that (i) HU ⊆ U, (ii) U is isotropic; i.e., x^HJy = 0, for all x, y ∈ U; and (iii) Re(λ(H|U)) ≤ 0. Here λ(H|U) denotes an eigenvalue of H restricted to U.

Definition 1.2 A subspace V ⊆ C²ⁿ with dimension n is called an (M, L)- stable Lagrangian subspace, if (i) V is invariant under (M, L); i.e., there is a subspace W such that MV, LV ⊆ W (44, pp. 303–305); (ii) V is isotropic and (iii) |λ((M, L)|_V)| ≤ 1. Here λ((M, L)|_V) denotes an eigenvalue of (M, L) restricted to V.

Unfortunately, an H-stable Lagrangian subspace and an (M, L)-stable La- grangian subspace do not always exist, when some purely imaginary eigen- values of H and some unimodular eigenvalues of (M, L) have odd partial multiplicities, respectively. Counterexamples can be found in (41). To guar- antee the existence of the H- and (M, L)-stable Lagrangian subspaces and the weakly stabilizing Hermitian solutions of CAREs and DAREs, we assume that H in (1.3) and (M, L) in (1.5), respectively, satisfy the conditions:

(4)

(A1) The partial multiplicities (the sizes of Jordan blocks) of H associated with the purely imaginary eigenvalues (if any) are all even.

(A2) The partial multiplicities of (M, L) associated with the unimodular eigen- values (if any) are all even.

Under these assumptions, equivalence statements for the existence of weakly stabilizing Hermitian solutions of CAREs and DAREs have first been given by (18; 29) and (28), respectively. In order to enhance the uniqueness of the weakly stabilizing Hermitian solution in some sense, we, further, give the following definitions.

Definition 1.3 Assume that (A1) holds. The CARE (1.1) is said to have a weakly stabilizing Hermitian solution with property (P), if the matrix Φ in (1.7) (i.e., Φ ≡ A − GX) satisfies that σ(Φ) ⊆ R− and each purely imaginary eigenvalue has a half of the partial multiplicity of H corresponding to the same eigenvalue. (For example, H has a purely imaginary eigenvalue with partial multiplicity (2, 4, 6) i.e., the jordan blocks are of size 2, 4 and 6, respectively, then Φ has the same eigenvalue with the partial multiplicity (1, 2, 3).)

Definition 1.4 Assume that (A2) holds. The DARE (1.2) is said to have a weakly stabilizing Hermitian solution with property (P), if the matrix Φ in (1.8) (i.e., Φ ≡ (I + GX)⁻¹A) satisfies that σ(Φ) ⊆ °₁ and each unimodular eigenvalue has a half of the partial multiplicity of (M, L) corresponding to the same eigenvalue.

For the continuous-time case, a well-known backward stable algorithm care (30) computes a stabilizing Hermitian solution X for the CARE by applying the QR algorithm with reordering (43) to H. Unfortunately, the QR algo- rithm preserves neither the Hamiltonian structure nor the associated split- ting of eigenvalues. When H in (1.3) has no purely imaginary eigenvalues, a strongly stable method has been proposed by (10) for computing the Hamil- tonian Schur form of H, and therefore, the H-stable Lagrangian subspace.

Efficient structured doubling algorithms (incorporating an appropriate Cay- ley transform) (11; 27) and the matrix sign function methods (5; 8; 14; 16) have been developed to compute the unique positive semidefinite solution of CARE (1.1). When H in (1.3) satisfies Assumption (A1), an eigenvector de- flation technique proposed by (13) guarantees that the eigenvalues appear with the correct pairing. This is certainly an advantage over the QR algo- rithm, but the method ignores most of the structure of the problem during computation. A structured algorithm proposed by (1) only using symplectic orthogonal transformations, computes the H-stable Lagrangian subspace. But there are numerical difficulties in the convergence of the deflation steps when purely imaginary eigenvalues occur (38, p. 143). To avoid the numerical difficulties mentioned above, another stable and structured algorithm has been developed in (33), preprocessing to deflate all purely imaginary eigenvalues.

(5)

When H satisfies Assumption (A1) with partial multiplicities equal to two, an efficient Newton’s method has been developed in (21) for solving the CAREs with global and linear convergence.

For the discrete-time case, a well-known backward stable algorithm dare (37; 42; 48) computes a stabilizing Hermitian solution X for DAREs by apply- ing the QZ algorithm with reordering to (M, L). Unfortunately, this algorithm does not take into account the symplectic structure of (M, L). Non-structure- preserving iterative processes spoil the symplectic structure, causing the al- gorithms to fail or lose accuracy in adverse circumstances. When (M, L) in (1.5) has no unimodular eigenvalues, an efficient doubling algorithm was firstly derived in (2) based on an acceleration scheme of the fixed point iteration for (1.2). Using different approaches, quadratic convergence of doubling algorithms has been shown in (26; 35). On the other hand, based on the viewpoint of the inverse-free iteration (4; 36), a matrix disk function method (MDFM) (6; 7) and a structure-preserving doubling algorithm (SDA) (12; 24) have been developed for solving DAREs. The symplectic structure in the MDFM and the SSF form in the SDA are preserved at each iterative step. However, the symplectic structure in the MDFM is preserved only in exact arithmetic.

When (M, L) in (1.5) satisfies Assumption (A2), a structured algorithm has been developed in (33), preprocessing to deflate all unimodular eigenvalues by the determining the isotropic Jordan subbasis using the S + S⁻¹-transform of M − λL (31). When (M, L) satisfies Assumption (A2) with partial multiplic- ities two, an efficient Newton-type method has been proposed by (20) to solve the DAREs with global and linear convergence.

As mentioned above, the MDFM and SDA have been proposed for solving DARE (1.2) with M − λL possessing no unimodular eigenvalues. To solve CARE (1.1) with H with no purely imaginary eigenvalues, the Hamiltonian matrix H is converted to a symplectic pencilM−λ^c L in SSF by an appropriate^b Cayley transform and then the MDFM or the SDA algorithm can be applied.

The main purpose of this paper is to apply the MDFM or SDA to solve CAREs and DAREs, where the associated H in (1.3) and M − λL in (1.5) satisfy Assumption (A1) and (A2), respectively. Under these assumptions, we prove the globally linear convergence of the MDFM and SDA.

This paper is organized as follows. In Sections 2 and 3, we describe structured doubling algorithms, the D-MDFM/D-SDA and C-MDFM/C-SDA, for solving DAREs and CAREs, respectively. In Section 4, we prove that under Assumption (A1) and (A2) structured doubling algorithms converge globally and linearly to the weakly stabilizing Hermitian solutions with property (P) of DAREs and CAREs, respectively. In Section 5, we test several numerical examples for illustrating the convergence behavior of the MDFM, SDA and Newton-type methods. Concluding remarks are given in Section 6.

(6)

Throughout this paper, we denote A^H = ¯A^T the conjugate transpose of A ∈ C^n×n, ı = √

−1, I ≡ In and 0 ≡ 0n (the identity and zero matrices of order n, respectively). The vector e_j the jth column of I_n, k · k a matrix norm, and σ(A) and ρ(A) the spectrum and the spectral radius of A, respectively. R₋ and °1, respectively, denote the closed left half plane and the closed unit disk.

2 Structured doubling algorithms for DAREs

The matrix disk function method (MDFM) in (6; 7) is developed to solve DARE (1.2) by using a swapping technique built on the QR factorization. We refer to this step as a QR-swap.

For a given symplectic pair (M, L), the QR-swap computes the QR-factorization of [L^T, −M^T]^T from

Q



 L

−M



≡



Q₁₁ Q₁₂ M∗ L∗







 L

−M



=



R 0



, (2.1)

where Q ∈ C^4n×4n is unitary and R ∈ C^2n×2n is upper triangular. Let

L := Lb _∗L, M := M^c _∗M. (2.2)

It is easily seen that (M,^c L) is symplectic. From (2.1)–(2.2), (^b M,^c L) satisfies^b the doubling property:

Mx = Mc _∗Mx = λM_∗Lx = λL_∗Mx = λ²L_∗Lx = λ²Lx,^b (2.3)

assuming that Mx = λLx.

Algorithm 2.1 (D-MDFM for DAREs)

(7)

Input: A, G, H; τ (a small tolerance);

Output: a weakly stabilizing Hermitian solution X to DARE.

Initialize: R ← 02n, M ←



 A 0

−H I



, L ←



I G 0 A^H



;

Repeat: Compute the QR-factorization:



Q₁₁ Q₁₂ M_∗ L_∗







 L

−M



=





Rc

0



;

If kR − Rk ≤ τ k^c Rk, Then solve the least squares problem for X:^c

−M(:, 1 : n) = M(:, n + 1 : 2n)X;

Stop Else

Set L ← L_∗L, M ← M_∗M, R ←R;^c Go To Repeat

End If

The sequence {(M_k, L_k)} generated by Algorithm 2.1 satisfies the recursive formula

Mk+1= M∗,kMk, Lk+1 = L∗,kLk, (2.4) where



Q_11,k Q_12,k M_∗,k L_∗,k







 L_k

−M_k



=





Rc_k

0





is the QR-factorization.

On the other hand, the SDA in (12) is developed to solve DARE (1.2) under conditions of stabilizability and detectability, by using a structured LU fac- torization instead of the QR factorization in (2.1). We refer to this step as a SLU-swap. As derived in (12), for (M, L) in SSF (1.5), we construct

M_∗ =



 A(I + GH)⁻¹ 0

−A^H(I + HG)⁻¹H I



, L_∗ =



I AG(I + HG)⁻¹ 0 A^H(I + HG)⁻¹



 (2.5)

(8)

and consequently deduce that

M∗L = L∗M. (2.6)

With L ≡ L^b _∗L andL ≡ M^b _∗M, and apply the Sherman-Morrison-Woodbury formula, we obtain

L =b



I G^b 0 A^b^H



, M =^c





A 0b

−H I^c



, (2.7)

where

A = A(I + GH)b ⁻¹A, (2.8a)

G = G + AG(I + HG)b ⁻¹A^H, (2.8b)

H = H + Ac ^H(I + HG)⁻¹HA. (2.8c)

Equations in (2.7) show that the newly derived matrix pair (M,^c L) is again^b in SSF form. From (2.6)–(2.7), (M,^c L) also satisfies the doubling property^b Mx = λc ²Lx.^b

Remark 2.1

(i) Equations in (2.8) have exactly the same form as the doubling algo- rithm (which has been first proposed and investigated in Anderson (2) and Kimura (26)). However, the original doubling algorithm was derived as an acceleration scheme for the fixed-point iteration from (1.2). Instead of producing the sequence {Xk}, the doubling algorithm produces {X₂^k}.

Furthermore, the convergence of the doubling algorithm was proven when A is nonsingular (2), and for (A, G, H) which is reachable and detectable, or stabilizable and observable (26). A stronger convergent result of the SDA algorithm under weaker conditions (stabilizability and detectability) can be found in (12; 35).

(ii) The matrix (I + GH) in (2.8) can possibly be singular in some step of SDA, thus A,^b G and^b H in (2.8) do not exist and the SDA may break^c down. In our numerical experiments in Section 6, this happens only in the limiting case.

Algorithm 2.2 (D-SDA for DAREs)

(9)

Input: A, G, H; τ (a small tolerance);

Output: a weakly stabilizing Hermitian solution X to DARE.

Repeat W ← I + GH;

If W is singular, then break down.

Solve W V₁ = A, V₂W = A for V₁ and V₂; Set G ← G + V₂GA^H,

H ← H + Ac ^HHV₁, A ← AV₁

If kH − Hk ≤ τ k^c Hk, then X ←^c H, Stop.^c Set H ←H;^c

Goto Repeat

Remark 2.2 The linear systems W V₁ = A, V₂W = A for V₁ and V₂ in Algorithm2.2 can be solved by the LU factorization of W or the GSVD (gen- eralized singular value decomposition) of (W, A) and (W^H, A^H), respectively.

For the latter case, let







T₁W U_a= C₁, T₁AV_a= S₁,







T₂W^HU_b = C₂, T₂A^HV_b = S₂

be the GSVD of (W, A) and (W^H, A^H), respectively, where U_a, V_a, U_b, V_b are unitary, T₁, T₂ are nonsingular, C₁, S₁, C₂, S₂ are positive diagonal (19, p.466).

Then V₁ and V₂ can be solved by V₁ = U_aC₁⁻¹S₁V_a^H, V₂ = V_bS₂C₂⁻¹U_b^H. A detailed error analysis of D-SDA will be given in Appendix.

The sequence {(A_k, G_k, H_k)} generated by Algorithm 2.2 satisfies the following recursive formula

A_k+1= A_k(I + G_kH_k)⁻¹A_k, (2.9a)

G_k+1= G_k+ A_kG_k(I + H_kG_k)⁻¹A^H_k, (2.9b) H_k+1= H_k+ A^H_kH_k(I + G_kH_k)⁻¹A_k. (2.9c)

(10)

3 Structured doubling algorithm for CAREs

To solve CARE (1.1), a structured doubling algorithm was first proposed by Kimura (27) using a Cayley transformation. With an appropriate parameter γ > 0, the Hamiltonian matrix H in (1.3) can be transformed to a symplectic pair (M, L) ≡ (H + γI, H − γI) (38; 39), and then simplifies to a symplectic pair (M₀, L₀) in the SSF form. Here

M₀ =



 A₀ 0

−H0 I



, L₀ =



I G₀ 0 A^H₀



, (3.1)

with

A0= I + 2γ(Aγ+ GA^−H_γ H)⁻¹, (3.2a)

G₀= 2γA⁻¹_r G(A^H_γ + HA⁻¹_γ G)⁻¹, (3.2b) H₀= 2γ(A^H_γ + HA⁻¹_γ G)⁻¹HA⁻¹_γ , (3.2c) and A_γ ≡ A − γI. The DARE associated with the symplectic pair (M₀, L₀) is

X = A^H₀ X(I + G₀X)⁻¹A₀+ H₀

on which Algorithm 2.1 or 2.2 can then be applied. For details on how a suitable γ is chosen for the Cayley transformation, see (11).

Algorithm 3.1 (C-MDFM/C-SDA for CAREs) Input: A, G, H; τ (a small tolerance);

Output: a stabilizing Hermitian solution X to CARE.

(I) Find an appropriate value γ > 0 so that Aγ and

A_γ+ GA^−H_γ H are well-conditioned (see (11) for details).

(II) Initialize A0 ← I + 2γ(Aγ+ GA^−H_γ H)⁻¹, G₀ ← 2γA⁻¹_r G(A^H_γ + HA⁻¹_γ G)⁻¹, H0 ← 2γ(A^H_γ + HA⁻¹_γ G)⁻¹HA⁻¹_γ ;

(III) Call Algorithm 2.1 (D-MDFM) or 2.2 (D-SDA).

(11)

4 Convergence of structured doubling algorithms

Let

M =



 A 0

−H I



, L =



I G 0 A^H



, (4.1)

where G = G^H and H = H^H. In the light of Definition 1.4, we assume that the matrix pair (M, L) is regular (i.e., det(M − λL) 6≡ 0) and satisfies Assumption (A2). In this section, we shall show that under these assumptions Algorithm 2.1 or 2.2 converges to a weakly stabilizing Hermitian solution with property (P) for DARE (1.2). A similar proof can be applied to the convergence of Algorithm 3.1 to a weakly stabilizing Hermitian solution with property (P) for the CARE when H satisfies Assumption (A1).

Denote the Jordan block of size p corresponding to a unimodular eigenvalue ω ≡ e^ıθ by

Jω,p =







ω 1 0 · · · 0 0 ω 1 . .. ...

... ... ... ... 0 ... . .. ... 1 0 · · · 0 ω







p×p

. (4.2)

We now quote or prove some useful lemmas. For example, the Jordan block J_ω,p to the power of 2^k can be explicitly evaluated.

Lemma 4.1 (19, pp. 557) Let Jω,p be given by (4.2). Then

J_ω,p²^k =







γ_1,k γ_2,k · · · γ_p,k 0 γ_1,k . .. ...

... ... ... γ_2,k 0 · · · 0 γ_1,k







, (4.3a)

where

(12)

γ_i,k = 1 (i − 1)!

dⁱ⁻¹ dxⁱ⁻¹x²^k

¯¯

¯x=ω

= 2^k(2^k− 1) · · · (2^k− i + 2)

(i − 1)! ω²^k⁻ⁱ⁺¹, (4.3b) for i = 1, . . . , p.

With p = 2m, let

Γ_k,m=







γ_m+1,k γ_m+2,k · · · γ_2m−1,k γ_2m,k γm,k . .. ... γ2m−1,k

... . .. ... ... ...

γ3,k . .. ... γ_m+2,k γ_2,k γ_3,k · · · γ_m,k γ_m+1,k







m×m

(4.4)

≡ J_ω,2m²^k (1 : m, m + 1 : 2m),

in which γ_i,k are defined in (4.3b), for i = 2, . . . , 2m.

To show that Γ_k,m is invertible, we first prove the following lemma.

Lemma 4.2 Let 2 ≤ r ≤ m and

F_r(m) =







f₁₁ f₁₂ · · · f_1r f₂₁ f₂₂ · · · f_2r ... ... ...

f_r1 f_r2 · · · f_rr







∈ R^r×r (4.5)

where

f_ij =







1, if j = 1,

Q_i+j−2

ν=i (m + r − ν), if 2 ≤ j ≤ r, for i = 1, 2, . . . , r. Then

|det(F_r(m))| =

Yr

ν=1

(ν − 1)!. (4.6)

Proof. Since

(13)

det(F₂(m)) =

¯¯

¯

1 m + 1

1 m

¯¯

¯

= −1,

(4.6) is true for r = 2. Suppose that

|det(Fr−1(m))| =

r−1Y

ν=1

(ν − 1)!.

Eliminating the first to (r − 1)th entries in the first column of F_r(m) by elementary row operations, we obtain







Ir−2 0 0 0 1 −1 0 0 1







· · ·







1 −1 0 0 1 0 0 0 Ir−2







F_r(m)

=



0 F^br−1(m)

1 m, m(m − 1), · · · , m(m − 1) · · · (m − r + 2)



,

where F^br−1(m) ≡ [f^bij] ∈ R(r−1)×(r−1) with

fb_ij =







1, if j = 1,

j^Q^i+j−2_ν=i [m + (r − 1) − ν], if 2 ≤ j ≤ r − 1, for i = 1, 2, . . . , r − 1. Using the factorization

Fb_r−1(m) = F_r−1(m)diag {1, 2, · · · , r − 1} , we have

|det(F_r(m))| = (r − 1)! |det(F_r−1(m))| =

Yr

ν=1

(ν − 1)!.

The proof is completed by mathematical induction.

Definition 4.1 Given ` ∈ Z. An m × m Toeplitz matrix T = [tpq]^m_p,q=1 can be written in the form

T = [t_i]^2m+`−2_i=`

(14)

with t_i = t_pq, where i = (m − 1) − (p − q) + ` and p, q = 1, . . . , m.

Since |ω| = 1, we assume w.l.o.g. that ω = 1 in the following discussion for convenience.

Lemma 4.3 The Toeplitz matrix Γ_k,m in (4.4) is nonsingular for sufficiently large k and satisfies

kΓ⁻¹_k,mJ_ω²^kk = O(2^−k), kJ_ω²^kΓ⁻¹_k,mJ_ω²^kk = O(2^−k), (4.7) where J_ω ≡ J_ω,m is defined in (4.2).

Proof. The Toeplitz matrix T =^h_i!¹ⁱ^2m−1

i=1 can be factorized by

T = diag

( 1

(2m − 1)!, · · · , 1 m!

)

F_m(m)Π, (4.8)

where Π = [en, · · · , e1] and Fm(m) is given by (4.5) with r = m. From Lemma 4.2, T is nonsingular. We write Γ_k,m of (4.4) in the form

Γ_k,m= 2^kD₁T D˜ ₂, (4.9)

where D₁ = diag^³2^(m−1)k, · · · , 2^k, 1^´, D₂ = diag^³1, 2^k, · · · , 2^(m−1)k^´ and ˜T = T + O(2^−k). Similarly, J_ω²^k = D₂⁻¹J˜₁D₂ = D₁J˜₂D⁻¹₁ , where ˜J₁ = I + O(2^−k) and ˜J₂ = I + O(2^−k). With these formulae, (4.7) obviously holds.

For the unimodular eigenvalues ωj = e^ıθ^j of (M, L) with an even partial multiplicity p = 2m_j, we have

Jωj,2mj =



J_ω_j_,m_j Γ_1,m_j 0_m_j J_ω_j_,m_j



, Γ1,mj = emje^T₁, (4.10)

for j = 1, . . . , r. From the symplectic Kronecker’s Theorem for (M, L) (see (32)), there exist a symplectic matrixZ (i.e.,^b Z^b^HJZ = J) and a nonsingular^b Q such thatb

(15)

QMb Z =^b



J_s⊕ J₁ 0_`⊕Γ^b₁ 0_n I_`⊕ J₁^−H



, (4.11a)

QLb Z =^b



I_`⊕ I_µ 0_n 0_n J_s^H ⊕ I_µ



, (4.11b)

where J_s ∈ C^`×` consists of asymptotically stable Jordan blocks (i.e., ρ(J_s) <

1),

J₁= J_ω₁_,m₁ ⊕ · · · ⊕ J_ω_r_,m_r, (4.12) Γb1=Γ^b1,m1 ⊕ · · · ⊕Γ^b1,mr

with Γ^b_1,m_j = e_m_je^T_m_j(j = 1, . . . , r), ` = n − (m₁ + · · · + m_r) ≡ n − µ and ⊕ denotes the direct sum of matrices.

Based on the standard Weierstrass form and the special eigen-structure shown in (4.11), there exists a suitable nonsingular W_j ∈ C^m^j^×m^j, for j = 1, . . . , r, such that



I_m_j 0 0 W_j⁻¹







J_ω_j_,m_j Γ^b_1,m_j 0 J_ω^−H_j_,m_j







I_m_j 0 0 W_j



=



J_ω_j_,m_j Γ_1,m_j 0 J_ω_j_,m_j



. (4.13)

Let

Z =Z(I^b _n+`⊕ W₁ ⊕ · · · ⊕ W_r), Q = (I_n+`⊕ W₁⁻¹⊕ · · · ⊕ W_r⁻¹)Q. (4.14)^b Then from (4.13)–(4.14), the equations (4.11a) and (4.11b), respectively, be- come

QMZ =



Js⊕ J1 0`⊕ Γ1

0_n I_`⊕ J₁



≡ J_M, (4.15a)

QLZ =



In 0n

0_n J_s^H ⊕ I_µ



≡ J_L, (4.15b)

where Γ₁ ≡ Γ_1,m₁⊕ · · · ⊕ Γ_1,m_r with Γ_1,m_j being given in (4.10). Since J_M and J_L in (4.15) commute with each other and from (4.15), one can derive

MZJ_L= Q⁻¹J_LJ_M = LZJ_M. (4.16)

From (4.11) and (4.14), it follows that span{Z(:, 1 : n)} forms the unique

(16)

stable Lagrangian subspace of (M, L) corresponding to J_s⊕ J₁. On the other hand, if we interchange the roles of M and L in (4.15) and consider the symplectic pair (L, M), there are nonsingular matrices P and Y such that

PLY = ¯J_M, PMY = ¯J_L. (4.17)

Similar arguments also produce

LY ¯J_L= MY ¯J_M, (4.18)

where span{Y(:, 1 : n)} forms the unique stable Lagrangian subspace of (L, M) corresponding to ¯J_s⊕ ¯J₁.

Let {(M_k, L_k)}^∞_k=0be the sequence of symplectic pairs generated by Algorithm 2.1 (see (2.4)), or {(M_k, L_k)}^∞_k=0 be the sequence of symplectic pairs in SSF form with

M_k=



 A_k 0

−Hk I



, L_k=



I G_k 0 A^H_k



 (4.19)

generated by Algorithm 2.2 (see (2.9)). With M₀ = M, L₀ = L, from (4.16) as well as (2.1)–(2.2) or (2.6)–(2.7), it follows that

M₁ZJ_L²= M_∗M₀ZJ_L² = M_∗L₀ZJ_MJ_L = L_∗M₀ZJ_LJ_M (4.20)

= L∗L0ZJ_M² = L1ZJ_M² . By induction, we have

M_kZJ_L²^k = L_kZJ_M²^k. (4.21)

By the result of Lemma 4.1 with p = 2m_j and the definitions of Γ_k,m, J_M and J_L in (4.4), (4.15a) and (4.15b), (4.21) can be rewritten as

M_kZ



I_n 0_n 0_n (J_s^H)²^k ⊕ I_µ



= L_kZ



J_s²^k⊕ J₁²^k 0_`⊕ Γ_k 0_n I_`⊕ J₁²^k



, (4.22)

where

Γ_k ≡ Γ_k,m₁ ⊕ · · · ⊕ Γ_k,m_r (4.23)

(17)

with Γ_k,m_j being defined as in (4.4), for j = 1, . . . , r. Similarly, from (4.18) it also holds

LkY ¯J_L²^k = MkY ¯J_M²^k. (4.24)

Lemma 4.4 Let J₁ and Γ_k be defined in (4.12) and (4.23), respectively. Then Γ_k is invertible and satisfies

kΓ⁻¹_k J₁²^kk = O(2^−k), kJ₁²^kΓ⁻¹_k J₁²^kk = O(2^−k). (4.25)

Proof. By Lemma 4.3.

We now partition Z and Y in (4.16) and (4.18):

Z =



Z₁ Z₃ Z2 Z4



, Y =



Y₁ Y₃ Y2 Y4



, (4.26)

where Z_i, Y_i ∈ C^n×n, for i = 1, 2, 3, 4.

Theorem 4.1 Let (M, L) be given in (1.5) satisfying (A2). Suppose the cor- responding DARE (1.2) has a weakly stabilizing Hermitian solution X with property (P). Then, the sequence {(Mk, Lk)} (see (2.4)) generated by Algo- rithm 2.1 satisfies

kM_k(:, 1 : n) + M_k(:, n + 1 : 2n)Xk (4.27)

≤ O^³ρ(J_s)²^k^´+ O^³2^−k^´→ 0, as k → ∞.

Proof. By the assumption it holds that (M, L) has a unique stable Lagrangian subspace of the form span^³[I, X^T]^T^´ satisfying (1.8), where Φ is similar to J_s⊕ J₁ as in (4.11). From (4.16) and (4.26), we have [I, X^T]^T and [Z₁^T, Z₂^T]^T spanning the same Lagrangian subspace corresponding to J_s⊕ J₁. So, Z₁ is invertible and X = Z₂Z₁⁻¹.

To show (4.27), we partition M_k, L_k conformally with (4.26) into

M_k=



M_k,1 M_k,3 Mk,2 Mk,4



, L_k=



L_k,1 L_k,3 Lk,2 Lk,4



. (4.28)

Substituting (4.26) and (4.28) into (4.22), we have

(18)

M_k,1Z₁+ M_k,3Z₂= L_k,1Z₁^³J_s²^k ⊕ J₁²^k^´+ L_k,3Z₂^³J_s²^k ⊕ J₁²^k^´, (4.29a) (M_k,1Z₃+ M_k,3Z₄)^³(J_s^H)²^k ⊕ I_µ^´= L_k,1^hZ₁(0_`⊕ Γ_k) + Z₃^³I_`⊕ J₁²^k^´i

+ Lk,3

hZ2(0`⊕ Γk) + Z4

³I`⊕ J₁²^k^´i,(4.29b) M_k,2Z₁+ M_k,4Z₂= L_k,2Z₁^³J_s²^k ⊕ J₁²^k^´+ L_k,4Z₂^³J_s²^k ⊕ J₁²^k^´,

(4.29c) (M_k,2Z₃+ M_k,4Z₄)^³(J_s^H)²^k ⊕ I_µ^´= L_k,2^hZ₁(0_`⊕ Γ_k) + Z₃^³I_`⊕ J₁²^k^´i

+ L_k,4^hZ₂(0_`⊕ Γ_k) + Z₄^³I_`⊕ J₁²^k^´i.(4.29d) Postmultiplying (4.29b) by ^³0_`⊕ Γ⁻¹_k J₁²^k^´Z₁⁻¹ to eliminate L_k,1Z₁^³0_`⊕ J₁²^k^´ and L_k,3Z₂^³0_`⊕ J₁²^k^´in (4.29a), we get

M_k,1+ M_k,3X = (M_k,1Z₃+ M_k,3Z₄)^³0_`⊕ Γ⁻¹_k J₁²^k^´Z₁⁻¹ (4.30)

− (L_k,1Z₃+ L_k,3Z₄)^³0_`⊕ J₁²^kΓ⁻¹_k J₁²^k^´Z₁⁻¹ + (L_k,1Z₁+ L_k,3Z₂)^³J_s²^k ⊕ 0_µ^´Z₁⁻¹. Similarly, from (4.29c) and (4.29d), we obtain

M_k,2+ M_k,4X = (M_k,2Z₃+ M_k,4Z₄)^³0_`⊕ Γ⁻¹_k J₁²^k^´Z₁⁻¹ (4.31)

− (L_k,2Z₃+ L_k,4Z₄)^³0_`⊕ J₁²^kΓ⁻¹_k J₁²^k^´Z₁⁻¹ + (L_k,2Z₁+ L_k,4Z₂)^³J_s²^k ⊕ 0_µ^´Z₁⁻¹.

By Algorithm 2.1 or 2.4, we deduce that L_k = L_∗,k−1L_k−1, M_k= M_∗,k−1M_k−1 with kL∗,k−1k ≤ 1 and kM∗,k−1k ≤ 1 for all k. So we have kLkk ≤ kL0k, kM_kk ≤ kM₀k and therefore kM_k,ik and kL_k,ik are bounded, for i = 1, . . . , 4.

From (4.30)–(4.31) and Lemma 4.4, assertion (4.27) follows.

Theorem 4.2 Let (M, L) be given in (1.5) satisfying (A2). Suppose the cor- responding DARE (1.2) and the dual DARE

Y = AY (I + HY )⁻¹A^H + G (4.32)

have weakly stabilizing Hermitian solutions X and Y with property (P), re- spectivly. If the sequence {(A_k, G_k, H_k)} generated by Algorithm 2.2 (see (2.9)) is well-defined, then

(i) kA_kk ≤ O(2^−k) → 0, as k → ∞,

(ii) kX − H_kk ≤ O^³ρ(J_s)²^k^´+ O(2^−k) → 0, as k → ∞,