Added a new minimization package with an implementation of the Brent algorithm

contributed by Gilles Sadowski.
The implementation needs some testing as it seems to never use the parabola fitting
and only relying on golden section. This may be due to the refactoring I did on
the original patch.
Jira: MATH-177

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@735475 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2009-01-18 15:39:35 +00:00
parent e7c3207b05
commit e60e3de474
9 changed files with 583 additions and 6 deletions

View File

@ -147,6 +147,9 @@
<contributor>
<name>Andreas Rieger</name>
</contributor>
<contributor>
<name>Gilles Sadowski</name>
</contributor>
<contributor>
<name>Joni Salonen</name>
</contributor>

View File

@ -0,0 +1,204 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.analysis.minimization;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.MaxIterationsExceededException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
/**
* Implements Richard Brent's algorithm (from his book "Algorithms for
* Minimization without Derivatives", p. 79) for finding minima of real
* univariate functions.
*
* @version $Revision$ $Date$
* @since 2.0
*/
public class BrentMinimizer extends UnivariateRealMinimizerImpl {
/** Serializable version identifier */
private static final long serialVersionUID = 7185472920191999565L;
/**
* Golden section.
*/
private static final double c = 0.5 * (3 - Math.sqrt(5));
/**
* Construct a solver.
*/
public BrentMinimizer() {
super(100, 1E-10);
}
/**
* Find a minimum in the given interval, start at startValue.
* <p>
* A minimizer may require that the interval brackets a single minimum.
* </p>
* @param f the function to minimize.
* @param min the lower bound for the interval.
* @param max the upper bound for the interval.
* @param startValue this parameter is <em>not</em> used at all
* @return a value where the function is minimum
* @throws ConvergenceException if the maximum iteration count is exceeded
* or the minimizer detects convergence problems otherwise.
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if min > max or the arguments do not
* satisfy the requirements specified by the minimizer
*/
public double minimize(final UnivariateRealFunction f,
final double min, final double max, final double startValue)
throws MaxIterationsExceededException, FunctionEvaluationException {
return minimize(f, min, max);
}
/** {@inheritDoc} */
public double minimize(final UnivariateRealFunction f,
final double min, final double max)
throws MaxIterationsExceededException,
FunctionEvaluationException {
clearResult();
return localmin(min, max, relativeAccuracy, absoluteAccuracy, f);
}
/**
* Find the minimum of the function {@code f} within the interval {@code (a, b)}.
*
* If the function {@code f} is defined on the interval {@code (a, b)}, then
* this method finds an approximation {@code x} to the point at which {@code f}
* attains its minimum.<br/>
* {@code t} and {@code eps} define a tolerance {@code tol = eps |x| + t} and
* {@code f} is never evaluated at two points closer together than {@code tol}.
* {@code eps} should be no smaller than <em>2 macheps</em> and preferable not
* much less than <em>sqrt(macheps)</em>, where <em>macheps</em> is the relative
* machine precision. {@code t} should be positive.
*
* @param f the function to solve
* @param a Lower bound of the interval.
* @param b Higher bound of the interval.
* @param eps Relative accuracy.
* @param t Absolute accuracy.
* @return the point at which the function is minimal.
* @throws MaxIterationsExceededException if the maximum iteration count
* is exceeded.
* @throws FunctionEvaluationException if an error occurs evaluating
* the function.
*/
private double localmin(double a, double b, final double eps,
final double t, final UnivariateRealFunction f)
throws MaxIterationsExceededException, FunctionEvaluationException {
double x = a + c * (b - a);
double v = x;
double w = x;
double e = 0;
double fx = f.value(x);
double fv = fx;
double fw = fx;
int count = 0;
while (count < maximalIterationCount) {
double m = 0.5 * (a + b);
double tol = eps * Math.abs(x) + t;
double t2 = 2 * tol;
// Check stopping criterion.
if (Math.abs(x - m) > t2 - 0.5 * (b - a)) {
double p = 0;
double q = 0;
double r = 0;
double d = 0;
double u = 0;
if (Math.abs(e) > tol) { // Fit parabola.
r = (x - w) * (fx - fv);
q = (x - v) * (fx - fw);
p = (x - v) * q - (x - w) * r;
q = 2 * (q - r);
if (q > 0) {
p = -p;
} else {
q = -q;
}
r = e;
e = d;
}
if (Math.abs(p) < Math.abs(0.5 * q * r) &&
(p < q * (a - x)) && (p < q * (b - x))) { // Parabolic interpolation step.
d = p / q;
u = x + d;
// f must not be evaluated too close to a or b.
if (((u - a) < t2) || ((b - u) < t2)) {
d = (x < m) ? tol : -tol;
}
} else { // Golden section step.
e = ((x < m) ? b : a) - x;
d = c * e;
}
// f must not be evaluated too close to a or b.
u = x + ((Math.abs(d) > tol) ? d : ((d > 0) ? tol : -tol));
double fu = f.value(u);
// Update a, b, v, w and x.
if (fu <= fx) {
if (u < x) {
b = x;
} else {
a = x;
}
v = w;
fv = fw;
w = x;
fw = fx;
x = u;
fx = fu;
} else {
if (u < x) {
a = u;
} else {
b = u;
}
if ((fu <= fw) || (w == x)) {
v = w;
fv = fw;
w = u;
fw = fu;
} else if ((fu <= fv) || (v == x) || (v == w)) {
v = u;
fv = fu;
}
}
} else { // Termination.
setResult(x, fx, count);
return x;
}
++count;
}
throw new MaxIterationsExceededException(maximalIterationCount);
}
}

