Conjugate Gradient Method - Stanford University

1 Conjugate Gradient Method direct and indirect methods positive definite linear systems Krylov sequence derivation of the Conjugate Gradient Method spectral analysis of Krylov sequence preconditioningEE364b, Stanford UniversityProf. Mert Pilanciupdated: May 5, 2022 Three classes of methods for linear equationsmethods to solve linear systemAx=b,A Rn n dense direct(factor-solve methods) runtime depends only on size; independent of data, structure, orsparsity work well fornup to a few thousand sparse direct(factor-solve methods) runtime depends on size, sparsity pattern; (almost) independent ofdata can work well fornup to104or105(or more) requires good heuristic for orderingEE364b, Stanford University1 indirect(iterative methods) runtime depends on data, size, sparsity, required accuracy requires tuning, preconditioning.

2 Good choice in many cases; only choice forn= 106or largerEE364b, Stanford University2 Symmetric positive definite linear systemsSPD system of equationsAx=b, A Rn n, A=AT 0examples Newton/interior-point search direction: 2 (x) x= (x) least-squares normal equations:(ATA)x=ATb regularized least-squares:(ATA+ I)x=ATb minimization of convex quadratic function(1/2)xTAx bTx solving (discretized) elliptic PDE ( , Poisson equation)EE364b, Stanford University3 analysis of resistor circuit:Gv=i vis node voltage (vector),iis (given) source current Gis circuit conductance matrixGij={total conductance incident on nodei i=j (conductance between nodesiandj)i6=jEE364b, Stanford University4CG overview proposed by Hestenes and Stiefel in 1952 (as direct Method ) solves SPD systemAx=b in theory ( , exact arithmetic)}

3 Inniterations each iteration requires a few inner products inRn, and onematrix-vector multiplyz Az forAdense, matrix-vector multiplyz Azcostsn2, so total cost isn3, same as direct methods get advantage over dense if matrix-vector multiply is cheaper thann2 with roundoff error, CG can work poorly (or not at all) but for someA(andb), can get good approximate solution in niterationsEE364b, Stanford University5 Solution and error x?=A 1bis solution x?minimizes (convex function)f(x) = (1/2)xTAx bTx f(x) =Ax bis Gradient off withf?=f(x?), we havef(x) f?= (1/2)xTAx bTx (1/2)x?TAx?+bTx?= (1/2)(x x?)

4 TA(x x?)= (1/2) x x? ,f(x) f?is half of squaredA-norm of errorx x?EE364b, Stanford University6 a relative measure (comparingxto0): =f(x) f?f(0) f?= x x? 2A x? 2A(fraction of maximum possible reduction inf, compared tox= 0)EE364b, Stanford University7 Residual r=b Axis called theresidualatx r= f(x) =A(x? x) in terms ofr, we havef(x) f?= (1/2)(x x?)TA(x x?)= (1/2)rTA 1r= (1/2) r 2A 1 a commonly used measure of relative accuracy: = r / b (A) 2( is easily computable fromx; is not)EE364b, Stanford University8 Krylov subspace( controllability subspace)Kk= span{b,Ab,..,Ak 1b}={p(A)b|ppolynomial,degp < k}we define theKrylov sequencex(1),x(2).

5 Asx(k)= argminx Kkf(x) = argminx Kk x x? 2 Athe CG algorithm (among others) generates the Krylov sequenceEE364b, Stanford University9 Properties of Krylov sequence f(x(k+1)) f(x(k))(but r can increase) x(n)=x?( ,x? Kneven whenKn6=Rn) x(k)=pk(A)b, wherepkis a polynomial withdegpk< k less obvious: there is atwo-term recurrencex(k+1)=x(k)+ kr(k)+ k(x(k) x(k 1))for some k, k(basis of CG algorithm)EE364b, Stanford University10 Cayley-Hamilton theoremcharacteristic polynomial ofA: (s) = det(sI A) =sn+ 1sn 1+ + nby Caley-Hamilton theorem (A) =An+ 1An 1+ + nI= 0and soA 1= (1/ n)An 1 ( 1/ n)An 2 ( n 1/ n)Iin particular, we see thatx?

6 =A 1b KnEE364b, Stanford University11 Deriving the Conjugate Gradient Method suppose{d0,..,dk 1}is an orthogonal basis forKkunder theA innerproduct, ,(di)TAdj= 0 i6=j letx(k)= k 1i=0 idif(x(k)) =12(x(k))TAx(k) bTx(k)=12(k 1 i=0 idi)A(k 1 i=0 idi) bT(k 1 i=0 idi)=12k 1 i=0 2i(di)TAdi ibTdi decomposable problem, optimal i=bTdi/ di 2A iEE364b, Stanford University12 Constructing the basis (slow) suppose{d0,..,dk 1}is an orthogonal basis forKk extend{d0,..,dk 1}to a basis ofKk+1using Gram-Schmidt proceduredk=g k 1 j=0dj(dj)TAg(dj)TAdj where we can pick anyg Kk+1such thatg / Kk, ,g=Akb after constructing the basis, setx(k)= k 1j=0 jdj= kj=0bTdj(dj)TAdjdj Conjugate Gradient Method constructs the basis online by pickingg.

