## Structured doubling algorithms for weakly stabilizing Hermitian solutions of algebraic

## Riccati equations

### Tsung-Ming Huang

^{a,1}### , Wen-Wei Lin

^{b,}*∗*

a*Department of Mathematics, National Taiwan Normal University, Taipei 11677,*
*Taiwan.*

b*Department 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 alge-
braic 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.

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*^{H}*X + XA + H = 0,* (1.1)

or

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

*X = A*^{H}*X(I + GX)*^{−1}*A + H,* (1.2)

*where A, G, H ∈ C*^{n×n}*with G = G*^{H}*, H = H*^{H}*and I ≡ I**n* 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

*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 σ(Φ) ⊆ °*_{1} *(the closed unit disk) and (I + GX) is invertible. The*
particular invariant subspaces spanned by ^{h}*I, X*^{T}^{i}* ^{T}* in (1.7) and (1.8) are
usually referred to as stable Lagrangian subspaces.

*Definition 1.1 A subspace U ⊆ C*^{2n}*with dimension n is called an H-stable*
*Lagrangian subspace, if U satisfies that (i) HU ⊆ U, (ii) U is isotropic; i.e.,*
*x*^{H}*Jy = 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*^{2n}*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:*

*(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 fol- lowing 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)*^{−1}*A) satisfies that σ(Φ) ⊆ °*_{1} *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 dif-
ficulties mentioned above, another stable and structured algorithm has been
developed in (33), preprocessing to deflate all purely imaginary eigenvalues.

*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 algo-
rithms 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** ^{−1}*-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 struc-
tured 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.

*Throughout this paper, we denote A** ^{H}* = ¯

*A*

^{T}*the conjugate transpose of A ∈*C

^{n×n}*, ı =*

*√*

*−1, I ≡ I**n* *and 0 ≡ 0**n* (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}*]

*from*

^{T}*Q*

*L*

*−M*

*≡*

*Q*_{11} *Q*_{12}
*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 := L*b _{∗}*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 = M*c _{∗}*Mx = λM*_{∗}*Lx = λL*_{∗}*Mx = λ*^{2}*L*_{∗}*Lx = λ*^{2}*Lx,*^{b} (2.3)

*assuming that Mx = λLx.*

Algorithm 2.1 (D-MDFM for DAREs)

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

*Output: a weakly stabilizing Hermitian solution X to DARE.*

*Initialize: R ← 0**2n**, M ←*

*A 0*

*−H I*

*, L ←*

*I G*
*0 A*^{H}

*;*

*Repeat: Compute the QR-factorization:*

*Q*_{11} *Q*_{12}
*M*_{∗}*L*_{∗}

*L*

*−M*

=

*R*c

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

*M**k+1**= M**∗,k**M**k**,* *L**k+1* *= L**∗,k**L**k**,* (2.4)
where

*Q*_{11,k}*Q*_{12,k}*M*_{∗,k}*L*_{∗,k}

*L*_{k}

*−M*_{k}

=

*R*c_{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)** ^{−1}* 0

*−A*^{H}*(I + HG)*^{−1}*H I*

*,* *L** _{∗}* =

*I AG(I + HG)*^{−1}*0 A*^{H}*(I + HG)*^{−1}

(2.5)

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 0*b

*−H I*^{c}

*,* (2.7)

where

*A = A(I + GH)*b ^{−1}*A,* (2.8a)

*G = G + AG(I + HG)*b ^{−1}*A*^{H}*,* (2.8b)

