f10-3

10.3One-DimensionalSearchwithFirstDerivatives 399Sample 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.3 One-Dimensional Search with First Derivatives Here we want to accomplish precisely the same goal as in the previous section, namely to isolate a functional minimum that is bracketed by the triplet ofabscissas (a, b, c ), but utilizing an additional capability to compute the function’s first derivative as well as its value. In principle, we might simply search for a zero of the derivative, ignoring the function value information, using a root finder like rtflsporzbrent(§§9.2–9.3). Itdoesn’ttakelongtoreject thatidea: Howdowedistinguishmaximafromminima? Where do we go from initial conditions where the derivatives on one or both of the outer bracketing points indicate that “downhill” is in the direction outof the bracketed interval? We don’twant to give up our strategy of maintaininga rigorousbracket on the minimum at all times. The only way to keep such a bracket is to update it using function (not derivative)information,with the central point in the bracketingtripletalways that with the lowest functionvalue. Thereforethe role of the derivativescan only be to help us choose new trial points within the bracket. Oneschoolofthoughtisto“useeverythingyou’vegot”: Computeapolynomial of relatively high order (cubic or above) that agrees with some number of previous functionandderivativeevaluations. Forexample,thereis a uniquecubicthat agreeswith function and derivative at two points, and one can jump to the interpolated minimum of that cubic (if there is a minimum within the bracket). Suggested by Davidon and others, formulas for this tactic are given in [1]. We like to be more conservative than this. Once superlinear convergence sets in, it hardly matters whether its order is moderately lower or higher. In practical problems that we have met, most function evaluations are spent in getting globally close enoughto the minimumfor superlinearconvergenceto commence. So we are more worried about all the funny “stiff” things that high-order polynomials can do(cf. Figure 3.0.1b), and about their sensitivities to roundoff error. This leads us to use derivative information only as follows: The sign of the derivative at the central point of the bracketing triplet (a, b, c )indicates uniquely whether the next test point should be taken in the interval (a, b )or in the interval (b, c). The value of this derivative and of the derivative at the second-best-so-far point are extrapolated to zero by the secant method (inverse linear interpolation), whichbyitselfissuperlinearoforder1.618. (Thegoldenmeanagain: see [1],p.57.) We imposethe same sortof restrictionsonthis newtrial pointas in Brent’s method.If the trial point must be rejected, we bisectthe interval under scrutiny. Yes,wearefuddy-duddieswhenitcomestomakingflamboyantuseofderivative informationin one-dimensionalminimization. But we havemet toomanyfunctionswhose computed “derivatives” don’tintegrate up to the function value and don’t accurately point the way to the minimum, usually because of roundoff errors, sometimes because of truncationerror in the method of derivativeevaluation. You will see that the following routine is closely modeled on brentin the previous section. 400 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).FUNCTION dbrent(ax,bx,cx,f,df,tol,xmin) INTEGER ITMAX REAL dbrent,ax,bx,cx,tol,xmin,df,f,ZEPS EXTERNAL df,fPARAMETER (ITMAX=100,ZEPS=1.0e-10) Given a function fand its derivative function df, and given a bracketing triplet of abscissas ax,bx,cx[such that bxis between axandcx,a n d f(bx)is less than both f(ax)and f(cx)], this routine isolates the minimum to a fractional precision of about tolusing a modification of Brent’s method that uses derivatives. The abscissa of the minimum is returned as xmin, and the minimum function value is returned as dbrent, the returned function value. INTEGER iterREAL a,b,d,d1,d2,du,dv,dw,dx,e,fu,fv,fw,fx,olde,tol1,tol2, * u,u1,u2,v,w,x,xm Comments following will point out onlydifferences from the routine brent.R e a dt h a t routine first. LOGICAL ok1,ok2 Will be used as flags for whether proposed steps are accept- able or not. a=min(ax,cx) b=max(ax,cx)v=bx w=v x=ve=0. fx=f(x) fv=fxfw=fxdx=df(x) All our housekeeping chores are doubled bythe necessityof moving derivative values around as well as function val- ues.dv=dx dw=dxdo 11iter=1,ITMAX xm=0.5*(a+b) tol1=tol*abs(x)+ZEPStol2=2.*tol1if(abs(x-xm).le.(tol2-.5*(b-a))) goto 3 if(abs(e).gt.tol1) then d1=2.*(b-a) Initialize these d’s to an out-of-bracket value. d2=d1if(dw.ne.dx) d1=(w-x)*dx/(dx-dw) Secant method with one point. if(dv.ne.dx) d2=(v-x)*dx/(dx-dv) And the other. Which of these two estimates of dshall we take? We will insist that theybe within the bracket, and on the side pointed to bythe derivative at x: u1=x+d1 u2=x+d2ok1=((a-u1)*(u1-b).gt.0.).and.(dx*d1.le.0.)ok2=((a-u2)*(u2-b).gt.0.).and.(dx*d2.le.0.) olde=e Movement on the step before last. e=dif(.not.(ok1.or.ok2))then Take onlyan acceptable d,a n di fb o t h are acceptable, then take the small- est one.goto 1 else if (ok1.and.ok2)then if(abs(d1).lt.abs(d2))then d=d1 else d=d2 endif else if (ok1)then d=d1 else d=d2 endif if(abs(d).gt.abs(0.5*olde))goto 1u=x+dif(u-a.lt.tol2 .or. b-u.lt.tol2) d=sign(tol1,xm-x) goto 2 10.3One-DimensionalSearchwithFirstDerivatives 401Sample 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 1 if(dx.ge.0.) then Decide which segment bythe sign of the derivative. e=a-x else e=b-x endif d=0.5*e Bisect, not golden section. 2 if(abs(d).ge.tol1) then u=x+d fu=f(u) else u=x+sign(tol1,d)fu=f(u) if(fu.gt.fx)goto 3 If the minimum step in the downhill direction takes us uphill, then we are done. endif du=df(u) Now all the housekeeping, sigh. if(fu.le.fx) then if(u.ge.x) then a=x else b=x endifv=w fv=fw dv=dww=xfw=fx dw=dx x=ufx=fu dx=du else if(u.lt.x) then a=u else b=u endifif(fu.le.fw .or. w.eq.x) then v=w fv=fwdv=dw w=u fw=fudw=du else if(fu.le.fv .or. v.eq.x .or. v.eq.w) then v=u fv=fudv=du endif endif enddo 11 pause ’dbrent exceeded maximum iterations’ 3 xmin=x dbrent=fxreturn END CITED REFERENCES AND FURTHER READING: Acton, F.S. 1970, Numerical Methods That Work ; 1990, corrected edition (Washington: Mathe- matical Association of America), pp. 55; 454–458. [1] Brent,R.P.1973, AlgorithmsforMinimizationwithoutDerivatives (EnglewoodCliffs,NJ:Prentice- Hall), p. 78. 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 with N+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,