HERMITIAN SOLUTIONS OF ALGEBRAIC RICCATI EQUATIONS
TSUNG-MIN HWANG AND WEN-WEI LIN
Abstract. In this paper, we propose structured doubling algorithms for the computation of weak Hermitian solutions of continuous/discrete-time algebraic Riccati equations. Under the assumptions that partial multiplicities of purely imaginary and unimodular eigenvalues (if any) of the associated Hamiltonian and symplectic pencil, respectively, are all even, we prove that the developed structured doubling algorithms converge to the desired Hermitian solutions globally and linearly. Numerical experiments show that structured doubling algorithms perform efficiently and reliably.
1. Introduction
In this paper, we investigate structured doubling algorithms for the computation of Hermitian solutions X to
(I) the continuous-time algebraic Riccati equation (CARE): −XGX + AHX+ XA + H = 0, (1.1)
or
(II) the discrete-time algebraic Riccati equation (DARE): X = AHX(I + GX)−1A+ H,
(1.2)
where A ∈ Cn×n, H = HH and G = GH = BBH ≥ 0 ∈ Cn×n with B ∈ Cn×m
being of full column rank, and I ≡ In denotes the identity matrix of order n.
Furthermore, in CARE (1.1), the pair (A, B) is c-stablizable, i.e., if wHA= µwH
and wHB = 0 with w 6= 0, then Re(µ) < 0; in DARE (1.2), the pair (A, B) is
d-stablizable, i.e., if wHA= µwH and wHB= 0 with w 6= 0, then |µ| < 1.
Equations (1.1) and (1.2) arise frequently in solving the unique “weak” sta-blizing controllers of continuous- and discrete-time H∞-optimal control systems,
respectively [13, 15, 18, 25]. In addition, several applications in Wiener filtering theory [47], network synthesis [2] and Moser-Veselov equations [9, 40] also need the Hermitian solution of CAREs.
We consider the 2n × 2n Hamiltonian matrix H associated with the CARE: H = A −G −H −AH (1.3) Date: June 20, 2006.
2000 Mathematics Subject Classification. 93B36, 93B40, 93B52, 93B60.
Key words and phrases. algebraic Riccati equation, Hermitian solution, structured doubling algorithm, purely imaginary eigenvalue, unimodular eigenvalue, global and linear convergence.
This work is partially supported by the National Science Council and the National Center for Theoretical Sciences of Taiwan.
which satisfies HJ = −JHH, 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 AH (1.5) which satisfies MJMH = LJLH. (1.6)
Note that the special symplectic pair (M, L) of the forms in (1.5) is referred to as a standard symplectic form (SSF).
It is well-known that (e.g., [30, 32, 38, 41]) if A −G −H −AH Ω1 Ω2 = Ω1 Ω2 Φ, Ω1,Ω2,Φ ∈ Cn×n, (1.7)
where the spectrum of Φ is contained in the closed left half plane, then X = Ω2Ω−11
(if Ω−1
1 exists) solves the CARE (1.1). In addition, if
A 0 −H I Ω1 Ω2 = I G 0 AH Ω1 Ω2 Φ, Ω1,Ω2,Φ ∈ Cn×n, (1.8)
where the spectrum of Φ is contained in the closed unit disk, then X = Ω2Ω−11
(if Ω−11 exists) solves the DARE (1.2). The particular invariant subspaces spanned
byΩH 1, ΩH2
H
in (1.7) and (1.8) are usually referred to as stable Lagrangian sub-spaces.
Definition 1.1. A subspace U ⊆ C2n with dimension n is called a H-stable La-grangian subspace, if U satisfies that (i) HU ⊆ U, (ii) U is isotropic, i.e., xHJ y= 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 ⊆ C2n with dimension n is called a (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 Lagrangian subspace do not always exist, when some purely imaginary eigenvalues of H and some unimodular eigenvalues of (M, L) have odd partial multiplicities, respectively. Counterexamples can be found in [41]. To guarantee the existence of H- and (M, L)-stable Lagrangian subspaces, and the existence of Hermitian solutions of CARE and DARE, respectively, throughout this paper, we assume that H in (1.3) and (M, L) in (1.5) satisfy the following assumptions.
(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.
We now describe essential equivalence statements for the existence of Hermitian solutions of CARE and DARE, respectively.
Theorem 1.1. [19, 29] Assume (A, B) in CARE (1.1) is c-stable. Then the fol-lowing statements are equivalent :
(i) The CARE (1.1) has a Hermitian solution; (ii) Assumption (A1) holds;
(iii) There is an H-stable Lagrangian subspace.
The equivalence statements for the existence of Hermitian solutions of DARE have been firstly given by [28] under three conditions related to a given matrix K so that A−BK is stable and invertible, and ψK(η) > 0, for some |η| = 1 (see Theorem
1.2 of [28] for details). In Theorem 1.2 we give different equivalence statements for a Hermitian solution of DARE under condition (1.9). These equivalence statements have a correspondence to the statements in Theorem 1.1.
Theorem 1.2. Assume (A, B) in DARE (1.2) is d-stable. Suppose that there is an α ∈ C with |α| = 1 such that
λmin(BHA−Hα HA−1α B) > −1
(1.9)
where Aα= I − ¯αA. Then the following statements are equivalent:
(i) The DARE (1.2) has a Hermitian solution; (ii) Assumption (A2) holds;
(iii) There is an (M, L)-stable Lagrangian subspace. Proof. See Appendix.
For the continuous-time case, a well-known backward stable algorithm care, pro-posed by [30], computes a Hermitian solution X for CARE by applying QR algo-rithm to H with reordering [43]. Unfortunately, the QR algoalgo-rithm preserves neither the Hamiltonian structure nor the associated splitting eigenvalues. When H in (1.3) has no purely imaginary eigenvalues, a numerical strongly stable method has been proposed by [10] for computing the Hamiltonian Schur form of H, and therefore, the H-stable Lagrangian subspace. Efficient structured doubling algorithms using an appropriate Cayley transform [11, 27], as well as matrix sign function methods [5, 8, 14, 17] have been developed for solving the unique positive semidefinite so-lution of CARE (1.1). When H in (1.3) satisfies assumption (A1), an eigenvector deflation technique proposed by [13] guarantees that the eigenvalues appear with the correct pairing. This is certainly an advantages over the QR algorithm, but the method ignores most of the structure of the problem during computation. A struc-tured 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] as a preprocessing to deflate all purely imaginary eigenvalues and to get a reduced Hamiltonian matrix with no purely imaginary eigenvalues. When H satisfies (A1) with even=2, an efficient Newton’s method has been developed in [22] for solving CARE with globally linear convergence.
For the discrete-time case, a well-known backward stable algorithm dare, pro-posed by [37, 42, 48], computes a Hermitian solution X for DARE by applying QZ algorithm to (M, L) with reordering. Unfortunately, this algorithm does not take into account the symplectic structure of (M, L). Non-structure-preserving itera-tive processes loosen the symplectic structure, this may cause that the algorithms
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 [3] based on an acceleration scheme of the fixed point iteration for (1.2). Quadratic convergence of doubling algorithm has been shown in [26] and [35] with different approach. 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 DARE. Both of the symplectic structure in MDFM and SDA are preserved in SSF form at each iterative step. However, the symplectic structure in MDFM can only be preserved in exact arithmetic. When (M, L) in (1.5) satisfies (A2), a structured algorithm has been developed in [33] to deflate all unimodular eigenvalues as a preprocessing by the determination of their isotropic Jordan subbasis using the S +S−1-transform
of M − λL [31]. When (M, L) satisfies (A2) with even=2, an efficient Newton-type method has been proposed by [21] to solve DARE with globally linear convergence. As mentioned above, MDFM and SDA have been proposed for solving DARE (1.2) with M−λL having no unimodular eigenvalues. To solve CARE (1.1) with H having no purely imaginary eigenvalues, the Hamiltonian matrix H is converted to a symplectic pencil cM − λ bL in SSF form by an appropriate Cayley transform and then the MDFM method or the SDA algorithm can be applied. The main purpose of this paper is to apply MDFM or SDA for solving CAREs and DAREs, where the associated H in (1.3) and M − λL in (1.5) satisfy (A1) and (A2), respectively. Under these conditions, we prove the global linear convergence of MDFM and SDA algorithms.
This paper is organized as follows. In Sections 2 and 3, we describe structured doubling algorithms, D-MDFM/D-SDA and C-MDFM/C-SDA, for solving DAREs and CAREs, respectively. In Section 4, we prove that under assumptions (A1) and (A2), the structured doubling algorithms converge globally and linearly to Hermitian solutions of DAREs and CAREs. In Section 5, we test several numerical examples, illustrating the convergence behavior of MDFM, SDA and Newton-type methods. Concluding remarks are given in Section 6.
Throughout this paper, we denote AH = ¯AT the conjugate transpose of A ∈
Cn×n, ı =√−1, and I ≡ Inand 0 ≡ 0nthe identity and zero matrices, respectively,
of order n. The jth column of In is denoted by ej. k · k denotes any matrix norm,
σ(A) and ̺(A) denote the spectrum and the spectral radius of A, respectively.
2. Structured doubling algorithms for DAREs
The matrix disk function method (MDFM) in [6, 7] is developed to solve the 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 [LT,−MT]T by Q L −M ≡ Q11 Q12 M∗ L∗ L −M = R 0 , (2.1)
where Q ∈ C4n×4n is unitary and R ∈ C2n×2n is upper triangular. Let
b
L := L∗L, M := Mc ∗M.
It is easily seen that ( bL, cM) is symplectic. From (2.1)-(2.2) the pair ( cM, bL) satisfies the doubling property,
c
Mx = M∗Mx = λM∗Lx = λL∗Mx = λ2L∗Lx = λ2Lx,b
(2.3)
provided that Mx = λLx.
Algorithm 2.1 (D-MDFM for DAREs). Input: A, G, H; τ (a small tolerance); Output: a Hermitian solution X to DARE.
Initialize: R ← 02n, M ← A 0 −H I , L ← I G 0 AH ; Repeat: Compute the QR-factorization:
Q11 Q12 M∗ L∗ L −M = b R 0 ; Ifk bR − Rk ≤ τk bRk, Then solve the least
squared problem for X: −M(:, 1 : n) = M(:, n + 1 : 2n)X; Stop Else Set L ← L∗L, M ← M∗M, R ← bR; Go To Repeat End If
The SDA algorithm in [12] is developed to solve DARE (1.2) with conditions of stablizability and detectability by using a structured LU factorization instead of the QR factorization in (2.1). We refer to this step as a SLU -swap. As derived in [12], for a given pair (M, L) in SSF form (1.5), we construct
M∗= A(I + GH)−1 0 −AH(I + HG)−1H I , L∗= I AG(I + HG)−1 0 AH(I + HG)−1 (2.4)
and consequently deduce that
M∗L = L∗M.
(2.5)
We now compute L∗L and M∗M, and apply the Sherman-Morrison-Woodbury
formula to produce b L = " I Gb 0 AbH # = L∗L, M =c " b A 0 − bH I # = M∗M, (2.6) where b A = A(I + GH)−1A, (2.7a) b G = G+ AG(I + HG)−1AH, (2.7b) b H = H+ AH(I + HG)−1HA. (2.7c)
Equations in (2.6) show that the matrix pair ( cM, bL) is again in SSF form. From (2.5)–(2.6), the pair ( cM, bL) also satisfies the doubling property: cMx = λ2Lx.b
(i) Equations in (2.7) have exactly the same form as the doubling algorithm (which has been first proposed and investigated in Anderson [3] and Kimura [26]). However, the original doubling algorithm was derived as an acceler-ation scheme for the fixed-point iteracceler-ation from (1.2). Instead of producing the sequence {Xk}, the doubling algorithm produces {X2k}. Furthermore,
the convergence of the doubling algorithm was proven when A is nonsingu-lar [3], and for (A, G, H) which is reachable and detectable, or stablizable and observable [26]. A stronger convergent result of the SDA algorithm under weak conditions (stablizability and detectability) can be found in [12, 35].
(ii) Note that when the assumption (A2) in Theorem 1.2 holds, the matrix (I + GH) in (2.7) can be singular, thus bA, bG and bH in (2.7) do not exist and the SDA algorithm can break down. The numerical experience for this case shall be discussed in Section 5.
Algorithm 2.2 (D-SDA for DAREs).
Input: A, G, H; τ (a small tolerance); Output: a Hermitian solution X to DARE.
Repeat W ← I + GH;
If W is singular or very ill-conditioned, then break down; Else solve for V1, V2 from W V1= A, V2WH = G;
Set G← G + AV2AH, b H ← H + VH 1 HA, A← AV1 Ifk bH− Hk ≤ τk bHk, then X ← bH, Stop. Elseset H ← bH; End if End if Goto Repeat
3. Structured doubling algorithm for CAREs
To solve the CARE (1.1) a structured doubling algorithm was first proposed by Kimura [27] using Cayley transformation. With an appropriate γ > 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 equivalently simplifies to a symplectic pair (M0,L0) in the SSF form. Here
M0= A0 0 −H0 I , L0= I G0 0 AH 0 , (3.1) with A0 = I+ 2γ(Aγ+ GA−Hγ H)−1, (3.2a) G0 = 2γA−1r G(AHγ + HA−1γ G)−1, (3.2b) H0 = 2γ(AHγ + HA−1γ G)−1HA−1γ , (3.2c)
and Aγ ≡ A − γI. The DARE associated with the symplectic pair (M0,L0) is
X= AH
0 X(I + G0X)−1A0+ H0
on which Algorithm 2.1 or Algorithm 2.2 can then be applied. How to choose a suitable γ for Cayley transformation can be found in [11] for details.
Algorithm 3.1 (C-MDFM/C-SDA for CAREs). Input: A, G, H; τ (a small tolerance);
Output: a Hermitian solution X to CARE.
(I) Findan 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)−1,
G0← 2γA−1r G(AHγ + HA−1γ G)−1,
H0← 2γ(AHγ + HA−1γ G)−1HA−1γ , k← 0;
(III) CallAlgorithm 2.1/Algorithm 2.2.
4. Convergence of structured doubling algorithms Let M = A 0 −H I , L = I G 0 AH , (4.1)
where G = GH and H = HH. In light of the results in Theorem 1.2, we suppose
that the matrix pair (M, L) is regular , i.e., det(M−λL) 6= 0, and satisfies (A2). In this section, we shall show that under these assumptions Algorithm 2.1 or Algorithm 2.2 will converge to a weak Hermitian solution of DARE (1.2). A similar proof can be applied to the convergence of Algorithm 3.1 to a weak Hermitian solution of CARE when H satisfies (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)
The unimodular eigenvalues ωj = eıθj of (M, L) with even sizes can then be
ex-pressed by Jωj,2mj = Jωj,mj Γ0,mj 0mj Jωj,mj , Γ0,mj = emjeT1, (4.3)
for j = 1, . . . , r. By Kronecker’s Theorem [16, Chap. 12] there are nonsingular matrices Q and Z such that
QMZ = diag (J s, Jω1,m1,· · · , Jωr,mr) diag (0ℓ,Γ0,m1,· · · , Γ0,mr) 0 diag (Iℓ, Jω1,m1,· · · , Jωr,mr) (4.4)
and QLZ = " I 0 0 diag ¯Js, Im1,· · · , Imr # , (4.5)
in which Jsis the Jordan block of size ℓ associated with d-stable eigenvalues, where
ℓ= n − (m1+ · · · + mr).
In what follows, without loss of generality, we shall assume that (M, L) has only one unimodular eigenvalue with even partial multiplicities. Thus, (4.4) and (4.5) become QMZ = Js⊕ Jω,m 0ℓ⊕ Γ0,m 0 Iℓ⊕ Jω,m ≡ JM (4.6a) and QLZ = " I 0 0 J¯s⊕ Im # ≡ JL, (4.6b)
where ℓ = n − m and “⊕” denotes the direct sum. Since JMand JLcommute with
each other and from (4.6), one can derive
MZJL= Q−1JLJM= LZJM.
(4.7)
On the other hand, if we interchange the roles of M and L in (4.6) and consider the symplectic pair (L, M), there are nonsingular matrices P and Y such that
PLY = ¯JM, PMY = ¯JL. (4.8)
Similar arguments also produce
LY ¯JL= MY ¯JM. (4.9)
Let {(Mk,Lk)}∞k=0be the sequence of symplectic pairs generated by Algorithm 2.1,
or {(Mk,Lk)}∞k=0 be the sequence of symplectic pairs in SSF with
Mk = Ak 0 −Hk I , Lk= I Gk 0 AH k (4.10)
generated by Algorithm 2.2. With M0= M, L0 = L, from (4.7) as well as (2.1)–
(2.2) or (2.5)–(2.6) follows that
M1ZJL2 = M∗M0ZJL2 = M∗L0ZJMJL= L∗M0ZJLJM
(4.11)
= L∗L0ZJM2 = L1ZJM2 .
By inductive process we have MkZJ2 k L = LkZJ2 k M. (4.12)
Similarly, from (4.9) it also holds LkY ¯J2 k L = MkY ¯J2 k M. (4.13)
We now prove some useful lemmas for the convergence theorem. The Jordan block Jω,pto the power of 2k can be explicitly evaluated.
Lemma 4.1. [20, pp. 557] Let Jω,p be given by (4.2). Then Jω,p2k = γ1,k γ2,k · · · γp,k 0 γ1,k . .. ... .. . . .. ... γ2,k 0 · · · 0 γ1,k , (4.14a) where γj,k= 1 (j − 1)! dj−1 dxj−1x 2k x=ω= 2k(2k− 1) · · · (2k− j + 2) (j − 1)! ω 2k −j+1, (4.14b) for j = 1, . . . , p.
With p = 2m, Jω≡ Jω,m and by using (4.6), equation (4.12) can be expressed
by MkZ I 0 0 J¯2k s ⊕ Im = LkZ " Js2k⊕ Jω2k 0ℓ⊕ Γk 0 Iℓ⊕ J2 k ω # , (4.15) where Γk = γ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.16) ≡ Jω,2m2k (1 : m, m + 1 : 2m),
in which γj,k are defined in (4.14b), for j = 2, . . . , 2m.
To show that Γk is invertible, we first prove the following lemma.
Lemma 4.2. Let 2 ≤ r ≤ m and
Fr(m) = f11 f12 · · · f1r f21 f22 · · · f2r .. . ... ... fr1 fr2 · · · frr ∈ R r×r (4.17) where fij = 1, if j = 1, Qi+j−2 ν=i (m + r − ν), if 2 ≤ j ≤ r, for i = 1, 2, . . . , r. Then |det(Fr(m))| = r Y ν=1 (ν − 1)!. (4.18) Proof. Since det(F2(m)) = 11 mm+ 1 = −1,
(4.18) 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 Fr(m) by elementary
row operations, we obtain Ir−2 0 0 0 1 −1 0 0 1 · · · 10 −11 00 0 0 Ir−2 Fr(m) = 0 Fˆr−1(m) 1 m, m(m − 1), · · · , m(m − 1) · · · (m − r + 2) , where ˆFr−1(m) ≡ [ ˆfij] ∈ R(r−1)×(r−1)with ˆ fij= 1, if j = 1, jQiν=i+j−2[m + (r − 1) − ν], if 2 ≤ j ≤ r − 1, for i = 1, 2, . . . , r − 1. Using the factorization
ˆ Fr−1(m) = Fr−1(m)diag {1, 2, · · · , r − 1} , we have |det(Fr(m))| = (r − 1)! |det(Fr−1(m))| = r Y ν=1 (ν − 1)!. The proof is completed by mathematical induction.
Definition 4.1. Given ℓ ∈ Z. An m × m Toeplitz matrix T = [tpq]mp,q=1 can be
written in the form
T = [tj]2m+ℓ−2j=ℓ
with tj = tpq, where j = (m − 1) − (p − q) + ℓ and p, q = 1, . . . , m.
Since |ω| = 1, for convenience, in the following discussion we w.l.o.g. assume that ω = 1.
Lemma 4.3. The Toeplitz matrix Γk in (4.16) is nonsingular for k sufficiently
large and Γ−1k is Toeplitz in the big “O” sense having the form Γ−1 k = 1 O 2m2 k O 2jkm 2 −1 j=(m−1)2. (4.19)
Proof. The Toeplitz matrix T =h1 j! i2m−1 j=1 can be factorized by T = diag 1 (2m − 1)!,· · · , 1 m! Fm(m)Π, (4.20)
where Π = [en,· · · , e1] and Fm(m) is given by (4.17) with r = m. From Lemma
4.2, T is nonsingular. We write Γk of (4.16) in the form
Γk= 2jk j! + O 2(j−1)k 2m−1 j=1 . (4.21)
By (4.20)–(4.21) we see that the dominant term of det(Γk) is det(T )2m 2
k with
det(T ) 6= 0. Thus, for sufficiently large k it holds det(Γk) 6= 0. From (4.21) follows
that Γ−1k = adj(Γk) det(Γk) = 1 O 2m2k O 2jkm 2 −1 j=(m−1)2.
Lemma 4.4. Let Jω ≡ Jω,m and Γk be defined in (4.2) and (4.16), respectively.
Then kΓ−1 k J 2k ω k = O(2−k), kJ2 k ω Γ−1k J 2k ω k = O(2−k). (4.22)
Proof. From (4.14) we write the Toeplitz matrix J2k ω by Jω2k = [γjk]m+1j=−m+1, (4.23a) where γjk= 0, for j = −m, . . . , −1, O 2jk, for j = 0, 1, . . . , m − 1. (4.23b) We compute Γ−1k J2k ω and J2 k ω Γ−1k J2 k
ω using (4.19) and (4.23) by the big “O”
esti-mation. We then have Γ−1 k J 2k ω = 1 O 2m2 k O 2jkm 2 −1 j=(m−1)2 (4.24a) and Jω2kΓ−1k Jω2k= 1 O 2m2k O 2jkm 2 −1 j=(m−1)2. (4.24b) This implies (4.22).
Theorem 4.1. Let (M, L) be given in (4.1) and the sequence {(Mk,Lk)} be
gen-erated by Algorithm 2.1. Suppose that there are nonsingular matrices Z and Y such that (4.7) and (4.9) hold. Denote
Z = Z1 Z3 Z2 Z4 , (4.25)
where Zi∈ Cn×n for i = 1, . . . , 4. Suppose that Z1 is invertible, then
kMk(:, 1 : n) + Mk(:, n + 1 : 2n)Xk
(4.26)
≤ O̺(Js)2 k
+ O 2−k→ 0, as k → ∞, where X = Z2Z1−1 solves DARE (1.2).
Proof. Partition Mk,Lk conformally with (4.25) into
Mk = Mk,1 Mk,3 Mk,2 Mk,4 , Lk = Lk,1 Lk,3 Lk,2 Lk,4 . (4.27)
Substituting (4.25) and (4.27) into (4.15) we have Mk,1Z1+ Mk,3Z2 = Lk,1Z1 Js2k⊕ J2k ω + Lk,3Z2 Js2k⊕ J2k ω , (4.28a) (Mk,1Z3+ Mk,3Z4) ¯ Js2k⊕ Im = Lk,1 h Z1(0ℓ⊕ Γk) + Z3 Iℓ⊕ J2 k ω i +Lk,3 h Z2(0ℓ⊕ Γk) + Z4 Iℓ⊕ J2 k ω i , (4.28b) Mk,2Z1+ Mk,4Z2 = Lk,2Z1 Js2k⊕ Jω2k + Lk,4Z2 Js2k⊕ Jω2k , (4.28c) (Mk,2Z3+ Mk,4Z4) ¯ Js2k⊕ Im = Lk,2 h Z1(0ℓ⊕ Γk) + Z3 Iℓ⊕ J2 k ω i +Lk,4 h Z2(0ℓ⊕ Γk) + Z4 Iℓ⊕ J2 k ω i . (4.28d) Postmultiplying (4.28b) by0ℓ⊕ Γ−1k J2 k ω Z1−1to eliminate Lk,1Z1 0ℓ⊕ J2 k ω and Lk,3Z2 0ℓ⊕ J2 k ω in (4.28a) we get Mk,1+ Mk,3X = (Mk,1Z3+ Mk,3Z4) 0ℓ⊕ Γ−1k J 2k ω Z1−1 (4.29) − (Lk,1Z3+ Lk,3Z4) 0ℓ⊕ J2 k ω Γ−1k J 2k ω Z1−1 + (Lk,1Z1+ Lk,3Z2) Js2k⊕ 0m Z1−1. Similarly, from (4.28c) and (4.28d) follows that
Mk,2+ Mk,4X = (Mk,2Z3+ Mk,4Z4) 0ℓ⊕ Γ−1k J 2k ω Z1−1 (4.30) − (Lk,2Z3+ Lk,4Z4) 0ℓ⊕ J2 k ω Γ−1k J 2k ω Z1−1 + (Lk,2Z1+ Lk,4Z2) Js2k⊕ 0m Z1−1.
Since kMk,ik and kLk,ik are bounded, for i = 1, . . . , 4, from (4.29)–(4.30) and
Lemma 4.4 follows the assertion (4.26).
Theorem 4.2. Let (M, L) be given in (4.1) and the sequence {Ak, Gk, Hk} be
generated by Algorithm 2.2. Suppose Z and Y satisfy (4.7) and (4.9), respectively. Denote Z = Z1 Z3 Z2 Z4 , Y = Y1 Y3 Y2 Y4 , (4.31)
where Zi, Yi∈ Cn×n, for i = 1, . . . , 4. Suppose that Z1 is invertible and
kGkk ≤ O(2k), for sufficiently large k.
(4.32) Then
(a) if kGkk = O(2k) → ∞, as k → ∞, it holds
(i) kAkk ≤ O(1) bounded,
(ii) kX − Hkk ≤ O
̺(Js)2
k
+ O(2−k) → 0 as k → ∞,
where X = Z2Z1−1 solves the DARE (1.2);
(b) if kGkk ≤ O(1), for all k, and Y2 is invertible, it holds
(i) kAkk ≤ O(2−k) → 0, as k → ∞, (ii) kX − Hkk ≤ O ̺(Js)2 k + O(2−k) → 0, as k → ∞,
(iii) kY − Gkk ≤ O
̺(Js)2
k
+ O(2−k) → 0, as k → ∞,
where Y = −Y1Y2−1 solves the dual DARE
Y = AY (I + HY )−1AT + G.
Proof. . Substituting (Mk,Lk) of (4.10) and Z of (4.31) into (4.15), and comparing
both sides we obtain
AkZ1 = (Z1+ GkZ2) Js2k⊕ Jω2k, (4.33a) AkZ3 ¯ Js2k⊕ Im = (Z1+ GkZ2) (0ℓ⊕ Γk) (4.33b) +(Z3+ GkZ4) Iℓ⊕ J2 k ω , −HkZ1+ Z2 = ATkZ2 Js2k⊕ Jω2k, (4.33c) (−HkZ3+ Z4) ¯ Js2k⊕ Im = AT kZ2(0ℓ⊕ Γk) + ATkZ4 Iℓ⊕ J2 k ω . (4.33d)
Claim (a) (i) and (b) (i): Postmultiplying (4.33b) by 0ℓ⊕ Γ−1k J2 k ω
Z1−1 and using (4.33a) we have
Ak h I− Z3 0ℓ⊕ Γ−1k J 2k ω Z1−1i = (Z1+ GkZ2) Js2k⊕ 0m Z1−1 (4.34) − (Z3+ GkZ4) 0ℓ⊕ J2 k ω Γ−1k J 2k ω Z1−1. By Lemma 4.4 and assumptions in (a) or (b), the sequence {Ak} satisfies
kAkk ≤ O(2k) · O
̺(Js)2
k
+ O(2k) · O(2−k) = O(1)
(4.35)
for case (a); or
kAkk ≤ O(1) · O ̺(Js)2 k + O(1) · O(2−k) → 0, (4.36) as k → ∞, for case (b).
Claim (a) (ii) and (b) (ii): Postmultiplying (4.33d) by0ℓ⊕ Γ−1k J2 k ω Z1−1 and using (4.33c), we get −Hk h I− Z3 0ℓ⊕ Γ−1k J2 k ω Z1−1i+ X (4.37) = Z4 0ℓ⊕ Γ−1k J2 k ω Z1−1− AT kZ4 0ℓ⊕ J2 k ω Γ−1k J 2k ω Z1−1 +AT kZ2 Js2k⊕ 0m Z1−1.
By Lemma 4.4 and the boundedness of {Ak} in (4.35) or (4.36), the matrix X − Hk
in (4.37) can be bounded by kX − Hkk ≤ O ̺(Js)2 k + O(2−k),
for k sufficiently large.
Claim (b) (iii): Substituting (Mk,Lk) of (4.10) and Y of (4.31) into (4.13) we
have Y1+ GkY2 = AkY1 ¯ Js2k⊕ ¯Jω2k, (4.38a) (Y3+ GkY4) Js2k⊕ Im = AkY1 0ℓ⊕ ¯Γk+ AkY3 Iℓ⊕ ¯J2 k ω . (4.38b)
SDA [12] MDFM [7] NTM [23] Flops 23 3n 3 352 3 n 3 30n3
Table 1. Flop counts in each iteration.
As above, postmultiplying (4.38b) by0ℓ⊕ ¯Γ−1k J¯2 k ω
Y2−1and using (4.38a), we get −Y + Gk h I− Y4 0ℓ⊕ ¯Γ−1k J¯2 k ω Y2−1i (4.39) = Y3 0ℓ⊕ ¯Γ−1k J¯ 2k ω Y2−1− AkY3 0ℓ⊕ ¯J2 k ω Γ¯−1k J¯ 2k ω Y2−1 +AkY1 ¯ Js2k⊕ 0m Y2−1.
By Lemma 4.4 and boundedness of {Ak} we show that
kY − Gkk ≤ O
̺(Js)2
k
+ O(2−k),
for k sufficiently large.
Note that in Section 5 we show examples which satisfy either the assumption of kGkk in (a) or (b) of Theorem 4.2.
5. Numerical results
In this section, we test some numerical examples satisfying assumptions (A1) and (A2), respectively, on the SDA algorithms for DAREs and CAREs to illustrate the convergence behavior. All computations were performed in MATLAB R2006a on a PC with an Intel Pentium-IV 3.4 GHz processor and 2 GB main memory, using IEEE double-precision floating-point arithmetic (eps ≈ 2.22 × 10−16). The
operating system running the machine is Fedora Core 2 with Kernel 2.6.10-1.771-FC2.
The MATLAB commands “dare” and “care” fail in the initial step because one of the matrices G and H in (1.1) or (1.2) is indefinite. Furthermore, the strongly stable method [10] and matrix sign function methods [5, 8, 14, 17] fail to converge because the associated Hamiltonian matrix has purely imaginary eigen-values. Therefore, in our numerical testing, we compare SDA algorithms [11, 12] with MDFM methods [6, 7] and Newton’s methods [21, 22]. We summarize the flop counts for each iteration in SDAs, MDFMs and NTMs in Table 1.
For tables in the following examples, data for various methods are listed in columns with obvious headings. The heading SDA”, MDFM” and “D-NTM” stand for Algorithm 2.2, Algorithm 2.1 and Newton’s method [21] applied to DARE, respectively. “C-SDA” and “C-MDFM” stand for Algorithm 3.1 while call-ing Algorithm 2.2 and 2.1, respectively, in Step (III). “C-NTM” stands for Newton’s method [22] applied to CARE.
We report the numbers of iterations by “ITs”, the total flops (= Flops × ITs) by “TFs” for algorithms, and the maximal error between the accurate and the approximate stable eigenvalues of (I + G ˜X)−1A or A − G ˜X by “Err”, where ˜X
stable eigenvalue and ˜λi be the corresponding approximate eigenvalue. The “Err”
is defined by
Err = max
1≤i≤n{|λi− ˜λi|}.
The algorithm is terminated when the residual of DARE or CARE cannot be re-duced further.
5.1. Discrete-time algebraic Riccati equations. In this subsection, we shall report the numerical comparison of D-SDA, D-MDFM and D-NTM. The conver-gence of D-NTM is globally linear [21] when the size of Jω,p in (4.2) is two and
the matrix G is positive semidefinite. Furthermore, it needs to choose an initial matrix L0 so that A − BL0 is d-stable. If p > 2, the convergence of D-NTM can
not be guaranteed. In D-SDA and D-MDFM, G and H are only required to be symmetric and p is an even number. In the following, we shall show three examples to compare the performance of D-SDA, D-MDFM and D-NTM. The matrix G in the first example is negative definite. Thus, D-NTM can not be applied. The size of Jordan block Jω,p in the second example is two, so the linear convergence of
D-NTM is guaranteed. The other example with p > 2 the convergence of D-NTM is numerically linear.
For the residual of DARE, we use the “normalized” residual (DNRes) formula DNRes ≡ kA
TX(I + G ˜˜ X)−1A+ H − ˜Xk
k ˜Xk + kATX˜(I + G ˜X)−1Ak + kHk,
proposed in [7], where ˜X is an approximate solution to DARE. Example 5.1. [34, Example 2.2] Let
A= " 0 −3+√5 2 3−√5 2 0 # , G= −5 √ 5 − 1 2 !2 I2 and H = I2.
The symplectic pencil (M, L) has eigenvalues {ı, ı, −ı, −ı} and is equivalent to ı 1 0 ı ⊕ I2, I2⊕ ı 1 0 ı
,because rank(M ± ıL) = 3. On the other hand, we have M I 2 2 5−√5 I2 = L I 2 2 5−√5 I2 0 1 −1 0 and AT 0 −G I I 2 √ 5−5 2 I2 = I H 0 A I 2 √ 5−5 2 I2 0 −1 1 0 .
These imply that the unique Hermitian solutions of the DARE and the dual DARE are X∗= 2 5 −√5 I2, Y∗ = √ 5 − 5 2 ! I2, (5.1)
respectively. It holds that I + Y∗X∗= 0.
The numerical results by using D-SDA and D-MDFM for computing ˜X are reported in Table 2. The notation “∗” in the fourth column of Table 2 means that D-NTM cannot be applied, since G is not positive semidefinite. We see that the
D-SDA D-MDFM D-NTM DNRes 5.55 × 10−17 5.68 × 10−12 ∗
ITs 25 18 ∗
TFs(n3) 192 2112 ∗
Err 2.48 × 10−8 5.72 × 10−6 ∗
Table 2. Results for Example 5.1.
normalized residual, i.e., the backward error, for ˜X by D-SDA has 16 significant digits and attains the machine accuracy. It has only 12 significant digits by D-MDFM.
Compared to the DNRes which has the machine accuracy, the forward errors of stable eigenvalues have only eight and six significant digits for SDA and D-MDFM, respectively. This is due to the poor separation between the d-stable and d-unstable subspectra of (M, L) [45, 46]. In fact, one can easily check that the symplectic pencil 0 0 −H∞ I − λ I G∞ 0 0
is singular if and only if I+G∞H∞is singular. Then from the fact that I+Y∗X∗= 0, it follows that the pencil Mk− λLk generated by the D-SDA algorithm tends to a
singular pencil as k → ∞.
Example 5.2. [21, Example 6.2] Let M0= A0 0 0 I8 and L0= I8 G0 0 AT 0 with A0= −1 1 1 ⊕ " √ 3 2 12 −12 √ 3 2 # ⊕ 1 2 1 1 2 1 1 2 and G0= 1 1 1 . .. ... 1 1 1 1 1 . .. ... 1 1 T ∈ R8×8.
The symplectic pencil (M0,L0) has eigenvalues {−1, 1, √ 3−ı 2 , √ 3+ı 2 ,2, 1 2} with
par-tial multiplicities {2, 2, 2, 2, 3, 3}. For the corresponding DARE, all of the matrices Hk generated by D-SDA are equal to zero matrix which is the Hermitian solution
of DARE. That is, D-SDA only needs one iteration to obtain exact solution. There-fore, for such symplectic pencil (M0,L0), the performance of D-SDA is obviously
better than D-MDFM and D-NTM.
In order to get a fair comparison of these three methods, we take an equivalent transformation for (M0,L0). Set the “state” in the MATLAB command “randn”
to be [3616447672, 4257789242]T and the matrix H
D-SDA D-MDFM D-NTM DNRes 4.54 × 10−15 1.39 × 10−12 1.08 × 10−16
ITs 24 20 20
TFs(n3) 184 2347 600
Err 7.14 × 10−6 1.17 × 10−5 7.23 × 10−6
Table 3. Results for Example 5.2.
normally distributed random diagonal elements so that I − G0H0 is nonsingular.
Define a new equivalent symplectic pencil (M, L) as M = A 0 −H I ≡ (I − G0H0)−1 0 0 I I 0 AT 0H0(I − G0H0)−1 I × A0 0 0 I I 0 −H0 I (5.2a) and L = I G 0 AT ≡ (I − G0H0)−1 0 0 I I 0 AT 0H0(I − G0H0)−1 I × I G0 0 AT 0 I 0 −H0 I . (5.2b)
Note that G and H in (5.2) are positive definite and indefinite, respectively. The matrix L0 in D-NTM is taken as a normally distributed random diagonal
matrix with state = [3080845898, 338958098]T. The numerical results by using
D-SDA, D-MDFM and D-NTM for computing ˜X of DARE with matrices A, G and H in (5.2) are reported in Table 3.
Example 5.3. Let M0= A0 0 0 I7 and L0= I7 G0 0 AT 0 with A0= −1 −1 2 −1 ⊕ −1 −1 2 − 1 8 −1 −12 −1 ⊕ −13 I2 and G0= 1 32 1 8 1 8 1 2 ⊕ 1 512 1 128 1 32 1 128 1 32 1 8 1 32 1 8 1 2 ⊕ 2 9 I2 .
The symplectic pencil (M0,L0) has eigenvalues {−3, −3,13,13,−1, −1} and
eigen-value -1 has partial multiplicities {4, 6}.
Taking H0= −3.5I7and using the equivalent transformation as in (5.2), we get
a new DARE with G and H being positive semidefinite and indefinite, respectively. The numerical results by using D-SDA, D-MDFM and D-NTM for computing ˜X of
D-SDA D-MDFM D-NTM DNRes 3.76 × 10−13 2.03 × 10−10 1.52 × 10−10
ITs 18 18 31
TFs(n3) 138 2112 930
Err 1.80 × 10−3 1.49 × 10−3 1.78 × 10−3
Table 4. Results for Example 5.3.
this DARE are reported in Table 4. Here, the matrix L0 in D-NTM is a normally
distributed random diagonal matrix with state = [2355717396, 3700125409]T.
From Table 4, we see that the backward error for ˜X by D-SDA has 13 significant digits which is better than that by D-MDFM and D-NTM. The forward errors of stable eigenvalues by using these three methods are almost equal to √6eps.
5.2. Continuous-time algebraic Riccati equations. In this subsection, we shall report the numerical comparison of C-SDA, C-MDFM and C-NTM by three examples. Note that, similar to D-NTM, C-NTM converges globally and linearly when the size of Jω,p in (4.2) is two. If p > 2, the convergence of C-NTM cannot
be guaranteed. Consequently, we give three examples to illustrate the numerical behavior. In Example 5.4, C-NTM is guaranteed to converge linearly, where p = 2. In Example 5.5, where the maximal size of Jω,pis four, the convergence of C-NTM
is numerically linear. However, in Example 5.6, where the maximal size of Jω,p
is eight, the sequence generated by C-NTM is divergent. Note that in C-NTM it requires to choose an initial matrix X0 so that A − GX0is stable.
For the residuals of CARE, we use the “normalized” residual (CNRes) formula CNRes ≡ k − ˜XG ˜X+ A
TX˜+ ˜XA+ Hk
k ˜XG ˜Xk + kATX˜k + k ˜XAk + kHk,
proposed in [11], where ˜X is an approximate solution to CARE. Example 5.4. [22, Example 4.3] Let
H0= A0 −G0 0 −AT 0 where A0= 02⊕ 0 1 −1 0 ⊕ 0 2 −2 0 ⊕ −1 1 0 −1 and G0= 2 1 1 1 2 . .. . .. ... 1 1 1 2 ∈ R 8×8.
The eigenvalues of H0 are {1, −1, 0, 0, ı, −ı, 2ı, −2ı} and the partial multiplicity of
each eigenvalue is two. For this CARE, each of Hkin C-SDA is a zero matrix which
is the Hermitian solution of CARE. Thus, the exact solution of CARE is obtained after one iteration by C-SDA.
C-SDA C-MDFM C-NTM CNRes 3.78 × 10−16 5.39 × 10−14 2.08 × 10−16
ITs 26 23 30
TFs(n3) 199 2699 900
Err 2.27 × 10−8 3.58 × 10−7 3.75 × 10−8
Table 5. Results for Example 5.4.
Now, we use the following similarity transformation to get a new Hamiltonian matrix H: H = A −G −H −AT ≡ " 1 2I √ 3 2 I −√3 2 I 1 2I # A0 −G0 0 −AT 0 " 1 2I − √ 3 2 I √ 3 2 I 1 2I # . The numerical results by using C-SDA, C-MDFM and C-NTM for computing ˜X of the new CARE are reported in Table 5. Here, the initial matrix X0in C-NTM is a
normally distributed random diagonal matrix with state = [549776141, 3491173770]T
and we take γ = 2.2 in the Cayley transformation for C-SDA and C-MDFM. From Table 5, it shows that the normalized residuals for ˜X by C-SDA and C-NTM have 16 significant digits which are better than that by C-MDFM. The accuracy for stable eigenvalues by C-SDA and C-NTM are also better than that by C-MDFM.
Example 5.5. Let A0 and G0be 5 × 5 real matrices defined by
A0= U I2 0 0 U 0 0 0 a , G0= 00 I02 00 0 0 g
where a, g are parameters and U = 0 β −β 0 , β6= 0. The matrix A0 −G0 0 −AT 0
has nonzero eigenvalues a, −a and pure imaginary eigenvalues βı, −βı with partial multiplicity 4. Define A −G −H −AT = PT 0 0 PT C S −S C A0 −G0 0 −AT 0 C −S S C P 0 0 P , where P = I5− 2uuT is a Householder matrix and
C= I4 0 0 c , S= 04 0 0 s
with c > 0, s > 0 and c2+ s2= 1. Taking
a = 0.6781616521431886, β= 5.513985806849778, g = 43.14437852853182, c= 0.3559455724227920, u = 0.8210373788415466 0.4365444378807448 0.1651283156254210 0.3286639287858224 0.6263991871530291 × 10−2 ,
C-SDA C-MDFM C-NTM CNRes 1.00 × 10−16 3.25 × 10−9 8.40 × 10−13
ITs 9 17 34
TFs(n3) 69 1995 1020
Err 3.27 × 10−8 5.08 × 10−4 2.45 × 10−5
Table 6. Results for Example 5.5.
we have that G and H are positive and negative semidefinite, respectively.
Taking γ = 35 in the Cayley transformation and the initial X0 in C-NTM to be
a normally distributed random matrix with state = [2406827087, 4217562603]T,
numerical results by using C-SDA, C-MDFM and C-NTM for computing ˜X are reported in Table 6. It shows that the CNRes for ˜Xby C-SDA attains the machine accuracy. C-MDFM and C-NTM have only 9 and 13 significant digits, respectively. Furthermore, the accuracy of stable eigenvalues by C-SDA is better than that by C-MDFM and C-NTM. Example 5.6. [33] Let A0 = 0 ⊕ 0 0 1 0 ⊕ 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 ⊕ (−I2), H0 = (−1) ⊕ 0 0 0 −1 ⊕ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 −1 ⊕ (−I2).
Construct a Hamiltonian matrix H in exact arithmetic: H = I V2 0 I V1T 0 0 V1−1 A0 0 H0 −AT0 V1−T 0 0 V1 I −V2 0 I , where V1= 1 1 0 . .. ... . .. 1 0 1 , V2= 1 1 0 1 −1 2 2 1 3 . .. ... ... 7 −1 8 0 8 1 .
It is easily seen that H has nonzero eigenvalues −1, −1, 1, 1 and the zero eigenvalue has partial multiplicities {2, 4, 8}.
Taking γ = 10.5 in the Cayley transformation and the initial X0 in C-NTM to
be a normally distributed random matrix with state = [2902100237, 857523738]T,
numerical results by using C-SDA, C-MDFM and C-NTM are reported in Table 7. It shows that the normalized residual for ˜X by C-SDA and C-MDFM has 13 and 9 significant digits, respectively, while the C-NTM is divergent. The CNRes by C-SDA lost 3 significant digits compared to the machine accuracy, because the
C-SDA C-MDFM C-NTM CNRes 7.51 × 10−13 6.07 × 10−9 ∗
ITs 20 17 ∗
TFs(n3) 153 1995 ∗
Err 6.56 × 10−2 5.55 × 10−2 ∗
Table 7. Results for Example 5.6.
highest partial multiplicity of the zero eigenvalue of H is eight, making it very sensitive to perturbation [45, 46]. The forward errors of stable eigenvalues by using C-SDA and C-MDFM are attained to √8eps.
5.3. Comments.
(i) For SDAs, MDFMs and NTMs, it requires23 3n
3[12], 352 3 n
3[7] and 30n3[23]
flops, respectively, in each iterative step. The flops count of SDA is about 7% and 26% of that of MDFM and NTM, respectively. This is mainly due to the fact that the main step in MDFM involves the QR-factorization of
LT,−MTT and the formation of Q ∈ C4n×4n, all in higher dimensions.
The operations in SDA are all in Cn×n while keeping the SSF form.
(ii) Examples investigated here are all ill-conditioned in solving DARE or CARE because the associated symplectic pencils or Hamiltonian matrices have eigenvalues on the unit circle or on the purely imaginary axis. However, the SDA algorithms solve them efficiently and accurately without failure and with less flops count.
(iii) Tables 2, 5 and 6 show that the approximate solutions ˜X from SDA are as accurate as the machine accuracy, and more accurate than those from MDFM. This behavior illustrates the importance of the SSF property, and the consequent global and linear convergence.
(iv) For Examples 5.1, 5.2, 5.3, 5.4 and 5.6, the convergence behavior of Ak,
Gk and Hk coincides with the result of Theorem 4.2 (b) which converges
linearly and the matrix I + GkHktends to a singular matrix as k → ∞. We
use Example 5.1 to illustrate such convergence behavior. Let X∗and Y∗ be defined in (5.1), then the Frobenius norms of Ak, Gk− Y∗ and Hk− X∗in
each iteration of D-SDA are shown in Figure 1 (a), (b) and (c), respectively. (v) For Example 5.5, the Frobenius norms of Ak, Gk and Hk+1− Hk in each
iteration of C-SDA are shown in Figure 2 (a), (b) and (c), respectively. In Figure 2, the convergence behavior of Ak, Gk and Hk coincides with the
results of Theorem 4.2 (a). It shows that kA9kF ≈ 20.3, kG9kF ≈ 4.82×103
and kH9− H8kF ≈ 5.04 × 10−10. Furthermore, according to Figure 3 (a)
and (b), HkAkconverges to zero matrix as k → ∞ and the matrix I +GkHk
is well-conditioned in each iteration so that Hk converges to the solution of
CARE. At the end of the algorithm, we have kH9A9kF ≈ 5.04 × 10−10and
cond(I + H9G9) ≈ 11.5 with kI + H9G9kF ≈ 2.
(vi) In our numerical experience, SDAs do not break down. That is, the matrix I+GkHk is nonsingular in each iteration of SDAs. In fact, when the matrix
I+ GkHk is close to be ill-conditioned, the sequence {Hk} has converged
0 10 20 10−8 10−7 10−6 10−5 10−4 10−3 10−2 10−1 k || A k || F (a) 0 10 20 10−8 10−7 10−6 10−5 10−4 10−3 10−2 10−1 k || G k − Y * || F (b) 0 10 20 10−8 10−7 10−6 10−5 10−4 10−3 10−2 10−1 k || H k − X * || F (c)
Figure 1. The Frobenius norms of Ak, Gk− Y∗ and Hk− X∗for
Example 5.1 by D-SDA.
6. Concluding remarks
In this paper, we propose the structured doubling algorithms for finding weak Hermitian solutions of DAREs and CAREs. Under assumptions (A1) and (A2), we prove the global and linear convergence for SDA algorithms. The advantage of structured doubling algorithms is evident that the Hermitian solutions are ob-tained by the iterative process without any deflation preprocessing of unimodular eigenvalues. The MATLAB commands “care” and “dare” fail for the selected test examples, because the initial matrices G0and H0 are not both symmetric positive
semidefinite. Nevertheless the normalized residuals of Hermitian solutions of almost all tested examples computed by SDA algorithms are accurate to the machine ac-curacy. Numerical experiments show that SDA algorithms converge to the desired Hermitian solutions efficiently and reliably.
Appendix: Proof of Theorem 1.2 With α defined in (1.9), we have (I + GA−H
α HA−1α )−1 being positive definite,
and therefore (I + GA−H
α HA−1α )−1Gis symmetric and positive semidefinite. This
implies that M − αL, where M, L are given in (1.5), is invertible and can be factorized as M − αL = −αAα −αG −H AHα = I −αGA−H α 0 I −αKα 0 −H AHα (A.1) with Kα≡ Aα+ GA−Hα H.
0 5 10 2 4 6 8 10 12 14 16 18 20 22 k || A k || F (a) 0 5 10 10−1 100 101 102 103 104 k || G k || F (b) 0 5 10 10−10 10−8 10−6 10−4 10−2 100 k || H k+1 − H k || F (c)
Figure 2. The Frobenius norms of Ak, Gk and Hk+1− Hk for
Example 5.5 by C-SDA.
Applying Cayley transformation to the symplectic pair (M, L) with parameter α, we get the Hamiltonian matrix
b H = (M − αL)−1(M + αL) = −I + 2(M − αL)−1M ≡ " b A − bG − bH − bAH # . (A.2) With b G= 2A−1 α I+ GA−Hα HA−1α −1 GA−Hα (A.3) and b A= −αA¯ + I + bGHA−1α . (A.4)
Let bG = bB bBH ≥ 0 be a full rank decomposition. Since the Lagrangian subspace
is invariant via Cayley transformation, by applying the result of Theorem 1.1, it suffices to show that ( bA, bB) is c-stable. Let bw6= 0 ∈ Cn satisfy
b
wHBb= 0, wbHAb= µ bwH. (A.5)
0 5 10 10−10 10−8 10−6 10−4 10−2 100 k || H k A k || F (a) 0 5 10 0 2 4 6 8 10 12 k condition number of I + H k G k (b)
Figure 3. The Frobenius norms of HkAk and condition numbers
of I + HkGk for Example 5.5 by C-SDA.
From (A.5), we have b
wHAb = − bwHαA¯ + I + bGHA−1 α
= − bwH(¯αA+ I)(I − ¯αA)−1
= µ bwH. (A.6) This implies wHA= −¯α 1 + µ 1 − µ wH, (A.7)
where wH = bwH(I − ¯αA)−1. From (A.3) and (A.5), it holds
b wHGb = wbHA−1α G I+ A−Hα HA−1α G −1 A−Hα = wHG I+ A−Hα HA−1α G −1 A−Hα = 0. (A.8)
From (A.8), it follows that
wHG= wHB = 0.
(A.9)
In view of the d-stabilizability of the pair (A, B), we have 1 + µ 1 − µ < 1 (A.10)
Acknowledgements
The authors thank Prof. E. K.-W. Chu from Monash University for many valu-able comments and helpful suggestions.
References
1. G. Ammar and V. Mehrmann, On Hamiltonian and symplectic Hessenberg forms, Lin. Alg. Appl. 149 (1991), 55–72.
2. B. D. O. Anderson and S. Vongpanitlerd, Network analysis and synthesis, Lecture Notes in Control and Information Science, Prentice-Hall, Englewood Cliffs, NJ, 1973.
3. B.D.O. Anderson, Second-order convergent algorithms for the steady-state Riccati equation, Int. J. of Control 28 (1978), no. 2, 295–306.
4. Z. Bai, J. Demmel, and M. Gu, An inverse free parallel spectral divide and conquer algorithm for nonsymmetric eigenproblems, Numerische Mathematik 76 (1997), no. 3, 279–308. 5. L. Balzer, Accelerated convergence of the matrix sign function, Int. J. Control 21 (1980), pp.
1057–1078.
6. P. Benner, Contributions to the numerical solutions of algebraic Riccati equations and re-lated eigenvalue problems, Ph.D. thesis, Fakult¨at f¨ur Mathematik, TU Chemnitz-Zwickau, Chemnitz, Germany, 1997.
7. P. Benner and R. Byers, Evaluating products of matrix pencils and collapsing matrix products, Num. Lin. Alg. Appl. 8 (2001), 357–380.
8. R. Byers, Solving the algebraic Riccati equation with the matrix sign function, Lin. Alg. Appl. 85 (1987), pp. 267–279.
9. J. R. Cardoso and F. Silva Leite, The Moser-Veselov equation, Lin. Alg. Appl. 360 (2003), 237–248.
10. D. Chu, X. Liu, and V. Mehrmann, A numerically backwards stable method for computing the Hamiltonian Schur form, Tech. report, Preprint 24-2004 TU Berlin, Institut f¨ur Mathematik, 2004.
11. E. K.-W. Chu, H.-Y. Fan, and W.-W. Lin, A structure-preserving doubling algorithm for continuous-time algebraic Riccati equations, Linear Algebra Appl. 396 (2005), 55–80. 12. E. K.-W. Chu, H.-Y. Fan, W.-W. Lin, and C.-S. Wang, Structure-preserving algorithms for
periodic discrete-time algebraic Riccati equations, Int. J. Control 77 (2004), 767–788. 13. D. Clements and K. Glover, Spectral transformations via Hermitian pencils, Linear Algebra
Appl. 123 (1989), 797–846.
14. E. Denman and R. Beavers, The matrix sign function and computations in systems, Appl. Math. Comp. 2 (1976), pp. 63–94.
15. B. A. Francis and J. C. Doyle, Linear control theory with an H∞
-optimality criterion, SIAM J. Control Optim. 25 (1987), 815–844.
16. F. R. Gantmacher, The theory of matrices, vol. ii, Chelsea Publishing Company, New York, N. Y., 1977.
17. J. Gardiner and A. J. Laub, A generalization of the matrix-sign-function solution to the algebraic riccati equations, Int. J. of Control 44 (1986), 823–832.
18. K. Glover, D. J. N. Limebeer, J. C. Doyle, E. M. Kasenally, and M. G. Safonov, A charac-terization of all solutions to the fourblock general distance problem, SIAM J. Control Optim. 39 (1991), 283–324.
19. I Gohberg, P Lancaster, and L. Rodman, On the Hermitian solution of the symmetric alge-braic Riccati equation, SIAM J. Control and Optimization 24 (1986), 1323–1334.
20. G. H. Golub and C. F. Van Loan, Matrix computations, 3rd ed., The Johns Hopkins University Press, 1996.
21. C.-H. Guo, Newton’s method for discrete algebraic Riccati equations with the closed-loop matrix has eigenvalues on the unit circle, SIAM J. Matrix Anal. Appl. 20 (1998), 279–294. 22. C.-H. Guo and P. Lancaster, Analysis and modificaton of Newton’s method for algebraic
Riccati equations, Math. Comp. 67 (1998), 1089–1105.
23. X.-X. Guo, W.-W. Lin, and S.-F. Xu, A structure-preserving doubling algorithm for nonsym-metric algebraic Riccati equation, accepted for Num. Math. 2006.
24. T.-M. Hwang, E. K.-W. Chu, and W.-W. Lin, A generalized structure-preserving doubling algorithm for generalized discrete-time algebraic Riccati equations, Int. J. Control 78 (2005), 1063–1075.
25. V. Ionescu and M. Weiss, Two Riccati formulae for the discrete-time H∞
-control problem, Internat. J. Control 57 (1993), 141–195.
26. M. Kimura, Convergence of the doubling algorithm for the discrete-time algebraic Riccati equation, Int. J. Syst. Sci. 19 (1988), no. 5, 701–711.
27. , Doubling algorithm for continuous-time algebraic Riccati equation, Int. J. Syst. Sci. 20 (1989), no. 2, 191–202.
28. P. Lancaster, A. C. M. Ran, and L. Rodman, Hermitian solutions of the discrete algebraic Riccati equation, Int. J. Control 44 (1986), 777–802.
29. P. Lancaster and L. Rodman, Existence and uniqueness theorems for the algebraic Riccati equation, Int. J. Control 32 (1980), no. 2, 285–309.
30. A. J. Laub, A Schur method for solving algebraic Riccati equations, IEEE Transactions on Automatic Control 24 (1979), no. 6, 913–921.
31. W.-W. Lin, A new method for computing the closed-loop eigenvalues of a discrete-time alge-braic Riccatic equation, Lin. Alg. Appl. 96 (1987), 157–180.
32. W.-W. Lin, V. Mehrmann, and H. Xu, Canonical forms for Hamiltonian and Symplectic matrices and pencils, Lin. Alg. Appl. 302-303 (1999), 469–533.
33. W.-W. Lin and C.-S. Wang, On computing stable Largangian subspaces of Hamiltonian ma-trices and symplectic pencils, SIAM J. Matrix Anal. Appl. 18 (1997), no. 3, 590–614. 34. W.-W. Lin, C.-S. Wang, and Q.-F. Xu, Numerical computation of the minimum H∞norm of
the discrete-time output feedback control problem, SIAM J. Numer. Anal. 38 (2000), 515–547. 35. W.-W. Lin and S.-F. Xu, Convergence analysis of structure-preserving doubling algorithms
for Riccati-tye matrix equations, SIAM Matrix Anal. Appl. 28(1) (2006), 26–39.
36. A. N. Malyshev, Parallel algorithm for solving some spectral problems of linear algebra, Lin. Alg. Appl. 188/189 (1993), 489–520.
37. MathWorks, MATLAB user’s guide (for UNIX workstations), The Math Works, Inc., 2002. 38. V. Mehrmann, The autonomous linear quadratic control problem, Springer-Verlag, 1991. 39. , A step toward a unified treatment of continuous and discrete time control problems,
Lin. Alg. Appl. 241-243 (1996), 749–779.
40. J. Moser and A. P. Veselov, Discrete versions of some classical integrable systems and fac-torization of matrix polynomials, Commun. Math. Phys. 139 (1991), 217–243.
41. C. Paige and C. Van Loan, A Schur decomposition for Hamiltonian matrices, Linear Algebra and its Applications 41 (1981), 11–32.
42. T. Pappas, A. J. Laub, and N. R. Sandell, On the numerical solution of the discrete-time algebraic Riccati equation, IEEE Transactions on Automatic Control 25 (1980), 631–641. 43. G. W. Stewart, HQR3 and EXCHNG: Fortran subroutines for calulating and ordering the
eigenvalues of a real upper Hessenberg matrix, ACM Trans. Math. Software 2 (1976), pp. 275–280.
44. G. W. Stewart and J.-G. Sun, Matrix perturbation theory, Academic Press, New York, 1990. 45. J.-G. Sun, Perturbation theory for algebraic Riccati equations, SIAM J. Matrix Anal. Appl.
19 (1998), 39–65.
46. , Condition numbers of algebraic Riccati equations in the Frobenius norm, Linear Algebra Appl. 350 (2002), 237–261.
47. A. J. Van Der Schaft and J. C. Willems, A new procedure for stochastic realization of spectral density matrices, SIAM J. Control Optim. 22 (1984), 845–855.
48. P. Van Dooren, A generalized eigenvalue approach for solving riccati equations, SIAM J. Sci. Statist. Comput. 2 (1981), no. 2, 121–135.
Department of Mathematics, National Taiwan Normal University, Taipei 11677, Tai-wan.
E-mail address: [email protected]
Department of Mathematics, National Tsing Hua University, Hsinchu 30043, Taiwan. E-mail address: [email protected]
