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