f10-0

Sample 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).Chapter 10. Minimization or Maximization of Functions 10.0 Introduction In a nutshell: You are given a single function fthat depends on one or more independent variables. You want to find the value of those variables where ftakes on a maximum or a minimum value. You can then calculate what value of fis achievedatthemaximumorminimum. Thetasksofmaximizationandminimization are trivially related to each other, since one person’s function fcould just as well be another’s −f. The computational desiderata are the usual ones: Do it quickly, cheaply, and in small memory. Often the computational effort is dominated bythe cost of evaluating f(and also perhaps its partial derivatives with respect to all variables, if the chosen algorithm requires them). In such cases the desiderata are sometimes replacedby the simplesurrogate: Evaluate fas few times as possible. An extremum (maximum or minimum point) can be either global(truly the highest or lowest function value) or local(the highest or lowest in a finite neighborhoodand not on the boundary of that neighborhood). (See Figure 10.0.1.) Finding a global extremum is, in general, a very difficult problem. Two standard heuristics are widely used: (i) find local extrema starting from widely varyingstarting values of the independent variables (perhaps chosen quasi-randomly, as in §7.7), and then pick the most extreme of these (if they are not all the same); or (ii) perturb a local extremum by taking a finite amplitude step away from it, andthen see if your routine returns you to a better point, or “always” to the same one. Relatively recently, so-called “simulated annealing methods” ( §10.9) have demonstratedimportantsuccesses on a varietyof globalextremizationproblems. Our chapter title could just as well be optimization , which is the usual name for this very large field of numerical research. The importance ascribed to thevarious tasks in this field depends strongly on the particular interests of whom you talk to. Economists, and some engineers, are particularly concerned with constrainedoptimization , where there are a priorilimitations on the allowed values of independent variables. For example, the production of wheat in the U.S. must be a nonnegative number. One particularly well-developed area of constrained optimization is linear programming , where both the function to be optimized and the constraints happen to be linear functions of the independent variables. Section 10.8,whichisotherwisesomewhatdisconnectedfromtherestofthematerialthatwehavechosentoincludeinthischapter,implementstheso-called“simplexalgorithm” for linear programming problems. 387 388 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).G Z YFXE DBA C X1 X2⊗⊗⊗ Figure 10.0.1. Extrema of a function in an interval. Points A,C, and Eare local, but not global maxima. Points Band Fare local, but not global minima. The global maximum occurs at G, which is on the boundary of the interval so that the derivative of the function need not vanish there. Theglobal minimum is at D. At point E, derivatives higher than the first vanish, a situation which can cause difficulty for some algorithms. The points X,Y, and Zare said to “bracket”the minimum F, since Yis less than both Xand Z. Oneothersection, §10.9,alsolies outsideofourmainthrust,butfora different reason: so-called “annealing methods ”are relatively new, so we do not yet know where they will ultimately fit into the scheme of things. However, these methods have solved some problems previously thought to be practically insoluble; theyaddress directly the problem of finding global extrema in the presence of large numbers of undesired local extrema. The other sections in this chapter constitute a selection of the best established algorithms in unconstrained minimization. (For de finiteness, we will henceforth regard the optimization problem as that of minimization.) These sections are connected, with later ones depending on earlier ones. If you are just looking for the one“perfect”algorithm to solve your particular application, you may feel that we are telling you more than you want to know. Unfortunately, there is noperfect optimizationalgorithm. This is a case where we strongly urge you to try more than one method in comparative fashion. Your initial choice of method can be based on the following considerations: You must choose between methods that need only evaluations of the functionto be minimizedand methodsthat also requireevaluationsof the derivative of that function. In the multidimensional case, this derivative is the gradient, a vector quantity. Algorithms using the derivative are somewhat more powerful than those using only the function, but not always enough so as to compensate for the additional calculations ofderivatives. We can easily construct examples favoring one approach or favoring the other. However, if you cancompute derivatives, be prepared to try using them. For one-dimensional minimization (minimize a function of one variable) withoutcalculationofthederivative,brackettheminimumasdescribedin §10.1,andthenuse Brent’smethod asdescribedin §10.2. Ifyourfunction has a discontinuous second (or lower) derivative, then the parabolic 10.0Introduction 389Sample 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).interpolationsofBrent ’s methodareofnoadvantage,andyoumightwish to use the simplest formof goldensectionsearch , as describedin §10.1. Forone-dimensionalminimization withcalculationofthederivative, §10.3 supplies a variant of Brent ’s method which makes limited use of the first derivative information. We shy away from the alternative of using derivative information to construct high-order interpolating polynomials. In our experience the improvement in convergence very near a smooth, analytic minimum does not make up for the tendency of polynomials sometimes to give wildly wrong interpolations at early stages, especially for functions that may have sharp, “exponential ”features. We now turn to the multidimensionalcase, bothwith and without computation offirst derivatives. You must choose between methods that require storage of order N2and thosethat requireonlyoforder N,where Nis thenumberofdimensions. For moderate values of Nand reasonable memory sizes this is not a serious constraint. There will be, however, the occasional application where storage may be critical. We give in §10.4 a sometimes overlooked downhill simplex method due to Nelder and Mead. (This use of the word “simplex”is not to be confused with the simplex method of linear programming.) This methodjust crawls downhill in a straightforward fashion that makes almost no special assumptionsaboutyourfunction. This can be extremelyslow, but it can also, in some cases, be extremely robust. Not to be overlooked is the fact that the code is concise and completely self-contained: a general N-dimensional minimization program in under 100 program lines! This method is most useful when the minimization calculation is only an incidental part of your overall problem. The storage requirement is of order N 2, and derivative calculations are not required. Section 10.5 deals with direction-set methods , of which Powell’s method is the prototype. These arethe methodsofchoice whenyoucannoteasily calculate derivatives, and are not necessarily to be sneered at even if you can. Although derivatives are not needed, the method does require a one-dimensionalminimizationsub-algorithmsuchas Brent ’s method(see above). Storage is of order N2. There are two major families of algorithmsfor multidimensionalminimization withcalculation of first derivatives. Both families require a one-dimensional minimization sub-algorithm, which can itself either use, or not use, the derivative information,asyousee fit(dependingontherelativeeffortofcomputingthefunction andofits gradientvector). We donotthinkthateitherfamilydominatestheotherinall applications; you should think of them as available alternatives: Thefirstfamilygoesunderthename conjugategradientmethods ,astypi- fiedbythe Fletcher-Reevesalgorithm andthecloselyrelatedandprobably superiorPolak-Ribiere algorithm . Conjugate gradient methods require only of order a few times Nstorage, require derivative calculations and 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 typi fied 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, aand b, 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 (or c<b<a ), such that f(b)is less than both f(a)and f(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 aand bor between band c. Suppose, to be speci fic, that we make the latter choice. Then we evaluate f(x).I f f(b)<f (x), then the new bracketing triplet of points is (a, b, x );