*H = H + A*c ^{H}*(I + HG)*^{−1}*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 ^{2}*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 {X**k**}, the doubling algorithm produces {X*_{2}^{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)

*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*_{1} *= A, V*_{2}*W = A for V*_{1} *and V*_{2}*;*
*Set G ← G + V*_{2}*GA*^{H}*,*

*H ← H + A*c ^{H}*HV*_{1}*,*
*A ← AV*_{1}

*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*_{1} *= A, V*_{2}*W = A for V*_{1} *and V*_{2} *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*_{1}*W U*_{a}*= C*_{1}*,*
*T*_{1}*AV*_{a}*= S*_{1}*,*

*T*_{2}*W*^{H}*U*_{b}*= C*_{2}*,*
*T*_{2}*A*^{H}*V*_{b}*= S*_{2}

*be the GSVD of (W, A) and (W*^{H}*, A*^{H}*), respectively, where U*_{a}*, V*_{a}*, U*_{b}*, V*_{b}*are*
*unitary, T*_{1}*, T*_{2} *are nonsingular, C*_{1}*, S*_{1}*, C*_{2}*, S*_{2} *are positive diagonal (19, p.466).*

*Then V*_{1} *and V*_{2} *can be solved by V*_{1} *= U*_{a}*C*_{1}^{−1}*S*_{1}*V*_{a}^{H}*, V*_{2} *= V*_{b}*S*_{2}*C*_{2}^{−1}*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*_{k}*H** _{k}*)

^{−1}*A*

_{k}*,*(2.9a)

*G*_{k+1}*= G*_{k}*+ A*_{k}*G*_{k}*(I + H*_{k}*G** _{k}*)

^{−1}*A*

^{H}

_{k}*,*(2.9b)

*H*

_{k+1}*= H*

_{k}*+ A*

^{H}

_{k}*H*

_{k}*(I + G*

_{k}*H*

*)*

_{k}

^{−1}*A*

_{k}*.*(2.9c)

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*_{0}*, L*_{0}) in the SSF form. Here

*M*_{0} =

*A*_{0} 0

*−H*0 *I*

*,* *L*_{0} =

*I G*_{0}
*0 A*^{H}_{0}

*,* (3.1)

with

*A*0*= I + 2γ(A**γ**+ GA*^{−H}_{γ}*H)*^{−1}*,* (3.2a)

*G*_{0}*= 2γA*^{−1}_{r}*G(A*^{H}_{γ}*+ HA*^{−1}_{γ}*G)*^{−1}*,* (3.2b)
*H*_{0}*= 2γ(A*^{H}_{γ}*+ HA*^{−1}_{γ}*G)*^{−1}*HA*^{−1}_{γ}*,* (3.2c)
*and A*_{γ}*≡ A − γI. The DARE associated with the symplectic pair (M*_{0}*, L*_{0}) is

*X = A*^{H}_{0} *X(I + G*_{0}*X)*^{−1}*A*_{0}*+ H*_{0}

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 A*0 *← I + 2γ(A**γ**+ GA*^{−H}_{γ}*H)*^{−1}*,*
*G*_{0} *← 2γA*^{−1}_{r}*G(A*^{H}_{γ}*+ HA*^{−1}_{γ}*G)*^{−1}*,*
*H*0 *← 2γ(A*^{H}_{γ}*+ HA*^{−1}_{γ}*G)*^{−1}*HA*^{−1}_{γ}*;*

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

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

*can be explicitly evaluated.*

^{k}*Lemma 4.1 (19, pp. 557) Let J**ω,p* *be given by (4.2). Then*

*J*_{ω,p}^{2}* ^{k}* =

*γ*_{1,k}*γ*_{2,k}*· · · γ*_{p,k}*0 γ*_{1,k}*. .. ...*

*... ... ... γ*_{2,k}*0 · · · 0 γ*_{1,k}

*,* (4.3a)

*where*

*γ** _{i,k}* = 1

*(i − 1)!*

*d*^{i−1}*dx*^{i−1}*x*^{2}^{k}

¯¯

¯¯

¯*x=ω*

= 2* ^{k}*(2

^{k}*− 1) · · · (2*

^{k}*− i + 2)*

*(i − 1)!* *ω*^{2}^{k}^{−i+1}*,* (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}^{2}^{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*_{11} *f*_{12} *· · · f*_{1r}*f*_{21} *f*_{22} *· · · 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))| =*

