A Structure-Preserving Doubling Algorithm for Continuous-Time Algebraic Riccati Equations

algebrai Ri ati equations E.K.-W. Chu  H.-Y. Fan y W.-W.Lin z Abstra t

Continuous-timealgebrai Ri atiequations(CAREs) anbetransformed, ala Cayley,to

dis rete-timealgebrai Ri atiequations(DAREs).TheeÆ ientstru ture-preservingdoublingalgorithm(SDA)

forDAREs,fromChuetal(2003), anthenbeapplied.Inthispaper,wedevelopthestru ture-preserving

doublingalgorithmfromanewpointofviewandshowitsquadrati onvergen eunderassumptionswhi h

areweakerthanstabilizability anddete tability,aswellaspra ti alissuesinvolvedintheappli ationof

theSDAtoCAREs. A modiedversionoftheSDA,developedfor DAREswitha\doublysymmetri "

stru ture,isalsopresented. Extensivenumeri alresultsshowthatourapproa hiseÆ ientand

ompet-itive.

AMS lassi ation: 93B50,93B52,93C05, 93D15

Keywor ds:Cayleytransform, ontinuous-timealgebrai Ri atiequation,doublingalgorithm,matrixsign

fun tion,stru ture-preserving

1 Introdu tion

Inthis paperweinvestigateastru ture-preservingdoublingalgorithm [24, 37℄ for the omputation of the

symmetri positivesemi-denite (s.p.s.d.) solutionX (i.e. X 0) tothe ontinuous-timealgebrai Ri ati

equation(CARE): XGX+A T X+XA+H=0; (1) where A 2 R nn , X 2 R nn , R 2 R mm

is symmetri positive denite (or s.p.d.; i.e. R > 0), G =

BR 1 B T 0and H=C T C0withB 2R nm andC T 2R np

beingoffull olumnrank.

Equation(1)arisesfrequentlyin solvingthe ontinuous-timelinearoptimal ontrolproblem:

min u J = 1 2 Z 1 0 (x T C T Cx+u T R u)dt subje tto x_ =Ax+Bu: (2)

Theoptimalfeedba k ontrolu  for(2)isgivenby u  = R 1 B T Xx; (4)

whereX isthes.p.s.d.solutiontotheCARE(1). Weassumethat thepair(A;B)isstabilizable(S)(i.e. if

w T B=0andw T A=w T

forsome2C, thenRe()<0orw=0)andthat thepair(A;C)isdete table

(D)(i.e.(A T

;C T

)isstabilizable). Underassumptions(S)and(D),theCARE(1)hasbeenprovedtopossess

auniques.p.s.d.solution[39℄.

Considerthe2n2nHamiltonianmatrixHasso iatedwiththeCARE(1):

H=  A G H A T  (5) 

S hoolofMathemati alS ien es,Building28,MonashUniversity,VIC3800,Australia;eri . hus i.monash.edu .au

y

DepartmentofMathemati s,NationalTsingHuaUniversity,Hsin hu30043,Taiwan;d887206oz.nthu.edu.tw

z

H J= JH T ; J =  0 I n I n 0  withI n

denotingtheidentitymatrixofordern. By(5),theCARE(1) an bewrittenas

H  I X  =  I X  ; (6) where  2R nn

and the spe trum () is on the stable left half planeC . Under assumptions (S) and

(D), the Hamiltonian matrix H has exa tly n eigenvalues on C . If the olumns of [X T 1 ;X T 2 ℄ T span the

stableinvariantsubspa eofH ,thenX

1 isnonsingularandX=X 2 X 1 1

0solvestheCARE(1)(see,e.g.,

[39,44℄).

A numeri ally ba kwardstable algorithm are, proposed by Laub [39℄, omputes X by applying the

QR algorithm with reordering [4, 16, 48℄ to the eigenvalue problem H x = x. Unfortunately, the QR

algorithm preservesneither theHamiltonian stru ture of H northe asso iated splitting of eigenvalues. A

stru ture-preservingalgorithm hasbeenproposed byAmmar andMehrmann [1℄whi h utilizesorthogonal

symple ti transformationsin omputingabasisforthestableinvariantsubspa eofH . Astablesymple ti

orthogonalmethodhasbeensuggestedbyByers[19℄butappliedonlytosystemswithsingleinputoroutput.

Manyiterativemethods havebeensuggestedforsolvingCAREsoverthepast20years. Newton's method

has been applied in extensive literature[28, 31, 38, 41, 47℄. A defe t orre tion method for modifyingan

approximate solution hasalso beenproposed by Mehrmannand Tan [43℄. These methods requireagood

startingapproximatesolution,and anthereforeberegardedasiterativerenementmethods,tobe ombined

withotherdire t methods(see Bunse-Gerstneretal [17,18℄or Mehrmann[41℄fordetails). The

stru ture-preservingmatrixsignfun tion methods(MSFM) [7,11,12,13,14, 20,21,27, 33,46℄havebeenextended

byBarraud[8,9℄andCardinerandLaub [29℄.

A lassofmethods,referredtoasthedoublingalgorithms(DA),hasattra tedmu hinterestsin the70s

and80s(see[2℄andthereferen estherein). Thesemethodsoriginatefromthexed-pointiterationderived

fromthedis rete-timealgebrai Ri atiequation(DARE):

X k +1 = b A T X k (I+ b G X k ) 1 b A+ b H :

Insteadofprodu ing thesequen e fX

k

g,doublingalgorithms produ efX

2

kg. CAREs anbeta kledafter

being transformed to DAREs via the Cayleytransform. However, the onvergen e of the algorithm was

proven only when b

A is nonsingular [2℄, and for ( b

A; b

G ; b

H) whi h is stabilizable and dete table [36℄. DAs

were largelyforgottenin thepastde ade. Re ently,DAshavebeenrevivedfor(periodi )DAREs, be ause

of abetter theoreti al understanding. Stronger onvergen e resultshave been provedfor ( b A; b G; b H) under

weaker assumptions than stabilizability and dete tability [24℄. Superior numeri al results, in omparison

to state-of-the-artmethods ona widerange of test problems, havebeen obtainedbe auseof the stronger

stru ture-preservingpropertiesandthesuperioroperations ount.

Inthispaper,weproposeadoublingalgorithmforCAREs. TheCAREsaretransformedtoDAREs,with

the orrespondingHamiltonianmatrixtransformedinto asymple ti matrixpairbytheCayleytransform.

Ni e onvergen eproperties areinheritedfrom thestru ture-preservingdoublingalgorithm(SDA) [24℄

ap-pliedtothe orrespondingDARE.TheSDApreservesmatrixpairsinSSFwhi hisastrongerpropertythan

symple ti ity. IntheCAREsetting,thematrixsignfun tionmethodspreservetheHamiltonianstru turein

HwhiletheSDApreserves,inea hiterativestep,theasso iatedsymple ti matrixpair( b

N; b

L)inSSF.

Al-thoughunderthein uen eofnumeri alerrors,thematrixpairsthroughtheSDAretaintheirstabilizability,

dete tability aswellas eigenstru tures (withexa tlyhalf ofthespe trumbeingstable; seedetails in [24℄).

Thisstrongerstru ture-preservingpropertyisitsmainstrengthandthereasonofitsa ura y. InSe tion4,

amodiedversionoftheSDA(SDA m)isdeveloped,for\doublysymmetri "DAREs,where ^ A, ^ G= ^ H are

symmetri andpersymmetri . TheSDA mpreservesthesymple ti anddoublysymmetri stru turesofthe

DARE,resultinginbettera ura ythantheSDA. WehaveextensivelytestedtheSDAagainsttheMSFM

and are. Numeri alresultsshowedthatthedoublingalgorithmforCAREsis ompetitiveandpromising.

Finally,itisimportanttostress thatmatrixsignfun tions an beapplied tomoregeneralHamiltonian

matri es in other appli ations, su h as those from H

1

bytheCayleytransform,whi hrequires theestimationof theparameter (seex3below).

2 SDA and matrix sign fun tion method

Inthisse tionweproposeastru ture-preservingdoublingalgorithm(SDA)forsolvingtheCARE(1)based

onthedoublingalgorithmsin[24,37℄. Inaddition,thewell-knownstru ture-preservingmatrixsignfun tion

methods[7,11,12,13,14,20,21,27,33,46℄arealsoreviewedfromthepointofviewofpreservingHamiltonian

stru ture.

LetH bethesetof2n2nHamiltonianmatri es,i.e.,

H=  H H=  A G H A T  ; A;H;G2R nn ; H;G0  : (7)

NotethatifH2HthenH J = JH T

. We alla2n2nmatrixpair(N;L)symple ti ifNJN T

=LJL T

.

LetSbetheset of2n2nsymple ti matrixpairsinthestandardsymple ti form(SSF):

S= ( ( b N; b L) b L= " I b G 0 b A T # ; b N = " b A 0 b H I # ; b A ; b H; b G2R nn ; b G; b H 0 ) : (8)

It iseasily seenthat symple ti ity isweakerthan symple ti ityin SSF. Our proposed algorithm preserves

thestrongerstru tureandgivesriseto betternumeri alperforman e.

WeshallshowhowtheCARE(1),asso iatedwiththe orrespondingHamiltonianmatrix

H  A G H A T  =  A BR 1 B T C T C A T  2R 2n2n ;

anbetransformedto anequivalentDARE.

Byusing the Cayleytransformwith someappropriate >0,theHamiltonian matrixH an be

trans-formed to asymple ti matrixpair (N;L) (H+ I;H I)[41, 42℄. In thefollowing,we onstru tan

equivalen etransformationfrom(N;L)toasymple ti matrixpair( b N; b L)2S. Let A A I;  A A+ I: Startingfrom N =   A G H A T  ; L=  A G H  A T  ;

we hoosea >0su hthatthematri esA

andA +GA T

H arewell- onditioned(seeSe tion3laterfor

details). Totransformthesymple ti matrixpair(N;L)to ( b N; b L)2S,let T 1   A 1 0 HA 1 I  ; T 2   I 0 0 ( HA 1 G A T ) 1  ; T 3   I A 1 G 0 I  :

Simple al ulationsprodu e

b N = " b A 0 b H I # =T 3 T 2 T 1 N =T 3 T 2  A 1  A A 1 G HA 1  A H HA 1 G A T  =T 3  A 1  A A 1 G ( HA 1 G A T ) 1 (HA 1  A H) I  and b L= " I b G 0 b A T # =T 3 T 2 T 1 L=T 3 T 2  I A 1 G 0 HA 1 G  A T 

=T 3 I A 1 G 0 ( HA 1 G A T ) 1 ( HA 1 G  A T ) ; where b A=(  A +GA T H)(A +GA T H) 1 ; b G= A 1 G+A 1 G(A T +HA 1 G) 1 (  A T +HA 1 G); b H =(A T +HA 1 G) 1 (HA 1  A H): Notethat L 1 N = b L 1 b N. Sin e  A =A

+2 I,itfollowsthat

b A=I+2 (A +GA T H) 1 ; (9) b G=2 A 1 G(A T +HA 1 G) 1 ; (10) b H =2 (A T +HA 1 G) 1 HA 1 : (11)

Thenweobtainthedesiredsymple ti matrixpairinSSF, i.e.,

( b N; b L) " b A 0 b H I # ; " I b G 0 b A T # ! 2S; where b A, b Gand b

H aregivenby(9){(11). TheDARE asso iatedwiththesymple ti matrixpair( b N; b L) in SSFis X = b A T X(I+ b GX) 1 b A+ b H (12)

on whi h the eÆ ient SDA [24℄ an be applied. Note that X is the unique s.p.s.d. solutionto the above

DARE as well as theCARE (1). Moreover, in Theorems 1and 2of [37℄, thepairs ( b A ; b B) and ( b A; b C) are

proventobestabilizable anddete table,respe tively,wherethematri es b G= b B b B T and b H = b C T b C arefull rankde ompositions(FRD).

Using (9){(11)to transformtheCARE(1) toan equivalentDARE (12)with theasso iatedsymple ti

matrix pair ( b

N; b

L ) in SSF, the SDA in [24℄ an then be modied to the followingalgorithm for CAREs:

(withImdenotingtheimaginaryaxis)

Stru ture-preserving doubling algorithm (SDA):

Input: H=  A G H A T  2H with(H )\Im=;;

Output: thestabilizingsolutionX=X T

0to theCARE(1).

Findanappropriatevalue^ >0.

Compute b A 0 I+2^ (A ^ +GA T ^ H) 1 , b G 0 2^ A 1 ^ G(A T ^ +HA 1 ^ G) 1 , b H 0 2^ (A T ^ +HA 1 ^ G) 1 HA 1 ^ ,j 0; Dountil onvergen e: Compute b A j+1 b A j (I+ b G j b H j ) 1 b A j , b G j+1 b G j + b A j b G j (I+ b H j b G j ) 1 b A T j , b H j+1 b H j + b A T j (I+ b H j b G j ) 1 b H j b A j ,j j+1; Ifk b H j b H j 1 kk b H j k,Stop; End SetX b H j .

Let b N = " b A 0 b H I # , b L= " I b G 0 b A T # , where b G= b G T , b H = b H T . Suppose b N  b L hasnoeigenvalueson

theunit ir leandthere existnonsingularQ,Z su h that

Q b NZ=  J s 0 0 I  ; Q b LZ=  I 0 0 J s  (13)

where thespe trum(J

s )2O

s

f: jj<1g. Inthefollowingwequotethe onvergen eresultsforthe

SDAalgorithm from[24,25℄.

Theorem 2.1. [24℄ Let b N = " b A 0 b H I # and b L = " I b G 0 b A T # , where b G = b G T , b H = b H T . Suppose b N  b

Lhas noeigenvaluesontheunit ir le andthereexistnonsingularQ,Z su hthat(13) holds. Denote

Z=  Z 1 Z 3 Z 2 Z 4  ,Z i 2R nn for i=1;2;3;4. IfZ 1 andZ 4

areinvertible, thenthesequen esf b A j ; b H j ; b G j g

omputedbythe SDA algorithm satisfy

(i) k b A j k=O(kJ 2 j s k)!0as j!1, (ii) b H j

!X, where X solvesthe DARE (12):

X = b A T X(I+ b GX) 1 b A+ b H; (iii) b G j

!Y, whereY solvesthe dualDARE

Y = b AY(I+ b HY) 1 b A T + b G: (14)

Moreover,the onvergen eratein(i){(iii)aboveisO  j n j 2 j  , wherej 1 j


j n j<1<j n j 1 


 j 1 j 1 with  i ;  1 i

being theeigenvalues of b

N 

b

L(in luding 0and1).

The following Lemma proves that the stabilizability and dete tability properties are preserved by the

SDAthroughout itsiterativepro ess.

Lemma 2.2. [24℄ The stabilizability of ( b A ; b B) implies that ( b A j ; b B j ) is stabilizable, where b G j = b B j b B T j 0 is a FRD of b G j

for ea h j  1. The dete tability of ( b A ; b C) implies that ( b A j ; b C j ) is dete table, where b H j = b C T j b C j 0is aFRDof b H j for ea h j1. Theorem 2.3. [24℄Let b N = " b A 0 b H I # and b L= " I b G 0 b A T #

, wherethe matri es b G= b B b B T 0(FRD) and b H = b C T b C0(FRD).Assumethat ( b A; b B)isstabilizable and ( b A; b

C)is dete table. Then thesequen es

f b A j ; b H j ; b G j

g omputedbythe SDA satisfy(i), (ii),(iii) asinTheorem 2.1.

Remark.Theorem2.1dire tlyproves,undertheassumptionsthat b

N 

b

Lhavenounitmoduloeigenvalues

and Z

1 ;Z

4

are invertible,that thesequen esf b A j ; b H j ; b G j

ggenerated bytheSDA onvergeto zeroandthe

uniques.p.s.d.solutionsoftheDAREsin (12)and(14),respe tively. Lemma2.2 showsthepreservationof

stabilizabilityanddete tabilityoftheiterates( b A j ; b G j ; b H j

)generatedbytheSDA.Furthermore,inTheorem

2.3,weseethattheassumptionsinTheorem2.2areweakerthan onditions(S)and(D).Thisdistin tionof

preservingthesymple ti stru tureinSSF, aswellas thedieren einoperation ounts,areresponsiblefor

thesuperiorperforman eoftheSDA.

On theother hand,foragivenH= 

A G

H A

T 

2H with(H )\Im=;,thematrixsign fun tion

mainMSFMalgorithm isdes ribedas follows. Othermodiedversions anbefoundin [5,8,9,22,29℄and

referen estherein.

Matrix signfun tion algorithm: [7,11,12,13,14,20,21,27, 33,46℄

Input: H=  A G H A T  2H with(H )\Im=;;.

Output: thestabilizingsolutionX=X T 0to theCARE(1). LetH 0 H ,j 0. Dountil onvergen e: ComputeH j+1 1 2 (H j +H 1 j ), j j+1; IfkH j H j 1 kkH j k,Stop; End sgn(H ) H j ; Solve(I sgn(H ))  X 1 X 2  =0; ComputeX X 2 X 1 1 . Remarks:

(i) Noti ethat  X 1 X 2 

spansthestableinvariantsubspa eoftheH .

(ii) Both theSDAandthematrixsignfun tionalgorithm require 32

3 n

3

opsforea h iterativestep.

(iii) WhenworkingwiththeHamiltonianmatrixH ,amoreeÆ ientandstru ture-preservingversionofthe

lassi almatrixsignfun tion iteration an bederivedby workingonly withsymmetri matri esJH .

Detailsmaybe onsultedin[21,41℄.

3 Pra ti al implementation of SDA

Sele tion of ^

Here we rst derived the forwarderror bounds of matri es b A 0  b A, b G 0  b G and b H 0  b H given in (9){

(11), respe tively. A ording to these forward errors, we an design a numeri al s heme to determine an

appropriatevalue^ >0. Inthefollowingroundoanalysis,weuse fl(
)todenote omputed oatingpoint

values. Thequantityuistheunitroundo(orma hinepre ision),whi histypi allyoforder10 8

or10 16

in singleand doublepre ision omputerarithmeti , respe tively. WhenA and B are mnreal matri es,

thematrixB :=jAj if b

ij =ja

ij

j foralli;j, andA B ifa

ij b

ij

foralli;j. The1-,1- andFrobenius

matrixnormsaredenoted byk
k

1 ,k
k 1 andk
k F , respe tively.

We assume that the LU fa torizations of A

and W  A +GA T

H are omputed by Gaussian

elimination with partial pivoting (GEPP).We write these omputed LUfa torsas L

A , U A , L W and U W ,

respe tively. Re allthat

A +
A =L A U A ; j
A j n jL A jjU A j; (15) W +
W =L W U W ; j
W j n jL W jjU W j; (16)

n fl(W 1 )=W 1 +E 1 ; jE 1 j n ujW 1 jjL W jjU W jjfl(W 1 )j; (17) where n

isamodest onstant. From(17),theforwarderrorbound inevaluating b Ain (9)is fl( b A)= b A+E 2 ; jE 2 j4 n ujW 1 jjL W jjU W jjfl(W 1 )j+uj b Aj+O(u 2 ): (18)

Furthermore,from(15),wehave

b G fl(2 A 1 G)=2 A 1 G+E 3 ; jE 3 j2 n ujA 1 jjL A jjU A jj b G j; (19)

hen etheforwarderrorbound inevaluating b Gin(10)is fl( b G)= b G+E 4 ; jE 4 j2 n ujA 1 jjL A jjU A jj b G jjW 1 j T + n ujfl( b G)jjU T W jjL T W jjW 1 j T : (20)

Finally,from(15),wehave

b H fl(2 HA 1 )=2 HA 1 +E 5 ; jE 5 j2 n uj b H jjL A jjU A jjA 1 j; (21)

andtheforwarderrorboundin evaluating b H in (11)is fl( b H)= b H+E 6 ; jE 6 j2 n ujW 1 j T j b H jjL A jjU A jjA 1 j+ n ujW 1 j T jU T W jjL T W jjfl( b H)j: (22)

ForGEPP,wehavein pra ti ekjL

A jjU A jk 1 kA k 1 andkjL W jjU W jk 1 kW k 1 ,and itfollows

from(18),(20)and(22)that

kfl( b A) b Ak 1 .4 n u  1 (W )kfl(W 1 )k 1 +uk b Ak 1 +O(u 2 ); (23) kfl( b G) b Gk 1 .2 n u  1 (A )kW 1 k 1 k b G k 1 + n u 1 (W )kfl( b G)k 1 ; (24) kfl( b H) b Hk 1 2 n u  1 (A )kW 1 k 1 k b H k 1 + n u 1 (W )kfl( b H)k 1 ; (25) where 1 (W )kW k 1 kW 1 k 1 , 1 (W )kW k 1 kW 1 k 1 and 1 (A )kA k 1 kA 1 k 1 .

In order to ontrol the forwarderror bounds of b

A , b

G and b

H, we onsider the following min-max

opti-mizationproblem, todetermineanoptimalvalue ^>0:

min >0 F( ) max i=1;2;3 ff i ( )g; (26) where f 1 ( ):=  1 (W ),f 2 ( ):=  1 (A )and f 3 ( ):= 1 (W

). Sin ethe onditionnumbers

1 (W ),  1 (A )and 1 (W

)approa h1as !1, itfollowsthatF( )be omesunboundedas !1. Extensive

numeri alexperiments onrandomly generatedmatri es indi ate that F( ) is astri tly onvexfun tion in

theneighborhoodoftheoptimal ^wheretheglobalminimumofF( )o urs. For illustration,wereporta

sample ofgraphsof f 1 ( ), f 2 ( ), f 3

( )and F( )in Figures 1and 2. FromTheorem 2.1, we knowthat if

approa hes0and1,thesymple ti matrixpair( b

N; b

L)haseigenvalues losetoone,leadingtoveryslow

onvergen eoftheSDA.This an beavoidedthroughthemin-maxoptimizationproblem(26).

We anapplytheFibona isear hmethodto omputeanapproximatevalueof ,^ see,e.g.,[10,p.272℄.

Ourexperien eindi atesthatthreetoveiterationsofFibona isear hareadequatetoobtainasuboptimal

0

5

10

15

20

25

30

35

40

45

10

0

10

1

10

2

10

3

r

f1(r)

f2(r)

f3(r)

Figure1: Thegraphsoffun tionsf

1 ,f 2 andf 3 .

10

0

10

1

10

2

10

3

10

0

10

1

10

2

10

3

10

4

r

F(r)

Symmetry of G

0

and H

0

If the matri es G and H are of low-rank, say G =gg T  0and H = h T h 0, then so are b G 0 and b H 0 .

Indeed,byusingtheSherman-Morrison-Woodburyformula(SMWF)twi e,it anbeseenthat

b G 0 =2 A 1 G(A T +HA 1 G) 1 =2 A 1 gg T (A T +hh T A 1 gg T ) 1 =2  A 1 gg T (A T A T hh T (I+A 1 gg T A T hh T ) 1 A 1 gg T A T )  =2 n (A 1 g) h I (A 1 g) T hh T I+(A 1 g)(A 1 g) T hh T
1 (A 1 g) i (A 1 g) T o =2 n (A 1 g) I+(A 1 g) T hh T (A 1 g)
1 (A 1 g) T o =2 n (A 1 g) K T g K g
1 (A 1 g) T o (Choleskyde omposition) =2 (A 1 gK 1 g )(A 1 gK 1 g ) T :

Similarly,byapplyingthesamete hniques,wealsohave

b H 0 =2 (A T +HA 1 G) 1 HA 1 =2 (A T +h T hA 1 gg T ) 1 h T hA 1 =2  A T A T h T hA 1 (I+gg T A T h T hA 1 ) 1 gg T A T  h T hA 1 =2 n (hA 1 ) T h I (hA 1 ) I+gg T A T h T hA 1
1 gg T (hA 1 ) T i (hA 1 ) o =2 n (hA 1 ) T I+(hA 1 )gg T (hA 1 ) T
1 (hA 1 ) o =2 n (hA 1 ) T K h K T h
1 (hA 1 ) o (Choleskyde omposition) =2 (K 1 h hA 1 ) T (K 1 h hA 1 ): Computation of b A j , b G j and b H j

We now propose a stru tured and eÆ ient pro edure for the omputation of b A j , b G j and b H j in theSDA

algorithm, respe tively, where b G 0 = b B 0 b B T 0  0, b H 0 = b C T 0 b C 0

 0 are FRDs. For j = 0;1;2;:::, we let

W j (I+ b G j b H j ) 1

. Itis easilyseenthat b H j W j =W T j b H j and b G j W T j =W j b G j

ares.p.s.d.for ea hj 1.

BytheSMWFwe an derivetheformulae

W j =(I+ b G j b H j ) 1 =I b B j (I+ b B T j b H j b B j ) 1 b B T j b H j ; (27) b G j W T j = b G j b G j b C T j (I+ b C j b G j b C T j ) 1 b C j b G j = b B j (I+ b B T j b H j b B j ) 1 b B T j ; (28) W T j b H j = b H j b H j b B j (I+ b B T j b H j b B j ) 1 b B T j b H j = b C T j (I+ b C j b G j b C T j ) 1 b C j : (29)

Whenthematri esBandCstartwithlowranksin(1),we animprovetheeÆ ien yofour omputation

further bythe following ompression pro ess. Computethe Choleskyde ompositionof the s.p.d.matri es

W G;j (I+ b B T j b H j b B j )=K T B;j K B;j andW H ;j (I+ b C j b G j b C T j )=K C;j K T C;j ,respe tively. Forj=0;1;2;:::, appli ationof(27){(29)leadsto b A j+1 = b A 2 j b A j b B j (I+ b B T j b H j b B j ) 1 b B T j b H j b A j ; (30) b G j+1 = b G j + b A j b B j (I+ b B T j b H j b B j ) 1 b B T j b A T j = h b B j ; b A j b B j K 1 B;j i " b B T j K T B;j b B T j b A T j #  b B j+1 b B T j+1 0 (FRD) (31)

b H j+1 = b H j + b A T j b C T j (I+ b C j b G j b C T j ) 1 b C j b A j = h b C T j ; b A T j b C T j K T C;j i " b C j K 1 C;j b C j b A j #  b C T j+1 b C j+1 0 (FRD) ; (32) where b B j+1 and b C T j+1

arethe full olumn rank ompressions of h b B j ; b A j b B j K 1 B;j i and h b C T j ; b A T j b C T j K T C;j i ,

respe tively. Ingeneral, rank( b B j+1 )>rank( b B j )and rank( b C j+1 )>rank( b C j

), andthe ompression pro ess

be omesunprotablewhentheranksof b B j+1 and b C j+1 approa hn.

Error analysis of SDA

We onsidertheforwarderrorboundsofthe omputedmatri es b A j+1 , b G j+1 and b H j+1

intheSDAalgorithm

foroneiterativestepj. Sin eK

B;j andK

C;j

arethe omputedCholeskyfa torsofmatri esW

G;j andW

H ;j ,

respe tively,itfollowsthat

b K B fl(K T B;j b B T j )=K T B;j b B T j +
E 1 ; j
E 1 j 1 ujK T B;j jjK T B;j jj b K B j; (33) and b K C fl(K 1 C;j b C j )=K 1 C;j b C j +
e E 1 ; j
e E 1 j~ 1 ujK 1 C;j jjK C;j jj b K C j; (34) where 1 and ~ 1

aremodest onstants. Therefore,wehave

fl(K T B;j b B T j b A T j )=fl( b K B b A T j )=K T B;j b B T j b A T j +
E 2 ; j
E 2 j 2 ujK T B;j jjK T B;j jj b K B jj b A T j j; (35) and fl(K 1 C;j b C j b A j )=fl( b K C b A j )=K 1 C;j b C j b A j +
e E 2 ; j
e E 2 j~ 2 ujK 1 C;j jjK C;j jj b K C jj b A j j; (36) where 2 and ~ 2

aremodest onstants.

Ifrank( b

B T

j

)=`,thenfromTheorem18.4of[32℄and(31),thereexistanorthogonalmatrixQ

B 2R

2`2`

anda omputeduppertriangularmatrixfl( b

B T

j+1

)withfullrowrank,su hthat

" b B T j fl( b K B b A T j ) # + 
B 1
B 2  =Q B  fl( b B T j+1 ) 0  ; (37) wherej
B j j 3 uG ` (j b B T j j+jK T B;j jj b B T j jj b A T j j)forj=1;2,with 3

beingamodest onstantandkG

` k F = 1 2 .

From(35)and(37),wehave

" b B T j K T B;j b B T j b A T j # + 
B 1
e B 2  =Q B  fl( b B T j+1 ) 0  ; (38) wherej
e B 2 j 2 ujK T B;j jjK T B;j jj b K B jj b A T j j+ 3 uG ` (j b B T j j+jK T B;j jj b B T j jj b A T j

j). Pre-multiplyingbothsides

of(38)byQ T B , itfollowsthat  b B T j+1 0  +Q T B 
B 1
e B 2  =  fl( b B T j+1 ) 0  ; (39)

andwededu ethat

kfl( b B T j+1 ) b B T j+1 k F  4 uk b B j k F + 5 u s (K T B;j )k b K b k F k b A j k F ; (40)

where

4 and

5

are modest onstants,and

s (K B;j )kjK B;j jjK B;j jk 1

istheSkeel onditionnumberof

K T

B;j

. Furthermore,applyingasimilarargumentwiththehelpof(36),we anderivethat

kfl( b C j+1 ) b C j+1 k F  6 uk b C j k F + 7 u s (K C;j )k b K C k F k b A j k F ; (41) where 6 and 7

aremodest onstants.

Ontheotherhand,itfollowsfrom(33)that

fl(K T B;j b B T j b H j b A j )=K T B;j b B T j b H j b A j +
E 3 ; j
E 3 j 8 ujK T B;j jjK T B;j jj b K B jj b H j jj b A j j; (42) where 8

isamodest onstant. From(35)and(42),theforwarderrorbound of omputing b A j+1 is kfl( b A j+1 ) b A j+1 k F  9 uk b A j k 2 F + 10 u s (K T B;j )k b B j k 2 F k b K 1 B;j k 2 F k b H j k F k b A j k 2 F ; (43) where 9 and 10

aremodest onstants.

WhentheSkeel onditionnumbers

s (K T B;j )and s (K C;j

)in (40)and(41)areboundedfrom aboveby

a eptablenumbers,thea umulationoferrorwillbedampenedbythefastrateof onvergen eatthenal

stageoftheiterativepro ess. Danger,ifany,liesintheearlystageofthepro essbeforethe 2

j

n

onvergen e

fa tordominates. It isunlikelyto haveanyill-ee t,asthe a umulatederrorin thematrixadditionsand

multipli ationsshouldbeofmagnitudearoundasmallmultipleofthema hine a ura y.

As the SSF propertiesare preserved in the SDA, any errorwill be astru tured one, only pushingthe

iterationtowardsasolutionofaneighboringSSFsystem. Thusthealgorithmisstableinthissense,whenthe

errorsarenottoolargeandwhenstabilizabilityand dete tabilityaremaintained. Forlargejs,as b A j !0, b G j and b H j

onvergeto theunique s.p.s.d.solutionsof (14)and (12),respe tively. Danger againwill only

omesattheinitialstage oftheiteration. Corresponding he ksmaybeprudentinthealgorithm.

4 SDA m

A matrix A is persymmetri when A is symmetri with respe t to the main anti-diagonal ([30, p. 193℄).

When the DARE transformed from the CARE (1) has the additional property that the initial data b A 0 , b G 0 = b H 0 2R 2`2`

aresymmetri andpersymmetri ,theadditionalstru ture anbepreservedinamodied

versionoftheSDA (SDAm). For simpli ity,we onsider onlywhen =1. This doublysymmetri typeof

DAREsappearintheExamples10and17ofSe tion5(originallyfrom [15℄).

For onvenien e,in theSDA,wedenoteforj=1;2;

A b A j ; G b G j = b H j ; A +  b A j+1 ; G +  b G j+1 = b H j+1 : (44)

Sin eA,G=H aresymmetri andpersymmetri ofevenorder,wewrite

A=  a 1 a 2  a 2 a 1   ; G=  g 1 g 2  g 2 g 1   ; (45) where a 1 , a 2 , g 1 and g 2 2R ``

are symmetri and  =[e

` ;


;e

1

℄ withe

j

beingthejth olumn ofI

` . In

theSDA,weshallshowthat b

A, b

Gand b

H arealsosymmetri andpersymmetri with b G= b H,with A + = A(I+G 2 ) 1 A=  b a 1 b a 2  ba 2 ba 1   ; (46) G + = G+AG(I+G 2 ) 1 A T =  b g 1 b g 2  bg 2 bg 1   : (47) Let q 1 g 1 +g 2 ; q 2 g 1 g 2 ;  1 a 1 +a 2 ;  2 a 1 a 2 : (48)

ba 1 = 1 2   1 (I+q 2 1 ) 1  1 + 2 (I+q 2 2 ) 1  2  ; (49) ba 2 = 1 2   1 (I+q 2 1 ) 1  1  2 (I+q 2 2 ) 1  2  ; (50) b g 1 = g 1 + 1 2   1 (I+q 2 1 ) 1 q 1  1 + 2 (I+q 2 2 ) 1 q 2  2  ; (51) b g 2 = g 2 + 1 2   1 (I+q 2 1 ) 1 q 1  1  2 (I+q 2 2 ) 1 q 2  2  : (52) Furthermore let q 1 =[U 1 ;V 1 ℄   1 0 0 1  U T 1 V T 1  ; q 2 =[U 2 ;V 2 ℄   2 0 0 2  U T 2 V T 2  (53)

bethespe tralde ompositionsofq

1 andq 2 ,respe tively,with 1 , 1 , 2 and 2

beingnonnegativediagonal

matri es. Thenba 1 ,ba 2 ,bg 1 andbg 2

in(49){(52) anbe omputedbythefollowingsymmetri forms:

 1   1 U 1 (I+ 2 1 ) 1 U T 1  1  1 V 1 (I+ 2 1 ) 1 V T 1  1 ;  2   2 U 2 (I+ 2 2 ) 1 U T 2  2  2 V 2 (I+ 2 2 ) 1 V T 2  2 ; b a 1 = 1 2 f 1 + 2 g; ba 2 = 1 2 f 1  2 g; (54)  1   1 U 1 (I+ 2 1 ) 1  1 U T 1  1  1 V 1 (I+ 2 1 ) 1 1 V T 1  1 ;  2   2 U 2 (I+ 2 2 ) 1  2 U T 2  2  2 V 2 (I+ 2 2 ) 1 2 V T 2  2 ; b g 1 = g 1 + 1 2 f 1 + 2 g; bg 2 =g 2 + 1 2 f 1  2 g: (55)

The SDA m omputes b

A , b

G in (46) and (47) using the symmetri forms (54) and (55) and onsiderably

improvesthea ura yofExamples10and17in thenextse tion.

5 Numeri al examples

FortheTablesinthefollowingexamples,dataforvariousmethodsarelistsin olumnswithobviousheadings.

Theheading \ are"is forthe are ommand in MATLAB [40℄, \MSGM" is forthe matrixsign fun tion

method [21℄, and\SDA" (or\SDA m")stands forour SDA (or SDA m) algorithm. There is noiteration

numberstoreportfor areandan`'intheTablesindi atesafailureof onvergen einobtainingasolution.

Inthegraphs,\ratio are"and\ratioMSGM"aretheratiooftheCPU-timesfor areandMSGMtothatof

theSDA,respe tively. Forthe omparisonofresiduals,the\normalized"residual(NRes)formulaisapplied

inthenumeri alexamples,i.e.,

NRes kA T ~ X+ ~ XA T ~ XG ~ XA+Hk kA T ~ Xk+k ~ XA T k+k ~ XG ~ Xk+kHk ; where ~

X isanapproximatesolutionandk
kdenotesthe2-normformatri es.

Some numeri al examplesfrom [15℄ involved verylarge data sets, whi h havenot been repeated here.

Twentyexampleswerepresentedin[23℄. Weretainthenumberingofexamplesin[23℄, ommentuponallof

thembut presentonlyverepresentativeonesin thispaper.

IntheMSFM, thes aling strategysuggestedin [21℄wasimplemented. For afairer omparison, similar

onvergen e riteriawereusedinallthemethodsandthesolutionswerenotrened.

All omputations were performed using MATLAB/Version 6.0 on a Compaq/DS20 workstation. The

ma hinepre isionis2:2210 16

isa9th-order ontinuousstatespa emodelofatubular ammoniarea tor. Thea tualsystemmatri esare A= 2 6 6 6 6 6 6 6 6 6 6 6 6 4 4:019 5:12 0 0 2:082 0 0 0 0:87 0:346 0:986 0 0 2:34 0 0 0 0:97 7:909 15:407 4:096 0 6:45 0 0 0 2:68 21:816 35:606 0:339 3:87 17:8 0 0 0 7:39 60:196 98:188 7:907 0:34 53:008 0 0 0 20:4 0 0 0 0 94:0 147:2 0 53:2 0 0 0 0 0 0 94:0 147:2 0 0 0 0 0 0 0 12:8 0 31:6 0 0 0 0 0 12:8 0 0 18:8 31:6 3 7 7 7 7 7 7 7 7 7 7 7 7 5 ; B T = 2 4 0:010 0:003 0:009 0:024 0:068 0 0 0 0 0:011 0:021 0:059 0:162 0:445 0 0 0 0 0:151 0 0 0 0 0 0 0 0 3 5 ; H =I 9 ; R=I 3 :

Thenumeri alresultsaregiveninTable1.

SDA MSGM are NRes 1:6810 15 1:7310 13 4:6410 14 Iter.no. 9 8

-Table1: ResultsforExample5.

Example 6. The example is identi al to Example 6 of [15℄, whi h has been presented originally in [26℄.

This ontrolproblemforaJ-100jetengineisaspe ial aseofamultivariableservome hanismproblem. To

savespa e,weshallnotlistthesystemmatri eshere. Wereportthenumeri alresultsinTable2.

SDA MSGM are NRes 5:7810 13 3:1110 8 1:9110 12 Iter.no. 10 9

-Table2: ResultsforExample6.

Example 10. Theexample isidenti al to Example 10 of[15℄, whi h has been presentedoriginally in [6℄.

Here,thesystemmatri esare

A=  "+1 1 1 "+1  ; G=I 2 ; H =  " 2 0 0 " 2  :

Theexa tstabilizingsolutionX isgivenby

x 11 =x 22 = 1 2 h 2("+1)+ p 2("+1) 2 +2+ p 2" i ; x 12 =x 21 =x 11 =[x 11 ("+1)℄ :

The orrespondingDARE isdoublysymmetri andtheSDA mwasapplied (seedetailsin Se tion 4). The

numeri alresultswith"=1;10 3

;10 5

and10 7

aregiveninTable3.

Example11. Theexampleis identi alto Example11of [15℄,whi h hasbeenpresentedoriginally in [35℄.

Thisexamplerepresentsanalgebrai Ri atiequationarisingfromaH

1

- ontrolproblem[49℄. Let

A=  3 " 1 4 2 "  ; B =  1 1  ; R=1; H =  4" 11 2" 5 2" 5 2" 2  :

"=1 NRes 0:0010 0 0:0010 0 4:6910 16 9:3610 17 Rel. err. 1:9610 16 1:9610 16 8:8010 16 3:8310 16 Iter.no. 4 4 2 -"=10 3 NRes 1:5810 14 1:4310 16 1:1110 13 9:2010 17 Rel. err. 1:8210 11 2:2210 16 2:2210 13 4:0810 16 Iter.no. 16 13 12 -"=10 5 NRes 2:2810 12 1:1110 16 1:0710 11 5:5310 17 Rel. err. 7:1610 7 1:7610 16 2:1410 11 2:6010 16 Iter.no. 22 19 18 -"=10 7 NRes 1:4910 10 1:3210 16 3:3110 9 2:0610 17 Rel. err. 6:0410 8 4:4410 16 6:6310 9 1:3610 16 Iter.no. 12 26 20

-Table3: ResultsforExample10.

Thematrix X =  2 1 1 1 

is the stabilizingsolutionfor " >0. For "= 0,the solutionX is obtainedby an H-invariant Lagrangian

subspa e,i.e.,asolutioninthesenseofH

1

- ontrol.Thenumeri alresultswith"=1;0aregiveninTable4.

SDA MSGM are "=1 NRes 0:0010 0 1:6910 16 1:9710 16 Rel. err. 1:2610 16 1:2510 15 9:6810 16 Iter.no. 5 2 -"=0 NRes 3:0610 16  5:0610 17 Rel. err. 2:6610 9  7:6810 9 Iter.no. 28 

-Table4: ResultsforExample11.

Example12. Theexampleis identi alto Example12of [15℄,whi h hasbeenpresentedoriginally in [34℄.

Let V =I 2 3 vv T ; v T =  1 1 1  ; A 0 ="diag(1;2;3); H 0 =diag(" 1 ;1;"); wehave A=VA 0 V; G=" 1 I 3 ; H=VH 0 V: Thesolutionis X =V diag(x 1 ;x 2 ;x 3 )V where x 1 =" 2 + p " 4 +1; x 2 =2" 2 + p 4" 4 +"; x 3 =3" 2 + p 9" 4 +" 2 :

Thenumeri alresultswith"=1;10 6

aregivenin Table5.

Example 15. Theexample isidenti al to Example15 of[15℄, whi h hasbeen presentedoriginally in [39,

Example 4℄ and [3℄. This example arises from a mathemati al model of position and velo ity ontrol for

"=1 NRes 2:0110 16 1:7810 15 3:0010 16 Rel. err. 4:3310 16 2:7810 15 5:0310 16 Iter.no. 6 4 -"=10 6 NRes 1:6210 15 2:2210 4 2:1910 15 Rel. err. 2:5810 15 6:3310 4 4:9210 15 Iter.no. 11 10

-Table5: ResultsforExample12.

n=2N 1, thenumberof ontrol inputswill bem=N, and thenumberof outputswill bep=N 1,

respe tively. The omparisonofnormalizedresidualsarereportedin Table6forN =5;20;60;100;140and

180. Figure3reportsthe omparisonofCPUtimesfor are,MSGNandtheSDA.

SDA MSGM are N =5 NRes 1:6110 16 8:7510 15 2:5310 15 Iter.no. 5 6 -N =20 NRes 3:8510 16 3:5510 14 6:1510 15 Iter.no. 5 6 -N =60 NRes 1:5310 15 2:3210 13 8:1410 15 Iter.no. 7 8 -N =100 NRes 2:1510 15 6:6210 13 2:5510 14 Iter.no. 8 9 -N =140 NRes 3:0510 15 6:5010 12 3:6010 14 Iter.no. 8 9 -N =180 NRes 1:2510 14 4:6410 12 2:0110 13 Iter.no. 9 9

-Table6: ComparisonofnormalizedresidualsforExample15.

Example 17. Theexample isidenti al to Example17 of[15℄, whi h hasbeen presentedoriginally in [39,

Example6℄. Thesystemmatri esare

A= 2 6 6 6 6 6 6 4 0 1 0


0 . . . . . . . . . . . . . . . . . . . . . . . . 0 0





0 1 0





0 0 3 7 7 7 7 7 7 5 ; B = 2 6 6 6 6 6 6 4 0 . . . . . . 0 1 3 7 7 7 7 7 7 5 ; R=r; C T = p q 2 6 6 6 6 6 6 4 1 0 . . . . . . 0 3 7 7 7 7 7 7 5 :

It isknownfrom [39℄ thatx

1n =

p

qr. Therefore,wemay usetherelativeerrorin x

1n ,i.e., RE (jx 1n p qrj)= p

qr,asanindi atorofthea ura yoftheresults. The orrespondingDAREisdoublysymmetri and

theSDAm was applied (seedetails in Se tion 4). Table7reports the omparison ofnormalizedresiduals

omputed by SDA, SDAm and arefor n = 6;12;18;24;30. We also report the omparison of relative

errorsin x

1n

omputedbyabovethreealgorithmsinTable8.

Comments on numeri al results

Wehavetestedtwentyexamplesin[23℄toillustratethea ura yandeÆ ien yoftheSDAappliedtoCAREs,

0

50

100

150

200

0

50

100

150

200

250

N

CPU time (in sec.)

care

MSGM

SDA

0

50

100

150

200

1

2

3

4

5

6

7

8

9

10

11

N

The ratio of CPU times w.r.t. t

SDA

ratio_care

ratio_MSGM

Figure3: ComparisonofCPUtimesforExample15.

n NResSDA NResSDA m NRes MSGM NRes are

q;r=1 6 4:5010 15 3:5610 16 8:8710 15 1:8010 14 12 3:6310 10 3:2210 14 9:6810 12 1:2310 11 18 9:4710 5 1:8310 11 4:6310 9 9:4610 9 24 2:4710 2 2:3410 8 9:8810 6 3:2510 7 30 4:8010 1 3:5210 5 3:4710 2 7:1710 4 q;r=100 6 2:5910 15 2:8210 16 1:2010 11 1:0210 15 12 4:8110 10 2:9410 14 4:1110 9 1:5810 11 18 4:3310 5 2:2610 11 1:7810 6 7:8310 9 24 7:2410 1 2:9010 8 1:3710 2 1:5010 5 30 3:0710 1 1:4510 5 2:9410 1 4:3810 3

Table7: ComparisonofnormalizedresidualsforExample17.

n RESDA RESDAm REMSGM RE are

q;r=1 6 1:9410 14 1:1110 15 2:2210 14 5:6810 14 12 3:9910 10 1:6810 13 4:9210 11 5:6110 11 18 8:7310 5 6:3710 11 2:3510 8 2:6810 8 24 3:0910 1 6:3910 8 3:2910 5 6:1610 6 30 6:4910 1 1:5710 4 1:4010 1 8:3710 3 q;r=100 6 2:1310 15 9:9510 16 3:2710 11 7:6710 15 12 6:0410 10 1:8310 13 2:3010 8 6:1310 11 18 5:8710 4 1:1610 10 2:9510 5 2:7010 8 24 2:0210 1 1:3210 7 2:3910 2 5:2010 5 30 4:6010 1 5:6710 5 2:9810 1 1:7110 2

Table8: Comparisonofrelativeerrorsin x

1n

retainingtheoldlabellingoftheexamples:

(1) Comparingwith areforalltheexamples,solutionswithbetteror omparablea ura ywereobtained

usingtheSDAin farlesstime. This omparisonhasbeendiÆ ultas areyieldsnoiterationnumbers

andtheCPUtimeinformationfrom MATLABisnotalwaysa urate.

(2) The best indi ation of the eÆ ien y of the SDA over are omes from Example 15 (with varying

dimensionn),where arerequiredtwotoeighttimesmoreCPUtimesthantheSDA.Thisis onsistent

withthendingsin[24℄forDAREs. KeepinmindthattheSDArequiresfarlessnumberof opsthan

are in ea h iteration, as theoperations in theSDA areperformedin < nn

whereasthose for are

are arriedoutin< 2n2n

.

(3) For exampleswith varying onditioning, su h as Examples 9-14,17 and 18, the SDA out-performed

areand onvergestomorea uratesolutionsin lesstime. Fortheill- onditionedExample20, are

failed whiletheSDAsu eededwithoutdiÆ ulty.

(4) InExample11(inH

1

ontrol),someeigenvalueswerenumeri allyontheimaginaryaxisand

assump-tions in thetheorywerepra ti allyviolated. The strongerstru ture-preservingpropertyoftheSDA

enabledittoprodu eana uratesolutionwhentheMSFMfailed. Somehow, areprodu edaslightly

lessa uratesolutionusingmu hmoreCPUtime.

(5) InExamples10and17,theCAREsgaverisetoDAREswhi hwere\doublysymmetri "(seeSe tion4

fordetails). TheSDA mimprovedtheeÆ ien yoftheSDAfortheseexamples,obtaining omparable

a ura yforExample10whileout-performing areforExample17.

(6) ComparingtotheMSFMforill- onditionedproblems,theSDAperformedbetterintermsofa ura ies

or number of iterations. This is onsistent with the fa t that while both the SDA and MSFM are

stru ture-preserving,the formerpreservesmorestru turethan thelatter. For somewell- onditioned

problems, theeÆ ien y anda ura yoftheSDAandMSFMare omparable. For afewsimplesmall

examples,theMSFM onvergedqui klyandwassuperiortotheSDA.Notethattheworkinvolvedin

aniterationforeithermethodissimilar.

(7) TheMSFM, withsimilaroperations ountto SDA, was generallymoreeÆ ientthan are,espe ially

for well- onditioned problems. For ill- onditioned problems (su h as Example 10), the MSFM was

sometimes lessa uratethan are.

6 Con lusions

SolvingCAREsasDAREs,afterapplyingtheCayleytransform,haspreviouslybeeninvestigatedbymany.

Re entdevelopmentsand better understanding ofdoublingalgorithms, espe iallythe stru ture-preserving

properties and eÆ ien y of the SDA [24℄, givethis old approa ha newlease of life. Inaddition, wehave

studied how the parameter in the Cayley transform an be hosen optimally. A Fibona i sear h for

hoosing was suggested in Se tion 3, together with the details of other issues involved in the pra ti al

implementation of the SDA. We have also developed the SDA m whi h preserves the stru ture of some

doublysymmetri DAREs. Extensivenumeri alresultsshowthatthisapproa hofsolvingCAREsusingthe

SDAis eÆ ientand ompetitive,espe iallyforill- onditionedproblems.

A knowledgements

Wewouldliketo thankProfessorC.-S. Wang (NationalChengKungUniversity, Tainan)formany helpful

[1℄ G. Ammar and V. Mehrmann, On Hamiltonian andsymple ti Hessenberg forms, Lin. Alg. Appl.,

149(1991),pp.55{72.

[2℄ B.D. O. Anderson, Se ond-order onvergent algorithms for the steady-state Ri ati equation, Int. J.

Control,28(1978),pp.295{306.

[3℄ M. Athans, W. Levine, and A. Levis, A system for the optimal and suboptimal position and velo ity

ontrol for a string of high-speed vehi les, in Pro . 5th Int. Analogue Computation Meetings,

Lausanne,Switzerland,1967.

[4℄ Z.BaiandJ.Demmel,Onswapping diagonal blo ksinreal S hurform,Lin.Alg. Appl.,186(1993),

pp.73{95.

[5℄ Z.BaiandJ.Demmel,Usingthematrixsignfun tionto omputeinvariantsubspa es,SIAMJ.Matrix

Analy. Appl.,19(1998),pp.205{225.

[6℄ Z.BaiandQ.Qian,Inversefreeparallelmethodforthenumeri alsolutionofalgebrai Ri atiequations,

inPro . FifthSIAMConf.Appl.Lin.Alg.,Snowbird,UT,June1994,J.G.Lewis,ed.,SIAM,

Philadelphia,PA, 1994,pp.167{171.

[7℄ L. Balzer,A elerated onvergen e of the matrixsign fun tion,Int. J. Control, 21(1980),pp. 1057{

1078.

[8℄ A. Y. Barraud, Investigation autour de la fon tion signe d'une matri e, appli ation a l'equation de

Ri ati, R.A.I.R.O.Automatique,13(1979),pp.335{368.

[9℄ A. Y. Barraud,Produitetoileet fon tion signe de matri e. appli ation al'equationde Ri ati dansle

asdis ret,R.A.I.R.O.Automatique,14(1980),pp.55{85.

[10℄ M.S.Bazaraa,H.D.Sheraii,andC.M.Shetty,NonlinearProgramming,JohnWiley&Sons,1993.

[11℄ A.N.BeaversandE.D.Denman,Asymptoti solutionstothematrixRi atiequation,Mathemati al

Bios ien es,20(1974),pp.339{344.

[12℄ A.N.BeaversandE.D.Denman,A omputationalmethodfor eigenvaluesandeigenve torsofamatrix

withreal eigenvalues,Numer. Math.,21(1974),pp.389{396.

[13℄ A.N. Beaversand E. D.Denman, Anew similarity transformationmethodfor eigenvalues and

eigen-ve tors,Mathemati alBios ien es,21(1974),pp.143{169.

[14℄ A.N.BeaversandE.D.Denman,Anewsolutionmethodformatrixquadrati equations,Mathemati al

Bios ien es,20(1974),pp.135{143.

[15℄ P. Benner, A. J. Laub, and V. Mehrmann, A olle tion of ben hmark examples for the

nu-meri al solution of algebrai Ri ati equations I: Continuous-time ase, Te h. Rep. SPC 9522,

Fakultat fur Mathematik, TU Chemnitz{Zwi kau, 09107 Chemnitz, FRG, 1995. Available from

http://www.tu- hemnitz.de/sfb393/sp 95pr.html.

[16℄ J.H.Brandts,Matlab odeforsortingrealS hurforms,Num.Lin.Alg.Appl.,9(2002),pp.249{261.

[17℄ A.Bunse-Gerstner,R.Byers,andV.Mehrmann,A hartofnumeri almethodsforstru turedeigenvalue

problems,SIAM J.Matrix Analy. Appl.,13(1992),pp.419{453.

[18℄ A.Bunse-Gerstner,V. Mehrmann,andWatkins,An SRalgorithm for Hamiltonianmatri es, basedon

Gaussianelimination,Methods ofOperationsResear h,58(1989),pp.15{26.

tational and Combinatorial Methods in System Theory, C. Byrnes and A. Lindquist, eds.,

North-Holland,1986, pp.185{200.

[21℄ R. Byers, Solving the algebrai Ri ati equation with the matrix sign fun tion, Lin. Alg. Appl., 85

(1987),pp.267{279.

[22℄ R.Byers,C.He,andV.Mehrmann,the matrixsignfun tion methodandthe omputationof invariant

subspa es, SIAM J.Matrix Analy. Appl.,18(1997),pp.615{632.

[23℄ E. K.-W. Chu, H.-Y. Fan, and W.-W. Lin, A stru ture-preservingdoubling algorithm for

ontinuous-timealgebrai Ri atiequations,preprint2002-28,NCTS,NationalTsingHuaUniversity,Hsin hu300,

Taiwan,2003.

[24℄ E. K.-W. Chu, H.-Y. Fan, W.-W. Lin, and C.-S. Wang, Stru ture-preserving algorithms for periodi

des rete-timealgebrai Ri ati equations. ToappearinInt. J.Control, 2004.

[25℄ E. K.-W. Chu, H.-Y. Fan, W.-W. Lin, and C.-S. Wang, A stru ture-preserving doubling algorithm

for periodi des rete-time algebrai Ri ati equations, preprint 2002-18, NCTS, National Tsing Hua

University,Hsin hu 300,Taiwan,2003.

[26℄ E. J. Davison and W. Gesing, The systemati design of ontrol systems for the multivariable

ser-vome henism problem, in Alternatives for Linear Multivariable Control, M. K. Sain and J. L.

Pe zkowsky,eds.,Nat.Eng.ConsortiumIn .,Chi ago,IL,1978.

[27℄ E. Denman and R. Beavers, The matrix sign fun tion and omputations in systems, Appl. Math.

Comp.,2(1976),pp.63{94.

[28℄ L. Die i, Some numeri al onsiderations and Newton's method revisited for solving algebrai Ri ati

equations,IEEE Trans.Auto.Control,36(1991),pp.608{616.

[29℄ J.GardinerandA.J.Laub,Ageneralizationofthematrix-sign-fun tionsolutiontothealgebrai Ri ati

equations,Int.J.Control, 44(1986),pp.823{832.

[30℄ G.H. Golub andC. F.Van Loan,Matrix Computations, 3rd ed., TheJohnsHopkins University

Press,1996.

[31℄ S. Hammarling, Newton's method for solving the algebrai Ri ati equation, NPL Rep. DITC 12/82,

Nat.Phys.Lab., Teddington,MiddlesexTW110LW, U.K.,1982.

[32℄ N. J. Higham, A ura y and Stability of Numeri al Algorithms, So iety for Industrial and

AppliedMathemati s,1996.

[33℄ J.L. Howland,The signmatrixandtheseparation ofmatrixeigenvalues,Lin.Alg. Appl.,49(1983),

pp.221{332.

[34℄ P. Hr. Petkov, N. D. Christov, and M. M. Konstantinov, On the numeri al properties of the S hur

approa h forsolving the matrixRi atiequation, Sys.Control Lett., 9(1987),pp.197{201.

[35℄ G.D. Ian ules u, J.Ly, A. J. Laub, andP. M. Papadopoulos,Spa e station freedom solar array H

1

ontrol. Talkat31stIEEEConf. onDe isionandControl,Tu son,AZ,De .1992.

[36℄ M.Kimura,Convergen eof thedoublingalgorithmfor thedis rete-timealgebrai Ri atiequation,Int.

J.Syst. S i.,19(1988),pp.701{711.

[37℄ M.Kimura,Doubling algorithm for ontinuous-timealgebrai Ri ati equation, Int. J. Syst. S i.,20

(1989),pp.191{202.

[38℄ D.Kleinman,Onaniterativete hniquefor Ri atiequation omputations,IEEE Trans.Auto.

(1979),pp.913{921.

[40℄ MathWorks,MATLAB user's guide (forUNIX Workstations),TheMathWorks, In .,1992.

[41℄ V.Mehrmann,The AutonomousLinear Quadrati ControlProblem,Springer-Verlag,1991.

[42℄ V.Mehrmann,Asteptowardauniedtreatmentof ontinuousanddis retetime ontrolproblems,Lin.

Alg.Appl.,241-243(1996),pp.749{779.

[43℄ V. Mehrmann and E. Tan, Defe t orre tion methods for the solution of algebrai Ri ati equations,

IEEE Trans. Auto.Control,AC-33(1988),pp.695{698.

[44℄ C.C.PaigeandC.F.VanLoan,A S hur de omposition for Hamiltonianmatri es, Lin.Alg. Appl.,

41(1981),pp.11{32.

[45℄ L.Patnaik,N.Viswanadham,andI.Sarma,Computer ontrolalgorithmsforatubularammoniarea tor,

IEEE Trans. Auto.Control,25(1980),pp.642{651.

[46℄ J.Roberts,Linear model redu tion andsolution ofthe algebrai Ri ati equationby theuse ofthe sign

fun tion,Int.J.Control, 32(1980),pp.667{687.

[47℄ N.Sandell, On Newton's method for Ri ati equation solution, IEEE Trans. Auto. Control,AC-19

(1974),pp.254{255.

[48℄ G.W. Stewart,HQR3andEXCHNG:Fortran subroutinesfor al ulating andorderingthe eigenvalues

ofareal upper Hessenbergmatrix,ACM Trans. Math.Software, 2(1976),pp.275{280.

[49℄ K.Zhou, J.C. Doyle,and K. Glover,Robust and Optimal Control, Prenti e-Hall, Upper Saddle