View File

@ -0,0 +1,118 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.analysis.minimization;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.ConvergingAlgorithm;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
/**
* Interface for (univariate real) minimization algorithms.
*
* @version $Revision$ $Date$
* @since 2.0
*/
public interface UnivariateRealMinimizer extends ConvergingAlgorithm {
// /**
// * Set the function value accuracy.
// * <p>
// * This is used to determine when an evaluated function value or some other
// * value which is used as divisor is zero.</p>
// * <p>
// * This is a safety guard and it shouldn't be necessary to change this in
// * general.</p>
// *
// * @param accuracy the accuracy.
// * @throws IllegalArgumentException if the accuracy can't be achieved by
// * the minimizer or is otherwise deemed unreasonable.
// */
// void setFunctionValueAccuracy(double accuracy);
//
// /**
// * Get the actual function value accuracy.
// * @return the accuracy
// */
// double getFunctionValueAccuracy();
//
// /**
// * Reset the actual function accuracy to the default.
// * The default value is provided by the minimizer implementation.
// */
// void resetFunctionValueAccuracy();
/**
* Find a minimum in the given interval.
* <p>
* A minimizer may require that the interval brackets a single minimum.
* </p>
* @param f the function to minimize.
* @param min the lower bound for the interval.
* @param max the upper bound for the interval.
* @return a value where the function is minimum
* @throws ConvergenceException if the maximum iteration count is exceeded
* or the minimizer detects convergence problems otherwise.
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if min > max or the endpoints do not
* satisfy the requirements specified by the minimizer
*/
double minimize(UnivariateRealFunction f, double min, double max)
throws ConvergenceException,
FunctionEvaluationException;
/**
* Find a minimum in the given interval, start at startValue.
* <p>
* A minimizer may require that the interval brackets a single minimum.
* </p>
* @param f the function to minimize.
* @param min the lower bound for the interval.
* @param max the upper bound for the interval.
* @param startValue the start value to use
* @return a value where the function is minimum
* @throws ConvergenceException if the maximum iteration count is exceeded
* or the minimizer detects convergence problems otherwise.
* @throws FunctionEvaluationException if an error occurs evaluating the
* function
* @throws IllegalArgumentException if min > max or the arguments do not
* satisfy the requirements specified by the minimizer
*/
double minimize(UnivariateRealFunction f, double min, double max, double startValue)
throws ConvergenceException, FunctionEvaluationException;
/**
* Get the result of the last run of the minimizer.
*
* @return the last result.
* @throws IllegalStateException if there is no result available, either
* because no result was yet computed or the last attempt failed.
*/
double getResult();
/**
* Get the result of the last run of the minimizer.
*
* @return the value of the function at the last result.
* @throws IllegalStateException if there is no result available, either
* because no result was yet computed or the last attempt failed.
*/
double getFunctionValue();
}

View File

