f10-1

390 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).one-dimensionalsub-minimization. Turn to §10.6 for detailed discussion and implementation. •Thesecondfamilygoesunderthenames quasi-Newton orvariablemetric methods, as typified by the Davidon-Fletcher-Powell (DFP) algorithm (sometimes referred to just as Fletcher-Powell ) or the closely related Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm. These methods require of order N2storage, require derivative calculations and one- dimensional sub-minimization. Details are in §10.7. You are now ready to proceed with scaling the peaks (and/or plumbing the depths) of practical optimization. CITED REFERENCES AND FURTHER READING: Dennis,J.E., andSchnabel,R.B. 1983, NumericalMethods forUnconstrained Optimizationand Nonlinear Equations (Englewood Cliffs, NJ: Prentice-Hall). Polak, E. 1971, Computational Methods in Optimization (New York: Academic Press). Gill,P.E.,Murray,W.,andWright,M.H.1981, PracticalOptimization (NewYork:AcademicPress). Acton, F.S. 1970, Numerical Methods That Work ; 1990, corrected edition (Washington: Mathe- matical Association of America), Chapter 17. Jacobs, D.A.H. (ed.) 1977, The State of the Art in Numerical Analysis (London: Academic Press), Chapter III.1. Brent,R.P.1973, AlgorithmsforMinimizationwithoutDerivatives (EnglewoodCliffs,NJ:Prentice- Hall). Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall), Chapter 10. 10.1 Golden Section Search in One Dimension Recall how the bisection method finds roots of functions in one dimension (§9.1): The root is supposed to have been bracketed in an interval (a, b ). One then evaluates the function at an intermediate point xand obtains a new, smaller bracketinginterval,either (a, x )or(x, b ). Theprocesscontinuesuntilthebracketing interval is acceptably small. It is optimal to choose xto be the midpoint of (a, b ) so that the decrease in the interval length is maximized when the function is asuncooperative as it can be, i.e., when the luck of the draw forces you to take the bigger bisected segment. There is a precise, thoughslightly subtle, translation of these considerationsto the minimization problem: What does it mean to bracketa minimum? A root of a function is known to be bracketed by a pair of points, aandb, when the function has opposite sign at those two points. A minimum, by contrast, is known to be bracketedonlywhen thereis a tripletof points, a<b<c (orc<b<a ), such that f(b)is less than both f(a)andf(c). In this case we know that the function (if it is nonsingular) has a minimum in the interval (a, c ). The analog of bisection is to choose a new point x, either between aandbor between bandc. Suppose, to be specific, that we make the latter choice. Then we evaluate f(x).I ff(b)<f (x), then the new bracketing triplet of points is (a, b, x ); 10.1GoldenSectionSearchinOneDimension 391Sample 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).644 6 351 52 Figure 10.1.1. Successive bracketing of a minimum. The minimum is originally bracketed by points 1,3,2. The function is evaluated at 4, which replaces 2; then at 5, which replaces 1; then at 6, which replaces 4. The rule at each stage is to keep a center point that is lower than the two outside points. After the steps shown, the minimum is bracketed by points 5,3,6. contrariwise,if f(b)>f (x), thenthe new bracketingtriplet is (b, x, c ). In all cases themiddlepointofthenewtripletistheabscissawhoseordinateisthebestminimumachieved so far; see Figure 10.1.1. We continue the process of bracketing until the distance between the two outer points of the triplet is tolerably small. How small is “tolerably”small? For a minimum located at a value b, you might naively think that you will be able to bracket it in as small a range as (1−/epsilon1)b<b< (1 + /epsilon1)b, where /epsilon1is your computer ’sfloating-point precision, a number like 3×10 −8(single precision) or 10−15(double precision). Not so! In general,the shapeof yourfunction f(x)nearbwill be givenby Taylor ’stheorem f(x)≈f(b)+1 2f/prime/prime(b)(x−b)2(10.1.1 ) The second term will be negligible compared to the first (that is, will be a factor /epsilon1 smaller and will act just like zero when added to it) whenever |x−b|<√/epsilon1|b|/radicalBigg 2|f(b)| b2f/prime/prime(b)(10.1.2 ) The reason for writing the right-hand side in this way is that, for most functions, thefinal square root is a number of order unity. Therefore, as a rule of thumb, it is hopeless to ask for a bracketing interval of width less than√/epsilon1times its central value, a fractional width of only about 10−4(single precision) or 3×10−8(double precision). Knowingthis inescapablefact will save youa lot of useless bisections! Theminimum- findingroutinesofthischapterwilloftencallforauser-supplied argument tol,andreturnwithanabscissawhosefractionalprecisionis about ±tol (bracketing interval of fractional size about 2×tol). Unless you have a better 392 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).estimate for the right-hand side of equation (10.1.2), you should set tolequal to (notmuchlessthan)thesquarerootofyourmachine ’sfloating-pointprecision,since smaller values will gain you nothing. It remains to decide on a strategy for choosingthe new point x,g i v e n (a, b, c ). Suppose that bis a fraction wof the way between aandc, i.e. b−a c−a=wc−b c−a=1−w (10.1.3 ) Also suppose that our next trial point xis an additional fraction zbeyond b, x−b c−a=z (10.1.4 ) Thenthenextbracketingsegmentwilleitherbeoflength w+zrelativetothecurrent one,orelse oflength 1−w. Ifwe want tominimizethe worstcase possibility,then we will choose zto make these equal, namely z=1−2w (10.1.5 ) We see atoncethatthenewpointis thesymmetricpointto bintheoriginalinterval, namely with |b−a|equal to |x−c|. This implies that the point xlies in the larger of the two segments ( zis positive only if w< 1/2). But where in the larger segment? Where did the value of witself come from? Presumably from the previous stage of applying our same strategy. Therefore, if z is chosen to be optimal, then so was wbefore it. This scale similarity implies that xshould be the same fraction of the way from btoc(if that is the bigger segment) as was bfrom atoc, in other words, z 1−w=w (10.1.6 ) Equations (10.1.5) and (10.1.6) give the quadratic equation w2−3w+1=0 yielding w=3−√ 5 2≈0.38197 ( 10.1.7 ) In other words, the optimal bracketinginterval (a, b, c )has its middle point ba fractional distance 0.38197 from one end (say, a), and 0.61803 from the other end (say, b). These fractions are those of the so-called golden mean orgolden section , whose supposedly aesthetic properties hark back to the ancient Pythagoreans. This optimal method of function minimization, the analog of the bisection method for finding zeros, is thus called the goldensection search , summarizedas follows: Given, at each stage, a bracketing triplet of points, the next point to be tried is that which is a fraction 0.38197 into the larger of the two intervals (measuring fromthe centralpoint of the triplet). If you start out with a bracketingtriplet whose segments are not in the golden ratios, the procedure of choosing successive pointsat the golden mean point of the larger segment will quickly converge you to the proper, self-replicating ratios. The golden section search guarantees that each new function evaluation will (afterself-replicatingratios havebeenachieved)brackettheminimumto aninterval 10.1GoldenSectionSearchinOneDimension 393Sample 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).just 0.61803times the size of the preceding interval. This is comparable to, but not quite as good as, the 0.50000 that holds when finding roots by bisection. Note that the convergence is linear(in the language of Chapter 9), meaning that successive significantfigures are won linearly with additional function evaluations. In the next section we will give a superlinear method, where the rate at which successivesignificantfiguresare liberatedincreases with eachsuccessivefunctionevaluation. Routinefor InitiallyBracketing a Minimum Theprecedingdiscussionhasassumedthatyouareabletobrackettheminimum in thefirst place. We consider this initial bracketing to be an essential part of any one-dimensional minimization. There are some one-dimensional algorithms that do not require a rigorous initial bracketing. However, we would nevertrade the secure feeling of knowingthat a minimum is “in there somewhere ”for the dubious reduction of function evaluations that these nonbracketing routines may promise. Please bracketyourminima(or,forthat matter,yourzeros)beforeisolating them! There is not much theory as to how to do this bracketing. Obviouslyyou want to step downhill. But how far? We like to take larger and largersteps, starting with some (wild?) initial guess and then increasing the stepsize at each step either by a constant factor, or else by the result of a parabolic extrapolation of the preceding points that is designed to take us to the extrapolated turning point. It doesn ’t much matter if the steps get big. After all, we are stepping downhill, so we already have theleftandmiddlepointsofthebracketingtriplet. Wejustneedtotakeabigenough step to stop the downhill trend and get a high third point. Our standard routine is this: SUBROUTINE mnbrak(ax,bx,cx,fa,fb,fc,func) REAL ax,bx,cx,fa,fb,fc,func,GOLD,GLIMIT,TINY EXTERNAL funcPARAMETER (GOLD=1.618034, GLIMIT=100., TINY=1.e-20) Given a function func , and given distinct initial points axandbx, this routine searches in the downhill direction (defined by the function as evaluated at the initial points) andreturns new points ax,bx,cxthat bracket a minimum of the function. Also returned are the function values at the three points, fa,fb,a n d fc. Parameters: GOLD is the default ratio by which successive intervals are magnified; GLIMIT is the maximum magnification allowed for a parabolic-fit step. REAL dum,fu,q,r,u,ulimfa=func(ax) fb=func(bx) if(fb.gt.fa)then Switch roles of aand bso that we can go downhill in the direction from ato b. dum=ax ax=bx bx=dumdum=fbfb=fa fa=dum endifcx=bx+GOLD*(bx-ax) First guess for c. fc=func(cx) 1 if(fb.ge.fc)then “do while ”: keep returning here until we bracket. r=(bx-ax)*(fb-fc) Compute uby parabolic extrapolation from a, b, c .TINY is used to prevent any possible division by zero. q=(bx-cx)*(fb-fa) u=bx-((bx-cx)*q-(bx-ax)*r)/(2.*sign(max(abs(q-r),TINY),q-r)) ulim=bx+GLIMIT*(cx-bx) We won’t go farther than this. Test various possibilities: if((bx-u)*(u-cx).gt.0.)then Parabolic uis between band c:t r y i t . fu=func(u) 394 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).if(fu.lt.fc)then Got a minimum between band c. ax=bx fa=fb bx=ufb=fu return else if(fu.gt.fb)then Got a minimum between between aand u. cx=ufc=fu return endifu=cx+GOLD*(cx-bx) Parabolic fit was no use. Use default magnification. fu=func(u) else if((cx-u)*(u-ulim).gt.0.)then Parabolic fit is between cand its allowed limit. fu=func(u) if(fu.lt.fc)then bx=cx cx=uu=cx+GOLD*(cx-bx)fb=fc fc=fu fu=func(u) endif else if((u-ulim)*(ulim-cx).ge.0.)then Limit parabolic uto maximum allowed value. u=ulim fu=func(u) else Reject parabolic u, use default magnification. u=cx+GOLD*(cx-bx) fu=func(u) endif ax=bx Eliminate oldest point and continue. bx=cxcx=ufa=fb fb=fc fc=fugoto 1 endif return END (Because of the housekeeping involved in moving around three or four points and their function values, the above program ends up looking deceptively formidable.That is true of several otherprogramsin this chapteras well. The underlyingideas, however, are quite simple.) Routinefor GoldenSection Search FUNCTION golden(ax,bx,cx,f,tol,xmin) REAL golden,ax,bx,cx,tol,xmin,f,R,C EXTERNAL fPARAMETER (R=.61803399,C=1.-R) Given a function f, 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) andf(cx) ), this routine performs a golden section search for the minimum, isolating it to a fractional precision of about tol. The abscissa of the minimum is returned as xmin , and the minimum function value is returned as golden , the returned function value. Parameters: The golden ratios. REAL f1,f2,x0,x1,x2,x3 x0=ax At any given time we will keep track of four points, x0,x1,x2,x3 . 10.2ParabolicInterpolationandBrent’sMethod 395Sample 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).x3=cx if(abs(cx-bx).gt.abs(bx-ax))then Make x0tox1the smaller segment, x1=bx x2=bx+C*(cx-bx) and fill in the new point to be tried. else x2=bx x1=bx-C*(bx-ax) endiff1=f(x1) The initial function evaluations. Note that we never need to evaluate the function at the original endpoints. f2=f(x2) 1 if(abs(x3-x0).gt.tol*(abs(x1)+abs(x2)))then Do-while loop: we keep returning here. if(f2.lt.f1)then One possible outcome, x0=x1 its housekeeping, x1=x2 x2=R*x1+C*x3f1=f2 f2=f(x2) and a new function evaluation. else The other outcome, x3=x2x2=x1 x1=R*x2+C*x0 f2=f1f1=f(x1) and its new function evaluation. endif goto 1 Back to see if we are done. endifif(f1.lt.f2)then We are done. Output the best of the two current values. golden=f1 xmin=x1 else golden=f2 xmin=x2 endifreturn END 10.2 ParabolicInterpolationandBrent’sMethod in One Dimension We already tipped our hand about the desirability of parabolic interpolation in the previous section ’smnbrakroutine, but it is now time to be more explicit. A golden section search is designed to handle, in effect, the worst possible case offunctionminimization,with the uncooperativeminimumhunteddownandcornered like a scared rabbit. But why assume the worst? If the function is nicely parabolic nearto the minimum —surelythe genericcase forsuf ficientlysmoothfunctions — then the parabola fitted through any three points ought to take us in a single leap to the minimum, or at least very near to it (see Figure 10.2.1). Since we want to find an abscissa rather than an ordinate, the procedure is technically called inverse parabolic interpolation . Theformulafortheabscissa xthat is the minimumofa parabolathroughthree points f(a),f(b), and f(c)is x=b−1 2(b−a)2[f(b)−f(c)]−(b−c)2[f(b)−f(a)] (b−a)[f(b)−f(c)]−(b−c)[f(b)−f(a)](10.2.1 )