Y*r*

*ν=1*

*(ν − 1)!.* (4.6)

Proof. Since

*det(F*_{2}*(m)) =*

¯¯

¯¯

¯¯

¯

*1 m + 1*

1 *m*

¯¯

¯¯

¯¯

¯

*= −1,*

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

*|det(F**r−1**(m))| =*

*r−1*Y

*ν=1*

*(ν − 1)!.*

*Eliminating the first to (r − 1)th entries in the first column of F*_{r}*(m) by*
elementary row operations, we obtain

*I**r−2* 0 0
*0 1 −1*
0 0 1

*· · ·*

*1 −1 0*
0 1 0
0 *0 I**r−2*

*F*_{r}*(m)*

=

0 *F*^{b}*r−1**(m)*

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

*,*

where *F*^{b}*r−1**(m) ≡ [f*^{b}*ij**] ∈ R**(r−1)×(r−1)* with

*f*b* _{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*

*F*b_{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))| =*

Y*r*

*ν=1*

*(ν − 1)!.*

The proof is completed by mathematical induction.

*Definition 4.1 Given ` ∈ Z. An m × m Toeplitz matrix T = [t**pq*]^{m}_{p,q=1}*can be*
*written in the form*

*T = [t** _{i}*]

^{2m+`−2}

_{i=`}*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Γ*^{−1}_{k,m}*J*_{ω}^{2}^{k}*k = O(2*^{−k}*),* *kJ*_{ω}^{2}* ^{k}*Γ

^{−1}

_{k,m}*J*

_{ω}^{2}

^{k}*k = O(2*

^{−k}*),*(4.7)

*where J*

_{ω}*≡ J*

_{ω,m}*is defined in (4.2).*

*Proof. The Toeplitz matrix T =*^{h}_{i!}^{1}^{i}^{2m−1}

*i=1* can be factorized by

*T = diag*

( 1

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

)

*F*_{m}*(m)Π,* (4.8)

*where Π = [e**n**, · · · , e*1*] and F**m**(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

^{k}*D*

_{1}

*T D*˜

_{2}

*,*(4.9)

*where D*_{1} = diag^{³}2^{(m−1)k}*, · · · , 2*^{k}*, 1*^{´}*, D*_{2} = diag^{³}*1, 2*^{k}*, · · · , 2*^{(m−1)k}^{´} and ˜*T =*
*T + O(2*^{−k}*). Similarly, J*_{ω}^{2}^{k}*= D*_{2}^{−1}*J*˜_{1}*D*_{2} *= D*_{1}*J*˜_{2}*D*^{−1}_{1} , where ˜*J*_{1} *= I + O(2** ^{−k}*)
and ˜

*J*

_{2}

*= I + O(2*

*). With these formulae, (4.7) obviously holds.*

^{−k}*For the unimodular eigenvalues ω**j* *= e*^{ıθ}^{j}*of (M, L) with an even partial*
*multiplicity p = 2m** _{j}*, we have

*J**ω**j**,2m**j* =

*J*_{ω}_{j}_{,m}* _{j}* Γ

_{1,m}*0*

_{j}

_{m}

_{j}*J*

_{ω}

_{j}

_{,m}

_{j}

*,* Γ*1,m**j* *= e**m**j**e*^{T}_{1}*,* (4.10)

*for j = 1, . . . , r. From the symplectic Kronecker’s Theorem for (M, L) (see*
(32)), there exist a symplectic matrix*Z (i.e.,*^{b} *Z*^{b}^{H}*JZ = J) and a nonsingular*^{b}
*Q such that*b

*QM*b *Z =*^{b}

*J*_{s}*⊕ J*_{1} 0_{`}*⊕*Γ^{b}_{1}
0_{n}*I*_{`}*⊕ J*_{1}^{−H}

*,* (4.11a)

*QL*b *Z =*^{b}

*I*_{`}*⊕ I** _{µ}* 0

*0*

_{n}

