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.math.optimization.univariate;
018
019 import org.apache.commons.math.FunctionEvaluationException;
020 import org.apache.commons.math.MaxIterationsExceededException;
021 import org.apache.commons.math.analysis.UnivariateRealFunction;
022 import org.apache.commons.math.optimization.GoalType;
023
024 /**
025 * Implements Richard Brent's algorithm (from his book "Algorithms for
026 * Minimization without Derivatives", p. 79) for finding minima of real
027 * univariate functions.
028 *
029 * @version $Revision: 811685 $ $Date: 2009-09-05 13:36:48 -0400 (Sat, 05 Sep 2009) $
030 * @since 2.0
031 */
032 public class BrentOptimizer extends AbstractUnivariateRealOptimizer {
033
034 /**
035 * Golden section.
036 */
037 private static final double GOLDEN_SECTION = 0.5 * (3 - Math.sqrt(5));
038
039 /**
040 * Construct a solver.
041 */
042 public BrentOptimizer() {
043 super(100, 1E-10);
044 }
045
046 /** {@inheritDoc} */
047 public double optimize(final UnivariateRealFunction f, final GoalType goalType,
048 final double min, final double max, final double startValue)
049 throws MaxIterationsExceededException, FunctionEvaluationException {
050 return optimize(f, goalType, min, max);
051 }
052
053 /** {@inheritDoc} */
054 public double optimize(final UnivariateRealFunction f, final GoalType goalType,
055 final double min, final double max)
056 throws MaxIterationsExceededException, FunctionEvaluationException {
057 clearResult();
058 return localMin(f, goalType, min, max, relativeAccuracy, absoluteAccuracy);
059 }
060
061 /**
062 * Find the minimum of the function {@code f} within the interval {@code (a, b)}.
063 *
064 * If the function {@code f} is defined on the interval {@code (a, b)}, then
065 * this method finds an approximation {@code x} to the point at which {@code f}
066 * attains its minimum.<br/>
067 * {@code t} and {@code eps} define a tolerance {@code tol = eps |x| + t} and
068 * {@code f} is never evaluated at two points closer together than {@code tol}.
069 * {@code eps} should be no smaller than <em>2 macheps</em> and preferable not
070 * much less than <em>sqrt(macheps)</em>, where <em>macheps</em> is the relative
071 * machine precision. {@code t} should be positive.
072 * @param f the function to solve
073 * @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE}
074 * or {@link GoalType#MINIMIZE}
075 * @param a Lower bound of the interval
076 * @param b Higher bound of the interval
077 * @param eps Relative accuracy
078 * @param t Absolute accuracy
079 * @return the point at which the function is minimal.
080 * @throws MaxIterationsExceededException if the maximum iteration count
081 * is exceeded.
082 * @throws FunctionEvaluationException if an error occurs evaluating
083 * the function.
084 */
085 private double localMin(final UnivariateRealFunction f, final GoalType goalType,
086 double a, double b, final double eps, final double t)
087 throws MaxIterationsExceededException, FunctionEvaluationException {
088 double x = a + GOLDEN_SECTION * (b - a);
089 double v = x;
090 double w = x;
091 double e = 0;
092 double fx = computeObjectiveValue(f, x);
093 if (goalType == GoalType.MAXIMIZE) {
094 fx = -fx;
095 }
096 double fv = fx;
097 double fw = fx;
098
099 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 = GOLDEN_SECTION * 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 }