A Bluffers ’ Guide to Iterative Methods for System

ABluffers’GuidetoIterativeMethodsforSystemsofLinearEquationsVictorEijkhoutMarch20031IntroductionIterativemethodsareapopularwayofsolvinglinearsystemsofequationAx=bwhereAissquareandnon-singular.Therearemanyiterativemethods,andthenthereisthemysterioustopicof‘preconditioning’.Thisarticlegoesunabashedlyforthebreadthratherthanthedepthindescribingthisuniverse.Therearenoproofsofanysort.1.1WheredolinearsystemscomefromThestatementthatlinearsystemscanbesolvedwithiterativemethodsisabitdeceptive:noteverylinearsystemcansuccessfullybesolvebyaniterativemethod.Seethenextsectionforafewremarksaboutdirectmethods,whichhaveamorecatholictasteinlinearsystems.Thelinearsystemsthataremostamenabletosolutionbyiterativemethodsarethosethatcomefromdiscretisationofpartialdifferentialequations,inparticularthosefrommathe-maticalphysics.Hereareafewofthepossiblesourcesoflinearsystems:•Ellipticequations.Thesealwaysgivealinearsystemtosolve.Onecharacteristic1ofellipticsystemsisthatallpointsofthesolutionare‘globallycoupled’:anypointdependsonthesolutioninanyotherpoint.Thishastheimmediateconsequencethatthemoreparallelasolutionmethodis,thesloweritwillconvergebecauseofthedecreasedglobalcoupling.•Parabolicequations.Thesegiveanellipticlinearsystemifanimplicitmethod,suchasbackwardEuler,isused.•Nonlinearproblems.TheJacobiansysteminaNewton-Raphsonmethodcanbesolvediteratively.Iterativemethodsareparticularlyattractiveherebecausefullpre-cisionisoftennotneeded,andonecanstopiteratingafterthedesiredprecisionhasbeenreached.•Eigenvaluecalculations.Inordertofindinterioreigenvaluesofamatrix’spectrum,oneoftenappliesaLanczosmethodtoashiftedandinvertedformofA:ineffectonecomputeseigenvaluesof(A−σ)−1.SinceaLanczosmethodrequiresmultiplicationwiththecoefficientmatrix,weneedsolveasystemwithA−σineachiteration.Thiscanbedoneiteratively,althoughthesesystemsarebytheirverynatureindefinite,whichmakesapplyingiterativemethodsharder.1.Sotospeak1Oftheabove,ellipticequationsarebyfarthemostnaturalenvironmentforiterativemeth-odstothrivein.Hence,somefurtherremarks.ThediscretisationofanellipticPDEonaregulardomainhasameshsizeparameterhassociatedwithit.Onecanthenshowthatthe‘spectralconditionnumber’κ(A)=kAk·kA−1kofamatrixisofO(h−2)forsecondor-derproblems,andO(h−4)forfourthorderproblems,independentofthenumberofspacedimensions.ThePDEoriginoflinearsystemshasimplicationsfortheirsparsitystructure.Discretisa-tionofPDEsbyfinitedifferences,finiteelements,orfinitevolumesestablishesconnectionofeachdiscretisationpointwithjustafewneighbours,thenumberofneighboursbeingindependentofthetotaldiscretisation.Thismeansthatthenumberofnonzerosperrowisrathersmall,andnotrelatedtothematrixsize.1.2IterativeversusdirectmethodsThealternativetoiterativemethodsisGaussianelimination,whichisequivalenttofindingaLUfactorizationofthematrixandusingthattosolvethesystem.Theuseofafactorizationiscalledadirectmethodsinceityieldsitsresultsinafixednumberofoperations.Directmethodshaveawiderapplicabilitythaniterativemethods;ontheotherhand,theycantakelargeamountsoftimeandspace,aselementsthatarezerointhesparsecoefficientmatrixwillfillinduringthefactorisation.Fill-incanbesomewhataleviatedbyfindinganappropriateorderingoftheequations.Forinstance,intwospacedimensions,therecursivebisection(ornesteddissection)method[25]reducesfill-intoaboutNlogN.However,itcanbeprovedthatthisdoesnotholdinhigherdimensions[35,36].Iterativemethods,unlikedirectmethods,dependonnumericalpropertiesofthelinearsys-temfortheirtimetosolution.Thismakesthembothpotentiallyfasterthandirectmethods,andpotentiallyunsuitableforcertainlinearsystems.Italsomeansthatapplyingthemisconsiderablylessstraightforwardthanusingdirectmethods.2Iterativemethods2.1StationaryiterativemethodsTheeasiestiterativemethodstoexplainaretheso-called‘stationaryiterativemethods’.Here’stheidea.SupposethatsolvingasystemAx=bisnotcomputationallyfeasible,butthatsolvingMx=bispossible,andthatpresumablyinsomesenseM≈A.Letx1(thefirstiterate)bethesolutionofMx1=b.Withthiswegettheerrorvectore1=x1−xandresidualvectorr1=Ax1−b.Thetwoarerelatedbyr1=Ae1.Sincex=A−1b=x1−A−1r1,wecandefinex2=x1−M−1r1,andx3=x2−M−1r2,etcetera.Thisprocedureisalsocalled‘iterativerefinement’.Thebasicschemeisthenxi+1=xi−M−1ri,(1)andthisiscalledstationarybecauseeveryiterationfollowsthesamerecipe,withoutitera-tiondependencies.2Theanalysisofthisprocedureisfairlysimple.Therelationbetweenthefirsttworesidualsisr2=Ax2−b=A(x1−M−1r1)−b=(I−AM−1)r1whichinductivelygivesrn=(I−AM−1)n−1r1(2)sothemethodwillconvergeifall|λ(I−AM−1)|1(3)Forfuturereferencewenotethatequation(2)canbewrittenmoreabstractlyasrn=P(n−1)(AM−1)ri(4)whereP(n−1)(x)isapolynomialofdegreen−1andP(n)(0)=1foralln.Becauseofthisform,theseiterativemethodsarecalled‘polynomialiterativemethods’.Another,equivalent,wayoflookingatpolynomialiterativemethodsistoconsiderthe‘Krylovsequence’ki+1=AM−1kiwherek1=r1.Theresidualsarethencombinationsofthissequence,andthemethodgeneratingtheresid-ualsiscalleda‘Krylovmethod’2.ThespanoftheKrylovvectors,andtheinitialpartsofit,arecalledaKrylovspace,orKrylovsubspace,respectively.2.1.1ExamplesHerear