Linear system solvers Sparse iterative methods

LINEARSYSTEMSOLVERS:SPARSEITERATIVEMETHODSHenkA.VanderVorstMathematicalInstituteUniversityofUtrechtUtrecht,theNetherlandse-mail:[e-mail]1DepartmentofMathematicsUniversityofCaliforniaatLosAngelesLosAngeles,CAe-mail:[e-mail]wewillpresentanoverviewofanumberofrelatediterativemethodsforthesolutionoflinearsystemsofequations.Thesemethodsareso-calledKrylovprojectiontypemethodsandtheyincludepopularmethodsasConjugateGradients,Bi-ConjugateGradients,LSQRandGMRES.Wewillsketchhowthesemethodscanbederivedfromsimplebasiciterationformulas,andhowtheyareinterrelated.Iterativeschemesareusuallyconsideredasanalternativeforthesolu-tionoflinearsparsesystems,likethosearisingin,e.g.,niteelementornitedierenceapproximationof(systemsof)partialdierentialequations.Thestructureoftheoperatorsplaysnoexplicitroleinanyoftheseschemes,andtheoperatormaybegivenevenasaruleorasubroutine.Althoughthesemethodsseemtobealmosttriviallyparallellizableatrstglance,thisissometimesapointofconcernbecauseoftheinnerproductsinvolved.Wewillconsiderthispointinsomedetail.Iterativemethodsareusuallyappliedincombinationwithso-calledpreconditioningoperatorsinordertofurtherimproveconvergenceproperties.Thisaspectwillreceivemoreattentioninaseparatechapterinthesamevolume.1TheworkofthisauthorwaspartiallysupportedbytheNationalScienceFoundationundercontractASC92-01266,theArmyResearchOceundercon-tractsDAAL03-91-G-0150andsubcontactfromU.Tenn.DAAL03-91-C-0047,andONRundercontractONR-N00014-92-J-1890.1IterativemethodsversusdirectmethodsThepresentstateoftheartinnumericalmethodsisthatdirectmethodscanbeusedasblackboxes.Thisisbyfarnotthecaseforiterativemethods,atleastnotifwedonotknowaboutspecicpropertiesofthematrixofthelinearsystemtobesolved.Andeventhenitisnotrivialmattertodecidewhentostoptheiterationprocessandtoobtainareasonableestimateoftheapproximationerrorintheresult.Therefore,westartwithsomeveryglobalanalysisofthecomplexityofGaussianeliminationandaniterativesolutionmethodforamodelproblem.Thismayserveasamotivationforthefurtherstudyofiterativemethods.Gaussianeliminationleadstoll-in,thismakesthemethodoftenexpensive.Usuallylargesparsematricesarerelatedtosomegridornetwork.Ina3Dsituationthisleadstypicallytoabandwidthn23(=m2andm3=n,1=mthegridsize).ThenumberofopsisthentypicallyO(nm4)n213(GolubandVanLoan,1989)(Note:weassumethatthesparsitystructureisnotparticularlyregularsothatspecialcomputationalvariantscannotbeemployed).Ifonehastosolvemanysystemswithdierentright-handsides,thenthecostsforsolvingeachsystemwillvaryliken53.IfweassumethatthematrixissymmetricpositivedenitethentheConjugateGradient(CG)iterationmethod(whichwewilldescribeinmoredetaillateron)canbeusedsafely.TheerrorreductionperiterationstepofCGisp1p+1,with=kAk2kA1k2(GolubandvanLoan,1989).Fordiscretizedsecondorderpde’s,overgridswithgridsize1mwetypicallyseem2.Hence,forsuch3Dproblemswehavethatn23.Foranerrorreductionofwemusthavethat0@11p1+1p1Aj(12p)je2jp;wherejisthenumberofiterationstepstobecarriedout.Hence,itfollowsforjthat)jlog2plog2n13:Ifweassumethenumberofopsperiterationtobefn(fstandsforthenumberofnonzerosperrowofthematrixandtheoverheadperunknownintroducedbytheiterativescheme)thenthenumberofopsnecessarytoachieveanerrorreductionbyisfn43log.Conclusion:Ifwehavetosolveonesystematatime,thenforlargen,orsmallf,ormodest,Iterativemethodsmaybepreferable.Ifwehavetosolvemanysimilarsystemswithdierentright-handside,andifweassumetheirnumbertobesolargethatthecostsforconstructingthedecompositionofAisrelativelysmallpersystem,thenitseemslikelythatfor2Dproblemsdirectmethodsmaybemoreecient,whereasfor3Dproblemsthisisstilldoubtful,sincetheopscountforadirectsolutionmethodvariesliken73,andthenumberofopsfortheiterativesolver(forthemodelsituation)variesliken43.Example:TheabovegivenargumentsarequitenicelyillustratedbyobservationsmadebySimon(Simon,1989).Heexpectsthatbytheendofthiscenturywewillhavetosolverepeatedlylinearproblemswithsome5109unknowns.Forwhathebelievestobeamodelproblematthattime,hehasestimatedtheCPUtimerequiredbythemosteconomicdirectmethod,availableatpresent,as520;040years,providedthatthecomputationcanbecarriedoutataspeedof1TFLOP.Ontheotherhand,heestimatestheCPUtimeforpreconditionedconjugategradients,assumingstillaprocessingspeedof1TFLOPS,as575seconds.Thoughweshouldnottakeitforgrantedthatinparticularaverysparsepreconditioningpartcanbecarriedoutatthathighprocessingspeed(forthedirectsolverthisismorelikely),weseethatthedierencesinCPUtimerequirementsaregigantic,indeed(wewillcometothispointinmoredetail).Alsotherequirementsformemoryspacefortheiterativemethodsaretypicallysmallerbyordersofmagnitude.Thisisoftentheargumentfortheusageofiterativemethodsin2Dsituations,whenopcountsforbothclassesofmethodsaremoreorlesscomparable.Finallyitshouldbenotedthatiterativemethodscanexploitgoodinitialguesses,e.g.,intimedependentproblems.Thepreconditionercanoftenbechosentoadapttothemachinearchitecture.Remarks:Forspecialproblemsspecialmethodsmaybefaster:multigrid,fastpoi