@ -0,0 +1,134 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.analysis.minimization;
import org.apache.commons.math.ConvergingAlgorithmImpl;
import org.apache.commons.math.MathRuntimeException;
/**
* Provide a default implementation for several functions useful to generic
* minimizers.
*
* @version $Revision$ $Date$
* @since 2.0
*/
public abstract class UnivariateRealMinimizerImpl
extends ConvergingAlgorithmImpl implements UnivariateRealMinimizer {
/** Serializable version identifier. */
private static final long serialVersionUID = 4543031162377070699L;
// /** Maximum error of function. */
// protected double functionValueAccuracy;
//
// /** Default maximum error of function. */
// protected double defaultFunctionValueAccuracy;
/** Indicates where a root has been computed. */
protected boolean resultComputed = false;
/** The last computed root. */
protected double result;
/** Value of the function at the last computed result. */
protected double functionValue;
/**
* Construct a solver with given iteration count and accuracy.
*
* @param defaultAbsoluteAccuracy maximum absolute error
* @param defaultMaximalIterationCount maximum number of iterations
* @throws IllegalArgumentException if f is null or the
* defaultAbsoluteAccuracy is not valid
*/
protected UnivariateRealMinimizerImpl(int defaultMaximalIterationCount,
double defaultAbsoluteAccuracy) {
super(defaultMaximalIterationCount, defaultAbsoluteAccuracy);
// this.functionValueAccuracy = defaultFunctionValueAccuracy;
}
/** Check if a result has been computed.
* @exception IllegalStateException if no result has been computed
*/
protected void checkResultComputed() throws IllegalArgumentException {
if (!resultComputed) {
throw MathRuntimeException.createIllegalStateException("no result available", null);
}
}
/** {@inheritDoc} */
public double getResult() {
checkResultComputed();
return result;
}
/** {@inheritDoc} */
public double getFunctionValue() {
checkResultComputed();
return functionValue;
}
// /** {@inheritDoc} */
// public void setFunctionValueAccuracy(double accuracy) {
// functionValueAccuracy = accuracy;
// }
//
// /** {@inheritDoc} */
// public double getFunctionValueAccuracy() {
// return functionValueAccuracy;
// }
//
// /** {@inheritDoc} */
// public void resetFunctionValueAccuracy() {
// functionValueAccuracy = defaultFunctionValueAccuracy;
// }
/**
* Convenience function for implementations.
*
* @param result the result to set
* @param iterationCount the iteration count to set
*/
protected final void setResult(double result, int iterationCount) {
this.result = result;
this.iterationCount = iterationCount;
this.resultComputed = true;
}
/**
* Convenience function for implementations.
*
* @param x the result to set
* @param fx the result to set
* @param iterationCount the iteration count to set
*/
protected final void setResult(double x, double fx, int iterationCount) {
this.result = x;
this.functionValue = fx;
this.iterationCount = iterationCount;
this.resultComputed = true;
}
/**
* Convenience function for implementations.
*/
protected final void clearResult() {
this.resultComputed = false;
}
}

View File

@ -0,0 +1,22 @@
<html>
<!--
Licensed to the Apache Software Foundation (ASF) under one or more
contributor license agreements. See the NOTICE file distributed with
this work for additional information regarding copyright ownership.
The ASF licenses this file to You under the Apache License, Version 2.0
(the "License"); you may not use this file except in compliance with
the License. You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
-->
<!-- $Revision$ $Date$ -->
<body>
Univariate real functions minimum finding algorithms.
</body>
</html>

View File

@ -39,6 +39,9 @@ The <action> type attribute can be add,update,fix,remove.
</properties>
<body>
<release version="2.0" date="TBD" description="TBD">
<action dev="luc" type="add" issue="MATH-177" due-to="Gilles Sadowski">
Added a new minimization sub-package below the analysis package.
</action>
<action dev="luc" type="update" >
The analysis package has been reorganized with several sub-packages.
</action>

View File

@ -237,7 +237,19 @@ double c = solver.solve(function, 1.0, 5.0);</source>
</table>
</p>
</subsection>
<subsection name="4.3 Interpolation" href="interpolation">
<subsection name="4.3 Minimization" href="minimization">
<p>
A <a href="../apidocs/org/apache/commons/math/analysis/minimization/UnivariateRealMinimizer.html">
org.apache.commons.math.analysis.minimization.UnivariateRealMinimizer</a>
is used to find the minimal values of a univariate real-valued function <code>f</code>.
</p>
<p>
Minimization algorithms usage is very similar to root-finding algorithms usage explained
above. The main difference is that the <code>solve</code> methods in root finding algorithms
is replaced by <code>minimize</code> methods.
</p>
</subsection>
<subsection name="4.4 Interpolation" href="interpolation">
<p>
A <a href="../apidocs/org/apache/commons/math/analysis/interpolation/UnivariateRealInterpolator.html">
org.apache.commons.math.analysis.interpolation.UnivariateRealInterpolator</a>
@ -288,7 +300,7 @@ System.out println("f(" + interpolationX + ") = " + interpolatedY);</source>
adding more points does not always lead to a better interpolation.
</p>
</subsection>
<subsection name="4.4 Integration" href="integration">
<subsection name="4.5 Integration" href="integration">
<p>
A <a href="../apidocs/org/apache/commons/math/analysis/integration/UnivariateRealIntegrator.html">
org.apache.commons.math.analysis.integration.UnivariateRealIntegrator.</a>
@ -303,7 +315,7 @@ System.out println("f(" + interpolationX + ") = " + interpolatedY);</source>
</ul>
</p>
</subsection>
<subsection name="4.5 Polynomials" href="polynomials">
<subsection name="4.6 Polynomials" href="polynomials">
<p>
The <a href="../apidocs/org/apache/commons/math/analysis/polynomials/package.html">
org.apache.commons.math.analysis.polynomials</a> package provides real coefficients

