408 Chapter 10. Minimization or Maximization of Functions du=(*df)(u); Now all the housekeeping,sigh. if (fu 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)asu;else bsu; if (fu <fw l w =x){ MOV3(v,fv,dv,w,fw,dw) MOV3(w,fw,dw,u,fu,du) else if (fu fv Il v =x l v==w){ MOV3(v,fv,dv,u,fu,du) 83 granted for (including this one) 19881992 nrerror("Too many iterations in routine dbrent"); return 0.0; Never get here. 111800.872 NUMERICAL RECIPES (Nort CITED REFERENCES AND FURTHER READING: America server computer, University Press. 令 THE Acton,F.S.1970,Numerica/Methods That Work,1990,corrected edition (Washington:Mathe- to make one paper matical Association of America),pp.55;454-458.[1] ART Brent,R.P.1973,Algorithms for Minimization without Derivatives(Englewood Cliffs,NJ:Prentice- Hall),p.78. Programs OF SCIENTIFIC 10.4 Downhill Simplex Method in to dir Multidimensions st COMPUTING (ISBN 19891892 With this section we begin consideration of multidimensional minimization, that is,finding the minimum of a function of more than one independent variable. 10-621 This section stands apart from those which follow.however:All of the algorithms after this section will make explicit use of a one-dimensional minimization algorithm Fuurggoglrion 43108 as a part of their computational strategy.This section implements an entirely Numerical Recipes self-contained strategy,in which one-dimensional minimization does not figure. The downhill simplex method is due to Nelder and Mead [1].The method (outside requires only function evaluations,not derivatives.It is not very efficient in terms North Software. of the number of function evaluations that it requires.Powell's method (810.5)is almost surely faster in all likely applications.However,the downhill simplex method may frequently be the best method to use if the figure of merit is "get something working quickly"for a problem whose computational burden is small. machine The method has a geometrical naturalness about it which makes it delightful to describe or work through: A simplex is the geometrical figure consisting,in N dimensions,of N+1 points(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
408 Chapter 10. Minimization or Maximization of Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). } } du=(*df)(u); Now all the housekeeping, sigh. if (fu = 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) } } } 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 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 after this section will make explicit use of a one-dimensional minimization algorithm as a part of their computational strategy. This section implements an entirely self-contained strategy, in which one-dimensional minimization does not figure. The downhill 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 almost surely faster in all likely applications. However, the downhill simplex method may frequently be the best method 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: A simplex is the geometrical figure consisting, in N dimensions, of N + 1 points (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
10.4 Downhill Simplex Method in Multidimensions 409 simplex at beginning of step high low http://www.nr.com or call 1-800-872- Permission is read able files reflection (a) (including this one) reflection and expansion 7423 (North America to any server computer,is strictly prohibited. t users to make one paper 1988-1992 by Cambridge University Press.Programs from NUMERICAL RECIPES IN C: (b) THE contraction only),or copy for their (c) multiple contraction (d) v@cambr personaluse.Further reproduction, To order Numerical Recipes books or 1988-1992 by Numerical Recipes ART OF SCIENTIFIC COMPUTING(ISBN 0-521-43108-5) Figure 10.4.1. Possible outcomes for a step in the downhill simplex method.The simplex at the beginning of the step,here a tetrahedron,is shown,top.The simplex at the end of the step can be any one of(a)a reflection away from the high point,(b)a reflection and expansion away from the high point,(c) a contraction along one dimension from the high point,or(d)a contraction along all dimensions towards the low point.An appropriate sequence of such steps will always converge to a minimum of the function. (outside described in 810.8.also makes use ofthe geometrical concept ofa simplex.Otherwise North Software. it is completely unrelated to the algorithm that we are describing in this section.)In general we are only interested in simplexes that are nondegenerate,i.e.,that enclose Ame a finite inner N-dimensional volume.If any point of a nondegenerate simplex is visit website taken as the origin,then the N other points define vector directions that span the machine N-dimensional vector space. In one-dimensional minimization,it was possible to bracket a minimum,so that the success of a subsequent isolation was guaranteed.Alas!There is no analogous procedure in multidimensional space.For multidimensional minimization,the best we can do is give our algorithm a starting guess,that is,an N-vector of independent variables as the first point to try.The algorithm is then supposed to make its own way
10.4 Downhill Simplex Method in Multidimensions 409 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). simplex at beginning of step reflection reflection and expansion contraction multiple contraction (a) (b) (c) (d) high low Figure 10.4.1. Possible outcomes for a step in the downhill simplex method. The simplex at the beginning of the step, here a tetrahedron, is shown, top. The simplex at the end of the step can be any one of (a) a reflection away from the high point, (b) a reflection and expansion away from the high point, (c) a contraction along one dimension from the high point, or (d) a contraction along all dimensions towards the 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. Otherwise it is completely unrelated to the algorithm that we are describing in this section.) In general we are only interested in simplexes that are nondegenerate, i.e., that enclose a finite inner N-dimensional volume. If any point of a nondegenerate simplex is taken as the origin, then the N other points define vector directions that span the N-dimensional vector space. In one-dimensional minimization, it was possible to bracket a minimum, so that the success of a subsequent isolation was guaranteed. Alas! There is no analogous procedure in multidimensional space. For multidimensional minimization, the best we can do is give our algorithm a starting guess, that is, an N-vector of independent variables as the first point to try. The algorithm is then supposed to make its own way
410 Chapter 10.Minimization or Maximization of Functions 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 these points (it matters not which)as being your initial starting point Po,then you can take the other N points to be Pi Po +Aei (10.4.1) where the ei's are N unit vectors,and where A is a constant which is your guess of the problem's characteristic length scale.(Or,you could have different Ai's for each vector direction. The downhill simplex method now takes a series of steps,most steps just moving 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, and they are constructed to conserve the volume of the simplex(hence maintain its nondegeneracy).When it can do so,the method expands the simplex in one or 9 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 the valley.If there is a situation where the simplex is trying to"pass through the eye of 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 of behavior,the basic 五是号台 moves are summarized in Figure 10.4.1. Termination criteria can be delicate in any multidimensional minimization routine.Without bracketing,and with more than one independent variable,we 6 no longer have the option of requiring a certain tolerance for a single independent variable.We typically can identify one "cycle"or "step"of our multidimensional algorithm.It is then possible to terminate when the vector distance moved in that step is fractionally smaller in magnitude than some tolerance tol.Alternatively, we could require that the decrease in the function value in the terminating step be fractionally smaller than some tolerance ftol.Note that while tol should not Numerica 10621 usually be smaller than the square root of the machine precision,it is perfectly appropriate to let ftol be of order the machine precision(or perhaps slightly larger E 4310 so as not to be diddled by roundoff). ecipes Note well that either of the above criteria might be fooled by a single anomalous 腿 step that,for one reason or another,failed to get anywhere.Therefore,it is frequently a good idea to restart a multidimensional minimization routine at a point where North it claims to have found a minimum.For this restart,you should reinitialize any ancillary input quantities.In the downhill simplex method,for example,you should reinitialize N of the N+1 vertices of the simplex again by equation(10.4.1),with Po being one of the vertices of the claimed minimum. Restarts should never be very expensive;your algorithm did,after all,converge to the restart point once,and now you are starting the algorithm already there. Consider,then,our N-dimensional amoeba:
410 Chapter 10. Minimization or Maximization of Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (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 these points (it matters not which) as being your initial starting point P0, then you can take the other N points to be Pi = P0 + λei (10.4.1) where the ei’s are N unit 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.) The downhill simplex method now takes a series of steps, most steps just moving 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, and they are constructed to conserve the volume of the simplex (hence maintain its nondegeneracy). When it can do so, the method expands the simplex in 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 the valley. If there is a situation where the simplex is trying to “pass through the eye of 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 of behavior; the basic moves are summarized in Figure 10.4.1. Termination criteria can be delicate in any multidimensional minimization routine. Without bracketing, and with more than one independent variable, we no longer have the option of requiring a certain tolerance for a single independent variable. We typically can identify one “cycle” or “step” of our multidimensional algorithm. It is then possible to terminate when the vector distance moved in that step is fractionally smaller in magnitude than some tolerance tol. Alternatively, we could require that the decrease in the function value in the terminating step be fractionally smaller than some tolerance ftol. Note that while tol should not usually be smaller than the square root of the machine precision, it is perfectly appropriate to let ftol be of order the machine precision (or perhaps slightly larger so as not to be diddled by roundoff). Note well that either of the above criteria might be fooled by a single anomalous step that, for one reason or another, failed to get anywhere. Therefore, it is frequently a good idea to restart a multidimensional minimization routine at a point where it claims to have found a minimum. For this restart, you should reinitialize any ancillary input quantities. In the downhill simplex method, for example, you should reinitialize N of the N + 1 vertices of the simplex again by equation (10.4.1), with P0 being one of the vertices of the claimed minimum. Restarts should never be very expensive; your algorithm did, after all, converge to the restart point once, and now you are starting the algorithm already there. Consider, then, our N-dimensional amoeba:
10.4 Downhill Simplex Method in Multidimensions 411 #include #include "nrutil.h" #define TINY 1.0e-10 A small number. #define NMAX 5000 Maximum allowed number of function evalua- #define GET_PSUM\ tions. for (j=1;jy[2]?(inhi=2,1):(inhi=1,2); SCIENTIFIC for (i=1;iNMAX)nrerror("NMAX exceeded"); (outside *nfunk +2; Begin a new iteration.First extrapolate by a factor-1 through the face of the simplex North Software. across from the high point,i.e.,reflect the simplex from the high point. ytry=amotry(p,y,psum,ndim,funk,ihi,-1.0); if (ytry y[inhi]){ The reflected point is worse than the second-highest,so look for an intermediate lower point,i.e.,do a one-dimensional contraction. ysave=y[ihi]; ytry=amotry(p,y,psum,ndim,funk,ihi,0.5); if (ytry >ysave){ Can't seem to get rid of that high point.Better for (i=1;i<mmpts;i++){ contract around the lowest(best)point
10.4 Downhill Simplex Method in Multidimensions 411 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). #include #include "nrutil.h" #define TINY 1.0e-10 A small number. #define NMAX 5000 Maximum allowed number of function evalua- #define GET_PSUM \ tions. for (j=1;jy[2] ? (inhi=2,1) : (inhi=1,2); for (i=1;i y[ihi]) { inhi=ihi; ihi=i; } else if (y[i] > y[inhi] && i != ihi) inhi=i; } rtol=2.0*fabs(y[ihi]-y[ilo])/(fabs(y[ihi])+fabs(y[ilo])+TINY); Compute the fractional range from highest to lowest and return if satisfactory. if (rtol = NMAX) nrerror("NMAX exceeded"); *nfunk += 2; Begin a new iteration. First extrapolate by a factor −1 through the face of the simplex across from the high point, i.e., reflect the simplex from the high point. ytry=amotry(p,y,psum,ndim,funk,ihi,-1.0); if (ytry = y[inhi]) { The reflected point is worse than the second-highest, so look for an intermediate lower point, i.e., do a one-dimensional contraction. ysave=y[ihi]; ytry=amotry(p,y,psum,ndim,funk,ihi,0.5); if (ytry >= ysave) { Can’t seem to get rid of that high point. Better for (i=1;i<=mpts;i++) { contract around the lowest (best) point
412 Chapter 10. Minimization or Maximization of Functions if(11=11o)[ for (j=1;j<=ndim;j++) p[i][j]=psum[j]=0.5*(p[i][j]+p[i1o][j]); y[i]=(*funk)(psum); *nfunk +ndim; Keep track of function evaluations GET_PSUM Recompute psum else --(*nfunk); Correct the evaluation count. Go back for the test of doneness and the next free_vector(psum,1,ndim); iteration granted for #include "nrutil.h" -200 float amotry(float *p,float y[],float psum[],int ndim, float (*funk)(float [])int ihi,float fac) Extrapolates by a factor fac through the face of the simplex across from the high point,tries from NUMERICAL RECIPES I 19881992 it,and replaces the high point if the new point is better. int ji (North University Press. 令 float faci,fac2,ytry,*ptryi ptry=vector(1,ndim); America computer, fac1=(1.0-fac)/ndim; one paper THE ART fac2=faci-fac; for (j=1;j<=ndim;j++)ptry[j]=psum[j]*fac1-p[ihi][j]*fac2; ytry=(*funk)(ptry); Evaluate the function at the trial point. Programs if (ytry y[ihi]){ If it's better than the highest,then replace the highest. ictly proh y[ihi]=ytry; for (j=1;j<=ndim;j++){ psum[j]ptry[j]-p[ihi][j]; p[ihi][i]=ptry[j]; to dir free_vector(ptry,1,ndim); rectcustsen OF SCIENTIFIC COMPUTING(ISBN return ytry; @cambridge.org 1988-1992 by Numerical Recipes 10-:6211 CITED REFERENCES AND FURTHER READING: 43108-5 Nelder,J.A.,and Mead,R.1965,Computer Journal,vol.7,pp.308-313.[1] Yarbro.L.A.,and Deming,S.N.1974,Ana/ytica Chimica Acta,vol.73,pp.391-398. (outside Jacoby,S.L.S,Kowalik,J.S.,and Pizzo,J.T.1972,Iterative Methods for Nonlinear Optimization Software. Problems (Englewood Cliffs,NJ:Prentice-Hall) ying of 10.5 Direction Set (Powell's)Methods in Multidimensions We know (810.1-810.3)how to minimize a function of one variable.If we start at a point P in N-dimensional space,and proceed from there in some vector
412 Chapter 10. Minimization or Maximization of Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). if (i != ilo) { for (j=1;j<=ndim;j++) p[i][j]=psum[j]=0.5*(p[i][j]+p[ilo][j]); y[i]=(*funk)(psum); } } *nfunk += ndim; Keep track of function evaluations. GET_PSUM Recompute psum. } } else --(*nfunk); Correct the evaluation count. } Go back for the test of doneness and the next free_vector(psum,1,ndim); iteration. } #include "nrutil.h" float amotry(float **p, float y[], float psum[], int ndim, float (*funk)(float []), int ihi, float fac) Extrapolates by a factor fac through the face of the simplex across from the high point, tries it, and replaces the high point if the new point is better. { int j; float fac1,fac2,ytry,*ptry; ptry=vector(1,ndim); fac1=(1.0-fac)/ndim; fac2=fac1-fac; for (j=1;j<=ndim;j++) ptry[j]=psum[j]*fac1-p[ihi][j]*fac2; ytry=(*funk)(ptry); Evaluate the function at the trial point. if (ytry < y[ihi]) { If it’s better than the highest, then replace the highest. y[ihi]=ytry; for (j=1;j<=ndim;j++) { psum[j] += ptry[j]-p[ihi][j]; p[ihi][j]=ptry[j]; } } free_vector(ptry,1,ndim); return ytry; } CITED REFERENCES AND FURTHER READING: Nelder, J.A., and Mead, R. 1965, Computer Journal, vol. 7, pp. 308–313. [1] Yarbro, L.A., and Deming, S.N. 1974, Analytica Chimica Acta, vol. 73, pp. 391–398. Jacoby, S.L.S, Kowalik, J.S., and Pizzo, J.T. 1972, Iterative Methods for Nonlinear Optimization Problems (Englewood Cliffs, NJ: Prentice-Hall). 10.5 Direction Set (Powell’s) Methods in Multidimensions We know (§10.1–§10.3) how to minimize a function of one variable. If we start at a point P in N-dimensional space, and proceed from there in some vector