00001 // This may look like C code, but it is really -*- C++ -*- 00002 /* 00003 ************************************************************************ 00004 * 00005 * Numerical Math Package 00006 * 00007 * Brent's one-dimensional minimizer 00008 * 00009 * finds a local minimum of a single argument function 00010 * over a given range 00011 * 00012 * Input 00013 * double fminbr(ax,bx,f,tol) 00014 * const double ax a and b, a < b, specify the interval 00015 * const double bx the minimum is to be sought in 00016 * UnivariateFunctor& f The function under consideration 00017 * const double tol Acceptable tolerance for the minimum 00018 * location. It is an optional parameter 00019 * with default value DBL_EPSILON 00020 * 00021 * Output 00022 * Fminbr returns an estimate to the location of the minimum 00023 * with accuracy 3*SQRT_EPSILON*abs(x) + tol. 00024 * The procedure can only determine a local minimum, which coincides with 00025 * the global one if and only if the function under investigation is 00026 * unimodular. 00027 * If a function being examined possesses no local minimum within 00028 * the given interval, Fminbr returns either the left or the right end 00029 * point of the interval, wherever the function value is smaller. 00030 * 00031 * Algorithm 00032 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 00033 * computations. M., Mir, 1980, p.202 of the Russian edition 00034 * 00035 * The function makes use of a "gold section" procedure combined with 00036 * a parabolic interpolation. 00037 * At each step the code operates three abscissae - x,v, and w. 00038 * x - the last and the best approximation to the minimum location, 00039 * i.e. f(x) <= f(a) or/and f(x) <= f(b) 00040 * (if the function f has a local minimum in (a,b), then both 00041 * conditions are met after one or two steps). 00042 * v,w are previous approximations to the location of the minimum. 00043 * They may coincide with a, b, or x (although the algorithm tries 00044 * to make all u, v, and w distinct). 00045 * Points x, v, and w are used to construct an interpolating parabola, 00046 * whose minimum is regarded as a new approximation to the minimum 00047 * of the function, provided the parabola's minimum falls within [a,b] 00048 * and reduces the current interval [a,b] to a larger extent than the 00049 * gold section procedure does. 00050 * When f(x) has a positive second derivative at the point of minimum 00051 * (which does not coincide with a or b) the procedure converges 00052 * superlinearly at a rate of about 1.324 00053 * 00054 * $Id: fminbr.cpp,v 1.1 2004/05/21 21:02:52 maxx Exp $ 00055 * 00056 ************************************************************************ 00057 */ 00058 00059 #include "math_num.h" 00060 00061 namespace linalg 00062 { 00063 using namespace linalg; 00064 00065 00066 double fminbr( // An estimate to the min location 00067 const double ax, // Specify the interval the minimum 00068 const double bx, // to be sought in 00069 UnivariateFunctor& f, // Function under investigation 00070 const double tol) // Acceptable tolerance 00071 { 00072 assure( tol > 0, "Tolerance must be positive"); 00073 assure( bx > ax, 00074 "Left end point of the interval should be strictly less than the " 00075 "right one" ); 00076 00077 static const double r = (3-sqrt(5.0))/2; // The golden section ratio 00078 static const double sqrt_eps = sqrt(DBL_EPSILON); 00079 00080 double a = ax, b = bx; // Current interval 00081 double v = a + r*(b-a); // First step - always gold section 00082 double fv = f(v); 00083 double x = v; // the last and the best approximation 00084 double fx = fv; 00085 double w = v; // a previous approx to the min 00086 double fw = fv; 00087 00088 for(;;) // Main iteration loop 00089 { 00090 const double range = b-a; // Interval where the minimum 00091 // is searched in 00092 const double midpoint = (a+b)/2; 00093 const double tol_act = // The effective tolerance 00094 sqrt_eps*fabs(x) + tol/3; 00095 00096 00097 00098 if( 2*fabs(x-midpoint) + range <= 4*tol_act ) 00099 return x; // Acceptable approximation is found 00100 00101 // Compute a new step with the gold 00102 // section 00103 double new_step = r * ( x < midpoint ? b-x : a-x ); 00104 00105 00106 // Decide on the interpolation 00107 if( fabs(x-w) >= tol_act ) // If x and w are distinct 00108 { // interpolatiom may be tried 00109 register double p; // Interpolation step is calcula- 00110 register double q; // ted as p/q; division operation 00111 // is delayed until last moment 00112 register double t; 00113 00114 t = (x-w) * (fx-fv); 00115 q = (x-v) * (fx-fw); 00116 p = (x-v)*q - (x-w)*t; 00117 q = 2*(q-t); 00118 00119 if( q > 0 ) // Formulas above computed new_step 00120 p = -p; // = p/q with a wrong sign (on purpose). 00121 else // Correct this, but in such a way so 00122 q = -q; // that q would be positive 00123 00124 if( fabs(p) < fabs(new_step*q) && // If x+p/q falls in [a,b] and is not 00125 p > q*(a-x+2*tol_act) && // too close to a and b, and isn't 00126 p < q*(b-x-2*tol_act) ) // too large, it is accepted 00127 new_step = p/q; 00128 // If p/q is too large then the 00129 // gold section procedure would 00130 // reduce [a,b] to larger extent 00131 } 00132 00133 if( fabs(new_step) < tol_act ) // Adjust the step to be not less 00134 new_step = new_step > 0 ? // than tolerance 00135 tol_act : -tol_act; 00136 00137 // Obtain the next approximation to min 00138 // and reduce the encompassing interval 00139 register const double t = x + new_step; // Tentative point for the min 00140 register const double ft = f(t); 00141 if( ft <= fx ) 00142 { // t is a better approximation 00143 ( t < x ? b : a ) = x; // Reduce the interval so that 00144 // t would fall within it 00145 v = w; w = x; x = t; // Assign the best approx to x 00146 fv=fw; fw=fx; fx=ft; 00147 } 00148 else // x remains the better approx 00149 { 00150 ( t < x ? a : b ) = t; // Reduce the interval encompassing x 00151 00152 if( ft <= fw || w==x ) 00153 { 00154 v = w; w = t; 00155 fv=fw; fw=ft; 00156 } 00157 else if( ft<=fv || v==x || v==w ) 00158 v = t, fv = ft; 00159 } 00160 } // ===== End of loop ===== 00161 00162 } 00163 }