c10-0 (Numerical Recipes in C)

10.0 IntroductionIn a nutshell: You are given a single function f that depends on one or moreindependent variables. You want to find the value of those variables where f takeson a maximum or a minimum value. You can then calculate what value of f isachieved at the maximum or minimum. The tasks of maximization and minimizationare trivially related to each other, since one person’s function f could just as wellbe 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 allvariables, if the chosen algorithm requires them). In such cases the desiderata aresometimes replaced by the simple surrogate: Evaluate f as few times as possible.An extremum (maximum or minimum point) can be either global (trulythe highest or lowest function value) or local (the highest or lowest in a finiteneighborhood and not on the boundary of that neighborhood).

(See Figure 10.0.1.)Finding a global extremum is, in general, a very difficult problem. Two standardheuristics 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 sameone.

Relatively recently, so-called “simulated annealing methods” (§10.9) havedemonstrated important successes on a variety of global extremization problems.Our chapter title could just as well be optimization, which is the usual namefor this very large field of numerical research. The importance ascribed to thevarious tasks in this field depends strongly on the particular interests of whomyou talk to. Economists, and some engineers, are particularly concerned withconstrained optimization, where there are a priori limitations on the allowed valuesof independent variables. For example, the production of wheat in the U.S.

mustbe a nonnegative number. One particularly well-developed area of constrainedoptimization is linear programming, where both the function to be optimized andthe constraints happen to be linear functions of the independent variables. Section10.8, which is otherwise somewhat disconnected from the rest of the material that wehave chosen to include in this chapter, implements the so-called “simplex algorithm”for linear programming problems.394Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).Chapter 10.

Minimization orMaximization of Functions39510.0 IntroductionGAC⊗ZE⊗X⊗FYDX1X2Figure 10.0.1. Extrema of a function in an interval. Points A, C, and E are local, but not globalmaxima. Points B and F are local, but not global minima. The global maximum occurs at G, whichis 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 cancause difficulty for some algorithms. The points X, Y , and Z are said to “bracket” the minimum F ,since Y is less than both X and Z.One other section, §10.9, also lies outside of our main thrust, but for a differentreason: so-called “annealing methods” are relatively new, so we do not yet knowwhere they will ultimately fit into the scheme of things. However, these methodshave solved some problems previously thought to be practically insoluble; theyaddress directly the problem of finding global extrema in the presence of largenumbers of undesired local extrema.The other sections in this chapter constitute a selection of the best establishedalgorithms in unconstrained minimization.

(For definiteness, we will henceforthregard the optimization problem as that of minimization.) These sections areconnected, with later ones depending on earlier ones. If you are just looking forthe one “perfect” algorithm to solve your particular application, you may feel thatwe are telling you more than you want to know. Unfortunately, there is no perfectoptimization algorithm. This is a case where we strongly urge you to try more thanone method in comparative fashion. Your initial choice of method can be basedon the following considerations:• You must choose between methods that need only evaluations of thefunction to be minimized and methods that also require evaluations of thederivative of that function. In the multidimensional case, this derivativeis the gradient, a vector quantity. Algorithms using the derivative aresomewhat more powerful than those using only the function, but notalways enough so as to compensate for the additional calculations ofderivatives.

We can easily construct examples favoring one approach orfavoring the other. However, if you can compute derivatives, be preparedto try using them.• For one-dimensional minimization (minimize a function of one variable)without calculation of the derivative, bracket the minimum as described in§10.1, and then use Brent’s method as described in §10.2. If your functionhas a discontinuous second (or lower) derivative, then the parabolicSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).B396Chapter 10.Minimization or Maximization of FunctionsWe now turn to the multidimensional case, both with and without computationof first derivatives.• You must choose between methods that require storage of order N 2 andthose that require only of order N , where N is the number of dimensions.For moderate values of N and reasonable memory sizes this is not aserious constraint.

There will be, however, the occasional applicationwhere storage may be critical.• We give in §10.4 a sometimes overlooked downhill simplex method dueto Nelder and Mead. (This use of the word “simplex” is not to beconfused with the simplex method of linear programming.) This methodjust crawls downhill in a straightforward fashion that makes almost nospecial assumptions about your function.

This can be extremely slow, butit can also, in some cases, be extremely robust. Not to be overlooked isthe fact that the code is concise and completely self-contained: a generalN -dimensional minimization program in under 100 program lines! Thismethod is most useful when the minimization calculation is only anincidental part of your overall problem. The storage requirement is oforder N 2 , and derivative calculations are not required.• Section 10.5 deals with direction-set methods, of which Powell’s methodis the prototype. These are the methods of choice when you cannot easilycalculate derivatives, and are not necessarily to be sneered at even if youcan. Although derivatives are not needed, the method does require aone-dimensional minimization sub-algorithm such as Brent’s method (seeabove). Storage is of order N 2 .There are two major families of algorithms for multidimensional minimizationwith calculation of first derivatives.

Both families require a one-dimensionalminimization sub-algorithm, which can itself either use, or not use, the derivativeinformation, as you see fit (depending on the relative effort of computing the functionand of its gradient vector). We do not think that either family dominates the other inall applications; you should think of them as available alternatives:• The first family goes under the name conjugate gradient methods, as typified by the Fletcher-Reeves algorithm and the closely related and probablysuperior Polak-Ribiere algorithm.

Conjugate gradient methods requireonly of order a few times N storage, require derivative calculations andSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).interpolations of Brent’s method are of no advantage, and you might wishto use the simplest form of golden section search, as described in §10.1.• For one-dimensional minimization with calculation of the derivative, §10.3supplies a variant of Brent’s method which makes limited use of the firstderivative information.

We shy away from the alternative of usingderivative 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 polynomialssometimes to give wildly wrong interpolations at early stages, especiallyfor functions that may have sharp, “exponential” features.10.1 Golden Section Search in One Dimension397You are now ready to proceed with scaling the peaks (and/or plumbing thedepths) of practical optimization.CITED REFERENCES AND FURTHER READING:Dennis, J.E., and Schnabel, R.B.