_{n}*J*

_{s}

^{H}*⊕ I*

_{µ}

*,* (4.11b)

*where J*_{s}*∈ C*^{`×`}*consists of asymptotically stable Jordan blocks (i.e., ρ(J*_{s}*) <*

1),

*J*_{1}*= J*_{ω}_{1}_{,m}_{1} *⊕ · · · ⊕ J*_{ω}_{r}_{,m}_{r}*,* (4.12)
Γb1=Γ^{b}*1,m*1 *⊕ · · · ⊕*Γ^{b}*1,m**r*

with Γ^{b}_{1,m}_{j}*= e*_{m}_{j}*e*^{T}_{m}_{j}*(j = 1, . . . , r), ` = n − (m*_{1} *+ · · · + 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}

^{−1}

*J*_{ω}_{j}_{,m}* _{j}* Γ

^{b}

_{1,m}*0*

_{j}*J*

_{ω}

^{−H}

_{j}

_{,m}

_{j}

*I*_{m}* _{j}* 0

*0 W*

_{j}

=

*J*_{ω}_{j}_{,m}* _{j}* Γ

_{1,m}*0*

_{j}*J*

_{ω}

_{j}

_{,m}

_{j}

*.* (4.13)

Let

*Z =Z(I*^{b} _{n+`}*⊕ W*_{1} *⊕ · · · ⊕ W*_{r}*),* *Q = (I*_{n+`}*⊕ W*_{1}^{−1}*⊕ · · · ⊕ W*_{r}* ^{−1}*)

*Q. (4.14)*

^{b}Then from (4.13)–(4.14), the equations (4.11a) and (4.11b), respectively, be- come

*QMZ =*

*J**s**⊕ J*1 0*`**⊕ Γ*1

0_{n}*I*_{`}*⊕ J*_{1}

*≡ J*_{M}*,* (4.15a)

*QLZ =*

*I**n* 0*n*

0_{n}*J*_{s}^{H}*⊕ I*_{µ}

*≡ J*_{L}*,* (4.15b)

where Γ_{1} *≡ Γ*_{1,m}_{1}*⊕ · · · ⊕ Γ*_{1,m}* _{r}* with Γ

_{1,m}

_{j}*being given in (4.10). Since J*

*and*

_{M}*J*

*in (4.15) commute with each other and from (4.15), one can derive*

_{L}*MZJ*_{L}*= Q*^{−1}*J*_{L}*J*_{M}*= LZJ*_{M}*.* (4.16)

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

*stable Lagrangian subspace of (M, L) corresponding to J*_{s}*⊕ J*_{1}. 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*_{1}.

*Let {(M*_{k}*, L*_{k}*)}*^{∞}* _{k=0}*be the sequence of symplectic pairs generated by Algorithm

*2.1 (see (2.4)), or {(M*

_{k}*, L*

_{k}*)}*

^{∞}*be the sequence of symplectic pairs in SSF form with*

_{k=0}*M** _{k}*=

*A** _{k}* 0

*−H**k* *I*

*,* *L** _{k}*=

*I G*_{k}*0 A*^{H}_{k}

(4.19)

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

*M*_{1}*ZJ*_{L}^{2}*= M*_{∗}*M*_{0}*ZJ*_{L}^{2} *= M*_{∗}*L*_{0}*ZJ*_{M}*J*_{L}*= L*_{∗}*M*_{0}*ZJ*_{L}*J** _{M}* (4.20)

*= L**∗**L*0*ZJ*_{M}^{2} *= L*1*ZJ*_{M}^{2} *.*
By induction, we have

*M*_{k}*ZJ*_{L}^{2}^{k}*= L*_{k}*ZJ*_{M}^{2}^{k}*.* (4.21)

*By the result of Lemma 4.1 with p = 2m** _{j}* and the definitions of Γ

_{k,m}*, J*

*and*

_{M}*J*

*in (4.4), (4.15a) and (4.15b), (4.21) can be rewritten as*

_{L}*M*_{k}*Z*

*I** _{n}* 0

*0*

_{n}

_{n}*(J*

_{s}*)*

^{H}^{2}

^{k}*⊕ I*

_{µ}

*= L*_{k}*Z*

*J*_{s}^{2}^{k}*⊕ J*_{1}^{2}* ^{k}* 0

_{`}*⊕ Γ*

*0*

_{k}

_{n}*I*

_{`}*⊕ J*

_{1}

^{2}

^{k}

*,* (4.22)

where

Γ_{k}*≡ Γ*_{k,m}_{1} *⊕ · · · ⊕ Γ*_{k,m}* _{r}* (4.23)

with Γ_{k,m}_{j}*being defined as in (4.4), for j = 1, . . . , r. Similarly, from (4.18) it*
also holds

*L**k**Y ¯J*_{L}^{2}^{k}*= M**k**Y ¯J*_{M}^{2}^{k}*.* (4.24)

*Lemma 4.4 Let J*_{1} *and Γ*_{k}*be defined in (4.12) and (4.23), respectively. Then*
Γ_{k}*is invertible and satisfies*

*kΓ*^{−1}_{k}*J*_{1}^{2}^{k}*k = O(2*^{−k}*),* *kJ*_{1}^{2}* ^{k}*Γ

^{−1}

_{k}*J*

_{1}

^{2}

^{k}*k = O(2*

^{−k}*).*(4.25)

Proof. By Lemma 4.3.

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

*Z =*

*Z*_{1} *Z*_{3}
*Z*2 *Z*4

*,* *Y =*

*Y*_{1} *Y*_{3}
*Y*2 *Y*4

*,* (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 {(M**k**, L**k**)} (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}*)

^{2}

^{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*

_{1}

*as in (4.11). From (4.16) and (4.26), we have [I, X*

*]*

^{T}

^{T}*and [Z*

_{1}

^{T}*, Z*

_{2}

*]*

^{T}

^{T}*spanning the same Lagrangian subspace corresponding to J*

_{s}*⊕ J*

_{1}

*. So, Z*

_{1}is

*invertible and X = Z*

_{2}

*Z*

_{1}

*.*

^{−1}*To show (4.27), we partition M*_{k}*, L** _{k}* conformally with (4.26) into

*M** _{k}*=

*M*_{k,1}*M*_{k,3}*M**k,2* *M**k,4*

*,* *L** _{k}*=

*L*_{k,1}*L*_{k,3}*L**k,2* *L**k,4*

*.* (4.28)

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

*M*_{k,1}*Z*_{1}*+ M*_{k,3}*Z*_{2}*= L*_{k,1}*Z*_{1}^{³}*J*_{s}^{2}^{k}*⊕ J*_{1}^{2}^{k}^{´}*+ L*_{k,3}*Z*_{2}^{³}*J*_{s}^{2}^{k}*⊕ J*_{1}^{2}^{k}^{´}*,*
(4.29a)
*(M*_{k,1}*Z*_{3}*+ M*_{k,3}*Z*_{4})^{³}*(J*_{s}* ^{H}*)

^{2}

^{k}*⊕ I*

_{µ}^{´}

*= L*

_{k,1}^{h}

*Z*

_{1}(0

_{`}*⊕ Γ*

_{k}*) + Z*

_{3}

^{³}

*I*

_{`}*⊕ J*

_{1}

^{2}

^{k}^{´i}

*+ L**k,3*

h*Z*2(0*`**⊕ Γ**k**) + Z*4

³*I**`**⊕ J*_{1}^{2}^{k}^{´i}*,(4.29b)*
*M*_{k,2}*Z*_{1}*+ M*_{k,4}*Z*_{2}*= L*_{k,2}*Z*_{1}^{³}*J*_{s}^{2}^{k}*⊕ J*_{1}^{2}^{k}^{´}*+ L*_{k,4}*Z*_{2}^{³}*J*_{s}^{2}^{k}*⊕ J*_{1}^{2}^{k}^{´}*,*

(4.29c)
*(M*_{k,2}*Z*_{3}*+ M*_{k,4}*Z*_{4})^{³}*(J*_{s}* ^{H}*)

^{2}

^{k}*⊕ I*

_{µ}^{´}

*= L*

_{k,2}^{h}

*Z*

_{1}(0

_{`}*⊕ Γ*

