Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 93
Текст из файла (страница 93)
In this case tryrestarting from a different initial guess.{void fdjac(int n, float x[], float fvec[], float **df,void (*vecfunc)(int, float [], float []));float fmin(float x[]);void lnsrch(int n, float xold[], float fold, float g[], float p[], float x[],float *f, float stpmax, int *check, float (*func)(float []));void qrdcmp(float **a, int n, float *c, float *d, int *sing);void qrupdt(float **r, float **qt, int n, float u[], float v[]);void rsolv(float **a, int n, float d[], float b[]);int i,its,j,k,restrt,sing,skip;float den,f,fold,stpmax,sum,temp,test,*c,*d,*fvcold;float *g,*p,**qt,**r,*s,*t,*w,*xold;c=vector(1,n);d=vector(1,n);fvcold=vector(1,n);g=vector(1,n);p=vector(1,n);9.7 Globally Convergent Methods for Nonlinear Systems of Equations391qt=matrix(1,n,1,n);r=matrix(1,n,1,n);s=vector(1,n);t=vector(1,n);w=vector(1,n);xold=vector(1,n);fvec=vector(1,n);Define global variables.nn=n;nrfuncv=vecfunc;f=fmin(x);The vector fvec is also computed by thistest=0.0;call.for (i=1;i<=n;i++)Test for initial guess being a root.
Use moreif (fabs(fvec[i]) > test)test=fabs(fvec[i]);stringent test than simif (test < 0.01*TOLF) {ply TOLF.*check=0;FREERETURN}for (sum=0.0,i=1;i<=n;i++) sum += SQR(x[i]);Calculate stpmax for line searches.stpmax=STPMX*FMAX(sqrt(sum),(float)n);restrt=1;Ensure initial Jacobian gets computed.for (its=1;its<=MAXITS;its++) {Start of iteration loop.if (restrt) {fdjac(n,x,fvec,r,vecfunc);Initialize or reinitialize Jacobian in r.qrdcmp(r,n,c,d,&sing);QR decomposition of Jacobian.if (sing) nrerror("singular Jacobian in broydn");for (i=1;i<=n;i++) {Form QT explicitly.for (j=1;j<=n;j++) qt[i][j]=0.0;qt[i][i]=1.0;}for (k=1;k<n;k++) {if (c[k]) {for (j=1;j<=n;j++) {sum=0.0;for (i=k;i<=n;i++)sum += r[i][k]*qt[i][j];sum /= c[k];for (i=k;i<=n;i++)qt[i][j] -= sum*r[i][k];}}}for (i=1;i<=n;i++) {Form R explicitly.r[i][i]=d[i];for (j=1;j<i;j++) r[i][j]=0.0;}} else {Carry out Broyden update.for (i=1;i<=n;i++) s[i]=x[i]-xold[i];s = δx.for (i=1;i<=n;i++) {t = R · s.for (sum=0.0,j=i;j<=n;j++) sum += r[i][j]*s[j];t[i]=sum;}skip=1;for (i=1;i<=n;i++) {w = δF − B · s.for (sum=0.0,j=1;j<=n;j++) sum += qt[j][i]*t[j];w[i]=fvec[i]-fvcold[i]-sum;if (fabs(w[i]) >= EPS*(fabs(fvec[i])+fabs(fvcold[i]))) skip=0;Don’t update with noisy components of w.else w[i]=0.0;}if (!skip) {for (i=1;i<=n;i++) {t = QT · w.for (sum=0.0,j=1;j<=n;j++) sum += qt[i][j]*w[j];t[i]=sum;}392Chapter 9.Root Finding and Nonlinear Sets of Equationsfor (den=0.0,i=1;i<=n;i++) den += SQR(s[i]);for (i=1;i<=n;i++) s[i] /= den;Store s/(s · s) in s.qrupdt(r,qt,n,t,s);Update R and QT .for (i=1;i<=n;i++) {if (r[i][i] == 0.0) nrerror("r singular in broydn");d[i]=r[i][i];Diagonal of R stored in d.}}}for (i=1;i<=n;i++) {Compute ∇f ≈ (Q · R)T · F for the line search.for (sum=0.0,j=1;j<=n;j++) sum += qt[i][j]*fvec[j];g[i]=sum;}for (i=n;i>=1;i--) {for (sum=0.0,j=1;j<=i;j++) sum += r[j][i]*g[j];g[i]=sum;}for (i=1;i<=n;i++) {Store x and F.xold[i]=x[i];fvcold[i]=fvec[i];}fold=f;Store f .for (i=1;i<=n;i++) {Right-hand side for linear equations is −QT · F.for (sum=0.0,j=1;j<=n;j++) sum += qt[i][j]*fvec[j];p[i] = -sum;}rsolv(r,n,d,p);Solve linear equations.lnsrch(n,xold,fold,g,p,x,&f,stpmax,check,fmin);lnsrch returns new x and f .
It also calculates fvec at the new x when it calls fmin.test=0.0;Test for convergence on function values.for (i=1;i<=n;i++)if (fabs(fvec[i]) > test) test=fabs(fvec[i]);if (test < TOLF) {*check=0;FREERETURN}if (*check) {True if line search failed to find a new x.if (restrt) FREERETURNFailure; already tried reinitializing the Jacoelse {bian.test=0.0;Check for gradient of f zero, i.e., spuriousden=FMAX(f,0.5*n);convergence.for (i=1;i<=n;i++) {temp=fabs(g[i])*FMAX(fabs(x[i]),1.0)/den;if (temp > test) test=temp;}if (test < TOLMIN) FREERETURNelse restrt=1;Try reinitializing the Jacobian.}} else {Successful step; will use Broyden update forrestrt=0;next step.test=0.0;Test for convergence on δx.for (i=1;i<=n;i++) {temp=(fabs(x[i]-xold[i]))/FMAX(fabs(x[i]),1.0);if (temp > test) test=temp;}if (test < TOLX) FREERETURN}}nrerror("MAXITS exceeded in broydn");FREERETURN}9.7 Globally Convergent Methods for Nonlinear Systems of Equations393More Advanced ImplementationsOne of the principal ways that the methods described so far can fail is if J (in Newton’smethod) or B in (Broyden’s method) becomes singular or nearly singular, so that δx cannotbe determined.
If you are lucky, this situation will not occur very often in practice. Methodsdeveloped so far to deal with this problem involve monitoring the condition number of J andperturbing J if singularity or near singularity is detected. This is most easily implementedif the QR decomposition is used instead of LU in Newton’s method (see [1] for details).Our personal experience is that, while such an algorithm can solve problems where J isexactly singular and the standard Newton’s method fails, it is occasionally less robust onother problems where LU decomposition succeeds. Clearly implementation details involvingroundoff, underflow, etc., are important here and the last word is yet to be written.Our global strategies both for minimization and zero finding have been based on linesearches.
Other global algorithms, such as the hook step and dogleg step methods, are basedinstead on the model-trust region approach, which is related to the Levenberg-Marquardtalgorithm for nonlinear least-squares (§15.5). While somewhat more complicated than linesearches, these methods have a reputation for robustness even when starting far from thedesired zero or minimum [1].CITED REFERENCES AND FURTHER READING:Dennis, J.E., and Schnabel, R.B. 1983, Numerical Methods for Unconstrained Optimization andNonlinear Equations (Englewood Cliffs, NJ: Prentice-Hall). [1]Broyden, C.G. 1965, Mathematics of Computation, vol.
19, pp. 577–593. [2]Chapter 10. Minimization orMaximization of Functions10.0 IntroductionIn a nutshell: You are given a single function f that depends on one or moreindependent variables. You want to find the value of those variables where f takeson a maximum or a minimum value. You can then calculate what value of f isachieved at the maximum or minimum. The tasks of maximization and minimizationare trivially related to each other, since one person’s function f could just as wellbe another’s −f. The computational desiderata are the usual ones: Do it quickly,cheaply, and in small memory.
Often the computational effort is dominated bythe cost of evaluating f (and also perhaps its partial derivatives with respect to allvariables, if the chosen algorithm requires them). In such cases the desiderata aresometimes replaced by the simple surrogate: Evaluate f as few times as possible.An extremum (maximum or minimum point) can be either global (trulythe highest or lowest function value) or local (the highest or lowest in a finiteneighborhood and not on the boundary of that neighborhood). (See Figure 10.0.1.)Finding a global extremum is, in general, a very difficult problem. Two standardheuristics are widely used: (i) find local extrema starting from widely varyingstarting values of the independent variables (perhaps chosen quasi-randomly, as in§7.7), and then pick the most extreme of these (if they are not all the same); or(ii) perturb a local extremum by taking a finite amplitude step away from it, andthen see if your routine returns you to a better point, or “always” to the sameone.
Relatively recently, so-called “simulated annealing methods” (§10.9) havedemonstrated important successes on a variety of global extremization problems.Our chapter title could just as well be optimization, which is the usual namefor this very large field of numerical research. The importance ascribed to thevarious tasks in this field depends strongly on the particular interests of whomyou talk to. Economists, and some engineers, are particularly concerned withconstrained optimization, where there are a priori limitations on the allowed valuesof independent variables.
For example, the production of wheat in the U.S. mustbe a nonnegative number. One particularly well-developed area of constrainedoptimization is linear programming, where both the function to be optimized andthe constraints happen to be linear functions of the independent variables. Section10.8, which is otherwise somewhat disconnected from the rest of the material that wehave chosen to include in this chapter, implements the so-called “simplex algorithm”for linear programming problems.39439510.0 IntroductionGACB⊗ZE⊗X⊗FYDX1X2Figure 10.0.1.
Extrema of a function in an interval. Points A, C, and E are local, but not globalmaxima. Points B and F are local, but not global minima. The global maximum occurs at G, whichis on the boundary of the interval so that the derivative of the function need not vanish there. Theglobal minimum is at D. At point E, derivatives higher than the first vanish, a situation which cancause difficulty for some algorithms.
The points X, Y , and Z are said to “bracket” the minimum F ,since Y is less than both X and Z.One other section, §10.9, also lies outside of our main thrust, but for a differentreason: so-called “annealing methods” are relatively new, so we do not yet knowwhere they will ultimately fit into the scheme of things. However, these methodshave solved some problems previously thought to be practically insoluble; theyaddress directly the problem of finding global extrema in the presence of largenumbers of undesired local extrema.The other sections in this chapter constitute a selection of the best establishedalgorithms in unconstrained minimization.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
