• 沒有找到結果。

On A Doubling Algorithm for The Nonlinear Matrix Equation when

N/A
N/A
Protected

Academic year: 2021

Share "On A Doubling Algorithm for The Nonlinear Matrix Equation when"

Copied!
16
0
0

加載中.... (立即查看全文)

全文

(1)

On A Doubling Algorithm for The Nonlinear Matrix Equation

X

+ A

T

X

−1

A

= Q when |λ(X

−1

A

)| ≤ 1

Eric King-wah Chu

Tsung-Min Hwang

Wen-Wei Lin

§

Shu-Fang Xu

Abstract

In [15], a structure-preserving doubling algorithm (SDA) was proposed for the solution of the nonlinear matrix equation X + ATX−1

A= Q (with Q > 0). The SDA was proven to be convergent quadratically to the maximal solution X+ when all eigenvalues of X+−1Alie inside the unit circle. In this paper, we

prove that SDA converges linearly to a symmetric solution X+ when all eigenvalues of X+−1Aare inside

or on the unit circle and the partial multiplicity of each unimodular eigenvalue is half of the partial multiplicity of the corresponding unimodular eigenvalue of the associated symplectic pencil.

Keywords. doubling, nonlinear matrix equations, unimodular eigenvalue, structure-preserving AMS subject classifications. 15A18, 65F15

1

Introduction

In this paper, we are interest in the study of the nonlinear matrix equation (NME)

X + ATX−1A = Q, (1)

where A, Q ∈Rn×nwith Q being symmetric positive definite. NMEs occur frequently in many applications.

Notable examples include algebraic Riccati equations [1, 2, 10, 11, 15], quadratic matrix equations [6, 7, 8] AX2

+ BX + C = 0,

where A, B, C are given coefficient matrices. Various aspects of NMEs, like solvability, numerical solution, perturbation and applications, can be found in [3, 4, 8, 18, 19, 20, 21, 22] and the references therein. One main application of the NME is related to the solution of quadratic eigenvalue problems [7, 9].

For symmetric matrices X and Y , we write X ≥ Y (X > Y ) if X − Y is positive semidefinite (definite). A symmetric solution X+ of an NME (1) is called maximal if X+ ≥ X for any symmetric solution X of (1).

In [15], a structure-preserving doubling algorithm (SDA) was proposed for finding the maximal X+ of

the NME (1), and moreover, it was proven that this algorithm converges quadratically when all eigenvalues

Version November 9, 2006

School of Mathematical Sciences, Building 28, Monash University, VIC 3800, Australia; eric.chu@sci.monash.edu.auDepartment of Mathematics, National Taiwan Normal University, Taipei 11677, Taiwan; min@math.nthu.edu.tw §Department of Mathematics, National Tsinghua University, Hsinchu 300, Taiwan; wwlin@am.nthu.edu.twLMAM, School of Mathematical Sciences, Peking University, Beijing 100871, China; xsf@pku.edu.cn

(2)

of X+−1A lie inside the unit circle. The main purpose of this paper is to prove that the SDA converges

linearly when the eigenvalues of X+−1A lie inside or on the unit circle and the partial multiplicity of each

unimodular eigenvalue is half of the partial multiplicity of the corresponding unimodular eigenvalues of the associated symplectic pencil. Note that other iterative solution processes, by Newton’s iteration or cyclic reduction, have been proposed in [6, 8, 18] and linear convergence for problems with semi-simple unimodular eigenvalues has been observed and proved in [6]. Consequently, because of the more general convergence result for unimodular eigenvalues, our doubling algorithm is more applicable than that in [6].

2

Preliminaries

Consider the NME (1) and define M ≡  −A 0 Q I  , L ≡  0 I AT 0  . (2)

It is easy to verify that the pencil M − λL is symplectic, i.e., it satisfies MJMT = LJLT , with J ≡

 0 I −I 0



and λ ∈ σ(M, L) if and only if 1/λ ∈ σ(M, L). Let R = X−1A. Then (1) can be rewritten as

M  I X  = L  I X  R. (3)

The existence of a symmetric positive definite solution and a maximal symmetric positive definite solution of (1) has been established in [4].

Theorem 2.1 [4] The NME (1) has a symmetric positive definite solution if and only if ψ(λ) ≡ λA + Q + λ−1AT is regular, and ψ(λ) ≥ 0, for all |λ| = 1.

Theorem 2.2 [4] If (1) has a symmetric positive definite solution, then it has a maximal symmetric pos-itive definite solution X+. Moreover, for the maximal solution X+, we have ρ(X+−1A) ≤ 1; for any other

symmetric positive definite solution X, we have ρ(X−1A) > 1. Here ρ( · ) denotes the spectral radius.

It follows from [6, Theorem 2.4] that the eigenvalues of the matrix X−1

+ A have the following

characteri-zation.

Theorem 2.3 [6] For (1), the eigenvalues of the matrix X−1

+ A are precisely the eigenvalues of the

ma-trix pencil M − λL inside or on the unit circle, with half of the partial multiplicities for each unimodular eigenvalue.