_{k}*) + Z*

_{3}

^{³}

*I*

_{`}*⊕ J*

_{1}

^{2}

^{k}^{´i}

*+ L*_{k,4}^{h}*Z*_{2}(0_{`}*⊕ Γ*_{k}*) + Z*_{4}^{³}*I*_{`}*⊕ J*_{1}^{2}^{k}^{´i}*.(4.29d)*
Postmultiplying (4.29b) by ^{³}0_{`}*⊕ Γ*^{−1}_{k}*J*_{1}^{2}^{k}^{´}*Z*_{1}^{−1}*to eliminate L*_{k,1}*Z*_{1}^{³}0_{`}*⊕ J*_{1}^{2}^{k}^{´}
*and L*_{k,3}*Z*_{2}^{³}0_{`}*⊕ J*_{1}^{2}^{k}^{´}in (4.29a), we get

*M*_{k,1}*+ M*_{k,3}*X = (M*_{k,1}*Z*_{3}*+ M*_{k,3}*Z*_{4})^{³}0_{`}*⊕ Γ*^{−1}_{k}*J*_{1}^{2}^{k}^{´}*Z*_{1}* ^{−1}* (4.30)

*− (L*_{k,1}*Z*_{3}*+ L*_{k,3}*Z*_{4})^{³}0_{`}*⊕ J*_{1}^{2}* ^{k}*Γ

