diff --git a/src/main/java/org/apache/commons/math3/optimization/univariate/BrentOptimizer.java b/src/main/java/org/apache/commons/math3/optimization/univariate/BrentOptimizer.java index cff5bfd33..25f2f5049 100644 --- a/src/main/java/org/apache/commons/math3/optimization/univariate/BrentOptimizer.java +++ b/src/main/java/org/apache/commons/math3/optimization/univariate/BrentOptimizer.java @@ -24,13 +24,19 @@ import org.apache.commons.math3.optimization.ConvergenceChecker; import org.apache.commons.math3.optimization.GoalType; /** - * Implements Richard Brent's algorithm (from his book "Algorithms for + * For a function defined on some interval {@code (lo, hi)}, this class + * finds an approximation {@code x} to the point at which the function + * attains its minimum. + * It implements Richard Brent's algorithm (from his book "Algorithms for * Minimization without Derivatives", p. 79) for finding minima of real - * univariate functions. This implementation is an adaptation partly - * based on the Python code from SciPy (module "optimize.py" v0.5). - * If the function is defined on some interval {@code (lo, hi)}, then - * this method finds an approximation {@code x} to the point at which - * the function attains its minimum. + * univariate functions. + *
+ * This code is an adaptation, partly based on the Python code from SciPy + * (module "optimize.py" v0.5); the original algorithm is also modified + * * * @version $Id$ * @since 2.0 @@ -141,6 +147,8 @@ public class BrentOptimizer extends BaseAbstractUnivariateOptimizer { UnivariatePointValuePair previous = null; UnivariatePointValuePair current = new UnivariatePointValuePair(x, isMinim ? fx : -fx); + // Best point encountered so far (which is the initial guess). + UnivariatePointValuePair best = current; int iter = 0; while (true) { @@ -224,10 +232,15 @@ public class BrentOptimizer extends BaseAbstractUnivariateOptimizer { // User-defined convergence checker. previous = current; current = new UnivariatePointValuePair(u, isMinim ? fu : -fu); + best = best(best, + best(current, + previous, + isMinim), + isMinim); if (checker != null) { if (checker.converged(iter, previous, current)) { - return best(current, previous, isMinim); + return best; } } @@ -264,7 +277,11 @@ public class BrentOptimizer extends BaseAbstractUnivariateOptimizer { } } } else { // Default termination (Brent's criterion). - return best(current, previous, isMinim); + return best(best, + best(current, + previous, + isMinim), + isMinim); } ++iter; } @@ -278,7 +295,8 @@ public class BrentOptimizer extends BaseAbstractUnivariateOptimizer { * @param isMinim {@code true} if the selected point must be the one with * the lowest value. * @return the best point, or {@code null} if {@code a} and {@code b} are - * both {@code null}. + * both {@code null}. When {@code a} and {@code b} have the same function + * value, {@code a} is returned. */ private UnivariatePointValuePair best(UnivariatePointValuePair a, UnivariatePointValuePair b, @@ -291,9 +309,9 @@ public class BrentOptimizer extends BaseAbstractUnivariateOptimizer { } if (isMinim) { - return a.getValue() < b.getValue() ? a : b; + return a.getValue() <= b.getValue() ? a : b; } else { - return a.getValue() > b.getValue() ? a : b; + return a.getValue() >= b.getValue() ? a : b; } } } diff --git a/src/test/java/org/apache/commons/math3/optimization/univariate/BrentOptimizerTest.java b/src/test/java/org/apache/commons/math3/optimization/univariate/BrentOptimizerTest.java index 2151c48d2..55b62cac3 100644 --- a/src/test/java/org/apache/commons/math3/optimization/univariate/BrentOptimizerTest.java +++ b/src/test/java/org/apache/commons/math3/optimization/univariate/BrentOptimizerTest.java @@ -184,6 +184,43 @@ public final class BrentOptimizerTest { Assert.assertEquals(804.9355825, result, 1e-6); } + /** + * Contrived example showing that prior to the resolution of MATH-855 + * (second revision), the algorithm would not return the best point if + * it happened to be the initial guess. + */ + @Test + public void testKeepInitIfBest() { + final double minSin = 3 * Math.PI / 2; + final double offset = 1e-8; + final double delta = 1e-7; + final UnivariateFunction f1 = new Sin(); + final UnivariateFunction f2 = new StepFunction(new double[] { minSin, minSin + offset, minSin + 2 * offset}, + new double[] { 0, -1, 0 }); + final UnivariateFunction f = FunctionUtils.add(f1, f2); + // A slightly less stringent tolerance would make the test pass + // even with the previous implementation. + final double relTol = 1e-8; + final UnivariateOptimizer optimizer = new BrentOptimizer(relTol, 1e-100); + final double init = minSin + 1.5 * offset; + final UnivariatePointValuePair result + = optimizer.optimize(200, f, GoalType.MINIMIZE, + minSin - 6.789 * delta, + minSin + 9.876 * delta, + init); + final int numEval = optimizer.getEvaluations(); + + final double sol = result.getPoint(); + final double expected = init; + +// System.out.println("numEval=" + numEval); +// System.out.println("min=" + init + " f=" + f.value(init)); +// System.out.println("sol=" + sol + " f=" + f.value(sol)); +// System.out.println("exp=" + expected + " f=" + f.value(expected)); + + Assert.assertTrue("Best point not reported", f.value(sol) <= f.value(expected)); + } + /** * Contrived example showing that prior to the resolution of MATH-855, * the algorithm, by always returning the last evaluated point, would @@ -200,7 +237,9 @@ public final class BrentOptimizerTest { final UnivariateFunction f = FunctionUtils.add(f1, f2); final UnivariateOptimizer optimizer = new BrentOptimizer(1e-8, 1e-100); final UnivariatePointValuePair result - = optimizer.optimize(200, f, GoalType.MINIMIZE, minSin - 6.789 * delta, minSin + 9.876 * delta); + = optimizer.optimize(200, f, GoalType.MINIMIZE, + minSin - 6.789 * delta, + minSin + 9.876 * delta); final int numEval = optimizer.getEvaluations(); final double sol = result.getPoint();