Te hi a

(1)

Fa ultyofMe hatroni sandInterdis iplinary EngineeringStudies

A ademy of S ien esof the Cze h Republi ,Prague

InstituteofComputerS ien e

Limiting A ura y of Iterative Methods

Pavel Jiranek

PhDThesis November,2007

(2)

Author: PavelJiranek

Study programme: P2612Ele trote hni sandInformati s

Field of study: 3901V025S ien eEngineering

Department: Fa ultyofMe hatroni s

andInterdis iplinaryEngineeringStudies

Te hni alUniversityofLibere

Halkova6,46117Libere

Cze hRepubli

InstituteofComputerS ien e

A ademyofS ien esoftheCze hRepubli

Pod Vodarenskouvez2,18207Prague8

Cze hRepubli

Supervisor: MiroslavRozloznk

2007PavelJiranek.

TypesetbyL A

T

E X .

(3)

zejsemtutopra ivypra ovalsamostatneauvedljsemveskereprameny,kter y hjsempouzil.

Datum: Podpis:

(4)

(5)

Jak je znamo,zaokrouhlova hyby a nepresneresen vnitrn h uloh maj vliv na numeri ke

hovanitera n hmetodvaritmeti eskone noupresnost;obe nesnizujjeji hry hlostkonver-

gen eaovlivnujkone noupresnostspo tenehoresen.Vpra isezab yvameanal yzoumaximaln

dosazitelnepresnostinekter y hitera n hmetodproresensoustavlinearn halgebrai k y hrov-

ni .

Dizerta eje rozdelenanadve asti.Prvnzni hobsahujeanal yzulimitnpresnostimetodkry-

lovovsk y h podprostoruproresenrozsahl y huloh sedlov y hbodu.Uvazujemedvatypysegre-

govan y h metod: metodu reduk e na S huruv doplneka metodu projek e na nulov y prostor

mimodiagonalnho bloku. Ukazuje se, ze v yber vzor e pro zpetnou substitu ima vliv na ma-

ximalndosazitelnoupresnostpribliznehoresenspo tenehovaritmeti eskone noupresnost.

Druha astjevenovanaanal yzenumeri keho hovannekter y hmetodminimaln hrezidu,ktere

jsoumatemati kyekvivalentnmetodezobe nen y hminimaln hreziduGMRES.Srovnavame

dvahlavnpostupy:jeden,kdepriblizneresenjevypo tenoze soustavshorntrojuhelnkovou

mati , a jeden,kde je priblizneresen upravovano pomo jednodu hehorekurentnho vzor e.

Ukazujese,zev yberbazemavlivnanumeri ke hovanv ysledneimplementa e.Zatm ometody

SimplerGMRESaORTHODIRjsoumenestabilndkyspatnepodmnenostizvolenebaze,baze

zkonstruovanazrezidumuzeb ytdobrepodmnena,jestlizejsounormyrezidudostate nekle-

saj .Tytov ysledkyvedouknoveimplementa i,kterajepodmnenezpetnestabiln,avjistem

smyslui vysvetluj experimentalneoveren yfakt,zemetoda GCR (ORTHOMIN)davav prak-

ti k y haplika hvelmipresneaproxima eresen.

Kl ova slova. Rozsahle linearn soustavy, metody krylovovsk y h podprostoru, ulohy sed-

lovehobodu,metoda reduk e na S huruv doplnek,metoda projek ena nulov y prostormimo-

diagonalnhobloku,metodyminimaln hrezidu,numeri kastabilita,anal yzazaokrouhlova h

hyb.

(6)

(7)

It is known that inexa t solution of inner systems and rounding errors ae t the numeri al

behaviorof iterativemethodsinnitepre isionarithmeti . Inparti ular,theyslowdowntheir

onvergen erateand have anee t on theultimate a ura y of the omputed solution. Here

wefo uson the analysisof the maximum attainablea ura y of several iterativemethods for

solvingsystemsoflinearalgebrai equations.

Thethesisisdividedintotwoparts. TherstpartisdevotedtotheanalysisofKrylovsubspa e

solversappliedtothelarge-s alesaddlepointproblems. Twomainrepresentativesofsegregated

solution approa hesareanalyzed: theS hur omplementredu tionmethod andthenull-spa e

proje tionmethod. Weshowthatthe hoi eoftheba k-substitutionformula an onsiderably

inuen ethemaximumattainablea ura yofapproximatesolutions omputedinnitepre ision

arithmeti .

Inthese ondpartweanalyzenumeri albehaviorofseveral minimumresidualmethods,whi h

aremathemati allyequivalenttotheGMRESmethod. Twomainapproa hesare ompared:the

approa h,whi h omputestheapproximatesolutionfromanuppertriangularsystem, andthe

approa hwheretheapproximatesolutionsareupdatedwithasimplere ursionformula. Weshow

thatadierent hoi eofthebasis ansigni antlyinuen ethenumeri albehaviorofresulting

implementation. WhileSimplerGMRESandORTHODIRarelessstableduetoill- onditioning

of hosenbasis,theresidualbasisremainswell- onditionedwhenwehaveareasonableresidual

normde rease. Theseresults leadto a new implementation, whi h is onditionally ba kward

stable, and in a sense explainan experimentally observed fa t that the GCR (ORTHOMIN)

methoddeliversinpra ti al omputationsverya urateapproximatesolutionswhenit onverges

fastenoughwithoutstagnation.

Keywords. large-s alelinearsystems,Krylovsubspa emethods,saddlepointproblems,S hur

omplement redu tion, null-spa e proje tion method, minimum residual methods, numeri al

stability,roundingerroranalysis.

(8)

(9)

Èçâåñòíî,÷òî íåàêêóðàòíûåðåøåíèÿâíóòðåííèõïðîáëåìè îøèáêèîêðóãëåíèÿ îòðàæà-

þòñÿ íà âû÷èñëèòåëüíîì ïîâåäåíèþ èòåðàöèîííûõ ìåòîäîâ. Îíè êîíêðåòíî çàòîðìîçÿò

èõ ñêîðîñòüñõîäèìîñòèè îêàçûâàþòâëèÿíèåíàèíàëüíóþàêêóðàòíîñòüâû÷èñëåííîãî

ðåøåíèÿ. Ìû çäåñü çàíèìàåìñÿ àíàëèçîììàêñèìàëüíîé äîñòèæèìîéàêêóðàòíîñòè íåêî-

òîðûõèòåðàöèîííûõìåòîäîâäëÿðåøåíèÿñèñòåìëèíåéíûõàëãåáðàè÷åñêèõóðàâíåíèé.

Ýòàäèññåðòàöèÿðàçäåëåíàíàäâå÷àñòè.Ïåðâàÿçàíèìàåòñÿàíàëèçîìëèìèòíîéàêêóðàò-

íîñòèìåòîäîâïðîñòðàíñòâ Êðûëîâà äëÿ ðåøåíèÿáîëüøèõ ñèñòåì ñåäåëüíûõòî÷åê. Ìû

ðàññìàòðèâàåì äâà òèïû ñåãðåãàöèîííûõ ìåòîäîâ: ìåòîäîì ïðåîáðàçîâàíèÿ íà äîïîëíå-

íèåØóðàèìåòîäîìïðîåêöèèíàÿäðîíåäèàãîíàëüíîãîáëîêà.Ìû óêàçûâàåì,÷òîâûáîð

îðìóëû îáðàòíîéïîäñòàíîâêèîòðàæàåòñÿ íàìàêñèìàëüíîé äîñòèæèìîé àêêóðàòíîñòè

ïðèáëèçèòåëüíîãîðåøåíèÿâû÷èñëåííîãîâàðèìåòèêåñêîíå÷íîéòî÷íîñòüþ.

Âòîðàÿ÷àñòüñîäåðæèòàíàëèçâû÷èñëèòåëüíîãîïîâåäåíèÿíåñêîëüêèõìåòîäîâìèíèìàëü-

íûõíåâÿçîê,êîòîðûåìàòåìàòè÷åñêèýêâèâàëåíòíûåìåòîäó¾GMRES¿.Ìûñðàâíèâàåìäâà

