View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.apache.commons.math3.optimization.univariate;
18  
19  import org.apache.commons.math3.util.Precision;
20  import org.apache.commons.math3.util.FastMath;
21  import org.apache.commons.math3.exception.NumberIsTooSmallException;
22  import org.apache.commons.math3.exception.NotStrictlyPositiveException;
23  import org.apache.commons.math3.optimization.ConvergenceChecker;
24  import org.apache.commons.math3.optimization.GoalType;
25  
26  /**
27   * For a function defined on some interval {@code (lo, hi)}, this class
28   * finds an approximation {@code x} to the point at which the function
29   * attains its minimum.
30   * It implements Richard Brent's algorithm (from his book "Algorithms for
31   * Minimization without Derivatives", p. 79) for finding minima of real
32   * univariate functions.
33   * <br/>
34   * This code is an adaptation, partly based on the Python code from SciPy
35   * (module "optimize.py" v0.5); the original algorithm is also modified
36   * <ul>
37   *  <li>to use an initial guess provided by the user,</li>
38   *  <li>to ensure that the best point encountered is the one returned.</li>
39   * </ul>
40   *
41   * @version $Id: BrentOptimizer.java 1591835 2014-05-02 09:04:01Z tn $
42   * @deprecated As of 3.1 (to be removed in 4.0).
43   * @since 2.0
44   */
45  @Deprecated
46  public class BrentOptimizer extends BaseAbstractUnivariateOptimizer {
47      /**
48       * Golden section.
49       */
50      private static final double GOLDEN_SECTION = 0.5 * (3 - FastMath.sqrt(5));
51      /**
52       * Minimum relative tolerance.
53       */
54      private static final double MIN_RELATIVE_TOLERANCE = 2 * FastMath.ulp(1d);
55      /**
56       * Relative threshold.
57       */
58      private final double relativeThreshold;
59      /**
60       * Absolute threshold.
61       */
62      private final double absoluteThreshold;
63  
64      /**
65       * The arguments are used implement the original stopping criterion
66       * of Brent's algorithm.
67       * {@code abs} and {@code rel} define a tolerance
68       * {@code tol = rel |x| + abs}. {@code rel} should be no smaller than
69       * <em>2 macheps</em> and preferably not much less than <em>sqrt(macheps)</em>,
70       * where <em>macheps</em> is the relative machine precision. {@code abs} must
71       * be positive.
72       *
73       * @param rel Relative threshold.
74       * @param abs Absolute threshold.
75       * @param checker Additional, user-defined, convergence checking
76       * procedure.
77       * @throws NotStrictlyPositiveException if {@code abs <= 0}.
78       * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
79       */
80      public BrentOptimizer(double rel,
81                            double abs,
82                            ConvergenceChecker<UnivariatePointValuePair> checker) {
83          super(checker);
84  
85          if (rel < MIN_RELATIVE_TOLERANCE) {
86              throw new NumberIsTooSmallException(rel, MIN_RELATIVE_TOLERANCE, true);
87          }
88          if (abs <= 0) {
89              throw new NotStrictlyPositiveException(abs);
90          }
91  
92          relativeThreshold = rel;
93          absoluteThreshold = abs;
94      }
95  
96      /**
97       * The arguments are used for implementing the original stopping criterion
98       * of Brent's algorithm.
99       * {@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 && checker.converged(iter, previous, current)) {
244                     return best;
245                 }
246 
247                 // Update a, b, v, w and x.
248                 if (fu <= fx) {
249                     if (u < x) {
250                         b = x;
251                     } else {
252                         a = x;
253                     }
254                     v = w;
255                     fv = fw;
256                     w = x;
257                     fw = fx;
258                     x = u;
259                     fx = fu;
260                 } else {
261                     if (u < x) {
262                         a = u;
263                     } else {
264                         b = u;
265                     }
266                     if (fu <= fw ||
267                         Precision.equals(w, x)) {
268                         v = w;
269                         fv = fw;
270                         w = u;
271                         fw = fu;
272                     } else if (fu <= fv ||
273                                Precision.equals(v, x) ||
274                                Precision.equals(v, w)) {
275                         v = u;
276                         fv = fu;
277                     }
278                 }
279             } else { // Default termination (Brent's criterion).
280                 return best(best,
281                             best(previous,
282                                  current,
283                                  isMinim),
284                             isMinim);
285             }
286             ++iter;
287         }
288     }
289 
290     /**
291      * Selects the best of two points.
292      *
293      * @param a Point and value.
294      * @param b Point and value.
295      * @param isMinim {@code true} if the selected point must be the one with
296      * the lowest value.
297      * @return the best point, or {@code null} if {@code a} and {@code b} are
298      * both {@code null}. When {@code a} and {@code b} have the same function
299      * value, {@code a} is returned.
300      */
301     private UnivariatePointValuePair best(UnivariatePointValuePair a,
302                                           UnivariatePointValuePair b,
303                                           boolean isMinim) {
304         if (a == null) {
305             return b;
306         }
307         if (b == null) {
308             return a;
309         }
310 
311         if (isMinim) {
312             return a.getValue() <= b.getValue() ? a : b;
313         } else {
314             return a.getValue() >= b.getValue() ? a : b;
315         }
316     }
317 }