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.math.optimization.univariate;
18  
19  import org.apache.commons.math.FunctionEvaluationException;
20  import org.apache.commons.math.MaxIterationsExceededException;
21  import org.apache.commons.math.analysis.UnivariateRealFunction;
22  import org.apache.commons.math.optimization.GoalType;
23  
24  /**
25   * Implements Richard Brent's algorithm (from his book "Algorithms for
26   * Minimization without Derivatives", p. 79) for finding minima of real
27   * univariate functions.
28   *  
29   * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
30   * @since 2.0
31   */
32  public class BrentOptimizer extends AbstractUnivariateRealOptimizer {
33      
34      /**
35       * Golden section.
36       */
37      private static final double c = 0.5 * (3 - Math.sqrt(5));
38  
39      /**
40       * Construct a solver.
41       */
42      public BrentOptimizer() {
43          super(100, 1E-10);
44      }
45  
46      /** {@inheritDoc} */
47      public double optimize(final UnivariateRealFunction f, final GoalType goalType,
48                             final double min, final double max, final double startValue)
49          throws MaxIterationsExceededException, FunctionEvaluationException {
50          return optimize(f, goalType, min, max);
51      }
52      
53      /** {@inheritDoc} */
54      public double optimize(final UnivariateRealFunction f, final GoalType goalType,
55                             final double min, final double max)
56          throws MaxIterationsExceededException, FunctionEvaluationException {
57          clearResult();
58          return localMin(f, goalType, min, max, relativeAccuracy, absoluteAccuracy);
59      }
60      
61      /**
62       * Find the minimum of the function {@code f} within the interval {@code (a, b)}.
63       *
64       * If the function {@code f} is defined on the interval {@code (a, b)}, then
65       * this method finds an approximation {@code x} to the point at which {@code f}
66       * attains its minimum.<br/>
67       * {@code t} and {@code eps} define a tolerance {@code tol = eps |x| + t} and
68       * {@code f} is never evaluated at two points closer together than {@code tol}.
69       * {@code eps} should be no smaller than <em>2 macheps</em> and preferable not
70       * much less than <em>sqrt(macheps)</em>, where <em>macheps</em> is the relative
71       * machine precision. {@code t} should be positive.
72       * @param f the function to solve
73       * @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE}
74       * or {@link GoalType#MINIMIZE}
75       * @param a Lower bound of the interval
76       * @param b Higher bound of the interval
77       * @param eps Relative accuracy
78       * @param t Absolute accuracy
79       * @return the point at which the function is minimal.
80       * @throws MaxIterationsExceededException if the maximum iteration count
81       * is exceeded.
82       * @throws FunctionEvaluationException if an error occurs evaluating
83       * the function. 
84       */
85      private double localMin(final UnivariateRealFunction f, final GoalType goalType,
86                              double a, double b, final double eps, final double t)
87          throws MaxIterationsExceededException, FunctionEvaluationException {
88          double x = a + c * (b - a);
89          double v = x;
90          double w = x;
91          double e = 0;
92          double fx = computeObjectiveValue(f, x);
93          if (goalType == GoalType.MAXIMIZE) {
94              fx = -fx;
95          }
96          double fv = fx;
97          double fw = fx;
98  
99          int count = 0;
100         while (count < maximalIterationCount) {
101             double m = 0.5 * (a + b);
102             double tol = eps * Math.abs(x) + t;
103             double t2 = 2 * tol;
104 
105             // Check stopping criterion.
106             if (Math.abs(x - m) > t2 - 0.5 * (b - a)) {
107                 double p = 0;
108                 double q = 0;
109                 double r = 0;
110                 double d = 0;
111                 double u = 0;
112 
113                 if (Math.abs(e) > tol) { // Fit parabola.
114                     r = (x - w) * (fx - fv);
115                     q = (x - v) * (fx - fw);
116                     p = (x - v) * q - (x - w) * r;
117                     q = 2 * (q - r);
118 
119                     if (q > 0) {
120                         p = -p;
121                     } else {
122                         q = -q;
123                     }
124 
125                     r = e;
126                     e = d;
127                 }
128 
129                 if (Math.abs(p) < Math.abs(0.5 * q * r) &&
130                     (p < q * (a - x)) && (p < q * (b - x))) { // Parabolic interpolation step.
131                     d = p / q;
132                     u = x + d;
133 
134                     // f must not be evaluated too close to a or b.
135                     if (((u - a) < t2) || ((b - u) < t2)) {
136                         d = (x < m) ? tol : -tol;
137                     }
138                 } else { // Golden section step.
139                     e = ((x < m) ? b : a) - x;
140                     d = c * e;
141                 }
142 
143                 // f must not be evaluated too close to a or b.
144                 u = x + ((Math.abs(d) > tol) ? d : ((d > 0) ? tol : -tol));
145                 double fu = computeObjectiveValue(f, u);
146                 if (goalType == GoalType.MAXIMIZE) {
147                     fu = -fu;
148                 }
149 
150                 // Update a, b, v, w and x.
151                 if (fu <= fx) {
152                     if (u < x) {
153                         b = x;
154                     } else {
155                         a = x;
156                     }
157                     v = w;
158                     fv = fw;
159                     w = x;
160                     fw = fx;
161                     x = u;
162                     fx = fu;
163                 } else {
164                     if (u < x) {
165                         a = u;
166                     } else {
167                         b = u;
168                     }
169                     if ((fu <= fw) || (w == x)) {
170                         v = w;
171                         fv = fw;
172                         w = u;
173                         fw = fu;
174                     } else if ((fu <= fv) || (v == x) || (v == w)) {
175                         v = u;
176                         fv = fu;
177                     }
178                 }
179             } else { // termination
180                 setResult(x, (goalType == GoalType.MAXIMIZE) ? -fx : fx, count);
181                 return x;
182             }
183 
184             ++count;
185         }
186 
187         throw new MaxIterationsExceededException(maximalIterationCount);
188 
189     }
190 
191 }