f10-5

406 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).ytry=funk(ptry) Evaluate the function at the trial point. if (ytry.lt.y(ihi)) then If it’s better than the highest, then replace the highest. y(ihi)=ytry do12j=1,ndim psum(j)=psum(j)-p(ihi,j)+ptry(j) p(ihi,j)=ptry(j) enddo 12 endifamotry=ytry return END CITED REFERENCES AND FURTHER READING: Nelder, J.A., and Mead, R. 1965, Computer Journal , vol. 7, pp. 308–313. [1] Yarbro, L.A., and Deming, S.N. 1974, Analytica Chimica Acta , vol. 73, pp. 391–398. Jacoby, S.L.S, Kowalik, J.S., and Pizzo, J.T. 1972, Iterative Methods for Nonlinear Optimization Problems (Englewood Cliffs, NJ: Prentice-Hall). 10.5 Direction Set (Powell’s) Methods in Multidimensions We know ( §10.1–§10.3) how to minimize a function of one variable. If we start at a point PinN-dimensional space, and proceed from there in some vector direction n, then any function of Nvariables f(P)can be minimized along the line nby our one-dimensional methods. One can dream up various multidimensional minimizationmethodsthatconsistofsequencesofsuchlineminimizations. Differentmethods will differ only by how, at each stage, they choose the next direction nto try. All such methodspresumethe existenceof a “black-box”sub-algorithm,which we mightcall linmin(givenas anexplicitroutineatthe endofthis section),whose definition can be taken for now as linmin: Given as input the vectors Pandn, and the function f,findthescalar λthatminimizes f(P+λn). ReplacePbyP+λn. Replace nbyλn. Done. All the minimization methods in this section and in the two sections following fall under this general schema of successive line minimizations. (The algorithm in§10.7 does not need very accurate line minimizations. Accordingly, it has its own approximate line minimization routine, lnsrch.) In this section we consider a class of methods whose choice of successive directions does not involve explicit computationofthefunction’sgradient;thenexttwosectionsdorequiresuchgradient calculations. You will note that we need not specify whether linminuses gradient information or not. That choice is up to you, and its optimization depends on your particular function. You would be crazy, however, to use gradients in linminand notuse them in the choice of directions, since in this latter role they can drastically reduce the total computational burden. 10.5DirectionSet(Powell’s)MethodsinMultidimensions 407Sample 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).starty x Figure 10.5.1. Successive minimizations along coordinate directions in a long, narrow “valley”(shown as contour lines). Unless the valley is optimally oriented, this method is extremely inef ficient, taking many tiny steps to get to the minimum, crossing and re-crossing the principal axis. Butwhatif,inyourapplication,calculationofthegradientisoutofthequestion. You might first think of this simple method: Take the unit vectors e1,e2,...eNas aset of directions . Using linmin, move along the first direction to its minimum, thenfrom there along the second direction to itsminimum, and so on, cycling through the whole set of directions as many times as necessary, until the function stops decreasing. This simple method is actually not too bad for many functions. Even more interesting is why it isbad, i.e. very inef ficient, for some other functions. Consider a function of two dimensions whose contour map (level lines) happens to de fine a long,narrowvalleyatsomeangletothecoordinatebasisvectors(seeFigure10.5.1). Then the only way “down the length of the valley ”going along the basis vectors at each stage is by a series of many tiny steps. More generally, in Ndimensions, if the function ’s second derivatives are much larger in magnitude in some directions than in others, then many cycles through all Nbasis vectors will be required in ordertoget anywhere. Thisconditionis notall thatunusual;accordingtoMurphy ’s Law, you should count on it. Obviouslywhatwe needis a betterset of directionsthanthe ei’s. Alldirection set methods consist of prescriptions for updatingthe set of directions as the method proceeds, attempting to come up with a set which either (i) includes some very 408 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).good directions that will take us far along narrow valleys, or else (more subtly) (ii) includes some number of “non-interfering ”directions with the special property that minimization along one is not “spoiled”by subsequent minimization along another,so that interminablecyclingthroughthe set ofdirectionscan be avoided. ConjugateDirections This concept of “non-interfering ”directions, more conventionally called con- jugate directions , is worth making mathematically explicit. First, note that if we minimize a function along some direction u, then the gradientofthefunctionmustbeperpendicularto uatthe lineminimum;ifnot,then there would still be a nonzero directional derivative along u. Next take some particular point Pas the origin of the coordinate system with coordinates x. Then any function fcan be approximatedby its Taylor series f(x)=f(P)+/summationdisplay i∂f ∂x ixi+1 2/summationdisplay i,j∂2f ∂x i∂x jxixj+··· ≈c−b·x+1 2x·A·x(10.5.1 ) where c≡f(P)b≡− ∇ f|P[A]ij≡∂2f ∂x i∂x j/vextendsingle/vextendsingle/vextendsingle/vextendsingleP(10.5.2 ) The matrix Awhose components are the second partial derivative matrix of the function is called the Hessian matrix of the function at P. In the approximationof (10.5.1),the gradientof fis easily calculated as ∇f=A·x−b (10.5.3 ) (Thisimpliesthatthegradientwillvanish —thefunctionwill beatanextremum — at a valueof xobtainedbysolving A·x=b. This ideawe will returnto in §10.7!) Howdoesthegradient ∇fchangeaswemovealongsomedirection? Evidently δ(∇f)=A·(δx)( 10.5.4 ) Suppose that we have moved along some direction uto a minimum and now proposetomovealongsomenewdirection v. Theconditionthatmotionalong vnot spoilour minimization along uis just that the gradient stay perpendicularto u, i.e., thatthechangeinthegradientbeperpendicularto u. Byequation(10.5.4)thisisjust 0=u·δ(∇f)=u·A·v (10.5.5 ) When (10.5.5) holds for two vectors uandv, they are said to be conjugate . When the relation holds pairwise for all members of a set of vectors, they are said to be a conjugate set. If you do successive line minimization of a function along a conjugate set of directions, then you don ’t need to redo any of those directions 10.5DirectionSet(Powell’s)MethodsinMultidimensions 409Sample 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).(unless, of course, you spoil things by minimizing along a direction that they are notconjugate to). A triumph for a direction set method is to come up with a set of Nlinearly independent,mutuallyconjugatedirections. Then,onepassof Nlineminimizations will put it exactly at the minimum of a quadratic form like (10.5.1). For functionsfthat are not exactly quadratic forms, it won ’t be exactly at the minimum; but repeated cycles of Nline minimizations will in due course converge quadratically to the minimum. Powell’sQuadraticallyConvergent Method Powellfirst discovered a direction set method that does produce Nmutually conjugate directions. Here is how it goes: Initialize the set of directions uito the basis vectors, ui=ei i=1,...,N (10.5.6 ) Now repeat the followingsequence of steps ( “basic procedure ”)until your function stops decreasing: •Save your starting position as P0. •Fori=1,...,N, movePi−1to the minimum along direction uiand call this point Pi. •Fori=1,...,N −1, setui←ui+1. •SetuN←PN−P0. •MovePNto the minimumalongdirection uNandcall this point P0. Powell, in 1964, showed that, for a quadratic form like (10.5.1), kiterations of the above basic procedure produce a set of directions uiwhose last kmembers are mutually conjugate. Therefore, Niterations of the basic procedure, amounting toN(N+1 )line minimizations in all, will exactly minimize a quadratic form. Brent[1]gives proofs of these statements in accessible form. Unfortunately, there is a problem with Powell ’s quadratically convergent al- gorithm. The procedure of throwing away, at each stage, u1in favor of PN−P0 tends to producesets of directions that “fold up on each other ”and becomelinearly dependent. Oncethishappens,thentheprocedure findstheminimumofthefunction fonly over a subspace of the full N-dimensional case; in other words, it gives the wronganswer. Therefore,the algorithmmustnot beusedin the formgivenabove. There are a number of ways to fix up the problem of linear dependence in Powell’s algorithm, among them: 1. Youcanreinitializethe set ofdirections uito thebasis vectors eiafterevery NorN+1iterations of the basic procedure. This produces a serviceable method, whichwecommendtoyouifquadraticconvergenceisimportantforyourapplication (i.e.,ifyourfunctionsare closetoquadraticformsandifyoudesirehighaccuracy). 2. Brent points out that the set of directions can equally well be reset to the columns of any orthogonal matrix. Rather than throw away the information on conjugate directions already built up, he resets the direction set to calculated principal directions of the matrix A(which he gives a procedure for determining). 410 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).The calculation is essentially a singular value decomposition algorithm (see §2.6). Brent has a number of other cute tricks up his sleeve, and his modi fication of Powell’s method is probably the best presently known. Consult [1]for a detailed description and listing of the program. Unfortunately it is rather too elaborate for us to include here. 3. You can give up the property of quadratic convergence in favor of a more heuristic scheme (due to Powell) which tries to find a few good directions along narrow valleys instead of Nnecessarily conjugate directions. This is the method thatwenowimplement. (ItisalsotheversionofPowell ’smethodgiveninActon [2], from which parts of the following discussion are drawn.) Discarding the Directionof Largest Decrease The fox and the grapes: Now that we are going to give up the property of quadratic convergence, was it so important after all? That depends on the functionthat you are minimizing. Some applications produce functions with long, twisty valleys. Quadratic convergence is of no particular advantage to a program which must slalom down the length of a valley floor that twists one way and another (and another, and another, ...–there are Ndimensions!). Along the long direction, a quadratically convergent method is trying to extrapolate to the minimum of aparabola which just isn ’t (yet) there; while the conjugacy of the N−1transverse directions keeps getting spoiled by the twists. Soonerorlater,however,wedoarriveatanapproximatelyellipsoidalminimum (cf. equation 10.5.1 when b, the gradient, is zero). Then, depending on how much accuracywerequire,amethodwithquadraticconvergencecansaveusseveraltimes N 2extra line minimizations, since quadratic convergence doublesthe number of significantfigures at each iteration. Thebasicideaofournow-modi fiedPowell ’smethodisstilltotake PN−P0as anewdirection;itis,afterall,theaveragedirectionmovedaftertryingall Npossible directions. For a valley whose long direction is twisting slowly, this direction is likely to give us a good run along the new long direction. The change is to discardthe old direction along which the function fmade itslargest decrease . This seems paradoxical, since that direction was the bestof the previous iteration. However, it is also likely to be a major component of the new direction that we are adding, so droppingit gives us the best chance of avoidinga buildupof linear dependence. There are a couple of exceptions to this basic idea. Sometimes it is better not to add a new direction at all. De fine f 0≡f(P0) fN≡f(PN) fE≡f(2PN−P0)( 10.5.7 ) Here fEis the function value at an “extrapolated ”point somewhat further along the proposed new direction. Also de fine∆fto be the magnitude of the largest decreasealongoneparticulardirectionofthepresentbasicprocedureiteration. ( ∆f is a positive number.) Then: 1. If fE≥f0, then keep the old set of directions for the next basic procedure, because the average direction PN−P0is all played out. 2. If 2(f0−2fN+fE)[ (f0−fN)−∆f]2≥(f0−fE)2∆f,thenkeeptheold set of directions for the next basic procedure, because either (i) the decrease along 10.5DirectionSet(Powell’s)MethodsinMultidimensions 411Sample 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).theaveragedirectionwas notprimarilyduetoanysingledirection ’sdecrease,or(ii) there is a substantial second derivative along the average direction and we seem tobe near to the bottom of its minimum. ThefollowingroutineimplementsPowell ’smethodintheversionjustdescribed. Intheroutine, xiisthematrixwhosecolumnsarethesetofdirections n i;otherwise the correspondence of notation should be self-evident. SUBROUTINE powell(p,xi,n,np,ftol,iter,fret) INTEGER iter,n,np,NMAX,ITMAX REAL fret,ftol,p(np),xi(np,np),func,TINYEXTERNAL func PARAMETER (NMAX=20,ITMAX=200,TINY=1.e-25) C USES func,linmin Minimization of a function func ofnvariables. ( func is not an argument, it is a fixed func- tion name.) Input consists of an initial starting point p(1:n) ; an initial matrix xi(1:n,1:n) with physical dimensions npbynp, and whose columns contain the initial set of directions (usually the nunit vectors); and ftol , the fractional tolerance in the function value such that failure to decrease by more than this amount on one iteration signals doneness. Onoutput, pis set to the best point found, xiis the then-current direction set, fret is the returned function value at p,a n d iter is the number of iterations taken. The routine linmin is used. Parameters: Maximum value of n, maximum allowed iterations, and a small number. INTEGER i,ibig,j REAL del,fp,fptt,t,pt(NMAX),ptt(NMAX),xit(NMAX)fret=func(p)do 11j=1,n Save the initial point. pt(j)=p(j) enddo 11 iter=0 1 iter=iter+1 fp=fretibig=0del=0. Will be the biggest function decrease. do 13i=1,n In each iteration, loop over all directions in the set. do12j=1,n Copy the direction, xit(j)=xi(j,i) enddo 12 fptt=fret call linmin(p,xit,n,fret) minimize along it, if(fptt-fret.gt.del)then and record it if it is the largest decrease so far. del=fptt-fret ibig=i endif enddo 13 if(2.*(fp-fret).le.ftol*(abs(fp)+abs(fret))+TINY)return Termination criterion. if(iter.eq.ITMAX) pause ’powell exceeding maximum iterations’do 14j=1,n Construct the extrapolated point and the average di- rection moved. Save the old starting point. ptt(j)=2.*p(j)-pt(j) xit(j)=p(j)-pt(j)pt(j)=p(j) enddo 14 fptt=func(ptt) Function value at extrapolated point. if(fptt.ge.fp)goto 1 One reason not to use new direction. t=2.*(fp-2.*fret+fptt)*(fp-fret-del)**2-del*(fp-fptt)**2 if(t.ge.0.)goto 1 Other reason not to use new direction. call linmin(p,xit,n,fret) Move to the minimum of the new direction, do15j=1,n and save the new direction. xi(j,ibig)=xi(j,n) xi(j,n)=xit(j) enddo 15 goto 1 Back for another iteration. END 412 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).ImplementationofLine Minimization Intheaboveroutine,youmighthavewonderedwhywedidn ’tmakethefunction name funcan argument of the routine. The reason is buried in a slightly dirty FORTRAN practicality in our implementation of linmin. Make no mistake, there is a rightway to implement linmin: I ti st ou s e themethods of one-dimensional minimization described in §10.1–§10.3, but to rewrite the programs of those sections so that their bookkeeping is done on vector- valued points P(all lying along a given direction n) rather than scalar-valued abscissas x. That straightforward task produces long routines densely populated with“do k=1,n ”loops. Wedonothavespacetoincludesuchroutinesinthisbook. Our linmin,which worksjust fine,isinsteadakindofbookkeepingswindle. Itconstructsan “artificial” function of one variable called f1dim, which is the value of your function func along the line going throughthe point pin the direction xi.linmincommunicates with f1dimthrough a common block. It then calls our familiar one-dimensional routines mnbrak(§10.1)and brent(§10.2)andinstructsthemto minimize f1dim. Stillfollowing? Thentrythis: brentreceivesthefunctionname f1dim,which itdutifullycalls. Butthereisnowaytosignalto f1dimthatitissupposedtouseyour functionname,whichcouldhavebeenpassedto linminasanargument. Therefore, we have to make f1dimuse afixedfunction name, namely func. The situation is reminiscentofHenryFord ’sblackautomobile: powellwill minimizeanyfunction, as long as it is named func. Needed to remedy this situation is a way to pass a function name through a common block; this is lacking in FORTRAN. Theonlythinginef ficientabout linministhis: Itsuseasaninterfacebetweena multidimensionalminimizationstrategyandaone-dimensionalminimizationroutine results in some unnecessary copying of vectors hither and yon. That should notnormallybeasigni ficantadditiontotheoverallcomputationalburden,butwecannot disguise its inelegance. SUBROUTINE linmin(p,xi,n,fret) INTEGER n,NMAX REAL fret,p(n),xi(n),TOL PARAMETER (NMAX=50,TOL=1.e-4) Maximum anticipated n,a n d TOLpassed to brent . C USES brent,f1dim,mnbrak Given an n-dimensional point p(1:n) and an n-dimensional direction xi(1:n) ,m o v e sa n d resets pto where the function func(p) takes on a minimum along the direction xifrom p, and replaces xiby the actual vector displacement that pwas moved. Also returns as fret the value of func at the returned location p. This is actually all accomplished by calling the routines mnbrak andbrent . INTEGER j,ncomREAL ax,bx,fa,fb,fx,xmin,xx,pcom(NMAX),xicom(NMAX),brentCOMMON /f1com/ pcom,xicom,ncom EXTERNAL f1dim ncom=n Set up the common block. do 11j=1,n pcom(j)=p(j) xicom(j)=xi(j) enddo 11 ax=0. Initial guess for brackets. xx=1. call mnbrak(ax,xx,bx,fa,fx,fb,f1dim)fret=brent(ax,xx,bx,f1dim,TOL,xmin) do 12j=1,n Construct the vector results to return. 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.