f10-4

402 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).10.4 Downhill Simplex Method in Multidimensions With this section we begin consideration of multidimensional minimization, that is, finding the minimum of a function of more than one independent variable.This section stands apart from those which follow, however: All of the algorithms afterthissectionwillmakeexplicituseofaone-dimensionalminimizationalgorithm as a part of their computational strategy. This section implements an entirely self-containedstrategy,in which one-dimensionalminimizationdoes not figure. Thedownhill simplex method is due to Nelder and Mead [1]. The method requires only function evaluations, not derivatives. It is not very efficient in terms of the number of function evaluations that it requires. Powell’s method ( §10.5) is almostsurelyfasterinalllikelyapplications. However,thedownhillsimplexmethodmay frequently be the bestmethod to use if the figure of merit is “get something working quickly” for a problem whose computational burden is small. The method has a geometrical naturalness about it which makes it delightful to describe or work through: Asimplexis the geometrical figure consisting, in Ndimensions, of N+1 points (orvertices) and all their interconnectingline segments, polygonalfaces, etc. In two dimensions, a simplex is a triangle. In three dimensions it is a tetrahedron, notnecessarilytheregulartetrahedron. (The simplexmethod oflinearprogramming, describedin §10.8,alsomakesuseofthegeometricalconceptofasimplex. Otherwise it is completelyunrelatedto the algorithmthat we are describingin this section.) In generalwe are onlyinterestedin simplexesthat are nondegenerate,i.e., that enclosea finite inner N-dimensional volume. If any point of a nondegenerate simplex is taken as the origin, then the Nother points define vector directions that span the N-dimensional vector space. Inone-dimensionalminimization,itwaspossibletobracketaminimum,sothat the success of a subsequent isolation was guaranteed. Alas! There is no analogousprocedure in multidimensional space. For multidimensional minimization, the best we candois giveouralgorithmastartingguess,thatis, an N-vectorofindependent variablesasthefirstpointtotry. Thealgorithmisthensupposedtomakeitsownwaydownhill through the unimaginable complexity of an N-dimensional topography, until it encounters a (local, at least) minimum. The downhill simplex method must be started not just with a single point, but withN+1points, defining an initial simplex. If you think of one of these points (it matters not which) as being your initial starting point P 0, then you can take the other Npoints to be Pi=P0+λei (10.4.1 ) where the ei’s are Nunit vectors, and where λis a constant which is your guess of the problem’s characteristic length scale. (Or, you could have different λi’s for each vector direction.) Thedownhillsimplexmethodnowtakesaseriesofsteps,moststepsjustmoving the point of the simplex where the function is largest (“highest point”) through the opposite face of the simplex to a lower point. These steps are called reflections, 10.4DownhillSimplexMethodinMultidimensions 403Sample 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).simplex at beginning of step reflection reflection and expansion contraction multiple contraction(a) (b) (c) (d)high low Figure 10.4.1. Possible outcomes for a step in the downhill simplex method. The simplex at the beginning ofthe step,here a tetrahedron, isshown,top. The simplex at the endofthe step can be any one of (a) a re flection away from the high point, (b) a re flection and expansion away from the high point, (c) a contraction along one dimension from the high point, or (d) a contraction along all dimensions towardsthelowpoint. Anappropriate sequence ofsuchstepswillalways converge toaminimumofthefunction. and they are constructed to conserve the volume of the simplex (hence maintain its nondegeneracy). When it can do so, the method expands the simplex in one or another direction to take larger steps. When it reaches a “valleyfloor,”the method contracts itself in the transversedirection and tries to ooze downthe valley. If thereis a situation where the simplex is trying to “pass through the eye of a needle, ”it contracts itself in all directions, pulling itself in around its lowest (best) point. The routinename amoebaisintendedtobedescriptiveofthiskindofbehavior;thebasic moves are summarized in Figure 10.4.1. Termination criteria can be delicate in any multidimensional minimization routine. Without bracketing, and with more than one independent variable, we no longer have the option of requiring a certain tolerance for a single independent 404 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).variable. We typically can identify one “cycle”or“step”of our multidimensional algorithm. It is then possible to terminate when the vector distance moved in thatstep is fractionally smaller in magnitude than some tolerance tol. Alternatively, we could require that the decrease in the function value in the terminating step be fractionally smaller than some tolerance ftol. Note that while tolshould not usually be smaller than the square root of the machine precision, it is perfectly appropriateto let ftolbe of orderthe machineprecision(orperhapsslightly larger so as not to be diddled by roundoff). Notewellthateitheroftheabovecriteriamightbefooledbyasingleanomalous stepthat,foronereasonoranother,failedtogetanywhere. Therefore,itisfrequentlya good idea to restarta multidimensional minimization routine at a point where it claims to have found a minimum. For this restart, you should reinitialize any ancillaryinputquantities. In the downhillsimplexmethod,forexample,youshouldreinitialize Nof the N+1vertices of the simplex again by equation (10.4.1),with P 0being one of the vertices of the claimed minimum. Restarts shouldneverbeveryexpensive;youralgorithmdid,afterall,converge to the restart point once, and now you are starting the algorithmalready there. Consider, then, our N-dimensional amoeba: SUBROUTINE amoeba(p,y,mp,np,ndim,ftol,funk,iter) INTEGER iter,mp,ndim,np,NMAX,ITMAXREAL ftol,p(mp,np),y(mp),funk,TINY PARAMETER (NMAX=20,ITMAX=5000,TINY=1.e-10) Maximum allowed dimensions and func- tion evaluations, and a small num-ber.EXTERNAL funk C USES amotry,funk Multidimensional minimization of the function funk(x) where x(1:ndim) is a vector inndim dimensions, by the downhill simplex method of Nelder and Mead. The matrix p(1:ndim+1,1:ndim) is input. Its ndim+1 rows are ndim -dimensional vectors which are the vertices of the starting simplex. Also input is the vector y(1:ndim+1) , whose compo- nents must be pre-initialized to the values of funk evaluated at the ndim+1 vertices (rows) ofp;a n d ftol the fractional convergence tolerance to be achieved in the function value (n.b.!). On output, pandywill have been reset to ndim+1 new points all within ftol of a minimum function value, and iter gives the number of function evaluations taken. INTEGER i,ihi,ilo,inhi,j,m,nREAL rtol,sum,swap,ysave,ytry,psum(NMAX),amotryiter=0 1d o 12n=1,ndim Enter here when starting or have just overall contracted. sum=0. Recompute psum . do11m=1,ndim+1 sum=sum+p(m,n) enddo 11 psum(n)=sum enddo 12 2 ilo=1 Enter here when have just changed a single point. if (y(1).gt.y(2)) then Determine which point is the highest (worst), next-highest, and lowest (best), ihi=1 inhi=2 else ihi=2inhi=1 endif do 13i=1,ndim+1 by looping over the points in the simplex. if(y(i).le.y(ilo)) ilo=iif(y(i).gt.y(ihi)) then inhi=ihi ihi=i else if(y(i).gt.y(inhi)) then if(i.ne.ihi) inhi=i 10.4DownhillSimplexMethodinMultidimensions 405Sample 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).endif enddo 13 rtol=2.*abs(y(ihi)-y(ilo))/(abs(y(ihi))+abs(y(ilo))+TINY) Compute the fractional range from highest to lowest and return if satisfactory. if (rtol.lt.ftol) then If returning, put best point and value in slot 1. swap=y(1) y(1)=y(ilo)y(ilo)=swapdo 14n=1,ndim swap=p(1,n) p(1,n)=p(ilo,n)p(ilo,n)=swap enddo 14 return endifif (iter.ge.ITMAX) pause ’ITMAX exceeded in amoeba’ iter=iter+2 Begin a new iteration. First extrapolate by a factor −1through the face of the simplex across from the high point, i.e., reflect the simplex from the high point. ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,-1.0) if (ytry.le.y(ilo)) then Gives a result better than the best point, so try an additional extrapolation by a factor 2. ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,2.0) else if (ytry.ge.y(inhi)) then The reflected point is worse than the second-highest, so look for an intermediate lower point,i.e., do a one-dimensional contraction. ysave=y(ihi) ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,0.5) if (ytry.ge.ysave) then Can’t seem to get rid of that high point. Better contract around the lowest (best) point. do 16i=1,ndim+1 if(i.ne.ilo)then do15j=1,ndim psum(j)=0.5*(p(i,j)+p(ilo,j))p(i,j)=psum(j) enddo 15 y(i)=funk(psum) endif enddo 16 iter=iter+ndim Keep track of function evaluations. goto 1 Go back for the test of doneness and the next iteration. endif else iter=iter-1 Correct the evaluation count. endifgoto 2 END FUNCTION amotry(p,y,psum,mp,np,ndim,funk,ihi,fac) INTEGER ihi,mp,ndim,np,NMAXREAL amotry,fac,p(mp,np),psum(np),y(mp),funkPARAMETER (NMAX=20) EXTERNAL funk C USES funk Extrapolates by a factor fac through the face of the simplex across from the high point, tries it, and replaces the high point if the new point is better. INTEGER jREAL fac1,fac2,ytry,ptry(NMAX)fac1=(1.-fac)/ndim fac2=fac1-fac do 11j=1,ndim ptry(j)=psum(j)*fac1-p(ihi,j)*fac2 enddo 11 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.