7 = f(x(k)) =b Ax(k)EE364b, Stanford University13 Simplifying Gram-Schmidt procedure initialize atx(0)= 0d0= f(x(k)) =b Ax(0)=bd1= (b Ax(1)) d0(d0)TA(b Ax(1))(d0)TAd0d2= (b Ax(2)) d1(d1)TA(b Ax(2))(d1)TAd1 d0(d0)TA(b Ax(2))(d0)TAd0d3= (b Ax(3)) d2(d2)TA(b Ax(3))(d2)TAd3 d1(d1)TA(b Ax(3))(d1)TAd1 d0(d0)TA(b Ax(3))(d0) turns out red terms are zero due to orthogonalityEE364b, Stanford University14 sincex(k)= kj=0 jdjwe havex(k+1)=x(k)+ kdk we pickd0= f(x(0))and fork= 1,..,n 1dk= f(x(k)) k 1 j=0djdTjA( f(x(k)))dTjAdj it holds thatdTj f(x(k)) = 0and f(x(j))T f(x(k)) = 0 j < kproof:we havedT0 f(x(1)) =dT0(Ax0 0 b) + 0dT0Ad0= 0anddTj f(x(k+1)) =dTj(Ax(k)+ kdk b) =dTj f(x(k)) kdTjAdk 0forj6=kin additionspan(d0.)

8 ,dk) =span( f(x(0),.., f(x(k)))EE364b, Stanford University15 Simplifying CG update finally, note that f(x(j+1)) f(x(j)) =A(x(j)+ jdj) b (Ax(j) b) = jAdj simplify the basis updatedk= f(x(k)) k 1 j=0djdTjA( f(x(k)))dTjAdj= f(x(k)) k 1 j=0dj( f(x(j+1)) f(x(j)))T( f(x(k)))( f(x(j+1)) f(x(j)))Tdj= f(x(k)) dk 1 f(x(k))T( f(x(k)))( f(x(k 1)))Tdk 1EE364b, Stanford University16 Final CG update rule CG algorithm simplifies todk= f(x(k)) +dk 1 f(x(k)) 22 f(x(k 1)) 22 x(k+1)=x(k)+ kdkwhere = arg min f(xk+ dk) =bTdkdTkAdksince we havef(x(k 1))Tdk 1= f(x(k 1)) : f(x(k))Tdk= f(x(k))T( f(x(k)) dk 1 f(x(k))T( f(x(k)))( f(x(k 1)))Tdk 1)= f(x(k))T f(x(k))EE364b, Stanford University17 Spectral analysis of Krylov sequence A=Q QT,Qorthogonal, =diag( 1.)

9 , n) definey=QTx, b=QTb,y?=QTx? in terms ofy, we havef(x) = f(y) = (1/2)xTQ QTx bTQQTx= (1/2)yT y bTy=n i=1((1/2) iy2i biyi)soy?i= bi/ i,f?= (1/2) ni=1 b2i/ iEE364b, Stanford University18 Krylov sequence in terms ofyy(k)= argminy Kk f(y), Kk= span{ b, b,.., k 1 b}y(k)i=pk( i) bi,degpk< kpk= argmindegp<kn i=1 b2i((1/2) ip( i)2 p( i))EE364b, Stanford University19f(x(k)) f?= f(y(k)) f?=mindegp<k(1/2)n i=1 b2i( ip( i) 1)2 i=mindegp<k(1/2)n i=1 y?2i i( ip( i) 1)2=mindegq k, q(0)=1(1/2)n i=1 y?2i iq( i)2=mindegq k, q(0)=1(1/2)n i=1 b2iq( i)2 iEE364b, Stanford University20 k=mindegq k, q(0)=1 ni=1 y?

10 2i iq( i)2 ni=1 y?2i i mindegq k, q(0)=1(maxi=1,..,nq( i)2) if there is a polynomialqof degreek, withq(0) = 1, that is small onthe spectrum ofA, thenf(x(k)) f?is small if eigenvalues are clustered inkgroups, theny(k)is a good approximatesolution if solutionx?is approximately a linear combination ofkeigenvectors ofA, theny(k)is a good approximate solutionEE364b, Stanford University21A bound on convergence rate takingqas Chebyshev polynomial of degreek, that is small on interval[ min, max], we get k ( 1 + 1)k, = max/ min convergence can be much faster than this, if spectrum ofAis spreadbut clusteredEE364b, Stanford University22 Small exampleA R7 7, spectrum shown as filled circles.