Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C, страница 98

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.!). On output, p andy will have been reset to ndim+1 new points all within ftol of a minimum function value, andnfunk gives the number of function evaluations taken.{float amotry(float **p, float y[], float psum[], int ndim,float (*funk)(float []), int ihi, float fac);int i,ihi,ilo,inhi,j,mpts=ndim+1;float rtol,sum,swap,ysave,ytry,*psum;psum=vector(1,ndim);*nfunk=0;GET_PSUMfor (;;) {ilo=1;First we must determine which point is the highest (worst), next-highest, and lowest(best), by looping over the points in the simplex.ihi = y[1]>y[2] ? (inhi=2,1) : (inhi=1,2);for (i=1;i<=mpts;i++) {if (y[i] <= y[ilo]) ilo=i;if (y[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 < ftol) {If returning, put best point and value in slot 1.SWAP(y[1],y[ilo])for (i=1;i<=ndim;i++) SWAP(p[1][i],p[ilo][i])break;}if (*nfunk >= NMAX) nrerror("NMAX exceeded");*nfunk += 2;Begin a new iteration.

First extrapolate by a factor −1 through the face of the simplexacross 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[ilo])Gives a result better than the best point, so try an additional extrapolation by afactor 2.ytry=amotry(p,y,psum,ndim,funk,ihi,2.0);else if (ytry >= y[inhi]) {The reflected point is worse than the second-highest, so look for an intermediatelower 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.

Betterfor (i=1;i<=mpts;i++) {contract around the lowest (best) point.412Chapter 10.Minimization or Maximization of Functionsif (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;GET_PSUM}} else --(*nfunk);}free_vector(psum,1,ndim);Keep track of function evaluations.Recompute psum.Correct the evaluation count.Go back for the test of doneness and the nextiteration.}#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, triesit, 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 OptimizationProblems (Englewood Cliffs, NJ: Prentice-Hall).10.5 Direction Set (Powell’s) Methods inMultidimensionsWe know (§10.1–§10.3) how to minimize a function of one variable. If westart at a point P in N -dimensional space, and proceed from there in some vector10.5 Direction Set (Powell’s) Methods in Multidimensions413direction n, then any function of N variables f(P) can be minimized along the linen by our one-dimensional methods. One can dream up various multidimensionalminimization methods that consist of sequences of such line minimizations.

Differentmethods will differ only by how, at each stage, they choose the next direction n totry. All such methods presume the existence of a “black-box” sub-algorithm, whichwe might call linmin (given as an explicit routine at the end of this section), whosedefinition can be taken for now aslinmin: Given as input the vectors P and n, and thefunction f, find the scalar λ that minimizes f(P + λn).Replace P by P + λn.

Replace n by λn. Done.All the minimization methods in this section and in the two sections followingfall under this general schema of successive line minimizations. (The algorithmin §10.7 does not need very accurate line minimizations. Accordingly, it has itsown approximate line minimization routine, lnsrch.) In this section we considera class of methods whose choice of successive directions does not involve explicitcomputation of the function’s gradient; the next two sections do require such gradientcalculations.

You will note that we need not specify whether linmin uses gradientinformation or not. That choice is up to you, and its optimization depends on yourparticular function. You would be crazy, however, to use gradients in linmin andnot use them in the choice of directions, since in this latter role they can drasticallyreduce the total computational burden.But what if, in your application, calculation of the gradient is out of the question.You might first think of this simple method: Take the unit vectors e1 , e2 , . . . eN as aset of directions. Using linmin, move along the first direction to its minimum, thenfrom there along the second direction to its minimum, and so on, cycling through thewhole set of directions as many times as necessary, until the function stops decreasing.This simple method is actually not too bad for many functions.

Even moreinteresting is why it is bad, i.e. very inefficient, for some other functions. Considera function of two dimensions whose contour map (level lines) happens to define along, narrow valley at some angle to the coordinate basis vectors (see Figure 10.5.1).Then the only way “down the length of the valley” going along the basis vectors ateach stage is by a series of many tiny steps.

More generally, in N dimensions, ifthe function’s second derivatives are much larger in magnitude in some directionsthan in others, then many cycles through all N basis vectors will be required inorder to get anywhere. This condition is not all that unusual; according to Murphy’sLaw, you should count on it.Obviously what we need is a better set of directions than the ei ’s.

All directionset methods consist of prescriptions for updating the set of directions as the methodproceeds, attempting to come up with a set which either (i) includes some verygood directions that will take us far along narrow valleys, or else (more subtly)(ii) includes some number of “non-interfering” directions with the special propertythat minimization along one is not “spoiled” by subsequent minimization alonganother, so that interminable cycling through the set of directions can be avoided.414Chapter 10.Minimization or Maximization of FunctionsystartxFigure 10.5.1.

Successive minimizations along coordinate directions in a long, narrow “valley” (shownas contour lines). Unless the valley is optimally oriented, this method is extremely inefficient, takingmany tiny steps to get to the minimum, crossing and re-crossing the principal axis.Conjugate DirectionsThis concept of “non-interfering” directions, more conventionally called conjugate directions, is worth making mathematically explicit.First, note that if we minimize a function along some direction u, then thegradient of the function must be perpendicular to u at the line minimum; if not, thenthere would still be a nonzero directional derivative along u.Next take some particular point P as the origin of the coordinate system withcoordinates x. Then any function f can be approximated by its Taylor seriesf(x) = f(P) + ∂f1 ∂2fxi +xi xj + · · ·∂xi2 i,j ∂xi ∂xji(10.5.1)1x·A·x≈ c − b·x +2wherec ≡ f(P)b ≡ −∇f|P[A]ij ≡∂ 2 f ∂xi ∂xj P(10.5.2)The matrix A whose components are the second partial derivative matrix of thefunction is called the Hessian matrix of the function at P.41510.5 Direction Set (Powell’s) Methods in MultidimensionsIn the approximation of (10.5.1), the gradient of f is easily calculated as∇f = A · x − b(10.5.3)(This implies that the gradient will vanish — the function will be at an extremum —at a value of x obtained by solving A · x = b.