ãëàâíûåìåòîäû: îäèí,êîòîðûé îïðåäåëÿåòïðèáëèæ¼ííîåðåøåíèå èçñèñòåìû ñâåðõíåé

òðåóãîëüíîéìàòðèöîé, è îäèí,ãäå ïðèáëèæ¼ííîåðåøåíèå êîððåêòèðîâàííîåñïîìîùüþ

ïðîñòîé ðåêóððåíòíîé îðìóëû. Ìû óêàçûâàåì, ÷òî âûáîð áàçû îòðàæàåòñÿ íà âû÷èñ-

ëèòåëüíîìïîâåäåíèèêîíå÷íîãîìåòîäà.Ïîêà ìåòîäû ¾SimplerGMRES¿è¾ORTHODIR¿

ìåíååñòàáèëüíûåçàñ÷åòïëîõîîáóñëîâëåííîéáàçû,áàçàíåâÿçîêìîæåòáûòüõîðîøîîáó-

ñëîâëåííàÿ,åñëèíîðìûíåâÿçîêäîñòàòî÷íîñíèæàþòñÿ.Ýòèðåçóëüòàòûâåäóòêíîâîìóìå-

òîäó,êîòîðûéóñëîâíîîáðàòíîñòàáèëüíûé,èâîïðåäåëåííîìñìûñëåîáúÿñíÿþòýêñïåðè-

ìåíòàëüíîóäîñòîâåðåííûéàêò,÷òîìåòîä ¾GCR¿(òàêæåèçâåñòíûéêàê¾ORTHOMIN¿)

äà¼òâïðàêòè÷åñêèõàïïëèêàöèÿõî÷åíüàêêóðàòíûåàïïðîêñèìàöèèðåøåíèÿ.

Êëþ÷åâûå ñëîâà. áîëüøèåëèíåéíûå óðàâíåíèÿ,ìåòîäû ïîäïðîñòðàíñòâÊðûëîâà,ìå-

òîäïðåîáðàçîâàíèÿíàäîïîëíåíèåØóðà,ìåòîäïðîåêöèèíàÿäðîíåäèàãîíàëüíîãîáëîêà,

ìåòîäûìèíèìàëüíûõíåâÿçîê,âû÷èñëèòåëüíàÿñòàáèëüíîñòü,àíàëèçîøèáîêîêðóãëåíèÿ.

(10)

(11)

Acknowledgements

Firstof all,I wouldliketo thankMiroslavRozloznkfor beinga patientsupervisor. His ideas

andourfruitfuldis ussionsshowedmethewaywheretogoandinspiredmetoexpressmyown

reativity. I am alsogratefulto him for introdu ingme to the fas inatingworld of numeri al

mathemati sandnumeri allinearalgebra.

