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

That straightforward task produces long routines densely populatedwith “for(k=1;k<=n;k++)” loops.We do not have space to include such routines in this book. Our linmin, whichworks just fine, is instead a kind of bookkeeping swindle. It constructs an “artificial”function of one variable called f1dim, which is the value of your function, say,func, along the line going through the point p in the direction xi. linmin calls ourfamiliar one-dimensional routines mnbrak (§10.1) and brent (§10.3) and instructsthem to minimize f1dim.

linmin communicates with f1dim “over the head” ofmnbrak and brent, through global (external) variables. That is also how it passesto f1dim a pointer to your user-supplied function.The only thing inefficient about linmin is this: Its use as an interface between amultidimensional minimization strategy and a one-dimensional minimization routineresults in some unnecessary copying of vectors hither and yon. That should not10.5 Direction Set (Powell’s) Methods in Multidimensions419normally be a significant addition to the overall computational burden, but we cannotdisguise its inelegance.#include "nrutil.h"#define TOL 2.0e-4Tolerance passed to brent.int ncom;Global variables communicate with f1dim.float *pcom,*xicom,(*nrfunc)(float []);void linmin(float p[], float xi[], int n, float *fret, float (*func)(float []))Given an n-dimensional point p[1..n] and an n-dimensional direction xi[1..n], moves andresets p to where the function func(p) takes on a minimum along the direction xi from p,and replaces xi by the actual vector displacement that p was moved.

Also returns as fretthe value of func at the returned location p. This is actually all accomplished by calling theroutines mnbrak and brent.{float brent(float ax, float bx, float cx,float (*f)(float), float tol, float *xmin);float f1dim(float x);void mnbrak(float *ax, float *bx, float *cx, float *fa, float *fb,float *fc, float (*func)(float));int j;float xx,xmin,fx,fb,fa,bx,ax;ncom=n;Define the global variables.pcom=vector(1,n);xicom=vector(1,n);nrfunc=func;for (j=1;j<=n;j++) {pcom[j]=p[j];xicom[j]=xi[j];}ax=0.0;Initial guess for brackets.xx=1.0;mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim);*fret=brent(ax,xx,bx,f1dim,TOL,&xmin);for (j=1;j<=n;j++) {Construct the vector results to return.xi[j] *= xmin;p[j] += xi[j];}free_vector(xicom,1,n);free_vector(pcom,1,n);}#include "nrutil.h"extern int ncom;Defined in linmin.extern float *pcom,*xicom,(*nrfunc)(float []);float f1dim(float x)Must accompany linmin.{int j;float f,*xt;xt=vector(1,ncom);for (j=1;j<=ncom;j++) xt[j]=pcom[j]+x*xicom[j];f=(*nrfunc)(xt);free_vector(xt,1,ncom);return f;}420Chapter 10.Minimization or Maximization of FunctionsCITED REFERENCES AND FURTHER READING:Brent, R.P.

1973, Algorithms for Minimization without Derivatives (Englewood Cliffs, NJ: PrenticeHall), Chapter 7. [1]Acton, F.S. 1970, Numerical Methods That Work; 1990, corrected edition (Washington: Mathematical Association of America), pp. 464–467. [2]Jacobs, D.A.H. (ed.) 1977, The State of the Art in Numerical Analysis (London: AcademicPress), pp. 259–262.10.6 Conjugate Gradient Methods inMultidimensionsWe consider now the case where you are able to calculate, at a given N dimensional point P, not just the value of a function f(P) but also the gradient(vector of first partial derivatives) ∇f(P).A rough counting argument will show how advantageous it is to use the gradientinformation: Suppose that the function f is roughly approximated as a quadraticform, as above in equation (10.5.1),f(x) ≈ c − b · x +1x·A·x2(10.6.1)Then the number of unknown parameters in f is equal to the number of freeparameters in A and b, which is 12 N (N + 1), which we see to be of order N 2 .Changing any one of these parameters can move the location of the minimum.Therefore, we should not expect to be able to find the minimum until we havecollected an equivalent information content, of order N 2 numbers.In the direction set methods of §10.5, we collected the necessary information bymaking on the order of N 2 separate line minimizations, each requiring “a few” (butsometimes a big few!) function evaluations.

Now, each evaluation of the gradientwill bring us N new components of information. If we use them wisely, we shouldneed to make only of order N separate line minimizations. That is in fact the casefor the algorithms in this section and the next.A factor of N improvement in computational speed is not necessarily implied.As a rough estimate, we might imagine that the calculation of each component ofthe gradient takes about as long as evaluating the function itself.

In that case therewill be of order N 2 equivalent function evaluations both with and without gradientinformation. Even if the advantage is not of order N , however, it is neverthelessquite substantial: (i) Each calculated component of the gradient will typically savenot just one function evaluation, but a number of them, equivalent to, say, a wholeline minimization.

(ii) There is often a high degree of redundancy in the formulasfor the various components of a function’s gradient; when this is so, especially whenthere is also redundancy with the calculation of the function, then the calculation ofthe gradient may cost significantly less than N function evaluations.10.6 Conjugate Gradient Methods in Multidimensions421(a)(b)Figure 10.6.1. (a) Steepest descent method in a long, narrow “valley.” While more efficient than thestrategy of Figure 10.5.1, steepest descent is nonetheless an inefficient strategy, taking many steps toreach the valley floor. (b) Magnified view of one step: A step starts off in the local gradient direction,perpendicular to the contour lines, and traverses a straight line until a local minimum is reached, wherethe traverse is parallel to the local contour lines.A common beginner’s error is to assume that any reasonable way of incorporating gradient information should be about as good as any other.

