Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 95
Текст из файла (страница 95)
If we want to minimize the worst case possibility, thenwe will choose z to make these equal, namelyz = 1 − 2w(10.1.5)We see at once that the new point is the symmetric point to b in the original interval,namely with |b − a| equal to |x − c|. This implies that the point x lies in the largerof the two segments (z is positive only if w < 1/2).But where in the larger segment? Where did the value of w itself come from?Presumably from the previous stage of applying our same strategy. Therefore, if zis chosen to be optimal, then so was w before it. This scale similarity implies thatx should be the same fraction of the way from b to c (if that is the bigger segment)as was b from a to c, in other words,z=w1−wEquations (10.1.5) and (10.1.6) give the quadratic equation√3− 52≈ 0.38197yieldingw=w − 3w + 1 = 02(10.1.6)(10.1.7)In other words, the optimal bracketing interval (a, b, c) has its middle point b afractional distance 0.38197 from one end (say, a), and 0.61803 from the other end(say, b).
These fractions are those of the so-called golden mean or golden section,whose supposedly aesthetic properties hark back to the ancient Pythagoreans. Thisoptimal method of function minimization, the analog of the bisection method forfinding zeros, is thus called the golden section search, summarized as follows:Given, at each stage, a bracketing triplet of points, the next point to be triedis that which is a fraction 0.38197 into the larger of the two intervals (measuringfrom the central point of the triplet). If you start out with a bracketing triplet whosesegments are not in the golden ratios, the procedure of choosing successive pointsat the golden mean point of the larger segment will quickly converge you to theproper, self-replicating ratios.The golden section search guarantees that each new function evaluation will(after self-replicating ratios have been achieved) bracket the minimum to an interval400Chapter 10.Minimization or Maximization of Functionsjust 0.61803 times the size of the preceding interval.
This is comparable to, but notquite as good as, the 0.50000 that holds when finding roots by bisection. Note thatthe convergence is linear (in the language of Chapter 9), meaning that successivesignificant figures are won linearly with additional function evaluations.
In thenext section we will give a superlinear method, where the rate at which successivesignificant figures are liberated increases with each successive function evaluation.Routine for Initially Bracketing a MinimumThe preceding discussion has assumed that you are able to bracket the minimumin the first place. We consider this initial bracketing to be an essential part of anyone-dimensional minimization. There are some one-dimensional algorithms thatdo not require a rigorous initial bracketing. However, we would never trade thesecure feeling of knowing that a minimum is “in there somewhere” for the dubiousreduction of function evaluations that these nonbracketing routines may promise.Please bracket your minima (or, for that matter, your zeros) before isolating them!There is not much theory as to how to do this bracketing.
Obviously you wantto step downhill. But how far? We like to take larger and larger steps, starting withsome (wild?) initial guess and then increasing the stepsize at each step either bya constant factor, or else by the result of a parabolic extrapolation of the precedingpoints that is designed to take us to the extrapolated turning point. It doesn’t muchmatter if the steps get big. After all, we are stepping downhill, so we already havethe left and middle points of the bracketing triplet.
We just need to take a big enoughstep to stop the downhill trend and get a high third point.Our standard routine is this:#include <math.h>#include "nrutil.h"#define GOLD 1.618034#define GLIMIT 100.0#define TINY 1.0e-20#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);Here GOLD is the default ratio by which successive intervals are magnified; GLIMIT is themaximum magnification allowed for a parabolic-fit step.void mnbrak(float *ax, float *bx, float *cx, float *fa, float *fb, float *fc,float (*func)(float))Given a function func, and given distinct initial points ax and bx, this routine searches inthe downhill direction (defined by the function as evaluated at the initial points) and returnsnew points ax, bx, cx that bracket a minimum of the function.
Also returned are the functionvalues at the three points, fa, fb, and fc.{float ulim,u,r,q,fu,dum;*fa=(*func)(*ax);*fb=(*func)(*bx);if (*fb > *fa) {Switch roles of a and b so that we can goSHFT(dum,*ax,*bx,dum)downhill in the direction from a to b.SHFT(dum,*fb,*fa,dum)}*cx=(*bx)+GOLD*(*bx-*ax);First guess for c.*fc=(*func)(*cx);while (*fb > *fc) {Keep returning here until we bracket.r=(*bx-*ax)*(*fb-*fc);Compute u by parabolic extrapolation fromq=(*bx-*cx)*(*fb-*fa);a, b, c.
TINY is used to prevent any posu=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/sible division by zero.10.1 Golden Section Search in One Dimension401(2.0*SIGN(FMAX(fabs(q-r),TINY),q-r));ulim=(*bx)+GLIMIT*(*cx-*bx);We won’t go farther than this. Test various possibilities:if ((*bx-u)*(u-*cx) > 0.0) {Parabolic u is between b and c: try it.fu=(*func)(u);if (fu < *fc) {Got a minimum between b and c.*ax=(*bx);*bx=u;*fa=(*fb);*fb=fu;return;} else if (fu > *fb) {Got a minimum between between a and u.*cx=u;*fc=fu;return;}u=(*cx)+GOLD*(*cx-*bx);Parabolic fit was no use. Use default magfu=(*func)(u);nification.} else if ((*cx-u)*(u-ulim) > 0.0) {Parabolic fit is between c and itsfu=(*func)(u);allowed limit.if (fu < *fc) {SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx))SHFT(*fb,*fc,fu,(*func)(u))}} else if ((u-ulim)*(ulim-*cx) >= 0.0) {Limit parabolic u to maximumu=ulim;allowed value.fu=(*func)(u);} else {Reject parabolic u, use default magnificau=(*cx)+GOLD*(*cx-*bx);tion.fu=(*func)(u);}SHFT(*ax,*bx,*cx,u)Eliminate oldest point and continue.SHFT(*fa,*fb,*fc,fu)}}(Because of the housekeeping involved in moving around three or four points andtheir function values, the above program ends up looking deceptively formidable.That is true of several other programs in this chapter as well.
The underlying ideas,however, are quite simple.)Routine for Golden Section Search#include <math.h>#define R 0.61803399The golden ratios.#define C (1.0-R)#define SHFT2(a,b,c) (a)=(b);(b)=(c);#define SHFT3(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);float golden(float ax, float bx, float cx, float (*f)(float), float tol,float *xmin)Given a function f, and given a bracketing triplet of abscissas ax, bx, cx (such that bx isbetween ax and cx, and f(bx) is less than both f(ax) and f(cx)), this routine performs agolden section search for the minimum, isolating it to a fractional precision of about tol.
Theabscissa of the minimum is returned as xmin, and the minimum function value is returned asgolden, the returned function value.{float f1,f2,x0,x1,x2,x3;402Chapter 10.Minimization or Maximization of Functionsx0=ax;At any given time we will keep track of fourx3=cx;points, x0,x1,x2,x3.if (fabs(cx-bx) > fabs(bx-ax)) {Make x0 to x1 the smaller segment,x1=bx;x2=bx+C*(cx-bx);and fill in the new point to be tried.} else {x2=bx;x1=bx-C*(bx-ax);}f1=(*f)(x1);The initial function evaluations. Note thatf2=(*f)(x2);we never need to evaluate the functionwhile (fabs(x3-x0) > tol*(fabs(x1)+fabs(x2))) {at the original endpoints.if (f2 < f1) {One possible outcome,SHFT3(x0,x1,x2,R*x1+C*x3)its housekeeping,SHFT2(f1,f2,(*f)(x2))and a new function evaluation.} else {The other outcome,SHFT3(x3,x2,x1,R*x2+C*x0)SHFT2(f2,f1,(*f)(x1))and its new function evaluation.}}Back to see if we are done.if (f1 < f2) {We are done.
Output the best of the two*xmin=x1;current values.return f1;} else {*xmin=x2;return f2;}}10.2 Parabolic Interpolation and Brent’sMethod in One DimensionWe already tipped our hand about the desirability of parabolic interpolation inthe previous section’s mnbrak routine, but it is now time to be more explicit. Agolden section search is designed to handle, in effect, the worst possible case offunction minimization, with the uncooperative minimum hunted down and corneredlike a scared rabbit.
But why assume the worst? If the function is nicely parabolicnear to the minimum — surely the generic case for sufficiently smooth functions —then the parabola fitted through any three points ought to take us in a single leapto the minimum, or at least very near to it (see Figure 10.2.1). Since we want tofind an abscissa rather than an ordinate, the procedure is technically called inverseparabolic interpolation.The formula for the abscissa x that is the minimum of a parabola through threepoints f(a), f(b), and f(c) isx=b−1 (b − a)2 [f(b) − f(c)] − (b − c)2 [f(b) − f(a)]2 (b − a)[f(b) − f(c)] − (b − c)[f(b) − f(a)](10.2.1)as you can easily derive.
This formula fails only if the three points are collinear,in which case the denominator is zero (minimum of the parabola is infinitely far10.2 Parabolic Interpolation and Brent’s Method403parabola through 1 2 3parabola through 1 2 431254Figure 10.2.1. Convergence to a minimum by inverse parabolic interpolation. A parabola (dashed line) isdrawn through the three original points 1,2,3 on the given function (solid line).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
