10.7 Variable Metric Methods in Multidimensions 429 for (its-1;its<=ITMAX;its++){ Main loop over the iterations. *itersits; Insrch(n,p,fp,g,xi,pnew,fret,stpmax,&check,func); The new function evaluation occurs in Insrch;save the function value in fp for the next line search.It is usually safe to ignore the value of check. fp *fret; for (i=1;i<=n;i++){ xi[i]=pnew[i]-p[i]; Update the line direction p[i]=pnew[i]; and the current point. test=0.0; Test for convergence on Ax for(1=1:i<=n:i++)[ http://www.nr Permission is read able files temp=fabs(xi [i])/FMAX(fabs(p[i]),1.0); if (temp test)test=temp; if (test TOLX){ (including this one) 18881992 FREEALL return; 11-600 for (i=1;i<=n;i++)dg[i]=g[i]; Save the old gradient, (*dfunc)(p,g); and get the new gradient. to any Cambridge test=0.0; Test for convergence on zero gradient from NUMERICAL RECIPES IN den=FMAX(*fret,1.0); for (i=1;i<=n;i++){ to make temp=fabs(g[i])*FMAX(fabs(p[i]),1.0)/den; if (temp test)test=temp; America server computer, (Nor one paper UnN电.t THE if (test gtol){ ART FREEALL return: Programs for (i=1;i<=n;i++)dg[i]-g[i]-dg[i]; Compute difference of gradients, for (i=1;i<=n;i++){ and difference times current matrix. hdg[i]=0.0; strictly prohibited for (j=1;j<=n;j++)hdg[i]+hessin[i][j]*dg[j]; to dir fac-fae=sumdg=sumxi=0.0; Calculate dot products for the denomi- for(i=1;i<=n;i++)[ nators. fac +dg[i]*xi[i]; 188819920 OF SCIENTIFIC COMPUTING(ISBN fae +dg[i]*hdg[i]; sumdg +SQR(dg[i]); sumxi +SQR(xi[i]); 2 v@cam 10-621 if (fac sqrt(EPS*sumdg*sumxi)){ Skip update if fac not sufficiently posi- fac=1.0/fac; tive. fad=1.0/fae; Further reproduction. Numerical Recipes -43108 The vector that makes BFGS different from DFP: for (i=1;i<=n;i++)dg[i]=fac*xi[i]-fad*hdg[i]; 10r (i=1;i<=n;i++){ The BFGS updating formula: (outside for(j=1;j<=n;j++){ 膜 hessin[i][i]fac*xi[i]*xi[i] Software. -fad*hdg[i]*hdg[j]+fae*dg[i]*dg[j]; hessin[j][i]-hessin[i][j]; for (i=1;i<=n;i++){ Now calculate the next direction to go. x1[i]=0.0; for (j=1;j<=n;j++)xi[i]-hessin[i][j]*g[j]; and go back for another iteration nrerror("too many iterations in dfpmin"); FREEALL10.7 Variable Metric Methods in Multidimensions 429 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 machine￾readable 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). for (its=1;its<=ITMAX;its++) { Main loop over the iterations. *iter=its; lnsrch(n,p,fp,g,xi,pnew,fret,stpmax,&check,func); The new function evaluation occurs in lnsrch; save the function value in fp for the next line search. It is usually safe to ignore the value of check. fp = *fret; for (i=1;i<=n;i++) { xi[i]=pnew[i]-p[i]; Update the line direction, p[i]=pnew[i]; and the current point. } test=0.0; Test for convergence on ∆x. for (i=1;i<=n;i++) { temp=fabs(xi[i])/FMAX(fabs(p[i]),1.0); if (temp > test) test=temp; } if (test < TOLX) { FREEALL return; } for (i=1;i<=n;i++) dg[i]=g[i]; Save the old gradient, (*dfunc)(p,g); and get the new gradient. test=0.0; Test for convergence on zero gradient. den=FMAX(*fret,1.0); for (i=1;i<=n;i++) { temp=fabs(g[i])*FMAX(fabs(p[i]),1.0)/den; if (temp > test) test=temp; } if (test < gtol) { FREEALL return; } for (i=1;i<=n;i++) dg[i]=g[i]-dg[i]; Compute difference of gradients, for (i=1;i<=n;i++) { and difference times current matrix. hdg[i]=0.0; for (j=1;j<=n;j++) hdg[i] += hessin[i][j]*dg[j]; } fac=fae=sumdg=sumxi=0.0; Calculate dot products for the denomi￾for (i=1;i<=n;i++) { nators. fac += dg[i]*xi[i]; fae += dg[i]*hdg[i]; sumdg += SQR(dg[i]); sumxi += SQR(xi[i]); } if (fac > sqrt(EPS*sumdg*sumxi)) { Skip update if fac not sufficiently posi￾fac=1.0/fac; tive. fad=1.0/fae; The vector that makes BFGS different from DFP: for (i=1;i<=n;i++) dg[i]=fac*xi[i]-fad*hdg[i]; for (i=1;i<=n;i++) { The BFGS updating formula: for (j=i;j<=n;j++) { hessin[i][j] += fac*xi[i]*xi[j] -fad*hdg[i]*hdg[j]+fae*dg[i]*dg[j]; hessin[j][i]=hessin[i][j]; } } } for (i=1;i<=n;i++) { Now calculate the next direction to go, xi[i]=0.0; for (j=1;j<=n;j++) xi[i] -= hessin[i][j]*g[j]; } } and go back for another iteration. nrerror("too many iterations in dfpmin"); FREEALL }