f10-6

10.6ConjugateGradientMethodsinMultidimensions 413Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).xi(j)=xmin*xi(j) p(j)=p(j)+xi(j) enddo 12 return END FUNCTION f1dim(x) INTEGER NMAXREAL f1dim,func,x PARAMETER (NMAX=50) C USES func Used by linmin as the function passed to mnbrak andbrent . INTEGER j,ncom REAL pcom(NMAX),xicom(NMAX),xt(NMAX) COMMON /f1com/ pcom,xicom,ncomdo 11j=1,ncom xt(j)=pcom(j)+x*xicom(j) enddo 11 f1dim=func(xt) return END CITED REFERENCES AND FURTHER READING: Brent,R.P.1973, AlgorithmsforMinimizationwithoutDerivatives (EnglewoodCliffs,NJ:Prentice- Hall), Chapter 7. [1] Acton, F.S. 1970, Numerical Methods That Work ; 1990, corrected edition (Washington: Mathe- matical Association of America), pp. 464–467. [2] Jacobs, D.A.H. (ed.) 1977, The State of the Art in Numerical Analysis (London: Academic Press), pp. 259–262. 10.6 Conjugate Gradient Methods in Multidimensions We consider now the case where you are able to calculate, at a given N- dimensional point P, not just the value of a function f(P)but also the gradient (vector of first partial derivatives) ∇f(P). Aroughcountingargumentwillshowhowadvantageousitistousethegradient information: Suppose that the function fis roughly approximated as a quadratic form, as above in equation (10.5.1), f(x)≈c−b·x+1 2x·A·x (10.6.1 ) Then the number of unknown parameters in fis equal to the number of free parameters in Aandb, which is1 2N(N+1 ), which we see to be of order N2. Changing any one of these parameters can move the location of the minimum. Therefore, we should not expect to be able to findthe minimum until we have collected an equivalent information content, of order N2numbers. 414 Chapter10. MinimizationorMaximizationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Inthedirectionsetmethodsof §10.5,wecollectedthenecessaryinformationby makingon the orderof N2separate line minimizations, eachrequiring“a few” (but sometimes a bigfew!) function evaluations. Now, each evaluation of the gradient will bring us Nnew components of information. If we use them wisely, we should need to make only of order Nseparate line minimizations. That is in fact the case for the algorithms in this section and the next. A factor of Nimprovementin computationalspeed is not necessarily implied. As a rough estimate, we might imagine that the calculation of each component of the gradient takes about as long as evaluating the function itself. In that case there will be of order N2equivalent function evaluations both with and without gradient information. Even if the advantage is not of order N, however, it is nevertheless quite substantial: (i) Each calculated component of the gradient will typically save not just one function evaluation, but a number of them, equivalent to, say, a wholeline minimization. (ii) There is often a high degree of redundancy in the formulas forthevariouscomponentsofafunction’sgradient;whenthisisso,especiallywhen there is also redundancywith the calculation of the function,then the calculationofthe gradient may cost significantly less than Nfunction evaluations. Acommonbeginner’serroristoassumethatanyreasonablewayofincorporating gradientinformationshouldbeaboutasgoodasanyother. Thislineofthoughtleads to the following not very good algorithm, the steepest descent method : Steepest Descent: Start at a point P0. As many times as needed, move from point Pito the point Pi+1by minimizing along the line from Piin the direction of the local downhill gradient −∇f(Pi). The problem with the steepest descent method (which, incidentally, goes back to Cauchy), is similar to the problem that was shown in Figure 10.5.1. The methodwillperformmanysmallstepsingoingdownalong,narrowvalley,evenifthevalley is a perfect quadratic form. You might have hoped that, say in two dimensions, your first step would take you to the valley floor, the second step directly downthe long axis; but remember that the new gradient at the minimum point of any line minimization is perpendicular to the direction just traversed. Therefore, with the steepest descent method, you mustmake a right angle turn, which does not,i n general, take you to the minimum. (See Figure 10.6.1.) Just as in the discussion that led up to equation (10.5.5), we really want a way of proceeding not down the new gradient, but rather in a direction that is somehow constructed to be conjugate to the old gradient, and, insofar as possible, to all previous directions traversed. Methods that accomplish this construction are calledconjugate gradient methods. In§2.7 we discussed the conjugate gradient method as a technique for solving linearalgebraicequationsbyminimizinga quadraticform. Thatformalismcan also be applied to the problem of minimizing a function approximated by the quadratic form (10.6.1). Recall that, starting with an arbitrary initial vector g 0and letting h0=g0, the conjugate gradient method constructs two sequences of vectors from the recurrence gi+1 =gi−λiA·hihi+1 =gi+1+γihi i=0,1,2,... (10.6.2 ) 10.6ConjugateGradientMethodsinMultidimensions 415Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).(a) (b) Figure 10.6.1. (a) Steepest descent method in a long, narrow “valley.”While more ef ficient than the strategy of Figure 10.5.1, steepest descent is nonetheless an inef ficient strategy, taking many steps to reach the valley floor. (b) Magni fied view of one step: A step starts off in the local gradient direction, perpendicular to the contour lines, and traverses a straight line until a local minimum is reached, wherethe traverse is parallel to the local contour lines. The vectors satisfy the orthogonality and conjugacy conditions gi·gj=0hi·A·hj=0gi·hj=0 j<i (10.6.3 ) The scalars λiand γiare given by λi=gi·gi hi·A·hi=gi·hi hi·A·hi(10.6.4 ) γi=gi+1·gi+1 gi·gi(10.6.5 ) Equations (10.6.2) –(10.6.5)are simply equations (2.7.32) –(2.7.35)for a symmetric Ain a new notation. (A self-contained derivation of these results in the context of function minimization is given by Polak [1].) Now suppose that we knew the Hessian matrix Ain equation (10.6.1). Then we could use the construction (10.6.2) to find successively conjugate directions hi along which to line-minimize. After Nsuch, we would ef ficiently have arrived at the minimum of the quadratic form. But we don ’t knowA. Here is a remarkable theorem to save the day: Suppose we happen to have gi=−∇f(Pi),forsomepoint Pi,where fis oftheform(10.6.1). Supposethatwe proceed from Pialong the direction hito the local minimum of flocated at some pointPi+1and then set gi+1 =−∇f(Pi+1). Then, this gi+1is the same vector as would have been constructed by equation (10.6.2). (And we have constructed it without knowledge of A!) Proof: By equation (10.5.3), gi=−A·Pi+b, and gi+1 =−A·(Pi+λhi)+b=gi−λA·hi (10.6.6 ) 416 Chapter10. MinimizationorMaximizationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).with λchosen to take us to the line minimum. But at the line minimum hi·∇f= −hi·gi+1 =0. This latter condition is easily combined with (10.6.6) to solve for λ. The result is exactly the expression (10.6.4). But with this value of λ, (10.6.6) is the same as (10.6.2), q.e.d. We have, then, the basis of an algorithmthat requiresneitherknowledgeof the Hessianmatrix A,noreventhestoragenecessarytostoresuchamatrix. Asequence of directions hiis constructed, using only line minimizations, evaluations of the gradientvector,and an auxiliaryvectorto store the latest in the sequenceof g’s. The algorithm described so far is the original Fletcher-Reeves version of the conjugate gradient algorithm. Later, Polak and Ribiere introduced one tiny, butsometimes signi ficant, change. They proposed using the form γ i=(gi+1−gi)·gi+1 gi·gi(10.6.7 ) insteadofequation(10.6.5). “Wait,”yousay,“aren’ttheyequalbytheorthogonality conditions (10.6.3)? ”They are equal for exact quadratic forms. In the real world, however, your function is not exactly a quadratic form. Arriving at the supposed minimum of the quadratic form, you may still need to proceed for another set ofiterations. There is some evidence [2]that the Polak-Ribiere formula accomplishes the transition to further iterations more gracefully: When it runs out of steam, it tends to reset hto be down the local gradient, which is equivalent to beginning the conjugate-gradient procedure anew. The followingroutine implements the Polak-Ribiere variant, which we recom- mend;butchangingoneprogramline,asshown,willgiveyouFletcher-Reeves. The routine presumesthe existence of a function func(p), where p(1:n)is a vectorof length n, andalso presumesthe existenceofa subroutine dfunc(p,df) thatreturns the vector gradient df(1:n) evaluated at the input point p. The routine calls linminto do the line minimizations. As already discussed, you may wish to use a modi fied version of linminthat uses dbrentinstead of brent, i.e.,that uses the gradientin doingthe line minimizations. See notebelow. SUBROUTINE frprmn(p,n,ftol,iter,fret) INTEGER iter,n,NMAX,ITMAXREAL fret,ftol,p(n),EPS,funcEXTERNAL func PARAMETER (NMAX=50,ITMAX=200,EPS=1.e-10) C USES dfunc,func,linmin Given a starting point pthat is a vector of length n, Fletcher-Reeves-Polak-Ribiere minimiza- t i o ni sp e r f o r m e do naf u n c t i o n func , using its gradient as calculated by a routine dfunc . The convergence tolerance on the function value is input as ftol . Returned quantities are p(the location of the minimum), iter (the number of iterations that were performed), andfret (the minimum value of the function). The routine linmin is called to perform line minimizations. Parameters: NMAX is the maximum anticipated value of n;ITMAX is the maximum allowed number of iterations; EPS is a small number to rectify special case of converging to exactly zero function value. INTEGER its,jREAL dgg,fp,gam,gg,g(NMAX),h(NMAX),xi(NMAX)fp=func(p) Initializations. call dfunc(p,xi) do 11j=1,n g(j)=-xi(j) h(j)=g(j) 10.6ConjugateGradientMethodsinMultidimensions 417Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).xi(j)=h(j) enddo 11 do14its=1,ITMAX Loop over iterations. iter=itscall linmin(p,xi,n,fret) Next statement is the normal return: if(2.*abs(fret-fp).le.ftol*(abs(fret)+abs(fp)+EPS))return fp=fretcall dfunc(p,xi)gg=0. dgg=0. do 12j=1,n gg=gg+g(j)**2 C dgg=dgg+xi(j)**2 This statement for Fletcher-Reeves. dgg=dgg+(xi(j)+g(j))*xi(j) This statement for Polak-Ribiere. enddo 12 if(gg.eq.0.)return Unlikely. If gradient is exactly zero then we are al- ready done. gam=dgg/gg do13j=1,n g(j)=-xi(j)h(j)=g(j)+gam*h(j) xi(j)=h(j) enddo 13 enddo 14 pause ’frprmn maximum iterations exceeded’returnEND Note on LineMinimizationUsing Derivatives Kindly reread the last part of §10.5. We here want to do the same thing, but using derivative information in performing the line minimization. Rather than reprint the whole routine linminjust to show one modi fied statement, let us just tell you what the change is: The statement fret=brent(ax,xx,bx,f1dim,tol,xmin) should be replaced by fret=dbrent(ax,xx,bx,f1dim,df1dim,tol,xmin) You must also include the following function, which is analogous to f1dimas discussed in §10.5. And remember, your function must be named func, and its gradient calculation must be named dfunc. FUNCTION df1dim(x) INTEGER NMAX REAL df1dim,x PARAMETER (NMAX=50) C USES dfunc INTEGER j,ncomREAL df(NMAX),pcom(NMAX),xicom(NMAX),xt(NMAX)COMMON /f1com/ pcom,xicom,ncomdo 11j=1,ncom xt(j)=pcom(j)+x*xicom(j) enddo 11 call dfunc(xt,df) df1dim=0. 418 Chapter10. MinimizationorMaximizationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).do12j=1,ncom df1dim=df1dim+df(j)*xicom(j) enddo 12 return END CITED REFERENCES AND FURTHER READING: Polak, E. 1971, Computational Methods inOptimization (NewYork: Academic Press), §2.3. [1] Jacobs,D.A.H.(ed.)1977, TheStateoftheArtinNumericalAnalysis (London:AcademicPress), Chapter III.1.7 (by K.W. Brodlie). [2] Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), §8.7. 10.7 Variable Metric Methods in Multidimensions Thegoalof variablemetric methods,whicharesometimescalled quasi-Newton methods,isnotdifferentfromthegoalofconjugategradientmethods: toaccumulate information from successive line minimizations so that Nsuch line minimizations lead to the exact minimum of a quadratic form in Ndimensions. In that case, the methodwill also be quadraticallyconvergentformoregeneralsmoothfunctions. Bothvariablemetricandconjugategradientmethodsrequirethatyouareableto computeyourfunction ’sgradient,or first partialderivatives,atarbitrarypoints. The variablemetricapproachdiffersfromtheconjugategradientinthewaythatit stores and updates the information that is accumulated. Instead of requiring intermediatestorage on the order of N, the number of dimensions, it requires a matrix of size N×N. Generally,for any moderate N, this is an entirely trivial disadvantage. Ontheotherhand,thereisnot,asfarasweknow,anyoverwhelmingadvantage thatthevariablemetricmethodsholdovertheconjugategradienttechniques,except perhapsa historicalone. Developedsomewhatearlier,andmorewidelypropagated,thevariablemetricmethodshavebynowdevelopedawiderconstituencyofsatis fied users. Likewise, some fancier implementations of variable metric methods (going beyondthe scope of this book, see below) have been developedto a greater level ofsophistication on issues like the minimization of roundofferror, handling of special conditions,andso on. Wetendtousevariablemetricratherthanconjugategradient, but we have no reason to urge this habit on you. Variablemetricmethodscomeintwomain flavors. Oneisthe Davidon-Fletcher- Powell (DFP) algorithm (sometimes referred to as simply Fletcher-Powell ). The othergoesbythename Broyden-Fletcher-Goldfarb-Shanno(BFGS) .TheBFGSand DFP schemes differ only in details of their roundoff error, convergence tolerances, and similar “dirty”issues which are outside of our scope [1,2]. However, it has becomegenerallyrecognizedthat,empirically,theBFGSschemeissuperiorinthese details. We will implement BFGS in this section. As before, we imagine that our arbitrary function f(x)can be locally approx- imated by the quadratic form of equation (10.6.1). We don ’t, however, have any