View File

@ -66,9 +66,10 @@
<ul>
<li><a href="analysis.html#4.1 Overview">4.1 Overview</a></li>
<li><a href="analysis.html#4.2 Root-finding">4.2 Root-finding</a></li>
<li><a href="analysis.html#4.3 Interpolation">4.3 Interpolation</a></li>
<li><a href="analysis.html#4.4 Integration">4.4 Integration</a></li>
<li><a href="analysis.html#4.5 Polynomials">4.5 Polynomials</a></li>
<li><a href="analysis.html#4.3 Minimization">4.3 Minimization</a></li>
<li><a href="analysis.html#4.4 Interpolation">4.4 Interpolation</a></li>
<li><a href="analysis.html#4.5 Integration">4.5 Integration</a></li>
<li><a href="analysis.html#4.6 Polynomials">4.6 Polynomials</a></li>
</ul></li>
<li><a href="special.html">5. Special Functions</a>
<ul>

View File

@ -0,0 +1,80 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.analysis.minimization;
import junit.framework.Test;
import junit.framework.TestCase;
import junit.framework.TestSuite;
import org.apache.commons.math.MathException;
import org.apache.commons.math.analysis.QuinticFunction;
import org.apache.commons.math.analysis.SinFunction;
import org.apache.commons.math.analysis.UnivariateRealFunction;
/**
* @version $Revision$ $Date$
*/
public final class BrentMinimizerTest extends TestCase {
public BrentMinimizerTest(String name) {
super(name);
}
public static Test suite() {
TestSuite suite = new TestSuite(BrentMinimizerTest.class);
suite.setName("BrentMinimizer Tests");
return suite;
}
public void testSinMin() throws MathException {
UnivariateRealFunction f = new SinFunction();
UnivariateRealMinimizer minimizer = new BrentMinimizer();
assertEquals(3 * Math.PI / 2, minimizer.minimize(f, 4, 5), 70 * minimizer.getAbsoluteAccuracy());
assertTrue(minimizer.getIterationCount() <= 50);
assertEquals(3 * Math.PI / 2, minimizer.minimize(f, 1, 5), 70 * minimizer.getAbsoluteAccuracy());
assertTrue(minimizer.getIterationCount() <= 50);
}
public void testQuinticMin() throws MathException {
// The quintic function has zeros at 0, +-0.5 and +-1.
// The function has extrema (first derivative is zero) at 0.27195613 and 0.82221643,
UnivariateRealFunction f = new QuinticFunction();
UnivariateRealMinimizer minimizer = new BrentMinimizer();
assertEquals(-0.27195613, minimizer.minimize(f, -0.3, -0.2), 1.0e-8);
assertEquals( 0.82221643, minimizer.minimize(f, 0.3, 0.9), 1.0e-8);
assertTrue(minimizer.getIterationCount() <= 50);
// search in a large interval
assertEquals(-0.27195613, minimizer.minimize(f, -1.0, 0.2), 1.0e-8);
assertTrue(minimizer.getIterationCount() <= 50);
}
public void testMinEndpoints() throws Exception {
UnivariateRealFunction f = new SinFunction();
UnivariateRealMinimizer solver = new BrentMinimizer();
// endpoint is minimum
double result = solver.minimize(f, 3 * Math.PI / 2, 5);
assertEquals(3 * Math.PI / 2, result, 70 * solver.getAbsoluteAccuracy());
result = solver.minimize(f, 4, 3 * Math.PI / 2);
assertEquals(3 * Math.PI / 2, result, 70 * solver.getAbsoluteAccuracy());
}
}