Convergence Analysis of Structure-Preserving Doubling
Algorithms for Riccati-Type Matrix Equations
Wen-Wei Lin
∗Shu-Fang Xu
†Abstract
In this paper, a structure-preserving transformation of a symplectic pencil is introduced, referred to as the doubling transformation, and its some basic properties are presented. Based on the nice properties of this kind of transformations, a unified convergence theory for the structure-preserving doubling algorithms for solving a class of Riccati-type matrix equations is established by using only the knowledge from elementary matrix theory.
Key words.matrix equation, structure-preserving doubling algorithm, convergence rate AMS subject classfications.15A24, 65H10, 93B50, 93D15
1
Introduction
In this paper we investigate the convergence of the structure-preserving doubling algorithms (SDA) for the computation of the symmetric positive definite or semi-definite solutions to the following Riccati-type matrix equations:
• Continuous-time algebraic Riccati equation (CARE) [20, 25]
−XGX + ATX + XA + H = 0, (1.1)
where A, H, G ∈ Rn×n with G and H symmetric positive semi-definite.
• Discrete-time algebraic Riccati equation (DARE) [20, 25]
X = ATX(I + GX)−1A + H, (1.2)
where A, H, G ∈ Rn×n with G and H symmetric positive semi-definite.
• Nonlinear matrix equation with the plus sign (NME-P) [3]
X + ATX−1A = Q, (1.3)
where A, Q ∈ Rn×n with Q symmetric positive definite.
• Nonlinear matrix equation with the minus sign (NME-M) [11]
X − ATX−1A = Q, (1.4)
where A, Q ∈ Rn×n with Q symmetric positive definite.
∗Department of Mathematics, National Tsing Hua University, Hsinchu, 300, Taiwan ([email protected]). †School of Mathematical Sciences, Peking University, Beijing, 100871, China ([email protected]). This research
The Riccati-type matrix equations occur in many important applications (see [3, 11, 20, 25] and references therein). The nonlinear matrix equations CARE and DARE have been studied extensively (see [1]–[2], [4]–[8], [13]–[14], [17]–[25], [27]–[29], and [31]); and the nonlinear matrix equations NME-P and NME-M has been studied recently by several authors (see [3], [9]–[11], [15]–[16], [26], [30], and [32]).
A class of methods, referred to as the doubling algorithms, has attracted much interests in 70s and 80s (see [2] and references therein). These methods originate from the fixed-point iteration derived from the DARE:
Xk+1= ATXk(I + GXk)−1A + H.
Instead of generating the sequence {Xk}, doubling algorithms generate {X2k}. Doubling algorithms
were largely forgotten in the past decade. Recently, they have been revived for DAREs and CAREs, because their nice numerical behavior: quadratical convergence rate, low cost computational cost per step, and good numerical stability (see [6, 7, 8]). Concerning the matrix equations NME-Ps and NME-Ms, in 2001 B. Meini proposed an iterative method with the same numerical behavior as the SDA algorithms for solving the DAREs and CAREs (see [26, 16]).
In this paper, by employing the same technique as in [8], we derive two SDA algorithms for solving the NME-Ps and NME-Ms, as a result we find they are essentially the same as Meini’s in [26]. In summary, we find all the SDA algorithms for solving the four Riccati-type matrix equa-tions can be viewed as applying some special structure-preserving transformaequa-tions to the associated symplectic pencils repeatedly. Therefore, we first introduce a structure-preserving transformation of a symplectic pencil, which is referred to as the doubling transformation, and then prove that an important feature of this kind of transformations is that it is structure preserving, eigenspace preserving, and eigenvalue doubling. Finally, based on the nice properties of the doubling trans-formations we develop a unified convergence theory for these SDA algorithms by using only the knowledge from elementary matrix theory.
Throughout this pater, the symbols k · k2 denote the matrix spectral norm. For a given n × n
matrix A, we use ρ(A) to denote the spectral radius of A. For real symmetric matrices X and Y , we write X > Y (X ≥ Y ) if X − Y is symmetric positive definite (semi-definite).
The paper is organized as follows. In Section 2, we introduce a structure-preserving transfor-mation of a symplectic pencil, and develop its basic properties. In Section 3, we do the convergence analysis of the SDA algorithm for solving the DAREs based on the theory established in Section 2. Section 4 derives a SDA algorithm for solving the NME-P (1.3) by using the doubling transfor-mations, and establishes the convergence theory of the alorithm. Concluding remarks are given in Section 5.
2
Doubling Transformation
In this section we introduce a structure-preserving transformation of a symplectic pencil and de-velop its some basic properties. We begin with the definition of the transformation.
Let M − λL ∈ R2n×2nbe a symplectic pencil, i.e.,
M JMT = LJLT, J = · 0 I −I 0 ¸ . (2.1) Define N (M, L) =n[M∗, L∗] : M∗, L∗∈ R 2n×2n, rank[M ∗, L∗] = 2n, [M∗, L∗] · L −M ¸ = 0o. (2.2) Since rank · L −M ¸
≤ 2n it follows that N (M, L) 6= ∅. For any given [M∗, L∗] ∈ N (M, L), define
c
M = M∗M, L = Lb ∗L. (2.3)
The transformation
is called a doubling transformation.
An important feature of this kind of transformations is that it is structure preserving, eigenspace preserving, and eigenvalue doubling, as shown in the following theorem.
Theorem 2.1. Assume that the pencil cM − λbL is a doubling transformation of a symplectic pencil M − λL. Then we have:
(a) The pencil cM − λbL is still symplectic. (b) If M · U V ¸ = L · U V ¸
S, where U, V ∈ Rn×m and S ∈ Rm×m, then cM
· U V ¸ = bL · U V ¸ S2.
(c) If the pencil M − λL has the Kronecker canonical form W M Z = · Jr 0 0 I2n−r ¸ , W LZ = · Ir 0 0 N2n−r ¸ , (2.4)
where W, Z are nonsingular, Jra Jordan matrix corresponding to the finite eigenvalues of M − λL,
and N2n−r a nilpotent Jordan matrix corresponding to the infinite eigenvalues of M − λL, then
there exists a nonsingular matrix cW such that c W cM Z = · J2 r 0 0 I2n−r ¸ , W bcLZ = · Ir 0 0 N2 2n−r ¸ . (2.5)
Proof. (a) Since [M∗, L∗] ∈ N (M, L) implies that M∗L = L∗M , it follows from (2.1) that
c
M J cMT = M
∗M JMTM∗T = M∗LJLTM∗T = (L∗M )J(MTL∗T) = L∗LJLTL∗T = bLJ bL T,
which means that the pencil cM − λbL is symplectic. (b) Using the equality M∗L = L∗M again, we have
M∗L · U V ¸ S = L∗M · U V ¸ S, and hence, c M · U V ¸ = M∗M · U V ¸ = M∗L · U V ¸ S = L∗M · U V ¸ S = L∗L · U V ¸ S2= bL · U V ¸ S2. (c) Let ˇ M∗= · Jr 0 0 I2n−r ¸ W, Lˇ∗= · Ir 0 0 N2n−r ¸ W. Then we have ˇ M∗LZ = · Jr 0 0 I2n−r ¸ W LZ = · Jr 0 0 I2n−r ¸ · Ir 0 0 N2n−r ¸ = · Jr 0 0 N2n−r ¸ , ˇ L∗M Z = · Ir 0 0 N2n−r ¸ W M Z = · Ir 0 0 N2n−r ¸ · Jr 0 0 I2n−r ¸ = · Jr 0 0 N2n−r ¸ .
This, together with Z nonsingular, implies that [ ˇM∗, ˇL∗]
· L −M
¸
= 0. It is clear that rank[ ˇM∗, ˇL∗] =
2n. Notice that (2.4) implies that the pencil M − λL is regular, it must have rank ·
M −L
¸ = 2n. Thus, the rows of [M∗, L∗] and [ ˇM∗, ˇL∗] form two different bases of the left null space of
· M −L
¸ . Consequently, there exists a nonsingular matrix cW such that cW [M∗, L∗] = [ ˇM∗, ˇL∗]. It can be
Remark 2.1.A subspace W of R2n is called a generalized eigenspace of a pencil M − λL ∈ R2n×2n
if W is spanned by the columns of W = ·
U V ¸
, where U, V ∈ Rn×m, W has full column rank and
satisfies M W = LW S with S ∈ Rm×m. Therefore, (b) of Theorem 2.1 tells us that if W is a
generalized eigenspace of a symplectic pencil M − λL, then it is also a generalized eigenspace of its doubling transformation, which is a corner stone for us to establish convergence theory of the SDA algorithms for solving the Riccati-type matrix equations in the next two sections.
Remark 2.2.A pencil M − λL is called regular if det(M − λL) does not vanish identically. It is well-known that a pencil is regular if and only if it has a Kronecker canonical form. Thus, (c) of Theorem 2.1 implies that if M − λL is regular, then its doubling transformation is also regular, and that λ is a eigenvalue of M − λL if and only if λ2 is an eigenvalue of its doubling transformation.
A symplectic pencil M − λL is said to be a first standard symplectic form (SSF-1) if it has the following form M = · A 0 −H I ¸ , L = · I G 0 AT ¸ , (2.6)
with H, G ≥ 0; A symplectic pencil M −λL is said to be a second standard symplectic form (SSF-2) if it has the following form
M = · A 0 Q −I ¸ , L = · −P I AT 0 ¸ , (2.7) with P, Q ≥ 0.
The following theorem shows that the two standard symplectic forms are preserved by an appropriate choice of the doubling transformations.
Theorem 2.2. (a) Let M − λL be a SSF-1 form. Then a matrix [M∗, L∗] ∈ N (M, L) can be
constructed such that its corresponding doubling transformation cM − λbL is still a SSF-1.
(b) Let M − λL be a SSF-2 form. If Q − P > 0 and Q − AT(Q − P )−1A ≥ 0, then a matrix
[M∗, L∗] ∈ N (M, L) can be constructed such that its corresponding doubling transformation cM −λbL
is still a SSF-2.
Proof. (a) Applying block Gaussian elimination and row permutation to the matrix · L −M ¸ , we can get M∗= · A(I + GH)−1 0 −AT(I + HG)−1H I ¸ , L∗= · I AG(I + HG)−1 0 AT(I + HG)−1 ¸ (2.8) such that M∗L = L∗M, (2.9)
i.e., [M∗, L∗] ∈ N (M, L). Here the Sherman-Morrison-Woodbury formula (see, e.g., [12, p. 50]) is
used. For more details see [8]. We then compute L∗L and M∗M to produce
c M = M∗M = " b A 0 − bH I # , L = Lb ∗L = " I Gb 0 AbT # , (2.10) where b A = A(I + GH)−1A, (2.11) b G = G + AG(I + HG)−1AT, (2.12) b H = H + AT(I + HG)−1HA. (2.13)
It is clear that the resulting pencil is still a SSF-1.
(b) Similarly, under the condition Q − P > 0 we can compute a matrix [M∗, L∗] ∈ N (M, L)
with M∗= · A(Q − P )−1 0 −AT(Q − P )−1 I ¸ , L∗= · I −A(Q − P )−1 0 AT(Q − P )−1 ¸ . (2.14)
Direct calculation gives rise to c M = M∗M = " b A 0 b Q −I # , L = Lb ∗L = " − bP I b AT 0 # , (2.15) where b A = A(Q − P )−1A, (2.16) b Q = Q − AT(Q − P )−1A, (2.17) b P = P + A(Q − P )−1AT. (2.18)
The assumption Q − AT(Q − P )−1A ≥ 0 implies that the resulting pencil is still a SSF-2.
Remark 2.3.The proof of Theorem 2.2 provided us with the well defined computation formulae for constructing the special structure preserving doubling transformations, which is the base for us to derive the SDA algorithms for solving the Riccati-type matrix equations.
3
SDA Algorithms for Preserving the SSF-1
In this section we shall use the theory established in the last section to develop the convergence theory of the SDA algorithms for for solving the DAREs and CAREs. The SDA algorithms pro-posed in [8] and [7] will present in the next two subsections by using the doubling transformations defined in the last section.
3.1
SDA Algorithm for Solving the DAREs
It is easy to verify that the DARE (1.2) has a symmetric positive semi-definite solution X (i.e., X ≥ 0) if and only if X satisfies that
M · I X ¸ = L · I X ¸ S (3.1)
for some matrix S ∈ Rn×n, where
M = · A 0 −H I ¸ , L = · I G 0 AT ¸ . (3.2)
Notice that the pencil M − λL is just a SSF-1 form. Therefore, applying the special doubling transformation defined by (2.11)–(2.13) repeatedly gives rise to the following structure-preserving doubling algorithm. Algorithm SDA-1. A0= A, G0= G, H0= H, Ak+1= Ak(I + GkHk)−1Ak Gk+1= Gk+ AkGk(I + HkGk) −1AT k, Hk+1= Hk+ ATk(I + HkGk)−1HkAk.
This is the SDA algorithm described in [8], in which extensive numerical experiments show that this algorithm is efficient and competitive.
3.2
SDA Algorithm for Solving the CAREs
Assume that X ≥ 0 solves the CARE (1.1). It is well-known that the CARE (1.1) can be rewritten as H · I X ¸ = · I X ¸ R, (3.3) where H = · A −G −H −AT ¸ , R = A − GX.
The matrix H is a Hamiltonian matrix, i.e., (HJ)T = HJ. Using Cayley transformation with some
appropriate γ > 0, we can transform (3.3) into the following form M · I X ¸ = L · I X ¸ S, (3.4) where
M = H + γI, L = H − γI, S = (R − γI)−1(R + γI).
Now assume that we have chosen a γ > 0 such that the matrices Aγ = A − γI and Wγ= ATγ + HA
−1
γ G (3.5)
are nonsingular. Chu et al. [7] proposed a method for computing γ such that both Aγ and Wγ are
well-conditioned. Let T1= · A−1 γ 0 HA−1 γ I ¸ , T2= · I −A−1 γ GWγ−1 0 −W−1 γ ¸ , (3.6)
which are obtained by alternatively applying block Gaussian elimination to the matrices L and M (see [7] for more details). Then, direct calculations gives rise to
c M = T2T1M = " b A 0 − bH I # , L = Tb 2T1L = " I Gb 0 AbT # , where b A = I + 2γW−T γ , G = 2γAb −1 γ GW −1 γ , H = 2γWb −1 γ HA −1 γ .
Here the Sherman-Morrison-Woodbury formula is used. Since γ > 0 and H, G ≥ 0 implies that b
G, bH ≥ 0, it follows that the resulting pencil cM − λbL is a SSF-1. In addition, it follows from (3.4) that c M · I X ¸ = bL · I X ¸ S, (3.7)
Thus, beginning with (3.7), following the same lines as Algorithm SDA-1 for solving the DARE we can construct a matrix sequence to approximate to the unique symmetric positive semi-definite solution X to the CARE (1.1). For more details see [7].
3.3
Convergence Analysis of Algorithm SDA-1
Now we establish the convergence theory of Algorithm SDA-1 based on Theorem 2.1. The main results are listed in the following theorem.
Theorem 3.1. Assume that X, Y ≥ 0 satisfy that
X = ATX(I + GX)−1A + H, (3.8)
where H, G ≥ 0, that is, assume that the DARE (3.8) and its dual DARE (3.9) have symmetric positive semi-definite solutions X and Y , respectively, and let
S = (I + GX)−1A, T = (I + HY )−1AT. (3.10)
Then the matrix sequences {Ak}, {Gk}, and {Hk} generated by Algorithm SDA-1 satisfy that
(a) Ak= (I + GkX)S2 k ; (b) H ≤ Hk ≤ Hk+1≤ X and X − Hk= (ST)2 k (X + XGkX)S2 k ≤ (ST)2k (X + XY X)S2k; (3.11) (c) G ≤ Gk ≤ Gk+1≤ Y and Y − Gk= (TT)2 k (Y + Y HkY )T2 k ≤ (TT)2k (Y + Y XY )T2k. (3.12) Proof. Notice that U, V ≥ 0 implies that I +U V is nonsingular and V (I +U V )−1, (I +U V )−1U ≥ 0,
it follows that Algorithm SDA-1 is well defined and
H = H0≤ Hk≤ Hk+1 and G = G0≤ Gk≤ Gk+1. (3.13) Define Mk= · Ak 0 −Hk I ¸ , Lk= · I Gk 0 AT k ¸ .
Then the pencil Mk+1− λLk+1 is a doubling transformation of the pencil Mk− λLk. Since (3.8)
implies that M0 · I X ¸ = L0 · I X ¸ S, (3.14)
where S is defined by (3.10), applying (b) of Theorem 2.1 repeatedly, we get Mk · I X ¸ = Lk · I X ¸ S2k. (3.15)
Equating the blocks of (3.15) yields that
Ak = (I + GkX)S2 k , (3.16) X − Hk= ATkXS 2k . (3.17)
Combining (3.16) with (3.17) gives rise to X − Hk= (ST)2
k
(X + XGkX)S2
k
. (3.18)
This, together with (I + XGk)X ≥ 0, implies that X − Hk ≥ 0, i.e., X ≥ Hk.
Similarly, (3.9) can be rewritten as M0 · −Y I ¸ T = L0 · −Y I ¸ , (3.19)
where T is defined by (3.10), and from (3.19) we can derive that Y − Gk= (TT)2
k
(Y + Y HkY )T2
k
, which implies that Y ≥ Gk. Thus, the theorem is proved.
Let W = " L · I X ¸ , M · −Y I ¸# , Z = · I −Y X I ¸ .
Noting that M0= M, L0= L and X, Y ≥ 0, it follows from (3.14) and (3.19) that W and Z are
nonsingular and satisfy
W−1M Z = · S 0 0 I ¸ , W−1LZ = · I 0 0 T ¸ .
Thus, by the spectral feature of a symplectic pencil, it follows that if ρ(S) < 1, then it must have ρ(T ) = ρ(S) < 1. In addition, it is well known that 0 ≤ U ≤ V implies that kU k2 ≤ kV k2.
Consequently, from Theorem 3.1 we immediately get the following convergence result of Algorithm SDA-1.
Corollary 3.2. Under the hypothesis of Theorem 3.1, ifρ(S) < 1, then we have (a) kAkk2≤ (1 + kXk2kY k2)kS2 k k2−→ 0 as k → 0; (b) kX − Hkk2≤ kX + XY Xk2kS2 k k2 2−→ 0 as k → 0; (c) kY − Gkk2≤ kY + Y XY k2kT2 k k2 2−→ 0 as k → 0.
Remark 3.1.The similar convergence result were obtained in [8]. In contrast to the work of [8], however, our analysis is much simpler and our convergence results are much stronger.
Remark 3.2. Let G = BR−1BT ≥ 0, with R > 0, and H = CTC ≥ 0 in the DARE (3.8),
and assume that (A, B) is stabilizable and (A, C) is detectable. Then it is well known that the DARE (3.8) and its dual DARE (3.9) have symmetric positive semi-definite solutions X and Y , respectively, and ρ(S) < 1 (see, for example, [23, 27], for details). Thus the conditions in Corollary 3.1 are satisfied. In fact, it is easy to verify that if the DARE (3.8) and its dual DARE (3.9) have symmetric positive semi-definite solutions X and Y , respectively, and ρ(S) < 1, with S = (I + GX)−1A, then (A, B) is stabilizable and (A, C) is detectable.
Remark 3.3.By Theorem 3.1 the matrix sequences {Hk} and {Gk} are monotonically increasing
and bounded above, and hence there exist symmetric positive semi-definite matrices ¯H and ¯G such that
lim
k→∞Hk= ¯H, k→∞lim Gk = ¯G.
Corollary 3.2 tells us that if ρ(S) < 1, then X = ¯H and Y = ¯G. Now it is natural to ask if this is still true without the condition ρ(S) < 1. This is a very interesting and worthwhile problem. See the following two simple examples.
Example 3.1.Let A = · α1 0 0 α2 ¸ , G = · ζ1 0 0 ζ2 ¸ , H = · θ1 0 0 0 ¸ , where α1, α2, ζ1, ζ2, θ1 ∈ R with ζ2 ≥ 0, ζ1, θ1 > 0 and α
2 1
1+ζ1 + θ1 = 1. Direct calculation verifies
that with these given data the DARE (1.2) has X = ·
1 0 0 0 ¸
as its symmetric positive semi-definite solution, and in this case
S = (I − GX)−1A = · α1 1+ζ1 0 0 α2 ¸ .
Since α2may be any real number, ρ(S) can be any positive real number. Applying Algorithm
SDA-1 to these given data we have Ak= " α1(k) 0 0 α(k)2 # , Gk = " ζ1(k) 0 0 ζ2(k) # , Hk= · θ(k)1 0 0 0 ¸ ,
where α(0)1 = α1, α(k+1)1 = (α(k)1 )2 1 + θ(k)1 ζ1(k), α(0)2 = α2, α(k+1)2 = (α (k) 2 )2, ζ1(0)= ζ1, ζ1(k+1)= ζ (k) 1 (1 + α (k+1) 1 ), ζ2(0)= ζ2, ζ2(k+1)= ζ (k) 2 (1 + α (k+1) 2 ), θ(0)1 = θ1, θ(k+1)1 = θ (k) 1 (1 + α (k+1) 1 ).
From the proof of Theorem 3.1 we can see that in this case the equality X − Hk= (ST)2
k
(X + XGkX)S2
k
still holds, and hence, letting λ1= 1+ζα11, we have
1 − θ(k)1 = (1 + ζ1(k))λ2k 1 → 0 (k → ∞), since |λ1| = ¯ ¯ ¯ α1 1+ζ1 ¯ ¯
¯ < 1 and 0 ≤ ζ1(k) ≤ ζθ11, which can be proved by taking ζ2 = 0 and using
the same method as the proof of (c) of Theorem 3.1. This shows that in this case we still have limk→∞Hk= X. However, when α2> 1, limk→∞α
(k) 2 = +∞. Example 3.2.Let A = · α1 0 0 1 ¸ , G = · ζ1 0 0 0 ¸ , H = · θ1 0 0 0 ¸ , where α1, ζ1, θ1 ∈ R with ζ1, θ1 > 0 and α
2 1
1+ζ1 + θ1 = 1. It easy to verify that with these given
data, for any ξ ≥ 0, the matrix Xξ =
· 1 0 0 ξ ¸
≥ 0 is a symmetric positive semi-definite solution of the DARE (1.2), which means that in this case the DARE (1.2) have infinitely many symmetric positive semi-definite solutions. Simple calculation gives rise to
S = (I − GXξ)−1A = · α1 1+ζ1 0 0 1 ¸ ,
and so we have that ρ(S) = 1 for any ξ ≥ 0. We can easily prove that in this case we still have lim k→∞Hk= X0= · 1 0 0 0 ¸ , lim k→∞Gk = Y = · ν 0 0 0 ¸ ,
where Hkand Gkare generated by applying Algorithm SDA-1 to these given data, Y is a symmetric
positive semi-definite solution of the dual DARE with ν = ζ1
θ1.
4
SDA Algorithms for Preserving the SSF-2
In this section, we shall first use the doubling transformations defined in the last section to de-rive two SDA algorithms for solving the NME-Ps and NME-Ms. Then, we shall use the theory established in the last section to develop the convergence theory of these SDA algorithms.
4.1
SDA Algorithm for Solving the NME-Ps
It is easy to verify that the NME-P (1.3) has a symmetric positive definite X (i.e., X > 0) if and only if X satisfies ·
for some matrix S ∈ Rn×n, where M = · A 0 Q −I ¸ , L = · 0 I AT 0 ¸ . (4.2)
Notice that the pencil M − λL is just a SSF-2. Therefore, applying the special doubling transfor-mation defined by (2.16)–(2.18) repeatedly gives rise to the following structure-preserving doubling algorithm. Algorithm SDA-2. A0= A, Q0= Q, P0= 0, Ak+1= Ak(Qk− Pk) −1A k, Qk+1= Qk− ATk(Qk− Pk)−1Ak, Pk+1= Pk+ Ak(Qk− Pk)−1ATk.
Remark 4.1.Of course, to ensure that this iteration is well defined, the matrix Qk− Pk must be
symmetric positive definite for all k. Below we shall prove that this condition can be guaranteed if the NME-P (1.3) has a symmetric positive solution (see Theorem 4.1).
Remark 4.2.It is interesting to note that this algorithm is essentially the same as proposed in [26]. In other words, Meini’s algorithm in [26] is a SDA algorithm. It has pointed out that this algorithm has very nice numerical behavior: having quadratical convergence rate, low computational cost per step, and good numerical stability. For more details see [26, 16].
4.2
SDA Algorithm for Solving the NEM-Ms
It is proved in [11] that there always exists a unique positive definite solution X to the NME-M
X − ATX−1A = Q, (4.3)
and moreover, the spectral radius of X−1A is strictly less than 1. The solution X is closely related
the generalized eigenspace of the following pencil M − λL = · A 0 −Q I ¸ − λ · 0 I AT 0 ¸ . (4.4)
In fact, it is easy to verify that a symmetric positive definite matrix X is a solution to the NME-M (4.3) if and only if X satisfies that
M · I X ¸ = L · I X ¸ S (4.5)
for some matrix S ∈ Rn×n.
Although the pencil (4.4) is not symplectic, we can use the same technique as described in Section 2 to transform it into a symplectic pencil. Take
M∗= · AQ−1 0 ATQ−1 −I ¸ , L∗= · I AQ−1 0 ATQ−1 ¸ , then we have M∗L = L∗M. (4.6)
Direct calculation leads to c M0= M∗M = " b A 0 b Q −I # , Lb0= L∗L = " b P I b AT 0 # , where b A = AQ−1A, Q = Q + Ab TQ−1A, P = AQb −1AT. (4.7)
Direct verification shows that cM0− λbL0 is a symplectic pencil, but neither a SSF-1 nor a SSF-2.
Assume that X > 0 is the unique symmetric positive solution to the NME-M (4.3). Then it satisfies the equality (4.5) with S = X−1A. Similar to the proof of (b) of Theorem 2.1 we can
show that c M0 · I X ¸ = bL0 · I X ¸ S2. (4.8) Now Let · I b X ¸ = ·I 0 b P I ¸ · I X ¸ , M =c " b A 0 b Q + bP −I # , L =b · 0 I b AT 0 ¸ . Then it follows from (4.8) that
c M · I b X ¸ = bL · I b X ¸ S2. (4.9)
Clearly, now the pencil cM − λbL is a SSF-2. Thus, beginning with (4.9), following the same lines as Algorithm SDA-2 for solving the NME-P (1.3) we can construct a matrix sequence to approximate to bX. Then, the unique symmetric positive definite solution X to the NME-M (4.3) can be obtained by computing X = bX − bP .
4.3
Convergence Analysis of Algorithm SDA-2
Now we establish the convergence theory of Algorithm SDA-2 based on Theorem 2.1. The main results are listed in the following theorem.
Theorem 4.1. Assume that X > 0 satisfy that
X + ATX−1A = Q, (4.10)
whereQ > 0, and let S = X−1A. Then the matrix sequences {A
k}, {Qk}, and {Pk} generated by
Algorithm SDA-2 satisfy that (a) Ak= (X − Pk)S2 k ; (b) 0 ≤ Pk ≤ Pk+1< X and Qk− Pk= (X − Pk) + AkT(X − Pk)−1Ak > 0; (4.11) (c) X ≤ Qk+1≤ Qk ≤ Q and Qk− X = (ST)2 k (X − Pk)S2 k ≤ (ST)2k XS2k. (4.12) Proof. (Apply Mathematical Induction.) Denote
Mk= · Ak 0 Qk −I ¸ , Lk = · −Pk I AT k 0 ¸ , where P0= 0.
For k = 1, since Q0− P0= Q > 0, it follows that A1, Q1, P1 are all well defined. Using the
equality (4.10) we have · X A AT Q ¸ = · I 0 ATX−1 I ¸ · X 0 0 X ¸ · I X−1A 0 I ¸ > 0. (4.13) Simple computation yields that
· I −AQ−1 0 I ¸ · X A AT Q ¸ · I 0 −Q−1AT I ¸ = · X − AQ−1AT 0 0 Q ¸ . (4.14) Combining (4.14) and (4.13) we get
On the other hand, From (4.10) it is easy to verify that X satisfies M0 · I X ¸ = L0 · I X ¸ S with S = X−1A. Since M
1− λL1 is a doubling transformation of M0− λL0, Applying (b) of
Theorem 2.1 we get M1 · I X ¸ = L1 · I X ¸ S2. (4.16)
Equating the blocks of (4.16) gives rise to
A1= (X − P1)S2, Q1− X = AT1S2.
This, together with (4.15), implies that
Q1− P1= (X − P1) + A1T(X − P1)−1A1> 0 (4.17)
Q1− X = (ST)2(X − P1)S2≥ 0. (4.18)
Obviously, the inequalities Q = Q0 ≥ Q1 and 0 = P0≤ P1 hold. Thus, we have proved that the
theorem is true for k = 1
Next, assume that the theorem is true for all positive integers less or equal to k. Consider the case of k + 1. Since Qk− Pk> 0, it follows that Ak+1, Qk+1, Pk+1are all well defined. Similar to
the proof of (4.15), using the equality (4.11) we can prove that
X − Pk+1= (X − Pk) − Ak(Qk− Pk)−1ATk > 0.
On the other hand, since Mj+1− λLj+1 is a doubling transformation of Mj − λLj for j =
0, 1, . . . k, By using (b) of Theorem 2.1 k + 1 times, we get Mk+1 · I X ¸ = Lk+1 · I X ¸ S2k+1 . (4.19)
From (4.19), following the same lines as the proof of (4.17) and (4.18) it can be proved that Qk+1− Pk+1= (X − Pk+1) + Ak+1T (X − Pk+1)−1Ak+1> 0, Qk+1− X = (ST)2 k+1 (X − Pk+1)S2 k+1 ≥ 0.
Clearly, Pk≤ Pk+1 and Qk ≥ Qk+1 are true. This shows that the theorem is also true for integer
k + 1. By induction principle the theorem is true for all positive integers.
Remark 4.3.The similar results were obtained in [26] by using some properties of cyclic reduction and the spectral properties of block Toeplitz matrices having nonnegative definite matrix-valued generating functions. In contrast to the work of [26], however, our analysis is much simpler and the results are much stronger.
It was proved in [10] that if the matrix equation (1.3) has a symmetric positive definite solution, then all symmetric solutions are positive definite and it has a maximal and minimal solution X+
and X−, respectively. Since Theorem 4.1 is true for any symmetric positive definite solution X to
the NME-P (1.3), the following result follows immediately. Corollary 4.2. Under the hypothesis of Theorem 4.1, we have
Qk− Pk> X+− X−≥ 0
for allk, where X+ andX− are the maximal and minimal solution of (1.3), respectively.
Corollary 4.3. Under the hypothesis of Theorem 4.1, ifρ(S) < 1, then we have (a) kAkk2≤ kXk2kS2 k k2−→ 0 as k → 0; (b) kX − Qkk2≤ kXk2kS2 k k2 2−→ 0 as k → 0.
Remark 4.4.By Theorem 4.1 the matrix sequence {Qk} is monotonically decreasing and bounded
below by X > 0, and hence there exists ¯Q > 0 such that limk→∞Qk = ¯Q. Corollary 4.3 tells us
that if ρ(S) < 1, then X = ¯Q. In fact, ρ(S) < 1 implies that X is the maximal solution of (1.3), and moreover, it has been proved that X is the maximal solution of (1.3) if and only if ρ(S) ≤ 1 (see [16]). Now assume that X = X+ is the maximal solution of (1.3), it is natural to ask whether
¯
Q = X+ still holds if ρ(S) = 1. In [16] Guo proved that if ρ(S) = 1 and all eigenvalues of S on
the unit circle are semisimple, then ¯Q = X+is sill true and in this case the convergence is at least
linear with rate 1/2. When S has nonsemisimple unimodular eigenvalues, it is still open whether ¯
Q = X+.
Remark 4.5.It is proved that the NME-P (1.3) has a symmetric positive definite solution X If and only if the nonlinear matrix equation
Y + AY−1AT = Q (4.20)
has a symmetric positive solution Y (see, for example, [26]). Assume that the maximal solution of (4.20) is Y+. Then it follows from (4.20) that
· A 0 Q −I ¸ · I Q − Y+ ¸ T = · 0 I AT 0 ¸ · I Q − Y+ ¸ , (4.21) where T = Y−1
+ AT. Similar to the proof of (4.12), from (4.21) we can show that
0 ≤ Q − Y+− Pk= (TT)2 k (Qk− Q + Y+)T2 k ≤ (TT)2k Y+T2 k ,
where Pk and Qk are generated by Algorithm SDA-2. Since ρ(T ) = ρ(Y+−1AT) = ρ(X+−1A) (see,
for example, [26]), where X+is the maximal solution of the NME-P (1.3), it follows that under the
conditions of Corollary 4.3 we have limk→∞Pk= Q − Y+. If A is nonsingular, then X−= Q − Y+
(see [26]), where X− is the minimal solution of the NME-P (1.3), and so in this case we have
limk→∞Pk = X−.
Remark 4.6.Since limk→∞(Qk − Pk) = X+− X− if A is nonsingular and ρ(S) < 1, the under
bound X+− X− in Corollary 4.2 is the best one. However, X+− X− may be singular, indeed it
can be zero matrix. For example, the NME-P (1.3) with Q = I and A = 1
2I has X+= X−= 12I.
5
Conclusions
In this paper, we have introduced a structure-preserving transformation of a symplectic pencil, referred to as the doubling transformation, and developed its basic properties. Based on these nice properties of the transformation, a unified convergence theory for the structure-preserving doubling algorithms for solving a class of Riccati-type matrix equations has been established by using only the knowledge from elementary matrix theory.
References
[1] G. Ammar and V. Mehrmann, On Hamiltonian and symplectic Hessenberg forms, Lin. Alg. Appl., 149 (1991), pp. 55–72.
[2] B. Anderson, Second-order convergent algorithms for the steady-state Riccati equation, Int. J. Control, 28 (1978), pp. 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), pp. 279–308.
[5] P. Benner and R. Byers, Evaluating products of matrix pencils and collapsing matrix products, Num. Lin. Alg. Appl., 8 (2001), pp. 357–380.
[6] E.K.-W. Chu, H.-Y. Fan, and W.-W. Lin, A generalized structure-preserving doubling al-gorithm for generalized discrete-time algebraic Riccati equations, preprint 2002-29, NCTS, National Tsing Hua University, Hsinchu 300, Taiwan, 2003.
[7] E.K.-W. Chu, H.-Y. Fan, and W.-W. Lin, A structure-preserving doubling algorithm for cotinuous-time algebraic Riccati equations, to appear in Lin. Alg. Appl..
[8] E.K.-W. Chu, H.-Y. Fan, W.-W. Lin, and C.-S. Wang, A structure-preserving doubling algo-rithm for periodic descrete-time algebraic Riccati equations, preprint 2002-18, NCTS, National Tsing Hua University, Hsinchu 300, Taiwan, 2003.
[9] J.C. Engwerda, On the existence of a positive definite solution of the matrix equation X + ATX−1A = I, Lin. Alg. Appl., 194 (1993), pp. 91–108.
[10] J.C. Engwerda, A.C.M. Ran, and A.L. Rijkeboer, Necessary and sufficient conditions for the existence of a positive definite solution of the matrix equation X + A∗
X−1A = Q, Lin. Alg.
Appl., 186 (1993), pp. 255–275.
[11] A. Ferrante and B.C. Levy, Hermitian solutions of the Equation X = Q + N X−1N∗, Lin. Alg.
Appl., 247 (1996), pp. 359–373.
[12] G.H. Golub and C.F. Van Loan, Matrix Computations, 3rd ed., The Johns Hopkins University Press, 1996.
[13] T. Gudmundsson, C. Kenney, and A. J. Laub, Scaling of the discrete-time algebraic Riccati equation to enhance stability of the Schur solution method, IEEE Trans. Auto. Control, 37 (1992), pp. 513–518.
[14] C.-H. Guo, Newton’s method for discrete algebraic Riccati equations when the closed-loop ma-trix has eigenvalues on the unit circle, SIAM J. Mama-trix Anal. Appl., 20 (1998), pp. 279–294. [15] C.-H. Guo and P. Lancaster, Iterative solution of two matrix equations, Math. Comp.,
68(1999), pp. 1589–1603.
[16] C.-H. Guo, Comvergence rate of an iterative method for a nonlinear matrix equation, SIAM J. Matrix Anal. Appl., 23 (2001), pp. 295–302.
[17] J.J. Hench and A.J. Laub, Numerical solution of the discrete-time periodic Riccati equation, IEEE Trans. Auto. Control, 39 (1994), pp. 1197–1210.
[18] M. Kimura, Convergence of the doubling algorithm for the discrete-time algebraic Riccati equation, Int. J. Systems Science, 19 (1988), pp. 701–711.
[19] D. Lainiotis, N. Assimakis, and S. Katsikas, New doubling algorithm for the discrete periodic Riccati equation, Applied Mathematics and Computation, 60 (1994), pp. 265–283.
[20] P. Lancaster and L. Rodman, Algebraic Riccati Equations, Oxford University Press, Oxford, 1995.
[21] A.J. Laub, A Schur method for solving algebraic Riccati equations, IEEE Trans. Auto. Control, 24 (1979), pp. 913–921.
[22] A.J. Laub, Invariant subspace methods for numerical solution of Riccati equation, in The Raccati Equation, S. Bittanti, A.J. Laub, and J.C. Willems eds, 163–196, Springer-Verlag, Berl
[23] L.-Z. Lu and W.-W. Lin, An iterative algorithm for the solution of the discrete time algebraic Riccati equations, Lin. Alg. Appl., 189 (1993), pp. 465–488.
[24] L.-Z. Lu, W.-W. Lin, and C. E. M. Pearce, An efficient algorithm for the discrete-time algebraic Riccati equation, IEEE Trans. Auto. Control, 44 (1999), pp. 1216–1220.
[25] V. Mehrmann, The Autonomous Linear Quadratic Control Problem, Lecture Notes in Control and Information Sciences, 163, Springer-Verlag, Berlin, 1991.
[26] B. Meini, Efficient computation of the extreme solutions of X + A∗X−1A = Q and X −
A∗X−1A = Q, Math. Comp., 71 (2001), pp. 1189–1204.
[27] C. Paige and C. Van Loan, A Schur decomposition for Hamiltonian matrices, Lin. Alg. Appl., 41 (1981), pp. 11–32.
[28] T. Pappas, A. J. Laub, and N. R. Sandell, On the numerical solution of the discrete-time algebraic Riccati equation, IEEE Trans. Auto. Control, 25 (1980), pp. 631–641.
[29] J.-G. Sun, Sensitivity analysis of the discrete-time algebraic Riccati equation, Lin. Alg. Appl., 275/276 (1998), pp. 595–615.
[30] J.-G. Sun and S.F. Xu, Perturbation analysis of the maximal solution of the matrix equation X + ATX−1A = P.II, Lin. Alg. Appl., 362 (2003), pp. 211–228.
[31] P. Van Dooren, A generalized eigenvalue approach for solving Riccati equations, SIAM J. Sci. & Statistical Computing, 2 (1981), pp. 121–135.
[32] S.F. Xu, On the maximal solution of the matrix equation X + ATX−1A = I, Acta Scientiarum
