c10-4 (Numerical Recipes in C)

408Chapter 10.Minimization or Maximization of Functions}nrerror("Too many iterations in routine dbrent");return 0.0;Never get here.}CITED REFERENCES AND FURTHER READING:Acton, F.S. 1970, Numerical Methods That Work; 1990, corrected edition (Washington: Mathematical Association of America), pp. 55; 454–458. [1]Brent, R.P. 1973, Algorithms for Minimization without Derivatives (Englewood Cliffs, NJ: PrenticeHall), p. 78.10.4 Downhill Simplex Method inMultidimensionsWith 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 algorithmsafter this section will make explicit use of a one-dimensional minimization algorithmas a part of their computational strategy.

This section implements an entirelyself-contained strategy, in which one-dimensional minimization does not figure.The downhill simplex method is due to Nelder and Mead [1]. The methodrequires only function evaluations, not derivatives. It is not very efficient in termsof the number of function evaluations that it requires. Powell’s method (§10.5) isalmost surely faster in all likely applications. However, the downhill simplex methodmay frequently be the best method to use if the figure of merit is “get somethingworking quickly” for a problem whose computational burden is small.The method has a geometrical naturalness about it which makes it delightfulto describe or work through:A simplex is the geometrical figure consisting, in N dimensions, of N + 1points (or vertices) and all their interconnecting line segments, polygonal faces, etc.In two dimensions, a simplex is a triangle.

In three dimensions it is a tetrahedron,not necessarily the regular tetrahedron. (The simplex method of linear programming,Sample 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).}}du=(*df)(u);Now all the housekeeping, sigh.if (fu <= fx) {if (u >= x) a=x; else b=x;MOV3(v,fv,dv, w,fw,dw)MOV3(w,fw,dw, x,fx,dx)MOV3(x,fx,dx, u,fu,du)} else {if (u < x) a=u; else b=u;if (fu <= fw || w == x) {MOV3(v,fv,dv, w,fw,dw)MOV3(w,fw,dw, u,fu,du)} else if (fu < fv || v == x || v == w) {MOV3(v,fv,dv, u,fu,du)}}10.4 Downhill Simplex Method in Multidimensions409simplex at beginning of stephighlow(a)reflection and expansion(b)contraction(c)multiplecontraction(d)Figure 10.4.1.Possible outcomes for a step in the downhill simplex method.

The simplex at thebeginning of the step, here a tetrahedron, is shown, top. The simplex at the end of the step can be any oneof (a) a reflection away from the high point, (b) a reflection and expansion away from the high point, (c) acontraction along one dimension from the high point, or (d) a contraction along all dimensions towardsthe low point. An appropriate sequence of such steps will always converge to a minimum of the function.described in §10.8, also makes use of the geometrical concept of a simplex. Otherwiseit is completely unrelated to the algorithm that we are describing in this section.) Ingeneral we are only interested in simplexes that are nondegenerate, i.e., that enclosea finite inner N -dimensional volume.

If any point of a nondegenerate simplex istaken as the origin, then the N other points define vector directions that span theN -dimensional vector space.In one-dimensional minimization, it was possible to bracket a minimum, so thatthe success of a subsequent isolation was guaranteed. Alas! There is no analogousprocedure in multidimensional space. For multidimensional minimization, the bestwe can do is give our algorithm a starting guess, that is, an N -vector of independentvariables as the first point to try. The algorithm is then supposed to make its own waySample 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).reflection410Chapter 10.Minimization or Maximization of FunctionsPi = P0 + λei(10.4.1)where the ei ’s are N unit vectors, and where λ is a constant which is your guessof the problem’s characteristic length scale. (Or, you could have different λi ’s foreach vector direction.)The downhill simplex method now takes a series of steps, most steps justmoving 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 calledreflections, and they are constructed to conserve the volume of the simplex (hencemaintain its nondegeneracy). When it can do so, the method expands the simplexin one or another direction to take larger steps. When it reaches a “valley floor,”the method contracts itself in the transverse direction and tries to ooze down thevalley.

If there is a situation where the simplex is trying to “pass through the eyeof a needle,” it contracts itself in all directions, pulling itself in around its lowest(best) point. The routine name amoeba is intended to be descriptive of this kind ofbehavior; the basic moves are summarized in Figure 10.4.1.Termination criteria can be delicate in any multidimensional minimizationroutine. Without bracketing, and with more than one independent variable, weno longer have the option of requiring a certain tolerance for a single independentvariable.

We typically can identify one “cycle” or “step” of our multidimensionalalgorithm. 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 befractionally smaller than some tolerance ftol. Note that while tol should notusually be smaller than the square root of the machine precision, it is perfectlyappropriate to let ftol be of order the machine precision (or perhaps slightly largerso as not to be diddled by roundoff).Note well that either of the above criteria might be fooled by a single anomalousstep that, for one reason or another, failed to get anywhere.

Therefore, it is frequentlya good idea to restart a multidimensional minimization routine at a point whereit claims to have found a minimum. For this restart, you should reinitialize anyancillary input quantities. In the downhill simplex method, for example, you shouldreinitialize N of the N + 1 vertices of the simplex again by equation (10.4.1), withP0 being one of the vertices of the claimed minimum.Restarts should never be very expensive; your algorithm did, after all, convergeto the restart point once, and now you are starting the algorithm already there.Consider, then, our N -dimensional amoeba:Sample 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).downhill 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 + 1 points, defining an initial simplex. If you think of one of thesepoints (it matters not which) as being your initial starting point P0 , then you cantake the other N points to be10.4 Downhill Simplex Method in Multidimensions411void amoeba(float **p, float y[], int ndim, float ftol,float (*funk)(float []), int *nfunk)Multidimensional minimization of the function funk(x) where x[1..ndim] is a vector in ndimdimensions, 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 ofthe starting simplex. Also input is the vector y[1..ndim+1], whose components must be preinitialized to the values of funk evaluated at the ndim+1 vertices (rows) of p; and ftol thefractional convergence tolerance to be achieved in the function value (n.b.!).