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 }