c10-2 (Numerical Recipes in C)

402Chapter 10.Minimization or Maximization of Functions}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 farSample 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).x0=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’s Method403parabola through 1 2 3parabola through 1 2 43254Figure 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). The function is evaluatedat the parabola’s minimum, 4, which replaces point 3. A new parabola (dotted line) is drawn throughpoints 1,4,2. The minimum of this parabola is at 5, which is close to the minimum of the function.away). Note, however, that (10.2.1) is as happy jumping to a parabolic maximumas to a minimum. No minimization scheme that depends solely on (10.2.1) is likelyto succeed in practice.The exacting task is to invent a scheme that relies on a sure-but-slow technique,like golden section search, when the function is not cooperative, but that switchesover to (10.2.1) when the function allows.

The task is nontrivial for severalreasons, including these: (i) The housekeeping needed to avoid unnecessary functionevaluations in switching between the two methods can be complicated. (ii) Carefulattention must be given to the “endgame,” where the function is being evaluatedvery near to the roundoff limit of equation (10.1.2). (iii) The scheme for detecting acooperative versus noncooperative function must be very robust.Brent’s method [1] is up to the task in all particulars.

At any particular stage,it is keeping track of six function points (not necessarily all distinct), a, b, u, v,w and x, defined as follows: the minimum is bracketed between a and b; x is thepoint with the very least function value found so far (or the most recent one incase of a tie); w is the point with the second least function value; v is the previousvalue of w; u is the point at which the function was evaluated most recently. Alsoappearing in the algorithm is the point xm , the midpoint between a and b; however,the function is not evaluated there.You can read the code below to understand the method’s logical organization.Mention of a few general principles here may, however, be helpful: Parabolicinterpolation is attempted, fitting through the points x, v, and w.

To be acceptable,the parabolic step must (i) fall within the bounding interval (a, b), and (ii) imply amovement from the best current value x that is less than half the movement of thestep before last. This second criterion insures that the parabolic steps are actuallyconverging to something, rather than, say, bouncing around in some nonconvergentlimit cycle. In the worst possible case, where the parabolic steps are acceptable butSample 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).1404Chapter 10.Minimization or Maximization of Functions#include <math.h>#include "nrutil.h"#define ITMAX 100#define CGOLD 0.3819660#define ZEPS 1.0e-10Here ITMAX is the maximum allowed number of iterations; CGOLD is the golden ratio; ZEPS isa small number that protects against trying to achieve fractional accuracy for a minimum thathappens to be exactly zero.#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);float brent(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 isolatesthe minimum to a fractional precision of about tol using Brent’s method.

The abscissa ofthe minimum is returned as xmin, and the minimum function value is returned as brent, thereturned function value.{int iter;float a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;float e=0.0;This will be the distance moved onthe step before last.a=(ax < cx ? ax : cx);a and b must be in ascending order,b=(ax > cx ? ax : cx);but input abscissas need not be.x=w=v=bx;Initializations...fw=fv=fx=(*f)(x);for (iter=1;iter<=ITMAX;iter++) {Main program loop.xm=0.5*(a+b);tol2=2.0*(tol1=tol*fabs(x)+ZEPS);if (fabs(x-xm) <= (tol2-0.5*(b-a))) {Test for done here.*xmin=x;return fx;}if (fabs(e) > tol1) {Construct a trial parabolic fit.r=(x-w)*(fx-fv);q=(x-v)*(fx-fw);p=(x-v)*q-(x-w)*r;q=2.0*(q-r);if (q > 0.0) p = -p;q=fabs(q);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).useless, the method will approximately alternate between parabolic steps and goldensections, converging in due course by virtue of the latter. The reason for comparingto the step before last seems essentially heuristic: Experience shows that it is betternot to “punish” the algorithm for a single bad step if it can make it up on the next one.Another principle exemplified in the code is never to evaluate the function lessthan a distance tol from a point already evaluated (or from a known bracketingpoint). The reason is that, as we saw in equation (10.1.2), there is simply noinformation content in doing so: the function will differ from the value alreadyevaluated only by an amount of order the roundoff error.