001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 package org.apache.commons.math3.optimization.univariate; 018 019 import org.apache.commons.math3.util.Precision; 020 import org.apache.commons.math3.util.FastMath; 021 import org.apache.commons.math3.exception.NumberIsTooSmallException; 022 import org.apache.commons.math3.exception.NotStrictlyPositiveException; 023 import org.apache.commons.math3.optimization.ConvergenceChecker; 024 import org.apache.commons.math3.optimization.GoalType; 025 026 /** 027 * For a function defined on some interval {@code (lo, hi)}, this class 028 * finds an approximation {@code x} to the point at which the function 029 * attains its minimum. 030 * It implements Richard Brent's algorithm (from his book "Algorithms for 031 * Minimization without Derivatives", p. 79) for finding minima of real 032 * univariate functions. 033 * <br/> 034 * This code is an adaptation, partly based on the Python code from SciPy 035 * (module "optimize.py" v0.5); the original algorithm is also modified 036 * <ul> 037 * <li>to use an initial guess provided by the user,</li> 038 * <li>to ensure that the best point encountered is the one returned.</li> 039 * </ul> 040 * 041 * @version $Id: BrentOptimizer.java 1422230 2012-12-15 12:11:13Z erans $ 042 * @deprecated As of 3.1 (to be removed in 4.0). 043 * @since 2.0 044 */ 045 @Deprecated 046 public class BrentOptimizer extends BaseAbstractUnivariateOptimizer { 047 /** 048 * Golden section. 049 */ 050 private static final double GOLDEN_SECTION = 0.5 * (3 - FastMath.sqrt(5)); 051 /** 052 * Minimum relative tolerance. 053 */ 054 private static final double MIN_RELATIVE_TOLERANCE = 2 * FastMath.ulp(1d); 055 /** 056 * Relative threshold. 057 */ 058 private final double relativeThreshold; 059 /** 060 * Absolute threshold. 061 */ 062 private final double absoluteThreshold; 063 064 /** 065 * The arguments are used implement the original stopping criterion 066 * of Brent's algorithm. 067 * {@code abs} and {@code rel} define a tolerance 068 * {@code tol = rel |x| + abs}. {@code rel} should be no smaller than 069 * <em>2 macheps</em> and preferably not much less than <em>sqrt(macheps)</em>, 070 * where <em>macheps</em> is the relative machine precision. {@code abs} must 071 * be positive. 072 * 073 * @param rel Relative threshold. 074 * @param abs Absolute threshold. 075 * @param checker Additional, user-defined, convergence checking 076 * procedure. 077 * @throws NotStrictlyPositiveException if {@code abs <= 0}. 078 * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}. 079 */ 080 public BrentOptimizer(double rel, 081 double abs, 082 ConvergenceChecker<UnivariatePointValuePair> checker) { 083 super(checker); 084 085 if (rel < MIN_RELATIVE_TOLERANCE) { 086 throw new NumberIsTooSmallException(rel, MIN_RELATIVE_TOLERANCE, true); 087 } 088 if (abs <= 0) { 089 throw new NotStrictlyPositiveException(abs); 090 } 091 092 relativeThreshold = rel; 093 absoluteThreshold = abs; 094 } 095 096 /** 097 * The arguments are used for implementing the original stopping criterion 098 * of Brent's algorithm. 099 * {@code abs} and {@code rel} define a tolerance 100 * {@code tol = rel |x| + abs}. {@code rel} should be no smaller than 101 * <em>2 macheps</em> and preferably not much less than <em>sqrt(macheps)</em>, 102 * where <em>macheps</em> is the relative machine precision. {@code abs} must 103 * be positive. 104 * 105 * @param rel Relative threshold. 106 * @param abs Absolute threshold. 107 * @throws NotStrictlyPositiveException if {@code abs <= 0}. 108 * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}. 109 */ 110 public BrentOptimizer(double rel, 111 double abs) { 112 this(rel, abs, null); 113 } 114 115 /** {@inheritDoc} */ 116 @Override 117 protected UnivariatePointValuePair doOptimize() { 118 final boolean isMinim = getGoalType() == GoalType.MINIMIZE; 119 final double lo = getMin(); 120 final double mid = getStartValue(); 121 final double hi = getMax(); 122 123 // Optional additional convergence criteria. 124 final ConvergenceChecker<UnivariatePointValuePair> checker 125 = getConvergenceChecker(); 126 127 double a; 128 double b; 129 if (lo < hi) { 130 a = lo; 131 b = hi; 132 } else { 133 a = hi; 134 b = lo; 135 } 136 137 double x = mid; 138 double v = x; 139 double w = x; 140 double d = 0; 141 double e = 0; 142 double fx = computeObjectiveValue(x); 143 if (!isMinim) { 144 fx = -fx; 145 } 146 double fv = fx; 147 double fw = fx; 148 149 UnivariatePointValuePair previous = null; 150 UnivariatePointValuePair current 151 = new UnivariatePointValuePair(x, isMinim ? fx : -fx); 152 // Best point encountered so far (which is the initial guess). 153 UnivariatePointValuePair best = current; 154 155 int iter = 0; 156 while (true) { 157 final double m = 0.5 * (a + b); 158 final double tol1 = relativeThreshold * FastMath.abs(x) + absoluteThreshold; 159 final double tol2 = 2 * tol1; 160 161 // Default stopping criterion. 162 final boolean stop = FastMath.abs(x - m) <= tol2 - 0.5 * (b - a); 163 if (!stop) { 164 double p = 0; 165 double q = 0; 166 double r = 0; 167 double u = 0; 168 169 if (FastMath.abs(e) > tol1) { // Fit parabola. 170 r = (x - w) * (fx - fv); 171 q = (x - v) * (fx - fw); 172 p = (x - v) * q - (x - w) * r; 173 q = 2 * (q - r); 174 175 if (q > 0) { 176 p = -p; 177 } else { 178 q = -q; 179 } 180 181 r = e; 182 e = d; 183 184 if (p > q * (a - x) && 185 p < q * (b - x) && 186 FastMath.abs(p) < FastMath.abs(0.5 * q * r)) { 187 // Parabolic interpolation step. 188 d = p / q; 189 u = x + d; 190 191 // f must not be evaluated too close to a or b. 192 if (u - a < tol2 || b - u < tol2) { 193 if (x <= m) { 194 d = tol1; 195 } else { 196 d = -tol1; 197 } 198 } 199 } else { 200 // Golden section step. 201 if (x < m) { 202 e = b - x; 203 } else { 204 e = a - x; 205 } 206 d = GOLDEN_SECTION * e; 207 } 208 } else { 209 // Golden section step. 210 if (x < m) { 211 e = b - x; 212 } else { 213 e = a - x; 214 } 215 d = GOLDEN_SECTION * e; 216 } 217 218 // Update by at least "tol1". 219 if (FastMath.abs(d) < tol1) { 220 if (d >= 0) { 221 u = x + tol1; 222 } else { 223 u = x - tol1; 224 } 225 } else { 226 u = x + d; 227 } 228 229 double fu = computeObjectiveValue(u); 230 if (!isMinim) { 231 fu = -fu; 232 } 233 234 // User-defined convergence checker. 235 previous = current; 236 current = new UnivariatePointValuePair(u, isMinim ? fu : -fu); 237 best = best(best, 238 best(previous, 239 current, 240 isMinim), 241 isMinim); 242 243 if (checker != null) { 244 if (checker.converged(iter, previous, current)) { 245 return best; 246 } 247 } 248 249 // Update a, b, v, w and x. 250 if (fu <= fx) { 251 if (u < x) { 252 b = x; 253 } else { 254 a = x; 255 } 256 v = w; 257 fv = fw; 258 w = x; 259 fw = fx; 260 x = u; 261 fx = fu; 262 } else { 263 if (u < x) { 264 a = u; 265 } else { 266 b = u; 267 } 268 if (fu <= fw || 269 Precision.equals(w, x)) { 270 v = w; 271 fv = fw; 272 w = u; 273 fw = fu; 274 } else if (fu <= fv || 275 Precision.equals(v, x) || 276 Precision.equals(v, w)) { 277 v = u; 278 fv = fu; 279 } 280 } 281 } else { // Default termination (Brent's criterion). 282 return best(best, 283 best(previous, 284 current, 285 isMinim), 286 isMinim); 287 } 288 ++iter; 289 } 290 } 291 292 /** 293 * Selects the best of two points. 294 * 295 * @param a Point and value. 296 * @param b Point and value. 297 * @param isMinim {@code true} if the selected point must be the one with 298 * the lowest value. 299 * @return the best point, or {@code null} if {@code a} and {@code b} are 300 * both {@code null}. When {@code a} and {@code b} have the same function 301 * value, {@code a} is returned. 302 */ 303 private UnivariatePointValuePair best(UnivariatePointValuePair a, 304 UnivariatePointValuePair b, 305 boolean isMinim) { 306 if (a == null) { 307 return b; 308 } 309 if (b == null) { 310 return a; 311 } 312 313 if (isMinim) { 314 return a.getValue() <= b.getValue() ? a : b; 315 } else { 316 return a.getValue() >= b.getValue() ? a : b; 317 } 318 } 319 }