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.optim.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.optim.ConvergenceChecker;
24  import org.apache.commons.math3.optim.nonlinear.scalar.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 1462503 2013-03-29 15:48:27Z luc $
42   * @since 2.0
43   */
44  public class BrentOptimizer extends UnivariateOptimizer {
45      /**
46       * Golden section.
47       */
48      private static final double GOLDEN_SECTION = 0.5 * (3 - FastMath.sqrt(5));
49      /**
50       * Minimum relative tolerance.
51       */
52      private static final double MIN_RELATIVE_TOLERANCE = 2 * FastMath.ulp(1d);
53      /**
54       * Relative threshold.
55       */
56      private final double relativeThreshold;
57      /**
58       * Absolute threshold.
59       */
60      private final double absoluteThreshold;
61  
62      /**
63       * The arguments are used implement the original stopping criterion
64       * of Brent's algorithm.
65       * {@code abs} and {@code rel} define a tolerance
66       * {@code tol = rel |x| + abs}. {@code rel} should be no smaller than
67       * <em>2 macheps</em> and preferably not much less than <em>sqrt(macheps)</em>,
68       * where <em>macheps</em> is the relative machine precision. {@code abs} must
69       * be positive.
70       *
71       * @param rel Relative threshold.
72       * @param abs Absolute threshold.
73       * @param checker Additional, user-defined, convergence checking
74       * procedure.
75       * @throws NotStrictlyPositiveException if {@code abs <= 0}.
76       * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
77       */
78      public BrentOptimizer(double rel,
79                            double abs,
80                            ConvergenceChecker<UnivariatePointValuePair> checker) {
81          super(checker);
82  
83          if (rel < MIN_RELATIVE_TOLERANCE) {
84              throw new NumberIsTooSmallException(rel, MIN_RELATIVE_TOLERANCE, true);
85          }
86          if (abs <= 0) {
87              throw new NotStrictlyPositiveException(abs);
88          }
89  
90          relativeThreshold = rel;
91          absoluteThreshold = abs;
92      }
93  
94      /**
95       * The arguments are used for implementing the original stopping criterion
96       * of Brent's algorithm.
97       * {@code abs} and {@code rel} define a tolerance
98       * {@code tol = rel |x| + abs}. {@code rel} should be no smaller than
99       * <em>2 macheps</em> and preferably not much less than <em>sqrt(macheps)</em>,
100      * where <em>macheps</em> is the relative machine precision. {@code abs} must
101      * be positive.
102      *
103      * @param rel Relative threshold.
104      * @param abs Absolute threshold.
105      * @throws NotStrictlyPositiveException if {@code abs <= 0}.
106      * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
107      */
108     public BrentOptimizer(double rel,
109                           double abs) {
110         this(rel, abs, null);
111     }
112 
113     /** {@inheritDoc} */
114     @Override
115     protected UnivariatePointValuePair doOptimize() {
116         final boolean isMinim = getGoalType() == GoalType.MINIMIZE;
117         final double lo = getMin();
118         final double mid = getStartValue();
119         final double hi = getMax();
120 
121         // Optional additional convergence criteria.
122         final ConvergenceChecker<UnivariatePointValuePair> checker
123             = getConvergenceChecker();
124 
125         double a;
126         double b;
127         if (lo < hi) {
128             a = lo;
129             b = hi;
130         } else {
131             a = hi;
132             b = lo;
133         }
134 
135         double x = mid;
136         double v = x;
137         double w = x;
138         double d = 0;
139         double e = 0;
140         double fx = computeObjectiveValue(x);
141         if (!isMinim) {
142             fx = -fx;
143         }
144         double fv = fx;
145         double fw = fx;
146 
147         UnivariatePointValuePair previous = null;
148         UnivariatePointValuePair current
149             = new UnivariatePointValuePair(x, isMinim ? fx : -fx);
150         // Best point encountered so far (which is the initial guess).
151         UnivariatePointValuePair best = current;
152 
153         int iter = 0;
154         while (true) {
155             final double m = 0.5 * (a + b);
156             final double tol1 = relativeThreshold * FastMath.abs(x) + absoluteThreshold;
157             final double tol2 = 2 * tol1;
158 
159             // Default stopping criterion.
160             final boolean stop = FastMath.abs(x - m) <= tol2 - 0.5 * (b - a);
161             if (!stop) {
162                 double p = 0;
163                 double q = 0;
164                 double r = 0;
165                 double u = 0;
166 
167                 if (FastMath.abs(e) > tol1) { // Fit parabola.
168                     r = (x - w) * (fx - fv);
169                     q = (x - v) * (fx - fw);
170                     p = (x - v) * q - (x - w) * r;
171                     q = 2 * (q - r);
172 
173                     if (q > 0) {
174                         p = -p;
175                     } else {
176                         q = -q;
177                     }
178 
179                     r = e;
180                     e = d;
181 
182                     if (p > q * (a - x) &&
183                         p < q * (b - x) &&
184                         FastMath.abs(p) < FastMath.abs(0.5 * q * r)) {
185                         // Parabolic interpolation step.
186                         d = p / q;
187                         u = x + d;
188 
189                         // f must not be evaluated too close to a or b.
190                         if (u - a < tol2 || b - u < tol2) {
191                             if (x <= m) {
192                                 d = tol1;
193                             } else {
194                                 d = -tol1;
195                             }
196                         }
197                     } else {
198                         // Golden section step.
199                         if (x < m) {
200                             e = b - x;
201                         } else {
202                             e = a - x;
203                         }
204                         d = GOLDEN_SECTION * e;
205                     }
206                 } else {
207                     // Golden section step.
208                     if (x < m) {
209                         e = b - x;
210                     } else {
211                         e = a - x;
212                     }
213                     d = GOLDEN_SECTION * e;
214                 }
215 
216                 // Update by at least "tol1".
217                 if (FastMath.abs(d) < tol1) {
218                     if (d >= 0) {
219                         u = x + tol1;
220                     } else {
221                         u = x - tol1;
222                     }
223                 } else {
224                     u = x + d;
225                 }
226 
227                 double fu = computeObjectiveValue(u);
228                 if (!isMinim) {
229                     fu = -fu;
230                 }
231 
232                 // User-defined convergence checker.
233                 previous = current;
234                 current = new UnivariatePointValuePair(u, isMinim ? fu : -fu);
235                 best = best(best,
236                             best(previous,
237                                  current,
238                                  isMinim),
239                             isMinim);
240 
241                 if (checker != null && checker.converged(iter, previous, current)) {
242                     return best;
243                 }
244 
245                 // Update a, b, v, w and x.
246                 if (fu <= fx) {
247                     if (u < x) {
248                         b = x;
249                     } else {
250                         a = x;
251                     }
252                     v = w;
253                     fv = fw;
254                     w = x;
255                     fw = fx;
256                     x = u;
257                     fx = fu;
258                 } else {
259                     if (u < x) {
260                         a = u;
261                     } else {
262                         b = u;
263                     }
264                     if (fu <= fw ||
265                         Precision.equals(w, x)) {
266                         v = w;
267                         fv = fw;
268                         w = u;
269                         fw = fu;
270                     } else if (fu <= fv ||
271                                Precision.equals(v, x) ||
272                                Precision.equals(v, w)) {
273                         v = u;
274                         fv = fu;
275                     }
276                 }
277             } else { // Default termination (Brent's criterion).
278                 return best(best,
279                             best(previous,
280                                  current,
281                                  isMinim),
282                             isMinim);
283             }
284             ++iter;
285         }
286     }
287 
288     /**
289      * Selects the best of two points.
290      *
291      * @param a Point and value.
292      * @param b Point and value.
293      * @param isMinim {@code true} if the selected point must be the one with
294      * the lowest value.
295      * @return the best point, or {@code null} if {@code a} and {@code b} are
296      * both {@code null}. When {@code a} and {@code b} have the same function
297      * value, {@code a} is returned.
298      */
299     private UnivariatePointValuePair best(UnivariatePointValuePair a,
300                                           UnivariatePointValuePair b,
301                                           boolean isMinim) {
302         if (a == null) {
303             return b;
304         }
305         if (b == null) {
306             return a;
307         }
308 
309         if (isMinim) {
310             return a.getValue() <= b.getValue() ? a : b;
311         } else {
312             return a.getValue() >= b.getValue() ? a : b;
313         }
314     }
315 }