This line of thoughtleads to the following not very good algorithm, the steepest descent method:Steepest Descent: Start at a point P0 . As many timesas needed, move from point Pi to the point Pi+1 byminimizing along the line from Pi in the direction ofthe local downhill gradient −∇f(Pi ).The problem with the steepest descent method (which, incidentally, goes backto Cauchy), is similar to the problem that was shown in Figure 10.5.1.

The methodwill perform many small steps in going down a long, narrow valley, even if the valleyis a perfect quadratic form. You might have hoped that, say in two dimensions,your first step would take you to the valley floor, the second step directly downthe long axis; but remember that the new gradient at the minimum point of anyline minimization is perpendicular to the direction just traversed. Therefore, withthe steepest descent method, you must make a right angle turn, which does not, ingeneral, take you to the minimum. (See Figure 10.6.1.)Just as in the discussion that led up to equation (10.5.5), we really want away of proceeding not down the new gradient, but rather in a direction that issomehow constructed to be conjugate to the old gradient, and, insofar as possible,to all previous directions traversed.

Methods that accomplish this construction arecalled conjugate gradient methods.In §2.7 we discussed the conjugate gradient method as a technique for solvinglinear algebraic equations by minimizing a quadratic form. That formalism can alsobe applied to the problem of minimizing a function approximated by the quadratic422Chapter 10.Minimization or Maximization of Functionsform (10.6.1). Recall that, starting with an arbitrary initial vector g0 and lettingh0 = g0 , the conjugate gradient method constructs two sequences of vectors fromthe recurrencegi+1 = gi − λi A · hihi+1 = gi+1 + γi hii = 0, 1, 2, .

. .(10.6.2)The vectors satisfy the orthogonality and conjugacy conditionsgi · gj = 0hi · A · hj = 0gi · h j = 0j<i(10.6.3)The scalars λi and γi are given byλi =gi · gigi · h i=hi · A · hihi · A · higi+1 · gi+1γi =gi · gi(10.6.4)(10.6.5)Equations (10.6.2)–(10.6.5) are simply equations (2.7.32)–(2.7.35) for a symmetricA in a new notation. (A self-contained derivation of these results in the context offunction minimization is given by Polak [1].)Now suppose that we knew the Hessian matrix A in equation (10.6.1). Thenwe could use the construction (10.6.2) to find successively conjugate directions hialong which to line-minimize. After N such, we would efficiently have arrived atthe minimum of the quadratic form.

But we don’t know A.Here is a remarkable theorem to save the day: Suppose we happen to havegi = −∇f(Pi ), for some point Pi , where f is of the form (10.6.1). Suppose that weproceed from Pi along the direction hi to the local minimum of f located at somepoint Pi+1 and then set gi+1 = −∇f(Pi+1 ). Then, this gi+1 is the same vectoras would have been constructed by equation (10.6.2). (And we have constructedit without knowledge of A!)Proof: By equation (10.5.3), gi = −A · Pi + b, andgi+1 = −A · (Pi + λhi ) + b = gi − λA · hi(10.6.6)with λ chosen to take us to the line minimum. But at the line minimum hi · ∇f =−hi · gi+1 = 0.

This latter condition is easily combined with (10.6.6) to solve forλ. The result is exactly the expression (10.6.4). But with this value of λ, (10.6.6)is the same as (10.6.2), q.e.d.We have, then, the basis of an algorithm that requires neither knowledge of theHessian matrix A, nor even the storage necessary to store such a matrix. A sequenceof directions hi is constructed, using only line minimizations, evaluations of thegradient vector, and an auxiliary vector to store the latest in the sequence of g’s.The algorithm described so far is the original Fletcher-Reeves version of theconjugate gradient algorithm.

Later, Polak and Ribiere introduced one tiny, butsometimes significant, change. They proposed using the formγi =(gi+1 − gi ) · gi+1gi · gi(10.6.7)10.6 Conjugate Gradient Methods in Multidimensions423instead of equation (10.6.5). “Wait,” you say, “aren’t they equal by the orthogonalityconditions (10.6.3)?” They are equal for exact quadratic forms. In the real world,however, your function is not exactly a quadratic form. Arriving at the supposedminimum of the quadratic form, you may still need to proceed for another set ofiterations. There is some evidence [2] that the Polak-Ribiere formula accomplishesthe transition to further iterations more gracefully: When it runs out of steam, ittends to reset h to be down the local gradient, which is equivalent to beginning theconjugate-gradient procedure anew.The following routine implements the Polak-Ribiere variant, which we recommend; but changing one program line, as shown, will give you Fletcher-Reeves. Theroutine presumes the existence of a function func(p), where p[1..n] is a vectorof length n, and also presumes the existence of a function dfunc(p,df) that setsthe vector gradient df[1..n] evaluated at the input point p.The routine calls linmin to do the line minimizations.

As already discussed,you may wish to use a modified version of linmin that uses dbrent instead ofbrent, i.e., that uses the gradient in doing the line minimizations. See note below.#include <math.h>#include "nrutil.h"#define ITMAX 200#define EPS 1.0e-10Here ITMAX is the maximum allowed number of iterations, while EPS is a small number torectify the special case of converging to exactly zero function value.#define FREEALL free_vector(xi,1,n);free_vector(h,1,n);free_vector(g,1,n);void frprmn(float p[], int n, float ftol, int *iter, float *fret,float (*func)(float []), void (*dfunc)(float [], float []))Given a starting point p[1..n], Fletcher-Reeves-Polak-Ribiere minimization is performed on afunction func, using its gradient as calculated by a routine dfunc.