(3)

3

SDA algorithm

We now consider a more general standard symplectic form (SSF) than that in (2) as follows M ≡  A 0 Q −I  , L ≡  −P I AT 0  , (4)

where Q − P > 0. Then we can compute a matrix [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  (5) which satisfies M?L = L?M. Direct calculations give rise to

ˆ M ≡ M?M =  ˆA 0 ˆ Q −I  , L ≡ Lˆ ?L =  − ˆP I ˆ AT 0  , (6) where ˆ A ≡ A(Q − P )−1A , Q ≡ Q − Aˆ T(Q − P )−1A , P ≡ P + A(Q − P )ˆ −1AT. (7) If ˆQ in (7) is again positive definite, then ( ˆM, ˆL) is again a SSF form.

We restate the structure-preserving doubling algorithm (SDA) developed in [15] as follows. Algorithm 3.1 (SDA for NME)

Input: A, Q; τ (a small tolerance); Output: a symmetric solution X to NME.

Set k = 0, Ak = A, Qk= Q and Pk = 0; Dountil convergence: ComputeAk+1= Ak(Qk− Pk)−1Ak, Qk+1= Qk− ATk(Qk− Pk)−1Ak, Pk+1= Pk+ Ak(Qk− Pk)−1ATk, k = k + 1; IfkQk− Qk−1k ≤ τ kQkk, Stop; End; Set X ← Qk+1.

Remark 3.1 To ensure that the iteration in Algorithm 3.1 is well-defined, the matrix Qk − Pk must be

positive definite for all k. Below we quote from [15, Theorem 2.1] to guarantee this condition is satisfied, provided that the NME (1) has a symmetric positive definite solution.

Theorem 3.1 Assume the NME (1) has a symmetric positive definite solution X. Then the matrix sequence {Ak, Qk, Pk} generated by Algorithm 3.1 is well-defined and satisfies

(i) Ak = (X − Pk)R2 k ; (ii) 0 ≤ Pk≤ Pk+1< X and Qk− Pk= (X − Pk) + ATk(X − Pk)−1Ak > 0; (iii) X ≤ Qk+1≤ Qk≤ Q and Qk− X = (RT)2 k (X − Pk)R2 k ≤ (RT)2kXR2k.

(4)

Remark 3.2 Algorithm 3.1 is essentially the same as Algorithm 3.1 proposed in [18] with Qk− Pk and Qk

in Algorithm 3.1 being replaced by Qk and Xk, respectively. In other words, Algorithm 3.1 in [18] is an SDA.

In [6], Guo proved that if ρ(R) = 1, where R is defined in (3), and all eigenvalues of R on the unit circle are semisimple, then the convergence of Algorithm 3.1 in [18] is at least linear with rate 1/2. But when R has nonsemisimple unimodular eigenvalues, the situation is unclear. In the following section, we shall prove that if Qk− Pk in each iteration of Algorithm 3.1 is invertible and {Pk} is bounded, then the convergence of

Algorithm 3.1 is globally linear with rate 1/2 when the size of Jordan block of each unimodular eigenvalue of (M, L) in (2) is even.

4

Convergence rate of SDA

Consider the matrix pair (M, L) defined in (2). Since Q > 0, it is easy to check that (M, L) is regular, i.e., det(M − λL) 6≡ 0. Suppose (M, L) satisfies the following assumption.

(H) The partial multiplicities of (M, L) associated with the unimodular eigenvalues (if any) are all even. (8) Let Jω,p be the p × p Jordan block with a unimodular eigenvalue ω = eiθ:

Jω,p≡          ω 1 0 · · · 0 0 . .. ... ... ... .. . . .. ... ... 0 .. . . .. ... 1 0 · · · 0 ω          . (9)

By the result in [5, p. 557], the 2k power of J

ω,psatisfies J2k ω,p=       γ1,k γ2,k · · · γp,k . .. ... ... . .. γ2,k 0 γ1,k      , (10) where γi,k≡ 2 k(2k− 1) · · · (2k− i + 2) (i − 1)! ω 2k−i+1, i = 1, · · · , p. (11)

For the unimodular eigenvalues ωj= eiθj of (M, L) with even partial multiplicity p = 2mj, we have

Jωj,2mj =  Jωj,mj Γ1,mj 0 Jωj,mj  , Γ1,mj ≡ emje T 1, (12)

for j = 1, · · · , r. From the theorem of symplectic Kronecker canonical form for (M, L) (see e.g., [12] or [13]) there exist a symplectic matrix bZ (i.e., bZHJ bZ = J) and a nonsingular bQ such that

b QM bZ =  Js⊕ J1 0`⊕ bΓ1 0n I`⊕ J1−H  , (13) b QL bZ =  I`⊕ Im 0n 0n JsH⊕ Im  , (14)

(5)

where Jsis the stable Jordan block of size ` (ρ(Js) < 1), J1= Jω1,m1⊕ · · · ⊕ Jωr,mr, bΓ1= bΓ1,m1⊕ · · · ⊕ bΓ1,mr

with bΓ1,mj = emje

T

mj (j = 1, · · · , r), l = n − (m1+ · · · + mr) ≡ n − m, ⊕ denotes the direct sum of matrices

and the supscription ’H’ denotes the conjugate transpose. It is easily seen that span{ bZ(:, 1 : n)} forms a

unique stable Lagrangian deflating subspace of (M, L) corresponding to Js⊕ J1([13]).

From (12), for each j = 1, · · · , r, there exists a suitable nonsingular Wj∈Cmj×mj such that

 Imj 0 0 Wj−1   Jωj,mj Γb1,mj 0 Jω−Hj,mj   Imj 0 0 Wj  =  Jωj,mj Γ1,mj 0 Jωj,mj  . (15) Set Z = bZ(In+`⊕ W1⊕ · · · ⊕ Wr), Q = (In+`⊕ W1−1⊕ · · · ⊕ W −1 r ) bQ. (16)

Thus, from (15)-(16) equations (13) and (14), respectively, become

QMZ =  Js⊕ J1 0`⊕ Γ1 0n I`⊕ J1  ≡ JM, (17) QLZ =  In 0n 0n JsH⊕ Im  ≡ JL, (18)

where Γ1≡ Γ1,m1⊕ · · · ⊕ Γ1,mr with Γ1,mj being given in (12). It also follows that span{Z(:, 1 : n)} forms

the unique stable Lagrangian deflating subspace of (M, L) corresponding to Js⊕ J1.

Since JMand JL commute with each other, it follows from (17) and (18) that

MZJL= Q−1JLJM= LZJM. (19)

On the other hand, if we interchange the roles of M and L in (17) and (18), and consider the symplectic pair (L, M), there are nonsingular matrices P and Y such that

PLY = JM, PMY = JL, (20)

where span{Y(:, 1 : n)} forms a unique stable Lagrangian deflating subspace of (M, L) corresponding to Js⊕ J1. Similar arguments also produce

LYJL= MYJM. (21)

Let {(Mk, Lk)}∞k=0 be the sequence of symplectic pairs in SSF with

Mk≡  Ak 0 Qk −I  , Lk≡  −Pk I AT k 0  (22) generated by Algorithm 3.1. By the construction of the algorithm we have

Mk+1= M?kMk, Lk+1= L?kLk, (23) where M?k =  Ak(Qk− Pk)−1 0 −AT k(Qk− Pk)−1 I  , L?k =  I −Ak(Qk− Pk)−1 0 AT k(Qk− Pk)−1  ,

(6)

with

M?kLk= L?kMk. (24)

For more details, see [15].

With M0= M and L0= L, it follows from (23), (24) and (19) that

M1ZJ 2 L= M?0M0ZJ 2 L= M?0L0ZJMJL= L?0M0ZJLJM= L?0L0ZJ 2 M= L1ZJ 2 M. (25)

By an inductive process, we have shown

MkZJ2

k

L = LkZJ2

k

M. (26)

Similarly, from (21), it also holds that

LkYJ2

k

L = MkYJ2

k

M. (27)

With ω = ωjand p = 2mj as in (10), for j = 1, · · · , r, using (17) and (18), equation (26) can be expressed

by MkZ  I 0 0 J2k s ⊕ Im  = LkZ " J2k s ⊕ J2 k 1 0`⊕ Γk 0 I`⊕ J2 k 1 # , (28) where Γk= Γk,m1⊕ · · · ⊕ Γk,mr with Γk,mj ≡          γmj+1,k γmj+2,k · · · γ2mj,k γmj,k . .. . .. ... .. . . .. . .. . .. ... γ3,k . .. . .. γmj+2,k γ2,k γ3,k · · · γmj,k γmj+1,k          ≡ J2k ωj,2mj(1 : mj, mj+ 1 : 2mj), (29)

in which γi,k are defined in (11) for i = 2, · · · , 2mj, j = 1, · · · , r.

We now state a useful Lemma in [11] for the convergence of the SDA algorithm.

Lemma 4.1 Let J1≡ Jω1,m1⊕ · · · ⊕ Jωr,mr and Γk ≡ Γk,m1⊕ · · · ⊕ Γk,mr, where Γk,mj is defined in (29),

for j = 1, · · · , r. Then Γk is invertible and satisfies

kΓ−1k J2k 1 k = O(2 −k), kJ2k 1 Γ −1 k J 2k 1 k = O(2 −k) as k → ∞. (30)

We now prove the convergence Theorem for Algorithm 3.1. Partition Z in (19) and Y in (21) by Z =  Z1 Z3 Z2 Z4  , Y =  Y1 Y3 Y2 Y4  , (31) where Zi, Yi∈Cn×n (i = 1, · · · , 4).

(7)

Theorem 4.2 Assume that the NME (1) has a symmetric positive definite solution and the matrix pair (M, L) of (2) satisfies the assumption (H) as in (8). Let X+ be its maximal solution. If ρ(X+−1A) = 1,

then the sequence {Ak, Qk, Pk} generated by Algorithm 3.1 satisfies

(i) kAkk = O  ρ(Js)2 k + O(2−k) → 0 as k → ∞; (ii) kQk− X+k = O  ρ(Js)2 k + O(2−k) → 0 as k → ∞ where X += Z2Z1−1 solves (1); (iii) kPk− X−k = O  ρ(Js)2 k + O(2−k) → 0 as k → ∞ where X −= Y2Y1−1 if Y1is invertible; moreover,

if in addition, A and Y2 are invertible, then X− solves (1);

(iv) {Qk− Pk} → singular matrix as k → ∞.

Proof. (i) From Theorem 3.1 the sequence {Ak, Qk, Pk} is well-defined. Substituting (Mk, Lk) of (22)

and Z of (31) into (28), we obtain

AkZ1 = (−PkZ1+ Z2)(J2 k s ⊕ J 2k 1 ), (32) AkZ3((JsH)2 k ⊕ Im) = (−PkZ1+ Z2)(0`⊕ Γk) + (−PkZ3+ Z4)(I`⊕ J2 k 1 ), (33) QkZ1− Z2 = ATkZ1(J 2k s ⊕ J 2k 1 ), (34) (QkZ3− Z4)((JsH)2 k ⊕ Im) = ATkZ1(0`⊕ Γk) + ATkZ3(I`⊕ J2 k 1 ). (35)

By assumption (H), Theorem 2.2 and 2.3 the maximal solution X+exists and the eigenvalues of R+≡ X+−1A

are inside or on the unit circle, with half of the partial multiplicities for each eigenvalue on the unit circle. Since the stable Lagrangian deflating subspace of (M, L) corresponding to Js⊕ J1 is unique and R+ is

similar to Js⊕ J1, it holds that span([I, X+T]T) = span([Z1T, Z2T]T). Then Z −1 1 exists and X+= Z2Z1−1. Post-multiplying (33) by (0`⊕ Γ−1k J 2k 1 )Z −1

1 , and using (32), we have

Ak h I − Z3(0`⊕ Γ−1k J 2k 1 )Z −1 1 i = (−PkZ1+ Z2)(J2 k s ⊕ 0m)Z1−1− (PkZ3+ Z4)(0`⊕ J2 k 1 Γ −1 k J 2k 1 )Z −1 1 . (36) Let Ek ≡ Z3(0`⊕ Γ−1k J 2k 1 )Z −1

1 . Then by Lemma 4.1 it holds that I − Ek = I + O(2−k) invertible, for k

sufficiently large. Consequently, (36) becomes Ak= h (−PkZ1+ Z2)(J2 k s ⊕ 0m)Z1−1− (−PkZ3+ Z4)(0`⊕ J2 k 1 Γ −1 k J 2k 1 )Z −1 1 i (I − Ek)−1.

By Lemma 4.1 and the boundedness of {Pk} (by Theorem 3.1) the sequence {Ak} satisfies

kAkk = O  ρ(Js)2 k + O(2−k) → 0 as k → ∞. (37) (ii) Post-multiplying (35) by (0`⊕ Γ−1k J 2k 1 )Z −1

1 and using (34), we get

Qk h I − Z3(0`⊕ Γ−1k J 2k 1 )Z −1 1 i − X+ = AT kZ1(J 2k s ⊕ 0m)Z1−1− A T kZ3(0`⊕ J2 k 1 Γ −1 k J 2k 1 )Z −1 1 − Z4(0`⊕ Γ−1k J 2k 1 )Z −1 1 . (38)

(8)

Let Ek= Z3(0`⊕ Γ−1k J 2k 1 )Z

−1

1 . Since I − Ek= I + O(2−k) as k → 0, the equation (38) can be rewritten by

Qk− X+(I − Ek) −1 = hAT kZ1(J 2k s ⊕ 0m)Z1−1− ATkZ3(0`⊕ J2 k 1 Γ −1 k J 2k 1 )Z −1 1 − Z4(0`⊕ Γ−1k J 2k 1 )Z −1 1 i (I − Ek)−1.

By Lemma 4.1 and (37), the matrix sequence {Qk− X+} satisfies

kQk− X+k = O

 ρ(Js)2

k

+ O(2−k) → 0 as k → 0

(iii) Substituting (Mk, Lk) of (22) and Y of (31) into (27), we have

Y2− PkY1 = AkY1  J2k s ⊕ J 2k 1  , (39) (Y4− PkY3)  (JsH) 2k ⊕ Im  = AkY1(0`⊕ Γk) + AkY3  I`⊕ J2 k 1  . (40)

Let X−= Y2Y1−1 .As above, post-multiplying (40) by



0`⊕ Γ−1k J2

k

1



Y1−1 and using (39), we get

X−− Pk h I − Y3  0`⊕ Γ−1k J 2k 1  Y−1 1 i = Y4  0`⊕ Γ−1k J 2k 1  Y1−1− AkY2  0`⊕ J2 k 1 Γ −1 k J 2k 1  Y1−1+ AkY1  J2k s ⊕ 0m  Y1−1. (41)

By Lemma 4.1 and the result of (i), we have kX−− Pkk = O

 ρ(Js)2

k

+ O(2−k) → 0 as k → ∞.

Moreover, substituting (L, M) of (2) and Y of (31) into (21), we get  0 I AT 0   I X−  =  A 0 Q −I   I X−  R−,

where R−≡ Y1(Js⊕ J1)Y1−1. It follows that

X−= AR−, AT = (Q − X−)R−.

Suppose A and Y2are invertible, then it holds that

X−+ ATX−−1A = Q

and

ρ(X−−1A) = ρ(R −1

− ) = ρ((Js⊕ J1)−1) > 1.

(iv) From (32) and (34), we get −PkZ1(J2 k s ⊕ J 2k 1 ) = AkZ1− Z2(J2 k s ⊕ J 2k 1 ), QkZ1(J2 k s ⊕ J 2k 1 ) = Z2(J2 k s ⊕ J 2k 1 ) + A T kZ1(J2 2k s ⊕ J 22k 1 ).

(9)

This implies that (Qk− Pk)Z1  0 Im  = AkZ1  0 J−2k 1  + ATkZ1  0 J2k 1  . (42)

Using the result of (37), the first column of (42) becomes (Qk− Pk)Z1  0 e1  = AkZ1  0 ω−2k e1  + ATkZ1  0 ω2ke 1  → 0, as k → ∞. Therefore, {Qk− Pk} converges to a singular matrix as k → ∞.

In the following, we prove a second convergence theorem for Algorithm 3.1 when the NME (1) only has a symmetric solution.

Theorem 4.3 Assume that the NME (1) has a symmetric solution X+ (not necessarily positive definite)

and (M, L) of (2) satisfies the assumption (H) as in (8), furthermore, the eigenvalues of X+−1A are inside

or on the unit circle, with half of the partial multiplicities for each eigenvalue on the unit circle. Suppose the sequence {Ak, Qk, Pk} generated by Algorithm 3.1 is well-defined for each k. If kPkk is bounded, then it

holds that (i)-(iv) as in Theorem 4.2.

Proof. By assumptions we have that the sequence {Ak, Qk, Pk} is well-defined and {Pk} is bounded.

Fur-thermore, σ(X−1

+ A) lies inside or on the unit circle, with half of the partial multiplicities for each eigenvalue

on the unit circle, and X−1

+ A is similar to Js⊕ J1, so the matrices [I, X+T]T and [Z1T, Z2T]T span the unique

stable Lagrangian deflating subspace of (M, L) corresponding to Js⊕ J1, and then X+= Z2Z1−1. Bases on

the above properties, assertions (i)-(iv) hold using the same argument as in Theorem 4.2.

Remark 4.1 (i) In Theorem 4.2 and 4.3, we see that the sequence {Qk} and {Pk} converges to X+ and

X−, respectively, with convergence rate 1/2

(ii) Note that the symmetric solution X+ in Theorem 4.2 is a maximal symmetric positive definite solution

for NME (1), but in Theorem 4.3 it is not necessarily a positive definite solution for NME (1).

5

Numerical Results

In this section, we test five numerical examples satisfying assumption (H) to illustrate the convergence behavior of Algorithm 3.1. The first three examples have symmetric positive definite solutions and the other two examples have not. All computations were performed in MATLAB R2006b 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.

In our numerical tests, we report the numbers of iterations by “ITs” for Algorithm 3.1, and the maxi-mal error between the accurate and the approximate stable eigenvalues of X+−1A by “Err+”, where X+ is

(10)

Example 5.1 Example 5.2 Example 5.3 Example 5.4 Example 5.5 NRes+ 9.49 × 10−17 1.29 × 10−14 9.52 × 10−10 3.99 × 10−14 3.67 × 10−9 Err+ 1.04 × 10−8 4.59 × 10−3 1.84 × 10−4 6.83 × 10−3 5.40 × 10−5 ITs 24 16 14 16 15 NRes− 1.38 × 10−6 1.73 × 10−14 7.88 × 10−10 1.40 × 10−14 3.29 × 10−9 Err− ∗ 4.61 × 10−3 1.84 × 10−4 6.87 × 10−3 5.41 × 10−5

Table 1: Numerical results for various examples.

an approximate solution to the NME. Let λi be the exact stable eigenvalue and ˜λi be the corresponding

approximate eigenvalue, we have

Err+= max

1≤i≤n{|λi− ˜λi|}.

Similarly, “Err−” represents the maximal error between the accurate and the approximate unstable

eigen-values of X−1 − A.

To measure accuracy, we use the “normalized” residual NRes±≡

kX±+ ATX±−1A − QkF

kX±kF + kATX±−1AkF+ kQkF

,

where X± is an approximate solution to the NME. The algorithm is terminated when NRes+ cannot be

reduced further.

The first example is provided by Meini [18].

Example 5.1 [18] Let n = 100, Q = In and A be symmetric whose entries are defined by the scheme:

For i = 1, . . . , n: Set aij = i2+ j for j = i, . . . , n; Compute s1= Pi−1 j=1aij, s2= Pn j=iaij; Set aij = aij(1/2 − s1)/s2, aji= aij for j = i, . . . , n.

The symplectic pair (M, L) in (2) has an eigenvalue 1 with partial multiplicity 2. The maximal solution equals X = 1

2(I + (I − 4A

TA)1/2) (see [22]).

The numerical results for X± are reported in Table 1. The notation “∗” in the Table 1 means that the

value of Err− cannot be computed since the exact unstable eigenvalues are not known. On the other hand,

the relative error kX+− XkF/kXkF equals 5.21 × 10−10 and the sequence {Pk} converges to X− which is

also symmetric positive definite.

In the following, we state the processes of constructing matrices A and Q in the NMEs for Examples 5.2– 5.5. Given matrices A0, G0 ∈Rn×n such that the partial multiplicities of the purely imaginary eigenvalues

of the Hamiltonian matrix



A0 −G0

0 −AT 0

(11)

are all even, we construct a Hamiltonian matrix H as follows: H ≡  ˜ A0 − ˜G0 − ˜H0 − ˜AT0  =  I V 0 I   A0 −G0 0 −AT 0   I −V 0 I  , (43) where V ∈Rn×n.

Using a Cayley transformation with an appropriate γ > 0, the Hamiltonian matrix H in (43) can be transformed to a symplectic pair (M, L) ≡ (H + γI, H − γI) [16, 17], and then equivalently simplifies to a symplectic pair (M0, L0). Here

M0=  ˆ A0 0 − ˆH0 I  , L0=  I Gˆ0 0 AˆT 0  , (44) with ˆ A0 = I + 2γ(Aγ+ ˜G0A−Tγ H˜0)−1, (45a) ˆ G0 = 2γA−1r G˜0(ATγ + ˜H0A−1γ G˜0)−1, (45b) ˆ H0 = 2γ(ATγ + ˜H0A−1γ G˜0)−1H˜0A−1γ , (45c)

and Aγ ≡ ˜A0− γI. Note that ˜H0 = ˆH0 = 0 and ˆG0 is symmetric. Choose V so that ˆG0 is nonsingular.

Then using following similarity transformation, the symplectic pair (M0, L0) in (44) can be transformed to

a symplectic pair ( ˆM, ˆL) with SSF form: ˆ M ≡  I 0 0 −I   I 0 − ˆAT 0 I   ˆ G−1 0 0 0 I  M0  I 0 − ˆG−1 0 I  ≡  A 0 Q −I  (46a) and ˆ L ≡  I 0 0 −I   I 0 − ˆAT 0 I   ˆ G−10 0 0 I  L0  I 0 − ˆG−10 I  =  0 I AT 0  , (46b) where Q = ˆH0+ ˆAT0Gˆ −1 0 Aˆ0+ ˆG−10 = ˆAT0Gˆ −1 0 Aˆ0+ ˆG−10 is symmetric.

Remark 5.1 In Examples 5.2–5.5, given A0, G0 and an appropriate γ > 0 in the Cayley transformation,

we shall choose V so that ˆG0is nonsingular and Q is symmetric positive definite. From Theorem 2.1, we can

numerically checked whether the associated NME has at least one or no symmetric positive definite solution, if ψ(λ) ≥ 0 for all |λ| = 1 or if ψ(λ) 0 for some |λ| = 1, where ψ(λ) is defined in Theorem 2.1.

Example 5.2 [14] Let A0 =  0 1 0 0  ⊕   00 10 01 0 0 0   ⊕ (−1 3I2), (47a) G0 =  0 0 0 1  ⊕   00 00 00 0 0 1   ⊕ I2. (47b)

(12)

Then H in (43) has nonzero eigenvalues {−1/3, −1/3, 1/3, 1/3} and the zero eigenvalue has partial multi-plicities {4, 6}. Taking V =           2 9 0 9 −1 −3 −3 9 3 3 1 1 1 −1 −1 −1 −4 3 0 3 −10           .

and γ = 4, then ˆG0 is nonsingular and Q is symmetric positive definite. The symplectic pair ( ˆM, ˆL) in

(46) has eigenvalues {−13/11, −13/11, −11/13, −11/13} and {−1, −1} with partial multiplicities {4, 6}. The matrix ψ(λ) is numerically checked to be symmetric positive semidefinite for λ = eı2π`/1000, ` = 0, 1, . . . , 1000. By Theorem 2.1, the NME (1) has a symmetric positive definite solution.

The numerical results for X±are reported in Table 1. In this example, all eigenvalues of X+are positive

and X+ is a maximal symmetric positive definite solution.

Example 5.3 Let A0 and G0 be 5 × 5 real matrices defined by

A0=   U0 IU2 00 0 0 2   , G0=   00 I02 00 0 0 1   (48) where U =  0 2 −2 0 

. Then H in (43) has nonzero eigenvalues {2, −2} and pure imaginary eigenvalues {2ı, −2ı} with partial multiplicity 4. Let B be a 5 × 5 uniformly distributed random matrix, generated with seed = 1307483728 in the MATLAB command rand, and consider the QR decomposition B = QR. Taking V = Q and γ = 1, then ˆG0is nonsingular and Q is symmetric positive definite. The symplectic pair ( ˆM, ˆL)

in (46) has eigenvalues {3, 1/3} and unimodular eigenvalues {0.6 + 0.8ı, 0.6 − 0.8ı} with partial multiplicities 4. Numerically, it can be checked that ψ(λ) ≥ 0.

The numerical results for X±are reported in Table 1. In this example, all eigenvalues of X+are positive

and X+ is a maximal symmetric positive definite solution.

Example 5.4 [14] Let A0 and G0 be defined in (47). Take

V =           8 5 0 5 2 −3 −3 6 8 8 4 −8 −8 −10 −1 −1 −6 −1 0 −1 −6          

and γ = 4, then ˆG0 is nonsingular and Q is symmetric positive definite. The symplectic pair ( ˆM, ˆL) in

(13)

X+ X− X+− X− λ7 −7.149 × 10−2 −8.073 × 10−2 1.879 × 10−10 λ6 9.923 × 10−2 8.720 × 10−2 1.867 × 10−09 λ5 1.828 × 10−1 1.798 × 10−1 1.543 × 10−05 λ4 4.117 × 10−1 2.948 × 10−1 1.253 × 10−03 λ3 5.363 × 10−1 3.839 × 10−1 5.955 × 10−02 λ2 3.020 × 10+0 3.009 × 10+0 1.179 × 10−01 λ1 4.724 × 10+2 4.723 × 10+2 1.542 × 10−01

Table 2: Eigenvalues of X± and X+− X− for Example 5.4.

X+ X− X+− X− λ5 −2.187 × 10−1 −2.188 × 10−1 −6.965 × 10−12 λ4 1.401 × 10−1 1.993 × 10−2 −6.963 × 10−12 λ3 1.813 × 10−1 1.415 × 10−1 5.397 × 10−04 λ2 2.870 × 10+0 2.870 × 10+0 5.397 × 10−04 λ1 5.516 × 10+0 5.516 × 10+0 1.600 × 10−01

Table 3: Eigenvalues of X± and X+− X− for Example 5.5.

matrix ψ(λ) with λ = −1 is symmetric indefinite. Therefore, by Theorem 2.1, the NME (1) with A and Q defined in (46a) has no symmetric positive definite solution.

The numerical results for X± are reported in Table 1. The eigenvalues of X± and X+− X− are list in

the second, third and fourth columns of Table 2. We see that X± are symmetric indefinite solutions, but

X+− X− is positive semidefinite.

Example 5.5 Let A0 and G0 be defined in (48). Take

V = 4 ⊕  8 2 2 6  ⊕  8 1 1 6  (49) and γ = 1, then ˆG0 is nonsingular and Q is symmetric positive definite. The symplectic pair ( ˆM, ˆL) in (46)

has eigenvalues {3, 1/3} and unimodular eigenvalues {0.6 + 0.8ı, 0.6 − 0.8ı} with partial multiplicity 4. Since ψ(λ) with λ = 0.6 − 0.8ı is symmetric indefinite, such NME has no symmetric positive definite solution.

The numerical results for computing X±are reported in Table 1. The eigenvalues of X±and X+−X−are

list in the second, third and fourth columns of Table 3. We see that X± are symmetric indefinite solutions,

but X+− X− is close to positive semidefinite.

Remark 5.2 (i) (Monotone properties for Qk and Pk) In Examples 5.1–5.3, there exist maximal symmetric

positive definite solutions X+ and Qs are symmetric positive definite. That is, the assumptions in

(14)

0 10 20 100 102 104 106 108 1010 1012 k κ 2 (Q k − P k ) (a) 0 10 20 −0.45 −0.4 −0.35 −0.3 −0.25 −0.2 −0.15 −0.1 −0.05 0 0.05 k λ min (Q k − Q k+1 ) (b) 0 10 20 −1.8 −1.6 −1.4 −1.2 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 k λ min (P k+1 − P k ) (c)

Figure 1: κ2(Qk− Pk), λmin(Qk− Qk+1) and λmin(Pk+1− Pk) for Example 5.5.

Theorem 3.1 are preserved. The convergence behaviors of Examples 5.1–5.3 comply with Theorem 4.2. However in Examples 5.4 and 5.5, they are no symmetric positive definite solutions, the assumptions of Theorem 3.1 are not satisfied, and the globally monotone properties do not hold. However, the numerical results show that Qk and Pk possess a locally monotone property; i.e. Qk and Pk are monotone when

k is large enough. Such numerical monotone properties for Example 5.5 are shown in (b) and (c) of Figure 1. The convergence behaviors of Examples 5.4 and 5.5 comply with Theorem 4.3.

(ii) (Condition numbers of Qk−Pk) From the above numerical results, we know that the sequences {Qk} and

{Pk} converge to X+ and X−, respectively. The absolute values of the eigenvalues for X+−1A are less

than or equal to one and those for X−1

− A are greater than or equal to one. The unimodular eigenvalues

of X−1

+ A and X −1

− A come from the same Jordan block Jω,p. This leads to ill-conditioned matrices

Qk− Pk when k is large. The condition numbers κ2(Qk − Pk) in each iteration for Example 5.5 are

shown in (a) of Figure 1. This behavior coincides with result (iv) of Theorem 4.3.

6

Conclusions

We have proved the convergence for the SDA in Algorithm 3.1, when the eigenvalues of X−1

+ A are inside

or on the unit circle. The results are confirmed by five numerical examples. Our paper fills in an existing gap in the solution of the NME (1), for equations with nonsimple eigenvalues on the unit circle for X+−1A.

(15)

Consequently, our results are more general than those in [6], which considered only semi-simple eigenvalues, and our doubling algorithm is more applicable than that in [6].

Acknowledgments

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

References

[1] 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:55–80, 2005.

[2] 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:767–788, 2004.

[3] J. C. Engwerda. On the existence of a positive definite solution of the matrix equation X +ATX−1A = I.

Lin. Alg. Appl., 194:91–108, 1993.

[4] 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:255–274,

1993.

[5] G. H. Golub and C. F. Van Loan. Matrix Computations, 3rd ed. The Johns Hopkins University Press, 1996.

[6] C. H. Guo. Convergence rate of an iterative method for a nonlinear matrix equation. SIAM J. Matrix Anal. Appl., 23, 2001.

[7] C. H. Guo. Numerical solution of a quadratic eigenvalue problem. Lin. Alg. Appl., 385:391–406, 2004. [8] C. H. Guo and P. Lancaster. Iterative solution of two matrix equations. Math. Comp., 68:1589–1603,

1999.

[9] C. H. Guo and P. Lancaster. Algorithms for hyperbolic quadratic eigenvalue problems. Math. Comp., 74:1777–1791, 2005.

[10] 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:1063–1075, 2005.

(16)

[11] T.-M. Hwang and W.-W. Lin. Structured doubling algorithms for weak Hermitian solutions of algebraic Riccati equations. Technical report, NCTS Preprints in Mathematics, National Tsing Hua University, Hsinchu, Taiwan, 2006-7-009, 2006.

[12] A. J. Laub and K. Meyer. Canonical forms for Symplectic and Hamiltonian matrices. Celestial Mech., 9:213–238, 1974.

[13] W.-W. Lin, V. Mehrmann, and H. Xu. Canonical forms for Hamiltonian and Symplectic matrices and pencils. Lin. Alg. Appl., 302-303:469–533, 1999.

[14] W.-W. Lin and C.-S. Wang. On computing stable Largangian subspaces of Hamiltonian matrices and symplectic pencils. SIAM J. Matrix Anal. Appl., 18(3):590–614, 1997.

[15] 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):26–39, 2006.

[16] V. Mehrmann. The Autonomous Linear Quadratic Control Problem. Springer-Verlag, 1991.

[17] V. Mehrmann. A step toward a unified treatment of continuous and discrete time control problems. Lin. Alg. Appl., 241-243:749–779, 1996.

[18] B. Meini. Efficient computation of the extreme solutions of X + A∗X−1A = Q and X − AX−1A = Q.

Math. Comp., 71:1189–1204, 2002.

[19] 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:211–228, 2003.

[20] S.F. Xu. Numerical methods for the maximal solution of the matrix equation X + ATX−1A = I. Acta

Scientiarum Naturalium Universitatis Pekinensis, 36:29–38, 2000.

[21] X. Zhan. Computing the extremal positive definite soluions of a matrix equation. SIAM J. Sci. Comput., 17:1167–1174, 1996.

數據

Table 3: Eigenvalues of X ± and X + − X − for Example 5.5.
Figure 1: κ 2 (Q k − P k ), λ min (Q k − Q k+1 ) and λ min (P k+1 − P k ) for Example 5.5.

參考文獻

相關文件

In section 4, based on the cases of circular cone eigenvalue optimization problems, we study the corresponding properties of the solutions for p-order cone eigenvalue

volume suppressed mass: (TeV) 2 /M P ∼ 10 −4 eV → mm range can be experimentally tested for any number of extra dimensions - Light U(1) gauge bosons: no derivative couplings. =&gt;

Define instead the imaginary.. potential, magnetic field, lattice…) Dirac-BdG Hamiltonian:. with small, and matrix

incapable to extract any quantities from QCD, nor to tackle the most interesting physics, namely, the spontaneously chiral symmetry breaking and the color confinement.. 

• Formation of massive primordial stars as origin of objects in the early universe. • Supernova explosions might be visible to the most

Theorem 5.6.1 The qd-algorithm converges for irreducible, symmetric positive definite tridiagonal matrices.. It is necessary to show that q i are in

The difference resulted from the co- existence of two kinds of words in Buddhist scriptures a foreign words in which di- syllabic words are dominant, and most of them are the

• elearning pilot scheme (Four True Light Schools): WIFI construction, iPad procurement, elearning school visit and teacher training, English starts the elearning lesson.. 2012 •