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 }