c10-6 (Numerical Recipes in C), страница 2

The convergence toleranceon the function value is input as ftol. Returned quantities are p (the location of the minimum),iter (the number of iterations that were performed), and fret (the minimum value of thefunction). The routine linmin is called to perform line minimizations.{void linmin(float p[], float xi[], int n, float *fret,float (*func)(float []));int j,its;float gg,gam,fp,dgg;float *g,*h,*xi;g=vector(1,n);h=vector(1,n);xi=vector(1,n);fp=(*func)(p);Initializations.(*dfunc)(p,xi);for (j=1;j<=n;j++) {g[j] = -xi[j];xi[j]=h[j]=g[j];}for (its=1;its<=ITMAX;its++) {Loop over iterations.*iter=its;linmin(p,xi,n,fret,func);Next statement is the normal return:if (2.0*fabs(*fret-fp) <= ftol*(fabs(*fret)+fabs(fp)+EPS)) {FREEALLreturn;}fp= *fret;(*dfunc)(p,xi);dgg=gg=0.0;for (j=1;j<=n;j++) {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).instead 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.424/*Chapter 10.Minimization or Maximization of Functionsgg += g[j]*g[j];dgg += xi[j]*xi[j];*/dgg += (xi[j]+g[j])*xi[j];Unlikely. If gradient is exactly zero thenwe are already done.}nrerror("Too many iterations in frprmn");}Note on Line Minimization Using DerivativesKindly reread the last part of §10.5. We here want to do the same thing, butusing derivative information in performing the line minimization.The modified version of linmin, called dlinmin, and its required companionroutine df1dim follow:#include "nrutil.h"#define TOL 2.0e-4Tolerance passed to dbrent.int ncom;Global variables communicate with df1dim.float *pcom,*xicom,(*nrfunc)(float []);void (*nrdfun)(float [], float []);void dlinmin(float p[], float xi[], int n, float *fret, float (*func)(float []),void (*dfunc)(float [], 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 dbrent.{float dbrent(float ax, float bx, float cx,float (*f)(float), float (*df)(float), float tol, float *xmin);float f1dim(float x);float df1dim(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;nrdfun=dfunc;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);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).}if (gg == 0.0) {FREEALLreturn;}gam=dgg/gg;for (j=1;j<=n;j++) {g[j] = -xi[j];xi[j]=h[j]=g[j]+gam*h[j];}This statement for Fletcher-Reeves.This statement for Polak-Ribiere.10.7 Variable Metric Methods in Multidimensions425*fret=dbrent(ax,xx,bx,f1dim,df1dim,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 dlinmin.extern float *pcom,*xicom,(*nrfunc)(float []);extern void (*nrdfun)(float [], float []);float df1dim(float x){int j;float df1=0.0;float *xt,*df;xt=vector(1,ncom);df=vector(1,ncom);for (j=1;j<=ncom;j++) xt[j]=pcom[j]+x*xicom[j];(*nrdfun)(xt,df);for (j=1;j<=ncom;j++) df1 += df[j]*xicom[j];free_vector(df,1,ncom);free_vector(xt,1,ncom);return df1;}CITED REFERENCES AND FURTHER READING:Polak, E.

1971, Computational Methods in Optimization (New York: Academic Press), §2.3. [1]Jacobs, D.A.H. (ed.) 1977, The State of the Art in Numerical Analysis (London: AcademicPress), Chapter III.1.7 (by K.W. Brodlie). [2]Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag),§8.7.10.7 Variable Metric Methods inMultidimensionsThe goal of variable metric methods, which are sometimes called quasi-Newtonmethods, is not different from the goal of conjugate gradient methods: to accumulateinformation from successive line minimizations so that N such line minimizationslead to the exact minimum of a quadratic form in N dimensions. In that case, themethod will also be quadratically convergent for more general smooth functions.Both variable metric and conjugate gradient methods require that you are able tocompute your function’s gradient, or first partial derivatives, at arbitrary points.

Thevariable metric approach differs from the conjugate gradient in the way that it storesSample 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).}.