Table Of Content

An Introduction to the Conjugate Gradient Method Without the Agonizing Pain Edition11 4 JonathanRichardShewchuk August4,1994 SchoolofComputerScience CarnegieMellonUniversity Pittsburgh,PA15213 Abstract TheConjugateGradientMethodisthemostprominentiterativemethodforsolvingsparsesystemsoflinearequations. Unfortunately, many textbook treatments of the topic are written with neither illustrations nor intuition, and their victims can be found to this day babbling senselessly in the corners of dusty libraries. For this reason, a deep, geometricunderstandingofthemethodhasbeenreservedfortheelitebrilliantfewwhohavepainstakinglydecoded themumblingsoftheirforebears. Nevertheless,theConjugateGradientMethodisacompositeofsimple,elegantideas thatalmostanyonecanunderstand. Ofcourse,areaderasintelligentasyourselfwilllearnthemalmosteffortlessly. TheideaofquadraticformsisintroducedandusedtoderivethemethodsofSteepestDescent,ConjugateDirections, and Conjugate Gradients. Eigenvectorsare explained and used to examine the convergenceof the Jacobi Method, SteepestDescent,andConjugateGradients. OthertopicsincludepreconditioningandthenonlinearConjugateGradient Method. I have taken pains to make this article easy to read. Sixty-six illustrations are provided. Dense prose is avoided. Conceptsareexplainedinseveraldifferentways. Mostequationsarecoupledwithanintuitiveinterpretation. SupportedinpartbytheNaturalSciencesandEngineeringResearchCouncilofCanadaundera1967ScienceandEngineering Scholarshipand by the National ScienceFoundationunder Grant ASC-9318163. The views and conclusionscontainedin this documentarethoseoftheauthorandshouldnotbeinterpretedasrepresentingtheofficialpolicies,eitherexpressorimplied,of NSERC,NSF,ortheU.S.Government. Keywords: conjugategradientmethod,preconditioning,convergenceanalysis,agonizingpain Contents 1. Introduction 1 2. Notation 1 3. TheQuadraticForm 2 4. TheMethodofSteepestDescent 6 5. ThinkingwithEigenvectorsandEigenvalues 9 (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:3)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0) 5.1. EigendoitifItry 9 (cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:3)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0) 5.2. Jacobiiterations 11 (cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:3)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0) 5.3. AConcreteExample 12 6. ConvergenceAnalysisofSteepestDescent 13 (cid:0)(cid:3)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:3)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0) 6.1. InstantResults 13 (cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:3)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0) 6.2. GeneralConvergence 17 7. TheMethodofConjugateDirections 21 (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:3)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:3)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0) 7.1. Conjugacy 21 (cid:0)(cid:3)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0) 7.2. Gram-SchmidtConjugation 25 (cid:0)(cid:3)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0) 7.3. OptimalityoftheErrorTerm 26 8. TheMethodofConjugateGradients 30 9. ConvergenceAnalysisofConjugateGradients 32 (cid:0)(cid:3)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0) 9.1. PickingPerfectPolynomials 33 (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:3)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0) 9.2. ChebyshevPolynomials 35 10.Complexity 37 11.StartingandStopping 38 (cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:3)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:3)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0) 11.1. Starting 38 (cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:3)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:3)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0) 11.2. Stopping 38 12.Preconditioning 39 13.ConjugateGradientsontheNormalEquations 41 14.TheNonlinearConjugateGradientMethod 42 (cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0) 14.1. OutlineoftheNonlinearConjugateGradientMethod 42 (cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:3)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0) 14.2. GeneralLineSearch 43 (cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:3)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0) 14.3. Preconditioning 47 A Notes 48 B CannedAlgorithms 49 (cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:3)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0) B1. SteepestDescent 49 (cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:3)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0) B2. ConjugateGradients 50 (cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0) B3. PreconditionedConjugateGradients 51 i (cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0) B4. NonlinearConjugateGradientswithNewton-RaphsonandFletcher-Reeves 52 (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0) B5. PreconditionedNonlinearConjugateGradientswithSecantandPolak-Ribie`re 53 (cid:4)(cid:6)(cid:5)(cid:8)(cid:7)(cid:10)(cid:9) C UglyProofs 54 (cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0) C1. TheSolutionto MinimizestheQuadraticForm 54 (cid:11) (cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0) C2. ASymmetricMatrixHas OrthogonalEigenvectors. 54 (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:3)(cid:0)(cid:1)(cid:0)(cid:2)(cid:0) C3. OptimalityofChebyshevPolynomials 55 D Homework 56 ii AboutthisArticle An electronic copy of this article is available by anonymous FTP to WARP.CS.CMU.EDU (IP address 128.2.209.103) under the filename quake-papers/painless-conjugate-gradient.ps. A PostScript file containing full-page copies of the figures herein, suitable for transparencies, is available as quake-papers/painless-conjugate-gradient-pics.ps. Most of the illustrations were createdusingMathematica. (cid:12) c1994byJonathanRichardShewchuk. Thisarticlemaybefreelyduplicatedanddistributedsolongas noconsiderationisreceivedinreturn,andthiscopyrightnoticeremainsintact. ThisguidewascreatedtohelpstudentslearnConjugateGradientMethodsaseasilyaspossible. Please mail me ([email protected]) comments, corrections, and any intuitions I might have missed; some of these will be incorporated into a second edition. I am particularly interested in hearing about use of this guideforclassroomteaching. Forthosewhowishtolearnmoreaboutiterativemethods,IrecommendWilliamL.Briggs’“AMultigrid Tutorial”[2],oneofthebest-writtenmathematicalbooksIhaveread. Special thanks to Omar Ghattas, who taught me much of what I know about numerical methods, and provided me with extensive comments on the first draft of this article. Thanks also to James Epperson, DavidO’Hallaron,JamesStichnoth,NickTrefethen,andDanielTunkelangfortheircomments. Tohelpyouskipchapters,here’sadependencegraphofthesections: 1 Introduction 10 Complexity 2 Notation 13 Normal Equations 3 Quadratic Form 4 Steepest Descent 5 Eigenvectors 7 Conjugate Directions 6 SD Convergence 8 Conjugate Gradients 9 CG Convergence 11 Stop & Start 12 Preconditioning 14 Nonlinear CG ThisarticleisdedicatedtoeverymathematicianwhousesfiguresasabundantlyasIhaveherein. iii iv 1. Introduction WhenIdecidedtolearntheConjugateGradientMethod(henceforth,CG),Ireadfourdifferentdescriptions, whichIshallpolitelynotidentify. Iunderstoodnoneofthem. Mostofthemsimplywrotedownthemethod, thenproveditspropertieswithoutanyintuitiveexplanationorhintofhowanybodymighthaveinventedCG inthefirstplace. Thisarticlewasbornofmyfrustration,withthewishthatfuturestudentsofCGwilllearn arichandelegantalgorithm,ratherthanaconfusingmassofequations. CG is the most popular iterative methodfor solving large systemsof linear equations. CG is effective (cid:4)(cid:13)(cid:5)(cid:14)(cid:7)(cid:15)(cid:9) forsystemsoftheform (1) (cid:5) (cid:9) (cid:4) where isanunknownvector, isaknownvector,and isaknown,square,symmetric,positive-definite (or positive-indefinite)matrix. (Don’t worry if you’ve forgotten what “positive-definite” means; we shall review it.) These systems arise in many important settings, such as finite difference and finite element methodsforsolvingpartialdifferentialequations,structuralanalysis,circuitanalysis,andmathhomework. (cid:4) (cid:4) Iterative methods like CG are suited for use with sparse matrices. If is dense, your best course of (cid:4) (cid:4) actionisprobablytofactor andsolvetheequationbybacksubstitution. Thetimespentfactoringadense (cid:9) isroughlyequivalenttothetimespentsolvingthesystemiteratively;andonce isfactored,thesystem (cid:4) can be backsolved quickly for multiple values of . Compare this dense matrix with a sparse matrix of (cid:4) larger size that fills the same amount of memory. The triangular factors of a sparse usually have many more nonzero elements than itself. Factoring may be impossible due to limited memory, and will be time-consuming as well; even the backsolving step may be slower than iterative solution. On the other hand,mostiterativemethodsarememory-efficientandrunquicklywithsparsematrices. I assume that you have taken a first course in linear algebra, and that you have a solid understanding of matrix multiplication and linear independence, although you probably don’t remember what those eigenthingieswereallabout. Fromthisfoundation,IshallbuildtheedificeofCGasclearlyasIcan. 2. Notation Letusbeginwithafewdefinitionsandnotesonnotation. (cid:4) (cid:5) (cid:9) Withafewexceptions,Ishallusecapitalletterstodenotematrices,lowercaseletterstodenotevectors, (cid:11)(cid:17)(cid:16)(cid:18)(cid:11) (cid:11)(cid:17)(cid:16) andGreekletterstodenotescalars. isan matrix,and and arevectors—thatis, 1matrices. Writtenoutfully,Equation1is (cid:19)(cid:20) (cid:4) (cid:4) (cid:4) (cid:25)(cid:27)(cid:26) (cid:19)(cid:20) (cid:5) (cid:25)(cid:27)(cid:26) (cid:19)(cid:20) (cid:9) (cid:25)(cid:27)(cid:26) (cid:20)(cid:20) (cid:0)(cid:22)(cid:0)(cid:22)(cid:0) (cid:23) (cid:26)(cid:26) (cid:20)(cid:20) (cid:26)(cid:26) (cid:20)(cid:20) (cid:26)(cid:26) (cid:20) (cid:4) (cid:4) (cid:4) (cid:26) (cid:20) (cid:5) (cid:26) (cid:20) (cid:9) (cid:26) 11 12 1(cid:23) 1 (cid:7) 1 (cid:0) 21 22 2 2 2 (cid:21) (cid:4) ... (cid:4) ... (cid:4) ... (cid:28) (cid:21) (cid:5) ... (cid:28) (cid:21) (cid:9) ... (cid:28) (cid:23) (cid:23) (cid:0)(cid:22)(cid:0)(cid:22)(cid:0) (cid:23)(cid:24)(cid:23) (cid:23) (cid:23) 1 2 (cid:5)(cid:30)(cid:29) (cid:31) (cid:23)"$# (cid:5) " (cid:31) " ! (cid:5)(cid:30)(cid:29) (cid:31)(cid:18)(cid:7)%(cid:31)&(cid:29) (cid:5) (cid:5) (cid:31) (cid:5)(cid:30)(cid:29) (cid:31)(cid:18)(cid:7) The inner product of two vectors is written , and represents the scalar sum . Note that 1 (cid:16) (cid:5)(cid:30)(cid:29) (cid:31) (cid:5)’(cid:29)((cid:4)(cid:6)(cid:5) . If and areorthogonal,then 0. Ingeneral,expressionsthatreduceto1 1matrices, suchas and ,aretreatedasscalarvalues. (cid:5) 2 JonathanRichardShewchuk 2 4 2 (cid:5) (cid:5) (cid:7) ) 3 2 2 (cid:5) 1 2 1 -4 -2 2 4 6 (cid:5) (cid:5) (cid:7)+* ) -2 2 1 6 2 8 -4 -6 Figure1: Sampletwo-dimensionallinearsystem. Thesolutionliesattheintersectionofthelines. (cid:4) (cid:5) Amatrix ispositive-definiteif,foreverynonzerovector , (cid:5) (cid:29) (cid:4)(cid:6)(cid:5)-, (cid:0) 0 (2) This may mean littleto you, butdon’t feel bad; it’s not a veryintuitive idea, and it’s hard to imaginehow a matrixthat is positive-definitemight look differently fromone that isn’t. We will get a feeling for what positive-definitenessisaboutwhenweseehowitaffectstheshapeofquadraticforms. (cid:4)0/(cid:18)12(cid:29)3(cid:7)4/5(cid:29) (cid:4)6(cid:29) (cid:4)0/7198 (cid:7)4/:8 (cid:4)(cid:3)8 . . Finally,don’tforgettheimportantbasicidentities and 1 1 1. 3. TheQuadraticForm Aquadraticformissimplyascalar,quadraticfunctionofavectorwiththeform ; (cid:5)<1=(cid:7) (cid:5) (cid:29) (cid:4)(cid:6)(cid:5)>*?(cid:9) (cid:29) (cid:5) . 1 )A@ (3) (cid:4) (cid:5) (cid:9) 2 (cid:4) @ ; (cid:5)<1 (cid:4)(cid:13)(cid:5)(cid:14)(cid:7)(cid:15)(cid:9) where isamatrix, and arevectors,and isascalarconstant. Ishallshowshortlythatif issymmetric . andpositive-definite, isminimizedbythesolutionto . Throughoutthispaper,Iwilldemonstrateideaswiththesimplesampleproblem (cid:4)4(cid:7)CB (cid:9)F(cid:7)CB (cid:7) 3 2 * 2 @ (cid:0) D(cid:2)E DGE 0 (4) 2 6 8 (cid:4)(cid:13)(cid:5)H(cid:7)I(cid:9) (cid:5) * (cid:5)%(cid:7)KJ * (cid:29) The system is illustrated in Figure 1. In general, the solution lies at the intersection point (cid:11) (cid:11) L of hyperplanes, each having; d(cid:5)Mim1 ension 1. For this problem, the so; lu(cid:5)<ti1on is 2E 2 . The . . correspondingquadratic form appears in Figure 2. A contour plot of is illustrated in Figure 3. TheQuadraticForm 3 150 ; (cid:5)<1 100 . 4 50 2 0 (cid:5) 0 66 2 -2 44 (cid:5)(cid:5) 22 -4 00 11 --22 -6 --44 N<OQPSR TUPWVYX Figure2: Graphofaquadraticform . Theminimumpointofthissurfaceisthesolutionto . (cid:5) 2 4 2 (cid:5) 1 -4 -2 2 4 6 -2 -4 -6 N<O$PZR Figure3: Contoursofthequadraticform. Eachellipsoidalcurvehasconstant . (cid:5) 4 JonathanRichardShewchuk 2 6 4 2 (cid:5) -4 -2 2 4 6 1 -2 -4 -6 -8 NZ[\O$PZR P Figure4: Gradient ofthequadraticform. Forevery ,thegradientpointsinthedirectionofsteepest N<O$PZR increaseof ,andisorthogonaltothecontourlines. (cid:4) ; (cid:5)M1 . Because ispositive-definite,thesurfacedefinedby isshapedlikeaparaboloidbowl. (I’llhavemore tosayaboutthisinamoment.) Thegradientofaquadraticformisdefinedtobe (cid:19)(cid:20) _ ; (cid:5)<1 (cid:25)(cid:26) (cid:20) _(cid:22)‘ . (cid:26) ;^] (cid:5)M1=(cid:7) (cid:20)(cid:20) _(cid:22)_‘ 1; . (cid:5)<1 (cid:26)(cid:26) . (cid:0) (cid:21) _a‘a_ b2;... . (cid:5)M1 (cid:28) (5) (cid:5) ; (cid:5)<1 . The gradient is a vectorfield that, for a given point , pointsin the directionof greatestincrease of . Figure4illustratesthegradientvectorsforEquation3withtheco; n(cid:5)<st1antsgivenin; E] q(cid:5)Mu1 ation4. Atthebottom . . oftheparaboloidbowl,thegradientiszero. Onecanminimize bysetting equaltozero. Withalittlebitoftediousmath,onecanapplyEquation5toEquation3,andderive ; ] (cid:5)M1c(cid:7) (cid:4) (cid:29) (cid:5) (cid:4)(cid:6)(cid:5):*?(cid:9) . 1 ) 1 (cid:0) (6) (cid:4) 2 2 If issymmetric,thisequationreducesto ; ] (cid:5)M1c(cid:7)4(cid:4)(cid:13)(cid:5)d*?(cid:9) . (cid:0) (7) (cid:4)(cid:13)(cid:5)%(cid:7)e(cid:9) ; (cid:5)M1 (cid:4) Setting the gradient to zero, we obtain Equation 1, the linear system we wish to solve. Therefore, the . solution to is a critical point of . If is positive-definite as well as symmetric, then this