c10-5 (Numerical Recipes in C)

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);}}} 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 vectorSample 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).}*nfunk += ndim;GET_PSUM10.5 Direction Set (Powell’s) Methods in Multidimensions413linmin: 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.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).direction 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 as414Chapter 10.Minimization or Maximization of FunctionsystartFigure 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) +X ∂f1 X ∂2fxi +xi xj + · · ·∂xi2 i,j ∂xi ∂xji(10.5.1)1x·A·x≈ c − b·x +2wherec ≡ f(P)b ≡ −∇f|P∂ 2 f [A]ij ≡∂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.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).x41510.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)δ(∇f) = A · (δx)(10.5.4)Suppose that we have moved along some direction u to a minimum and nowpropose to move along some new direction v.

The condition that motion along v notspoil our minimization along u is just that the gradient stay perpendicular to u, i.e.,that the change in the gradient be perpendicular to u. By equation (10.5.4) this is just0 = u · δ(∇f) = u · A · v(10.5.5)When (10.5.5) holds for two vectors u and v, they are said to be conjugate.When the relation holds pairwise for all members of a set of vectors, they are saidto be a conjugate set.

If you do successive line minimization of a function alonga conjugate set of directions, then you don’t need to redo any of those directions(unless, of course, you spoil things by minimizing along a direction that they arenot conjugate to).A triumph for a direction set method is to come up with a set of N linearlyindependent, mutually conjugate directions. Then, one pass of N line minimizationswill put it exactly at the minimum of a quadratic form like (10.5.1). For functionsf that are not exactly quadratic forms, it won’t be exactly at the minimum; butrepeated cycles of N line minimizations will in due course converge quadraticallyto the minimum.Powell’s Quadratically Convergent MethodPowell first discovered a direction set method that does produce N mutuallyconjugate directions. Here is how it goes: Initialize the set of directions ui tothe basis vectors,ui = eii = 1, .

. . , N(10.5.6)Now repeat the following sequence of steps (“basic procedure”) until your functionstops decreasing:• Save your starting position as P0 .• For i = 1, . . . , N , move Pi−1 to the minimum along direction ui andcall this point Pi .• For i = 1, .

. . , N − 1, set ui ← ui+1 .• Set uN ← PN − P0 .• Move PN to the minimum along direction uN and call this point P0 .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).(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.

This idea we will return to in §10.7!)How does the gradient ∇f change as we move along some direction? Evidently416Chapter 10.Minimization or Maximization of FunctionsDiscarding the Direction of Largest DecreaseThe fox and the grapes: Now that we are going to give up the property ofquadratic convergence, was it so important after all? That depends on the functionthat you are minimizing. Some applications produce functions with long, twistyvalleys.