My thanksalso go to Martin Gutkne ht, the oauthor of the paper [

62

^℄, ^Mario^Arioli, ^Julien

Langou,andYvanNotayfortheirsuggestionsanddis ussionswhi h onsiderablyimprovedthe

presentedresults.

Finan ial support of my work on the thesis was provided in part by the Ministry of Edu a-

tion, Youth and Sportsof theCze h Republi undertheproje t1M0554\Advan edRemedial

Te hnologies",andbytheproje t1ET400300415withintheNationalProgramofResear h\In-

formationSo iety".

Finally,Iwouldliketothankmyfamily,espe iallymyparentsandsister,fortheirattentionand

supportthroughouttheyearsofmy lifeandstudiesinLibere . Andmostof all,I amgrateful

tomygirlfriendEvawhohasbeenin rediblysupportive,understandinganden ouragingduring

thepastyearsandIthankherforloveandstandingbymeduringthegoodandbadtimes.

(12)

(13)

Introduction

Considerasystemoflinearalgebrai equationsintheform

Ax=b; (1.1)

whereAisanNNnonsingularmatrixandbisaright-handsideve tor. Usuallyweassumethat

Aislargeandsparseasitis,e.g.,whenAisadis reterepresentationofsomepartialdierential

operator. Wearelookingforthesolutionof(1.1)orforitssuÆ ientlya urateapproximation.

The methods for solving (1.1) are usually lassied as dire t and iterative. Dire t methods

aremostlybasedonthesu essiveeliminationofunknowns. Theyfa torizethesystemmatrix

(with suitablyordered rowsor olumns),e.g., into the produ t of lowerand uppertriangular

matri es as in the Gaussianelimination, or to the produ t of an orthogonal and a triangular

matrixas intheQRfa torization. Thesolution of(1.1) anbethenfound by solvingsystems

with these fa tors. In general, dire tmethods are well suitedfor dense and moderately large

systems. However, when solvinga large sparsesystem, the numberof newly reatednon-zero

elementsinbothfa tors anheavilyae tthe omputationaltimeandstoragerequirements. In

addition,even thoughdire tmethods deliverintheory theexa t solution,there is noneedfor

su hana ura yinpra ti eduetoun ertaindataordis retizationerrors.

Therefore, iterative methods be ame very popular when solving sparse systems. An iterative

method forthesolutionof(1.1)generatesa sequen eofapproximationsx

k

sothattheyideally

onverge to the exa t solution. The system matrix need not to be expli itly stored. In ea h

iterationweneedonlytoperformamatrix-ve tormultipli ation. Moreover,theapproximations

onvergeoftenmonotonously(oralmostmonotonously)insomexed normandsowe anstop

theiterationpro esswhentheapproximationisa urateenough. However,the onvergen erate

of iterativemethods an beslow ingeneral (dependingon properties of the system)and thus

hybridte hniques ombiningtheiterativeanddire tapproa h,su haspre onditionediterations,

arewidelyusedto makethepro essmoreeÆ ient.

In general, a solution method (no matter if a dire t or iterative one) an be interpreted as

a solution ofa sequen e of subproblemswhi h aresimpler to solve. In dire tmethods we an

identifyfollowingsubproblems:thefa torizationofthesystemmatrixandthesolutionofsystems

with omputedfa tors. Inea hstepofaniterativemethod,wemultiplyave torbythesystem

matrixandoptionallysolvethesystemwithapre onditionerwhi h anbealsoregardedasthe

subproblemssolvedrepeatedlyintheiterationloop. E.g.,thematrix-ve tormultipli ation an

involvethesolutionofaninnersystemasintheS hur omplementredu tionmethodwhi hwe

willdis usslater.

(16)

1. The state of the art

Fromnowonwerestri tourselvestoiterativemethods.Inpra ti e,the omputationsareae ted

byerrors. Theyareneverperformedexa tlyduetoroundingerrorsandsomeofthemaredone

inexa tly with a pres ribedlevel of a ura y, espe ially when omputations with the working

a ura y ould bea waste of time and resour es. E.g., matrix-ve torprodu ts may involvea

solutionofinnersystems, whi h (beinglargeandsparse) anbesolved inexa tlywith another

iterativemethod. Pre onditioning anbealsoappliedthroughsomeiterativepro ess. Usually,

a method is alled inexa t if some involved subproblems are solved only approximately even

thoughwe assumeexa tarithmeti . Roundingerrors analso onsiderablyae tthe behavior

of iterativemethods. Sin e thebehaviorof inexa t iterativemethods and \exa t"methods in

nitepre isionarithmeti issimilar,wewillnotstri tlydistinguishbetweenthesour esoferrors

andwewilltreatthem ommonlyinauniedapproa hinthefollowingdis ussion.

Whenaninexa tnessistakenintoa ount,thereareseveralimportantquestionswhi hneedtobe

answered. Inthefollowingwegiveabriefoverviewofthestateofartinthiseld(in ludingresults

innite pre isionarithmeti ). Generallytheinexa tnessintrodu edinaniterativemethodhas

twomainee ts:

The errors aused by inexa t omputations are propagatedthroughout the iterative

pro ess. Ideallytheerrorpropagationshouldberestrainedsothatthelo alerrorsare

not magnied. Thereis a limit in the a ura y whi h annot be ex eeded and it is

usually alledthemaximumattainable(orlimiting)a ura y.

The onvergen eofaninexa titerativemethod anbedelayedwithrespe ttothe on-

vergen eofthesamemethod,whereall omputationsareperformedexa tly. Wemay

askhowmanyadditionaliterationsshouldbeperformedsu hthatthesamea ura y

isattainedasintheideal(exa t) ase.

In this thesis we fo us on the limiting a ura y of inexa t iterative methods. The ee ts of

inexa tmatrix-ve tormultipli ationsiniterativemethods(alsoreferredasrelaxedmethods)on

themaximumattainablea ura ywere studiedsimultaneouslybyvan denEshofand Sleijpen

[

97

^℄, ^and ^by ^Simon
ini ^and ^Szyld ^[

90

^℄. ^Their ^analysis ^explains ^the experimental results of BourassandFraysse[

18

^℄^(the^report^with^an^extensiveexperimentalbasiswaspublishedin2000) who proposeda relaxation strategy for the a ura y of the omputed matrix-ve torprodu t.

They have shown that to a hieve the pres ribed a ura y of the omputed solution we need

to ompute the matrix-ve tor produ t with the a ura y (measured by the ba kward error)

inverselyproportionalto thea tualresidualnorm. Thepapers[

97, 90

^℄^provide^thetheoreti al support forthis strategy furtherdevelopedin[

98

^℄. ^This^topi^is ^losely ^related^to ^the^exible

pre onditioning,see,e.g.,[

11, 43, 76, 90, 39

^℄. ^Here^we^try^toâdopt^the^ba
kwardêrrorânalysis,

widelyusedinthestudyofroundingerrors,andweanalyzetheee tsofinexa t omputations

onthe limitinga ura y of ertain iterativemethods. The omputationsareperformed inthe

presen eof roundingerrorswhilesolutionsto ertainsubproblemsaredonewith morerelaxed

a ura y. Wewanttoknowhowtheinexa tnessoftheseinnersystemstogetherwith theerrors

aused by roundo ae t the behavior of the onsidered algorithms. It appears that some

measuresof the a ura y are ultimatelyon the level proportionalto the unit roundo, while

othermeasuresdependonthea ura yofinnersystems.

Theproblemofnumeri alstabilityof lassi aliterativemethodswasaddressedinseveralpapers.

Therstanalyzes arriedoutbyGolub[

40

^℄^and^Lynn^[

69

^℄^providestatisti alandnon-statisti al

(17)

results for the se ond order Ri hardson and SOR method. The statisti al error analysis of

lassi aliterativemethods was alsoperformedbyArioliand Romani[

5

^℄ ^larifying^the^relation

between the onditioningof thepre onditionedsystem matrixandthe onvergen erateofthe

iterativemethod. In[

56

^℄^Highamând^Knight^give^the^forwardând^ba
kwardêrrorânalysisôfâ

generalone-stepstationarymethod. Theiranalysisamongotherthingsshowsthatthea ura y

of the omputed solution strongly depends on the os illationsof norms of the iterates whi h

is a ommon observation not only in the ase of lassi al iterativemethods. Moreover, even

though the onvergen e is driven by the spe tral radius of the iteration matrix, the limiting

a ura y dependsrather on thenormof its powerswhi h anbearbitrarilylargein theearly

stage of the iterative pro ess. This was observed by Hammarling and Wilkinson [

53

^℄. ^The

stability of lassi al iterative methods was also analyzed by Wozniakovski in[

107, 108

^℄. ^He

provedtheforwardstabilityof lassi almethodslikeJa obi,Ri hardson,Gauss-SeidelandSOR

(for symmetri systemswith the PropertyA) and Chebyshevmethod (forsymmetri positive

denite systems). However, the Chebyshev method appeared to be not normwise ba kward

stable. In[

41

^℄ ^Golub^and^Overton^dis
uss^the onvergen erateofthese ondorderRi hardson and Chebyshev method. They onsider the inexa t solution of inner systems with uniformly

bounded relativeresiduals. The a ura yof the omputedsolution inthe Chebyshev method

is further analyzed by Giladi, Goluband Keller [

37

^℄ ^who ^show ^the ôptimality ôf ^the ûniform

toleran eusedin[

41

^℄. ^When^the^systemîs^solved^by^the^lassi
alîterative^methodⁱⁿêa
h^step

wemustsolvea simplersystemindu ed bythesplittingof thesystemmatrix. However,these

systems anbe alsosolved iteratively. Thesemethods,referredto as two-stagemethods, were

addressed,e.g.,in[

73, 64, 36

^℄.

One of the most important result in the study of Krylov subspa e methods is due to Paige

[

77

^℄. ^He ^provides ^the ânalysis ôf ^the ^behavior ôf ^the ^symmetri ^Lan
zos âlgorithm ^[

65

^℄ ⁱⁿ

the presen e of rounding errors. This algorithm is losely related to the onjugate gradient

method by Hestenes and Stiefel [

54

^℄. ^It ^was ^rst ^studied ^by ^Wo^zniakowski ^[

109

^℄ ^and ^Bollen

[

17

^℄. ^Wo^zniakowski ^shows ^that ^this ^method ônverges ⁱⁿ ^nite ^pre
ision ârithmeti ât ^least

linearlywiththe onvergen eratesimilartothesteepestdes entmethod. However,hisanalysis

does not re e t the realityverywell, sin e the onvergen eof the onjugate gradient method

annotbe hara terizedlo allybutitsa tualbehaviordependsonthewholeiterationpro ess;

see, e.g., [

99, 68

^℄ ând ^the ^referen
es^therein. ^The ^new însightînto ^this ^problem^was ^brought

by Greenbaum[

45

^℄ ^and^further^developed ^together^with ^Strako^s^[

95, 49

^℄. ^It ^appears^that ^the

nitepre isionLan zospro essaswellasthenitepre ision onjugategradientmethodbehave

as theirexa t ounterpartsappliedtothematrixof (possiblymu h) largerdimensionwiththe

eigenvalues lusteredneartheeigenvaluesoftheoriginalmatrix. Thisissuewasfurtherdis ussed

byNotayin[

75

^℄.

Theanalysisoflimitinga ura yofsome lassesofiterativemethods anbeperformedinrather

generalsettingwithoutreferringtoanyparti ularmethod. Themethodsbasedonthe oupled

two-term re urren eswere analyzed by Greenbaum in [

46, 47

^℄. ^The ^papers ^fo
us ^mainly ^on

the onjugategradientmethod buttheanalysisholdsforalargersetofmethods.Inparti ular,

theresultsofGreenbaumshowthatthehighlyirregular onvergen ebehavior(expressedbythe

os illationsofnormsofiterates)observedinthe aseofnon-optimaliterativemethods(su has

BiCG[

35

^℄^or^CGS^[

93

^℄)^an^have^anunfavorableee tonthelimitinga ura yofthe omputed solution. A similar phenomenon is mentioned also by van der Vorst in [

100

^℄, ^where^the ^loss

of a ura yis explainedbyos illationsof residual norms. On theother hand,su h os ilations

(18)

do not o ur (or an be a prioribounded) in the ase of optimal methods su h as onjugate

gradientsand onjugateresiduals[

94

^℄^applied^to^symmetri^positive^denite^problems,^orⁱⁿ^the

aseof residual minimizingmethods (Orthodir [

110

^℄, ^Orthomin ^[

102

^℄, ^GCR ^[

29

^℄) ^for ^general

nonsymmetri systems. The numeri al stability of various (equivalent) methods using short

re urren eswasfurtherstudiedbyGutkne htandStrakosin[

52

^℄^and^by^Sleijpen,^van^der^Vorst

andModersitzkiin[

92

^℄. ^In^[

51

^℄^Gutkne
htând^Rozlo^zn^k^dis
uss^theêe
tôf^residual^smoothing

onthelimitinga ura y.

Finallywesurveytheresultsforthenitepre isionbehaviorofnonsymmetri Krylovsubspa e

methods with the full-termre urren es su h as GMRES [

88

^℄. ^The Householder implementa- tionofthe underlyingArnoldipro ess[

6

^℄ ^is^quitestraightforwardto analyze,seethe paperby Drkosova,Greenbaum,RozloznkandStrakos[

27

^℄,ând^byÂrioliând^Fâssino^[

4

^℄. ^This^is^due^to

thealmostexa torthogonalityof the omputedKrylovsubspa ebasis. However,whenweuse

the heapermodiedGram-S hmidtimplementation,theorthogonalityis graduallylostduring

theiteration pro ess. The lossof orthogonalityhowevergoes handinhand with the de rease

of the ba kward error of the a tual omputedsolution as observed by Greenbaum, Rozloznk

andStrakosin[

48

^℄ând^furtherânalyzed^by^Paige,^Rozlo^zn^kând^Strako^sⁱⁿ^[

80, 78

^℄. ^F^or^more

detailssee[

67

^℄^and^the^referen
es^therein.

2. Organization of the thesis

Thisthesisisdividedintotwomainpartsandisorganizedasfollows. Chapter3,whi hisbased

onthe papers[

61, 60

^℄, îs ^devoted ^to ^theânalysis ôf înexa
t^methods ^for^solving ^saddle^point

problemsoftheform

A B

B T

0

x

y

=

f

0

:

AbriefoverviewonsaddlepointproblemsispresentedinChapter2. Weanalyzetwosegregated

methodsbasedonthetransformationofthewholeindeniteproblemtoaredu edsystemwith

morepreferableproperties(smallerdimension,positive(semi)deniteness). Theredu edsystem

issolvedbyasuitableiterativemethodgivingtheapproximationstooneoftheblo k omponents

ofthesolutionve tor(xory). Theremaining omponentis omputedviasomeba k-substitution

formula. We onsiderthreedierentbutmathemati allyequivalentformulas. Inea hiteration

wehaveto solveeitheranonsingularsystemwith A,orafullrankleastsquares problemwith

B. Sin e su h systems are not usually solved exa tly, we assume here that they are solved

witha pres ribedba kwarderrorand studytheee t onthemaximumattainablea ura yof

the solution method together with the ee ts of rounding errors. Su h inexa t methods have

beenalso onsideredin manypapersbutmostof them analyzed thedelayof onvergen e; see

thereferen esinChapter3. Here weprovideaqualitativeanalysisofthemaximumattainable

a ura yofthe omputedsolutionmeasuredbytrueresidualsinthesaddlepointsystem,bytrue

residualsinredu edsystemsand byforwarderrorsofthe omputedsolutions. Inaddition,we

showwhi h residuals(andhow) anbeae ted bythepossiblyirregular onvergen ebehavior

inthe aseof the nonsymmetri blo k A. Thetheoreti al results areillustrated on numeri al

experiments.

Chapter 4, based on the paper [

62

^℄, îs ^devoted ^to ^the ânalysisôf ^several ^residual ^minimizing

Krylov subspa e methods, whi h are mathemati ally equivalentto the GMRES method [

88

^℄.

In ontrast to GMRES, they, in the nth iteration, build an orthonormalbasis of AK

n (A;r

0 )

insteadofK (A;r ): K (A;r )denotesthenthKrylovsubspa egeneratedbythematrixAand

(19)

the ve torr

0

. Twoapproa hesare ompared: theapproa h,whi h omputesthe approximate

solution fromanuppertriangularsystem, and theapproa h,wheretheapproximatesolutions

areupdatedstepbystepwith a simplere ursionformula. We onsidera generalbasisto gen-

erate the orthonormal basis of AK

n (A;r

0

), and it appears that, while Simpler GMRES and

ORTHODIR are less stabledue to ill- onditioningof the hosen basis, the residual basis an

bewell- onditioned,when wehavea reasonableresidual normde rease. Theseresults leadto

a new implementation, whi h is onditionallyba kward stable, and to the well known GCR

(ORTHOMIN)method,andinasenseexplainanexperimentallyobservedfa t thatGCR(OR-

THOMIN) deliversverya urateapproximateapproximatesolutionsin pra ti alappli ations.

Thetheoreti alresultsareillustratedonnumeri alexperiments.

InChapter5wegive on lusionsanddire tionsofthefuturework.

3. List of related publications and conference talks Journal papers.

P.Jiranek,M.Rozloznk. Maximumattainablea ura yofinexa tsaddlepointsolvers.

A epted for publi ationin SIAM Journal on Matrix Analysis and Appli ations,

2007.

P. Jiranek,M.Rozloznk. Limitinga ura yof segregatedsolutionmethodsfor non-

symmetri saddlepointproblems. A eptedfor publi ationinJournal of Computa-

tional andApplied Mathemati s,2007.

P.Jiranek,M.Rozloznk,M.H.Gutkne ht. Howto makeSimplerGMRESandGCR

more stable. Submitted to SIAM Journal on Matrix Analysis and Appli ations,

2007.

Proceedings contributions.

P.Jiranek.Onamaximumattainablea ura yofsomesegregatedte hniquesforsaddle

pointproblems. Pro eedings of the XI. PhD. Conferen e,pages26{34,Institute of

ComputerS ien e,CAS,Matfyzpress,Prague,2006.

P. Jiranek,M.Rozloznk. Ona limitinga ura yof segregatedte hniques forsaddle

point problems, Pro eedings of the 3rd International Workshop on Simulation,

ModellingandNumeri alAnalysisSIMONA2006,pages62{69,Libere ,September

2006.

Conference talks.

P. Jiranek,M. Rozloznk. Numeri al behaviorof iterativemethods for saddle point

problems. GAMMAnnualMeeting2006,Berlin,Mar h27{31,2006.

P. Jiranek. On a maximum attainable a ura y of some segregated te hniques for

saddlepointsolvers. XI.PhD.Conferen e,InstituteofComputerS ien e,A ademyof

S ien esof theCze hRepubli ,Monne {Sedle -Pr i e,September18{20,2006.

P. Jiranek,M.Rozloznk. Ona limitinga ura yof segregatedte hniques forsaddle

pointsolvers. Simulation,ModellingandNumeri alAnalysis SIMONA2006,Libere ,

September18{20,2006.

P.Jiranek,M.Rozloznk. Numeri alsolutionofsaddlepointproblems. SNA'07,Sem-

(20)

P. Jiranek,M.Rozloznk. Onthelimitinga ura yofsegregatedsaddlepointsolvers.

MAT-TRIAD2007{threedaysfullofmatri es,B

,

edlewo,Poland,Mar h22{24,2007.

P. Jiranek,M.Rozloznk. Onthelimitinga ura yofsegregatedsaddlepointsolvers.

VIII. vede ka konferen ias medzinarodnou u ast '

ou, Te hni al University of Kosi e,

Slovakia,May28{30,2007.

P. Jiranek, M. Rozloznk. Limiting a ura y of inexa t saddle point solvers. 22nd

BiennialConferen eonNumeri alAnalysis,UniversityofDundee,S otland,UK,June

26{29,2007.

P. Jiranek, M. Rozloznk, M. H. Gutkne ht. On the stability of Simpler GMRES.

CEMRACS'07,Lumini,Fran e,Juny22{August31,2007.

P. Jiranek,M.Rozloznk,M.H.Gutkne ht. HowtomakeSimplerGMRESandGCR

morestable. IMAConferen eonNumeri alLinearAlgebraandOptimisation,Univer-

sityofBirmingham,UK,September13{15,2007.

3.1. Posters.

P.Jiranek,M.Rozloznk. Numeri alstabilityofinexa tsaddlepointsolvers.ICIAM'07,

6thInternationalCongress on Industrial and Applied Mathemati s, Zuri h,Switzer-

land,July16{20,2007.

(21)

Saddle point problems

Thesolutionoflarge-s alesystemsinthesaddlepointformattra tedalotofattentioninre ent

years. Theyappearinalargevarietyofappli ationsandmanysolutionmethodsweredeveloped

sofar. Thenext hapterisdevotedtothenumeri alstabilityanalysisof ertainiterativemethods

for saddlepointsystemsand, before we start,wegive ashortintrodu tionintothis eld. For

anexhaustiveoverview werefertothepaperbyBenzi,GolubandLiesen[

14

^℄.

We onsiderthelargesparsesystemoflinearalgebrai equationsintheblo kform

A

x

y

A B

B T

C

x

y

=

f

g

; (2.1)

where A 2

R

ⁿⁿ^, ^B ²

R

^nm ^and ^C ²

R

^mm^. ^The ^solution ^and ^right-hand^side ^ve
tors^are

partitioned onsistentlywithrespe ttothepartitioningofthesystemmatrix. LetAandBare

nonzeromatri es andfurthermoreweassumethat theright-handsideis always hosenso that

thesystemis onsistent.

Thepropertiesof blo ksA,B andC mayvary dependingontheappli ation. Inthe following

se tion wementionseveral importantexamples of problemsleading to a saddle pointsystem.

Notethat thesystem(2.1)hasasymmetri blo kstru turewhi h anberelaxedwhen solving

so alledgeneralizedsaddlepointproblems. However,wedonot onsiderthis asehere.

1. Applications leading to saddle point problems

Saddlepointproblemsariseina widesele tionof problemsof omputationals ien eand engi-

neering. WhenA issymmetri positivedenite,B hasa full olumnrankand C=0,wehave

themost ommonversionofthesaddlepointsystem

A B

B T

0

x

y

=

f

g

; (2.2)

whi happears,e.g.,whensolvingellipti se ondorderpartialdierentialequationsbythemixed

nite element method [

24

^℄ ^or ^quadrati programming problems with linear onstraints [

38, 74

^℄. ^The ômponent^x ôf ^the ^solution ^ve
tor^(x;^y)ôf ^(2.2) îs ^the ^solution ôf ^the onstrained minimizationproblem

min

u2

R

ⁿ

J(u)= 1

2 u

T

Au f T

u s.t. B T

u=g: (2.3)

The orrespondingLagrangianisdenedas

L(u;v)=J(u)+(B T

u g) T

v 8u2

R

ⁿ^; ^8v²

R

^m^;

(22)

wherevistheve torofLagrangemultipliers.Theve tor(x;y)isthesaddlepointofL,

L(x;v)L(x;y)L(u;y):

The nonsymmetri blo k A appears, e.g., when solving linearized Navier-Stokesequation via

thesequen eofStokesandOseenproblems. If,inthemixednite elements,theapproximation

spa esdonotfullltheLBB ondition,thestabilizationshouldbeappliedleadingtothenonzero

symetri positivesemidenitematrixC [

24, 32

^℄.

Another important appli ation of saddle point systems is the solution of linear least squares

problems. LetB beannmmatrixofafull olumnrankand onsider

nd y s.t. kf Byk= min

v2

R

^m

kf Bvk:

Itiswell-known[

16, 42

^℄^that^the^solutionôf^this^problemîsûniqueândîtîs hara terizedbythe orthogonality onditionx=f By?R (B)=N(B

T

)

?

forthe residualve torx(where R (B)

andN(B T

)denotes therangeandnull-spa eofthe matrixB and B T

,respe tively). Hen e we

havex+By=f,B T

x=0leadingtothesystem

I B

B T

0

x

y

=

f

0

:

Ingeneral,thesystemofthe form(2.2)(withg =0) orrespondstotheweightedleastsquares

problem, where A 1

-norm is minimized instead of the Eu lidean one (when A is symmetri

positivedenite).

2. Properties of saddle point matrices

Here we brie y re all the basi properties of saddle point matri es and relate their spe tral

and nonsingularityproperties with respe t to the properties of parti ularblo ks. We restri t

ourselvestothesymmetri asebutsomeresults anbeextendedtoamoregeneralsetting. For

amore ompletedis ussion,see [

14

^℄.

Theorem2.1. LetAbe a symmetri positivedenitematrixwith eigenvalues ontained in

theinterval[ ;℄andletB be ofafull olumn rankwithsingularvalues ontained in[;℄

with >0and >0andC is symmetri positive semidenite. Then

Ahas n positiveandm negativeeigenvalues;

ifC=0, the eigenvalues of Aare lo alized asfollows:

(A)I [I +

;

where

I

1

2

q

2

+4

2

; 1

2

q

2

+4

2

;

I +

; 1

2

+ q

2

+4

2

:

Proof. ThesaddlepointmatrixA anbefa torizedasfollows

A=

I 0

T 1

A 0

T 1

I A 1

B

:

(23)

TherststatementimmediatelyfollowsfromtheSylvester'slawofinertia[

57

^℄,^sin
e^the^S
hur

omplement B T

A 1

B Cissymmetri negativedenite. Fortheproofofthese ondstatement,

see[

85

^℄.

Thematrix Ais indenite,sin eit hasbothpositiveandnegativeeigenvalues. Solving highly

indenitematri es(withnm) anleadtotheslow onvergen ewhen usingKrylovsubspa e

methodslikeMINRES[

79

^℄,^see^[

34

^℄. ^A^simplemodi ationofthesystemmatrixintheform

^

A

A B

B T

C

;

asobserved,e.g.,in[

13, 34

^℄, ^leads^to^anonsingularsystemwithaspe trummovedtotheright half-planeofthe omplexplanebut,however,forthepri eoflosingthesymmetry.

Thenonsingularity onditionsaresummarizedinthefollowingtheorem(see[

13

^℄).

Theorem 2.2. Let A be symmetri nonnegative real (that is, 1

2 (A+A

T

) is positive semi-

denite), B has a full olumn rankandlet C be symmetri positive semidenite. Then

if Ais nonsingular,then N(A)\N(B T

)=0;

if N(

1

2 (A+A

T

))\N(B T

)=0, then Ais nonsingular.

Here 0 represents the null subspa e of

R

ⁿ^. ^In parti ular,if A is symmetri positivesemidenite, then Ais nonsingularif andonly ifN(A)\N(B

T

)=0.

3. Solution methods

Solutionmethodsforsystemsoftheform(2.1) anbedividedintotwo ategories alled oupled

and segregated methods. Coupled methods solve the system (2.1) as a whole and therefore

ompute both omponentsx and y of the solution ve tor at on e. They an be bothdire t,

e.g., usingLDL T

fa torizationwith 11and 22pivots, and iterative,e.g.,usingMINRES

[

79

^℄ ⁱⁿ ^the ^symmetri âse. Ôn ^the ôther ^hand, ^segregated ^methods ^transform ^the ^system

(2.1) of the dimension n+m to a redu ed system of a smaller dimensionsolving either for

the omponentx ory. Theremaining omponent isthen found by theba k-substitution into

(2.1). The redu edsystems an bealso solvedeitherdire tly oriteratively. They anbehard

to ompute expli itly, so the iterative approa h is more preferable in many ases. Moreover,

besidesthesmallerdimension,theredu edsystems anbeeasiertosolvethanthewholesaddle

pointsystem (e.g., the redu edsystem anbe positive(semi)denite). Sometimes the border

between oupledandsegregatedapproa hesisnotsharp,sin e oupledmethods anbetreated

assegregatedandvi eversa. Herewereviewtwomainrepresentativesofsegregatedapproa hes

whi h willbe analyzedin thenext hapter: theS hur omplementredu tionmethod and the

null-spa eproje tionmethod. Wewillnotdis ussotherissuesrelatedtothetopi andsolution

methods,espe iallypre onditioningofsaddlepointproblems;see[

14

^℄^for^moreinformation.

3.1. The Schur complement reduction method.

ÂssumeÂîs^symmetri^positive^def-

inite, B hasa full olumn rankand C is symmetri positivesemidenite. Then Theorem2.2

impliesthat the system (2.1) hasa unique solution. It anbe regarded as twomatrix-ve tor

equationsintheform

T

(24)

Sin eAisnonsingular,we antoeliminatexfromtherstequation,i.e.,x anbeexpressedas

x=A 1

(f By); (2.5)

andsubstitutedintothese ondequation. Thenweobtainthesystem

Sy=B T

A 1

f g; S B T

A 1

B+C (2.6)

withtheS hur omplementmatrixS (whi his,morepre isely,thenegativeS hur omplement

ofAinA). Thesolutionofan(n+m)-dimensionalindeniteproblem(2.1)isthustransformed

to the solution of two systems of orders m and n with symmetri positive denite matri es.

First,thesystem(2.6)is solvedfor y. It isnotalwayspreferableto omputeS dire tly, sin e,

eventhoughAissparse,Sneednottobe. Sometimestheeliminationpro ess anbeperformed

su hthatthesparsityispreserved[

71

^℄. ^When^(2.6)^is^solvediteratively,weneedto omputethe produ twithSwhi hinvolvesthesolutionofasystemwiththematrixA. Theiterativemethod

produ esthesequen eofapproximationsy

k

(k=0;1;2;:::) onvergingideallyto y. Whenthe

ve toryoraniteratey

k

isavailable,the orrespondingapproximationtox anbe omputedby

thesubstitutioninto(2.5).

One of the most popular methods for solving saddle point systems based on the S hur om-

plementredu tion is the Uzawa method [

7

^℄. ^The âlgorithmîs âs ^follows: ^hoose ^y0

, then for

k=0;1;2;:::do

(

solve Ax

k +1

=f By

k

;

y

k +1

=y

k

(g B T

x

k +1 +Cy

k ):

Here>0isa relaxationparameter. Hen ewe anwritetheiterationintheform

A 0

B T

1

I

x

k +1

y

k +1

=

0 B

0

1

I C

x

k

y

k

+

f

g

:

Thedire t omputationshowsthat theiterationmatrixoftheasso iatedstationarymethod is

A 0

B T

1

I

1

0 B

0

1

I C

=

0 A

1

B

0 I S

:

Thus the Uzawa method onvergesif and only if the spe tralradiusof I S is stri tlyless

thanone. Itiseasy tosee that theUzawamethod isbasedonthe S hur omplementmethod,

sin eitisnothingbuttheRi hardsoniterationappliedto theS hur omplementsystem(2.6).

Ontheotherhand,theUzawamethod anberegardedasablo kGauss-Seidelmethod(witha

regularizationintheblo k(2;2))appliedtothesaddlepointsystem(2.1).

3.2. The null-space projection method.

^The^S
hur^omplement^redu
tion^relies^on^the

ee tive solutionof systems with the matrix A. Sometimesthe appli ation of A 1

is hard to

omputeinwhi h asethenull-spa eproje tionmethod anbethemethod of hoi e. Assume

herethatAissymmetri positivedeniteonN(B T

),B hasafull olumnrankandC=0. The

system(2.2)isthusbyTheorem2.2uniquelysolvableand anbeexpressedastwomatrix-ve tor

equations

Ax+By=f; B T

x=g: (2.7)

Letx

0

beaparti ularsolutionofthese ondequationandZ2

R

ⁿ⁽ⁿ ^m) ^be^a^matrix^ontaining

a basis of the null-spa e of B T

. Everysu h solution lies in the aÆne spa e x

0

+N(B T

) and

hen e has the form x = x

0 +Zx

Z

, where x

Z

2

R

ⁿ ^m ^are ^the oordinates of x x 0

in the

(25)

null-spa ebasisZ. Substitutionintotherstequationof(2.4)andpremultiplyingbyZ T

gives

thesymmetri positivedenitesystem

Z T

AZx

Z

=Z T

(f Ax

0

) (2.8)

that is, the redu ed system of the ordern mfor the omponents of x x

0

in the basis of

N(B T

). The system Z T

AZ an be solved dire tly oriteratively. Whenwe haveZ expli itly

available(e.g., bythesparseQRfa torization)bothapproa hes anbeapplied. However,when

usinganiterativemethod, it anbeimplementedso that thematrix Z is keptonly impli itly

[

44

^℄. ^Weân^view^the^solutionôf^(2.8)âs ^the^solutionôfâ ^proje
ted^system

(I )A(I )x

1

=(I )f; (2.9)

wherex

1

=Zx

Z

and isthe orthogonalproje toronto R (B). Thesolution omponenty an

bethenfoundviathesolutionoftheleastsquaresproblem

kf Ax Byk= min

v2

R

^m

kf Ax Bvk: (2.10)

When (2.8) or (2.9) is solved iteratively produ ing the sequen e of approximations x

k (k =

0;1;2;:::),solving(2.10)givesanapproximationy

k

toywith xrepla edbyx

k .

(26)

(27)

Limiting accuracy of segregated saddle point solvers

Wewant to solvea saddle pointsystemwhi h is infa t the symmetri indenitesystem with

22blo kstru ture

A B

B T

0

x

y

=

f

0

; (3.1)

wherethediagonalnnblo kAissymmetri positivedeniteandthenmo-diagonalblo k

B hasfull olumnrank. Saddlepointproblemshave re entlyattra ted a lotof attention and

appearto bea time- riti al omponentinthesolutionoflarge-s aleproblemsinmanyappli a-

tionsof omputationals ien eandengineering. Alargeamountofwork hasbeendevotedto a

widesele tionofsolutionte hniquesvaryingfromthefullydire tapproa h,throughtheuseof

iterativestationaryorKrylovsubspa emethods,upto the ombinationof dire tanditerative

te hniquesin ludingpre onditionediteratives hemes. Foranex ellentsurveyonappli ations,

methods, andresults onnumeri al solutionof saddle point problems,we referto [

14

^℄ ^and ^nu-

merousreferen estherein (relevantreferen eswillbegivenlater inthetext). Signi antlyless

attention, however, has been paidso far to the numeri al stability aspe ts. Here we on en-

trate onthe numeri al behaviorof s hemes whi h omputeseparately the unknown ve tors x

and y: oneof them is rstobtained froma redu ed systemof a smaller dimension,and, on e

it has been omputed, the other unknown is obtained by ba k-substitution solvingexa tly or

inexa tly another redu ed problem. The main representativesof su h a segregated approa h

aretheS hur omplementredu tionmethodandthenull-spa eproje tionmethod. Weanalyze

su halgorithmswhi h an beinterpretedas iterationsfortheredu edsystembut omputethe

approximatesolutionsx

k andy

k

tobothunknownve torsxandysimultaneously.

TheS hur omplementredu tionmethodusestheblo kfa torizationintheform

A B

B T

0

=

I 0

B T

A 1

I

A B

0 B

T

A 1

B

;

wherethematrix B T

A 1

B istheS hur omplementof Ain(3.1). Su hde ompositionleads

tosolvingtheresultingblo ktriangularsystem

A B

0 B

T

A 1

B

x

y

=

f

B T

A 1

f

; (3.2)

whi h is nothing buta blo k Gaussian eliminationapplied to the original system (3.1). The

blo ktriangularsystem(3.2)issolvedby omputingtheunknownyfromthesymmetri positive

deniteS hur omplementsystem

B T

A 1

By=B T

A 1

f (3.3)

ofordermandthenby omputingtheunknownxfromasystemofordernwiththesymmetri

positivedenitematrixA. Thisapproa hleadsto theexpli itformulafortheunknownve tor

(28)

x = A 1

(f By). The null-spa e proje tion method is based on the proje tion of the rst

blo kequationin(3.1) onto the null-spa eN(B T

) andonto itsorthogonal omplementR (B),

respe tively. A ordingtothese ondblo kequationof(3.1)theunknownxbelongsto N(B T

)

andthereforewegettheblo ktriangularsystem

(I )A(I ) 0

B T

A B

T

B

x

y

=

(I )f

B T

f

; (3.4)

whereB(B T

B) 1

B T

denotes theorthogonalproje toronto R (B). Thistriangularsystem

is solvedby forward substitution, where we rst ompute the unknown x fromthe proje ted

system

(I )A(I )x=(I )f (3.5)

ofordernwith thesymmetri positivesemi-denitematrix(I )A(I ). On e ithasbeen

omputed,theunknowny isobtainedasy=B y

(f Ax)bysolvingtheleastsquaresproblem

kf Ax Byk= min

v2

R

^m

kf Ax Bvk; (3.6)

whereB y

denotestheMoore{PenrosepseudoinverseofB. Thesu essofalgorithmsforsolving

theblo k triangularsystems(3.2) or(3.4)dependsonthe availabilityof good approximations

to the inverse of the blo k A orto the pseudoinverse of B, respe tively. More pre isely, one

looks for a heap approximatesolution to the innersystems with the matrixA and/or to the

asso iatedleastsquaresproblemswiththematrixB. Numerousinexa ts hemeshavebeenused

andanalyzed,see,e.g.,theanalysisofinexa tUzawaalgorithms[

31, 22, 23, 12, 112

^℄,^inexa
t

null-spa e methods [

89, 105, 111

^℄, ^multilevel ^or ^multigrid ^methods ^[

21, 20, 111

^℄, ^domain

de omposition methods [

19

^℄, ^two-stage ^iterative ^pro
esses ^[

73, 36

^℄ ^or inner-outer iterations [

43

^℄. ^These^papersôntain^mainly ^the ânalysis ôf â onvergen e delay aused by the inexa t solutionofinnersystemsorleastsquaresproblems.

We on entrateonthe questionof what isthe best a ura ywe anget frominexa ts hemes

solving either (3.2) or (3.4) when implemented in nite pre ision arithmeti . The fa t that

the inner solution toleran e strongly in uen es the a ura y of omputed iterates is known

andwas studiedinseveral ontexts. The generalframework forunderstanding inexa tKrylov

subspa emethods hasbeendevelopedin[

90

^℄^and ^[

97

^℄. ^Assuming^exa
t arithmeti ,Simon ini andSzyld[

90

^℄ând ^van^denÊshofând^Sleijpen^[

97

^℄investigatedthe ee tofanapproximately omputedmatrix-ve torprodu tineveryiterationontheultimatea ura yofseveralsolversand

explainedthesu essofrelaxationstrategiesfortheinnera ura ytoleran efrom[

18, 19, 39

^℄.

The developed theory strongly exploits the parti ularproperties of an iterativemethod used

forsolvingtheasso iatedsystem. Inthe ontextofsaddlepointproblems,this requiresa deep

analysisof the outer iteration s heme for solvingthe redu ed S hur omplementor proje ted

system(inparti ular,wereferto[

90

^,^Se
tion^8℄).

The ee ts of roundingerrors in the S hur omplement redu tion method and the null-spa e

proje tionmethod have beenstudied, e.g., in[

2, 3, 26, 70

^℄, ^where^the ^maximum^attainable

a ura yof omputedapproximatesolutions bymeans ofresiduals anderrorsis estimatedde-

pendingonthe usertoleran espe ied intheouter iteration. We analyzethein uen eof the

inexa tsolutionofinnersystems/leastsquaresproblemsonthesamequantities. Ourapproa h

isbasedonastandardba kwardanalysiswhi hallowsustotakeintoa ountboththeinexa t-

nessoftheinneriterationloopsaswellasthea ompanyingroundingerrorsthato urinnite

(29)

The theory developed for the outer iteration pro ess is similar to the analysisof Greenbaum

in [

47, 46

^℄ ^who ^estimated ^the ^gap ^between ^the ^true ^and re ursively updated residual for a general lassofiterativemethodsusing oupledtwo-termre ursions. Thedieren ehereisthat

every omputedapproximatesolutionof innerproblemisinterpretedas anexa t solutionof a

perturbedproblemindu edbythea tualstopping riterion,whilethetheoryof[

47

^℄^onsidered

onlytheroundingerrorsasso iatedwithaxedmatrix-ve tormultipli ation. In ontrasttothe

theoryofinexa tKrylovmethods[

90, 97

^℄,^the^bounds^for^the^true^residualⁱⁿ^the^outer^iteration

loop are obtainedwithout spe ifying the solverused forsolving the S hur omplementorthe

proje tedHessiansystem. It appearsthat themaximumattainablea ura ylevelintheouter

pro ess is mainly given by the inexa tness of solvingthe inner problemsand it is not further

magniedby theasso iatedroundingerrors. Theseresults arethus similarto oneswhi h an

beobtainedinexa tarithmeti .

Thesituationisdierentwhenlookingatthenumeri albehaviorofresidualsasso iatedwiththe

originalsaddlepointsystem,whi h des ribehowa uratelytheblo kequations (3.1)aresatis-

ed. It isshown thattheattainablea ura yof omputedapproximatesolutionsthendepends

signi antlyontheba k-substitutionformulausedfor omputingtheremainingunknowns. Our

results show that, independent of the fa t that the inner systems are solved inexa tly, some

ba k-substitution s hemes leadultimately to residuals on the roundounit level. Indeed,our

results onrmthat dependingwhi h ba k-substitutionformulais usedthe omputediterates

maysatisfyeithertherstorthese ondblo kequationtotheworkinga ura y. Webelievethat

su hresults annotbeobtainedusingtheexa tarithmeti onsiderationsandareofimportan e

inappli ationsrequiringa urateapproximations(seee.g. [

44, 38, 24

^℄). ^On^the^other^hand,^we

agreethat inmanyappli ationsthesaddlepointsystem omesfroma dis retizationof ertain

partialdierentialequationsandmu hlowera ura yissuÆ ient. Inany ase,wegiveatheo-

reti alexplanationforthebehaviorwhi hwasprobablyobservedorisalreadyimpli itlyknown.

However,wehavenotfound anyexpli itreferen esto thisissue. Theimplementationsthatwe

pointoutasoptimalarea tuallythosewhi harewidelyusedandsuggestedinappli ations.

The hapterisorganizedasfollows. Se tions1and2aredevotedtotheroundingerroranalysis

oftheS hur omplementredu tionmethodandthenull-spa eproje tionmethod,respe tively.

Ea hse tionisdividedintovesubse tions. Insubse tions1.1and2.1weanalyzethein uen e

of inexa t solution of inner systems or least squares on the maximum attainable a ura y in

theouteriterationpro essforsolving(3.2)or(3.4),andweestimatetheultimatenormsofthe

trueresiduals B T

A 1

f+B T

A 1

By

k

and (I )f (I )A(I )x

k

. Inthe onsequent

threesubse tionsofSe tions1and2,wegiveboundsfortheultimatenormofthetrueresiduals

f Ax

k By

k

and B T

x

k

. As we will see in subse tions 1.2{1.4 and 2.2{2.4, the limiting

a ura y of these residuals may signi antly dier for various ba k-substitution formulas for

omputing x

k or y

k

, respe tively. Subse tions 1.5 and 2.5 ontain forward analysis with the

bounds forthe errors x x

k

and y y

k

. Throughoutthis hapterour theoreti al resultsare

illustratedonthemodelexampletakenfrom[

83

^℄: ^we^putⁿ⁼^100; ^m⁼^20,^and

A=tridiag(1;4;1)2

R

ⁿⁿ^; ^B⁼^rand(n;^m); ^f ⁼^rand(n;^1):

The spe trum ofA and singularvaluesof B liein theinterval[2:001;5:999℄and [2:173;7:170℄,

respe tively. Therefore the onditioning of A or B does not play an important role in our

(30)

Fordistin tion,wedenotequantities omputedinnitepre isionarithmeti bybars.Weassume

that theusualrules ofa well-designed oating-pointarithmeti hold, anduse o asionallythe

notation () for a omputed result of an expression. The roundo unit is denoted by u. In

parti ular,fora matrix-ve tormultipli ationthe boundk (Ax) Axk O(u)kAkkxk is used

andkxkdenotesthe2-normoftheve torx;forageneralmatrixAwemakeuseofthespe tral

normkAk and the orresponding ondition number (A) = kAk=

min

(A), where

min (A) is

theminimalsingularvalueofA. Forasymmetri positivedenitematrixA,kxk

A

denotesthe

A-normoftheve torx. Finally,weapplytheO-notationwhensuitable.

1. Schur complement reduction method

Inthisse tionwewilldis ussalgorithmswhi h omputesimultaneouslyapproximationsx

k and

y

k

totheunknownsxandy andideallyfullltherstblo kequationof(3.1)

Ax

k +By

k

=f: (3.7)

Ourgoalhereisnottosurveyallexistings hemesbasedon (3.7)buttoanalyzethenumeri al

behaviorofthreeimplementationswhi husedierentba k-substitutionformulasfor omputing

the approximate solution x

k

. More pre isely, without spe ifying any parti ular method, we

assume that we have omputed the approximate solution y

k +1

and the residual ve tor r (y)

k +1

usingthere ursions

y

k +1

=y

k +

k p

(y)

k

; (3.8)

r (y)

k +1

=r (y)

k +

k B

T

A 1

Bp (y)

k

(3.9)

withr (y)

0

= B

T

A 1

(f By

0

). Wewilldistinguishbetweenthefollowingthreemathemati ally

equivalentformulas:

x

k +1

=x

k +

k ( A

1

Bp (y)

k

); (3.10)

x

k +1

=A 1

(f By

k +1

); (3.11)

x

k +1

=x

k +A

1

(f Ax

k By

k +1

): (3.12)

Theresultings hemesaresummarizedinFigure3.1. Theses hemeshavebeenusedandstudied

in the ontext of many appli ations, in luding various lassi al Uzawa algorithms, two-level

pressure orre tion approa h, or inner-outer iteration method for solving (3.1); see, e.g., the

s hemes with (3.10) in [

82, 10

^℄, ^(3.11) ⁱⁿ ^[

31

^℄, ^or ^(3.12) ⁱⁿ ^[

22, 23, 12, 112

^℄, respe tively.

Be ausethesolveswith matrixAinformulas(3.10){(3.12)areexpensive,these systemsarein

pra ti esolvedonlyapproximately. Ouranalysisisbasedontheassumptionthateverysolution

ofasymmetri positivedenitesystemwiththematrixAisrepla edbyanapproximatesolution

produ edbyanarbitrarymethod. Theresultingve toristhen interpretedasanexa tsolution

ofthesystem with thesameright-handside ve torbutwith a perturbedmatrixA+A. We

alwaysrequirethat the relativenormof theperturbationis bounded as kAkkAk, where

represents a ba kwarderror asso iated with the omputedsolution ve tor. We will always

assumethattheperturbationAdoesnotex eed thelimitationgivenbythedistan eofA to

thenearestsingularmatrixandputrestri tionintheform(A)1. Itfollowsthenfromthe

standardperturbationanalysis(see,e.g.,[

55, 16

^℄)^that

k(A+A) 1

A 1

k

(A)

kA 1

k:

(31)

outer iteration

y

0

; solveAx

0

=f By

0

; r (y)

0

= B

T

x

0

for

k=0;1;2;::: y

k +1

=y

k +

k p

(y)

k

inner iteration / back-substitution

solveAp (x)

k

= Bp (y)

k

A

⁾ ^x^{k +1}⁼^x^k⁺^k^p^(x)k

B

⁾ ^solve^Ax^{k +1}⁼^f ^By^{k +1}

C

⁾ ^solve^Auk

=f Ax

k By

k +1

; x

k +1

=x

k +u

k

r (y)

k +1

=r (y)

k

k B

T

p (x)

k

Figure3.1. S hur omplementredu tion: Threedierents hemesfor omput-

ingtheapproximatesolutionx

k +1

( alledinthetexttheupdatedapproximate

solution(A), theapproximatesolution omputedbya dire tsubstitution(B),

andtheapproximatesolution omputedbya orre teddire tsubstitution(C),

respe tively).

Notethat if=O(u),thenwehaveaba kwardstablemethodforsolvingthepositivedenite

system with A. In our numeri al experiments, we solve the systems with A inexa tly using

the onjugategradientmethodorwiththe Choleskyfa torizationas indi atedbythenotation

=O(u).

1.1. The attainable accuracy in the Schur complement system.

^In^this^subse
tion

we look at the ultimate a ura y in the outer iteration pro ess by means of the true resid-

ual B T

A 1

f +B T

A 1

By

k

. It is lear that if we perturb the S hur omplement system

B T

A 1

By = B T

A 1

f to B T

(A+A) 1

By^ = B T

A 1

f, where kAk kAk, then

theresidualasso iatedwithy^ anbeboundedas

k B

T

A 1

f+B T

A 1

Byk^

(A)

1 (A) kA

1

kkBk 2

k^yk: (3.13)

We see from(3.13) that there is a limitationto the a ura y of theresidual obtaineddire tly

fromy^anditsboundisproportionalto . Notethatthese onsiderationsweremadeassuming

exa t arithmeti . The ee ts of rounding errors on the same quantity have been studied by

Greenbaum[

47

^℄,^whoônsideredâ^general^lassôf^methods^for^solving^the^xed^systemôf^linear

equations usingtwo-termre ursionsgivenby(3.8)and (3.9). Using asimilarapproa hwe an

(32)

Theorem3.1. The gapbetween the true residual B T

A 1

f+B T

A 1

By

k

andthe updated

residualr (y)

k

an be bounded as

k B

T

A 1

f+B T

A 1

By

k

r

(y)

k k

[(2k+1)+O(u)℄(A)

1 (A)

kA 1

kkBk(kfk+kBk

Y

k );

where

Y

k

is dened as a maximum norm over all omputed approximate solutions

Y

k

max

i=0;:::;k ky

i k.

Proof. The initialresidual r (y)

0

is omputedas r (y)

0

= (B T

x

0

),where(A+A

0 )x

0

=

(f By

0 ),kA

0

kkAk. Itiseasytoseethatthestatementholdsfork=0. The omputed

approximatesolutiony

k +1

andtheresidualr (y)

k +1 satisfy

y

k +1

=y

k +

k

p

(y)

k +y

k +1

;

ky

k +1

kuky

k

k+(2u+u 2

)k

k

p

(y)

k k;

(3.14)

r

(y)

k +1

=r (y)

k

k B

T

p

(x)

k +r

(y)

k +1

;

kr (y)

k +1

kukr (y)

k

k+O(u)kBkk

k

p

(x)

k k;

(3.15)

wherep (x)

k

is theexa t solutionoftheperturbedsystem

(A+A

k )p

(x)

k

= (Bp (y)

k

); kA

k

kkAk: (3.16)

Multiplying(3.14)by B T

A 1

B,substituting(3.16) intothe re urren e(3.15), and subtra ting

thesetwoequationswegetthere urren e

B T

A 1

f+B T

A 1

By

k +1

r

(y)

k +1

= B

T

A 1

f +B T

A 1

By

k

r

(y)

k

k (B

T

p

(x)

k +B

T

A 1

Bp (y)

k )+B

T

A 1

By

k

r (y)

k :

Thenormoftheve tor

k

p

(y)

k

anbeboundedas k

k

p

(y)

k

kky

k +1 k+ky

k

k+ky

k +1

k. This

boundin ombinationwith(3.14)givesky

k +1

kO(u)

Y

k +1

andk

k

p

(y)

k k3

Y

k +1

whi halso

implies

k

p

(x)

k k

3kA 1

k

1 (A) kBk

Y

k +1

: (3.17)

Using(3.16), the boundonk

k

p

(y)

k

k,and someelementarymanipulation,we anestimatethe

term

k (B

T

p

(x)

k +B

T

A 1

Bp (y)

k )

k

k (B

T

p

(x)

k +B

T

A 1

Bp (y)

k

)kk

k B

T

[(A+A

k )

1

A 1

℄ (Bp (y)

k )k

+k

k B

T

A 1

[ (Bp (y)

k

) Bp (y)

k

℄k

[+O(u)℄(A)

1 (A) kA

1

kkBk 2

Y

k +1 :

Considering (3.15), (3.17), and the indu tion assumption on the gap between B T

A 1

f +

B T

A 1

By

k andr

(y)

k

(similarto theoneusedin[

47

^℄),^we^obtain^the^bound^for^the^error^ve
tor

r (y)

k +1

intheform

kr (y)

k +1 k

O(u)(A)

kA 1

kkBk(kfk+kBk

Y

k +1 )

Te hi a

Acknowledgements

62

Contents

Introduction

1. The state of the art

97

90

18

97, 90

98

11, 43, 76, 90, 39

40

69

5

56

53

107, 108

41

37

41

73, 64, 36

77

65

54

109

17

99, 68

45

95, 49

75

46, 47

35

93

100

94

110

102

29

52

92

51

88

6

27

4

48

80, 78

67

2. Organization of the thesis

61, 60

62

88

3. List of related publications and conference talks Journal papers.

Proceedings contributions.

Conference talks.

3.1. Posters.

Saddle point problems

14

R

R

R

1. Applications leading to saddle point problems

24

38, 74

R

R

R

24, 32

R

16, 42

2. Properties of saddle point matrices

14

57

85

79

34

13, 34

13

R