^{−1}

_{k}*J*

_{1}

^{2}

^{k}^{´}

*Z*

_{1}

^{−1}*+ (L*

_{k,1}*Z*

_{1}

*+ L*

_{k,3}*Z*

_{2})

^{³}

*J*

_{s}^{2}

^{k}*⊕ 0*

_{µ}^{´}

*Z*

_{1}

^{−1}*.*Similarly, from (4.29c) and (4.29d), we obtain

*M*_{k,2}*+ M*_{k,4}*X = (M*_{k,2}*Z*_{3}*+ M*_{k,4}*Z*_{4})^{³}0_{`}*⊕ Γ*^{−1}_{k}*J*_{1}^{2}^{k}^{´}*Z*_{1}* ^{−1}* (4.31)

*− (L*_{k,2}*Z*_{3}*+ L*_{k,4}*Z*_{4})^{³}0_{`}*⊕ J*_{1}^{2}* ^{k}*Γ

^{−1}

_{k}*J*

_{1}

^{2}

^{k}^{´}

*Z*

_{1}

^{−1}*+ (L*

_{k,2}*Z*

_{1}

*+ L*

_{k,4}*Z*

_{2})

^{³}

*J*

_{s}^{2}

^{k}*⊕ 0*

_{µ}^{´}

*Z*

_{1}

^{−1}*.*

*By Algorithm 2.1 or 2.4, we deduce that L*_{k}*= L*_{∗,k−1}*L*_{k−1}*, M*_{k}*= M*_{∗,k−1}*M*_{k−1}*with kL**∗,k−1**k ≤ 1 and kM**∗,k−1**k ≤ 1 for all k. So we have kL**k**k ≤ kL*0*k,*
*kM*_{k}*k ≤ kM*_{0}*k and therefore kM*_{k,i}*k and kL*_{k,i}*k 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 )*^{−1}*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*_{k}*k ≤ O(2*^{−k}*) → 0, as k → ∞,*

*(ii) kX − H*_{k}*k ≤ O*^{³}*ρ(J** _{s}*)

^{2}

^{k}^{´}

*+ O(2*

^{−k}*) → 0, as k → ∞,*