diff --git a/src/main/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolver.java b/src/main/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolver.java new file mode 100644 index 000000000..dc5f16f0b --- /dev/null +++ b/src/main/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolver.java @@ -0,0 +1,397 @@ +/* + * 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.solvers; + + +import org.apache.commons.math.analysis.UnivariateRealFunction; +import org.apache.commons.math.exception.MathInternalError; +import org.apache.commons.math.exception.NoBracketingException; +import org.apache.commons.math.exception.NumberIsTooSmallException; +import org.apache.commons.math.util.FastMath; +import org.apache.commons.math.util.MathUtils; + +/** + * This class implements a modification of the Brent algorithm. + *

+ * The changes with respect to the original Brent algorithm are: + *

+ *

+ * The given interval must bracket the root. + * + * @version $Id$ + */ +public class BracketingNthOrderBrentSolver + extends AbstractUnivariateRealSolver + implements BracketedUnivariateRealSolver { + + /** Default absolute accuracy. */ + private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6; + + /** Default maximal order. */ + private static final int DEFAULT_MAXIMAL_ORDER = 5; + + /** Maximal aging triggering an attempt to balance the bracketing interval. */ + private static final int MAXIMAL_AGING = 2; + + /** Reduction factor for attempts to balance the bracketing interval. */ + private static final double REDUCTION_FACTOR = 1.0 / 16.0; + + /** Maximal order. */ + private final int maximalOrder; + + /** The kinds of solutions that the algorithm may accept. */ + private AllowedSolutions allowed; + + /** + * Construct a solver with default accuracy and maximal order (1e-6 and 5 respectively) + */ + public BracketingNthOrderBrentSolver() { + this(DEFAULT_ABSOLUTE_ACCURACY, DEFAULT_MAXIMAL_ORDER); + } + + /** + * Construct a solver. + * + * @param absoluteAccuracy Absolute accuracy. + * @param maximalOrder maximal order. + * @exception NumberIsTooSmallException if maximal order is lower than 2 + */ + public BracketingNthOrderBrentSolver(final double absoluteAccuracy, + final int maximalOrder) + throws NumberIsTooSmallException { + super(absoluteAccuracy); + if (maximalOrder < 2) { + throw new NumberIsTooSmallException(maximalOrder, 2, true); + } + this.maximalOrder = maximalOrder; + this.allowed = AllowedSolutions.ANY_SIDE; + } + + /** + * Construct a solver. + * + * @param relativeAccuracy Relative accuracy. + * @param absoluteAccuracy Absolute accuracy. + * @param maximalOrder maximal order. + * @exception NumberIsTooSmallException if maximal order is lower than 2 + */ + public BracketingNthOrderBrentSolver(final double relativeAccuracy, + final double absoluteAccuracy, + final int maximalOrder) + throws NumberIsTooSmallException { + super(relativeAccuracy, absoluteAccuracy); + if (maximalOrder < 2) { + throw new NumberIsTooSmallException(maximalOrder, 2, true); + } + this.maximalOrder = maximalOrder; + this.allowed = AllowedSolutions.ANY_SIDE; + } + + /** + * Construct a solver. + * + * @param relativeAccuracy Relative accuracy. + * @param absoluteAccuracy Absolute accuracy. + * @param functionValueAccuracy Function value accuracy. + * @param maximalOrder maximal order. + * @exception NumberIsTooSmallException if maximal order is lower than 2 + */ + public BracketingNthOrderBrentSolver(final double relativeAccuracy, + final double absoluteAccuracy, + final double functionValueAccuracy, + final int maximalOrder) + throws NumberIsTooSmallException { + super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy); + if (maximalOrder < 2) { + throw new NumberIsTooSmallException(maximalOrder, 2, true); + } + this.maximalOrder = maximalOrder; + this.allowed = AllowedSolutions.ANY_SIDE; + } + + /** Get the maximal order. + * @return maximal order + */ + public int getMaximalOrder() { + return maximalOrder; + } + + /** + * {@inheritDoc} + */ + @Override + protected double doSolve() { + + // prepare arrays with the first points + final double[] x = new double[maximalOrder + 1]; + final double[] y = new double[maximalOrder + 1]; + x[0] = getMin(); + x[1] = getStartValue(); + x[2] = getMax(); + verifySequence(x[0], x[1], x[2]); + + // evaluate initial guess + y[1] = computeObjectiveValue(x[1]); + if (MathUtils.equals(y[1], 0.0, 1)) { + // return the initial guess if it is a perfect root. + return x[1]; + } + + // evaluate first endpoint + y[0] = computeObjectiveValue(x[0]); + if (MathUtils.equals(y[0], 0.0, 1)) { + // return the first endpoint if it is a perfect root. + return x[0]; + } + + int nbPoints; + int signChangeIndex; + if (y[0] * y[1] < 0) { + + // reduce interval if it brackets the root + nbPoints = 2; + signChangeIndex = 1; + + } else { + + // evaluate second endpoint + y[2] = computeObjectiveValue(x[2]); + if (MathUtils.equals(y[2], 0.0, 1)) { + // return the second endpoint if it is a perfect root. + return x[2]; + } + + if (y[1] * y[2] < 0) { + // use all computed point as a start sampling array for solving + nbPoints = 3; + signChangeIndex = 2; + } else { + throw new NoBracketingException(x[0], x[2], y[0], y[2]); + } + + } + + // prepare a work array for inverse polynomial interpolation + final double[] tmpX = new double[x.length]; + + // current tightest bracketing of the root + double xA = x[signChangeIndex - 1]; + double yA = y[signChangeIndex - 1]; + double absYA = FastMath.abs(yA); + int agingA = 0; + double xB = x[signChangeIndex]; + double yB = y[signChangeIndex]; + double absYB = FastMath.abs(yB); + int agingB = 0; + + // search loop + while (true) { + + // check convergence of bracketing interval + final double xTol = getAbsoluteAccuracy() + + getRelativeAccuracy() * FastMath.max(FastMath.abs(xA), FastMath.abs(xB)); + if (((xB - xA) <= xTol) || (FastMath.max(absYA, absYB) < getFunctionValueAccuracy())) { + switch (allowed) { + case ANY_SIDE : + return absYA < absYB ? xA : xB; + case LEFT_SIDE : + return xA; + case RIGHT_SIDE : + return xB; + case BELOW_SIDE : + return (yA <= 0) ? xA : xB; + case ABOVE_SIDE : + return (yA < 0) ? xB : xA; + default : + // this should never happen + throw new MathInternalError(null); + } + } + + // target for the next evaluation point + double targetY; + if (agingA >= MAXIMAL_AGING) { + // we keep updating the high bracket, try to compensate this + targetY = -REDUCTION_FACTOR * yB; + } else if (agingB >= MAXIMAL_AGING) { + // we keep updating the low bracket, try to compensate this + targetY = -REDUCTION_FACTOR * yA; + } else { + // bracketing is balanced, try to find the root itself + targetY = 0; + } + + // make a few attempts to guess a root, + double nextX; + int start = 0; + int end = nbPoints; + do { + + // guess a value for current target, using inverse polynomial interpolation + System.arraycopy(x, start, tmpX, start, end - start); + nextX = guessX(targetY, tmpX, y, start, end); + + if (!((nextX > xA) && (nextX < xB))) { + // the guessed root is not strictly inside of the tightest bracketing interval + + // the guessed root is either not strictly inside the interval or it + // is a NaN (which occurs when some sampling points share the same y) + // we try again with a lower interpolation order + if (signChangeIndex - start >= end - signChangeIndex) { + // we have more points before the sign change, drop the lowest point + ++start; + } else { + // we have more points after sign change, drop the highest point + --end; + } + + // we need to do one more attempt + nextX = Double.NaN; + + } + + } while (Double.isNaN(nextX) && (end - start > 1)); + + if (Double.isNaN(nextX)) { + // fall back to bisection + nextX = xA + 0.5 * (xB - xA); + start = signChangeIndex - 1; + end = signChangeIndex; + } + + // evaluate the function at the guessed root + final double nextY = computeObjectiveValue(nextX); + if (MathUtils.equals(nextY, 0.0, 1)) { + // we have found an exact root, since it is not an approximation + // we don't need to bother about the allowed solutions setting + return nextX; + } + + if ((nbPoints > 2) && (end - start != nbPoints)) { + + // we have been forced to ignore some points to keep bracketing, + // they are probably too far from the root, drop them from now on + nbPoints = end - start; + System.arraycopy(x, start, x, 0, nbPoints); + System.arraycopy(y, start, y, 0, nbPoints); + signChangeIndex -= start; + + } else if (nbPoints == x.length) { + + // we have to drop one point in order to insert the new one + nbPoints--; + + // keep the tightest bracketing interval as centered as possible + if (signChangeIndex >= (x.length + 1) / 2) { + // we drop the lowest point, we have to shift the arrays and the index + System.arraycopy(x, 1, x, 0, nbPoints); + System.arraycopy(y, 1, y, 0, nbPoints); + --signChangeIndex; + } + + } + + // insert the last computed point + //(by construction, we know it lies inside the tightest bracketing interval) + System.arraycopy(x, signChangeIndex, x, signChangeIndex + 1, nbPoints - signChangeIndex); + x[signChangeIndex] = nextX; + System.arraycopy(y, signChangeIndex, y, signChangeIndex + 1, nbPoints - signChangeIndex); + y[signChangeIndex] = nextY; + ++nbPoints; + + // update the bracketing interval + if (nextY * yA <= 0) { + // the sign change occurs before the inserted point + xB = nextX; + yB = nextY; + absYB = FastMath.abs(yB); + ++agingA; + agingB = 0; + } else { + // the sign change occurs after the inserted point + xA = nextX; + yA = nextY; + absYA = FastMath.abs(yA); + agingA = 0; + ++agingB; + + // update the sign change index + signChangeIndex++; + + } + + } + + } + + /** Guess an x value by nth order inverse polynomial interpolation. + *

+ * The x value is guessed by evaluating polynomial Q(y) at y = targetY, where Q + * is built such that for all considered points (xi, yi), + * Q(yi) = xi. + *

+ * @param targetY target value for y + * @param x reference points abscissas for interpolation, + * note that this array is modified during computation + * @param y reference points ordinates for interpolation + * @param start start index of the points to consider (inclusive) + * @param end end index of the points to consider (exclusive) + * @return guessed root (will be a NaN if two points share the same y) + */ + private double guessX(final double targetY, final double[] x, final double[] y, + final int start, final int end) { + + // compute Q Newton coefficients by divided differences + for (int i = start; i < end - 1; ++i) { + final int delta = i + 1 - start; + for (int j = end - 1; j > i; --j) { + x[j] = (x[j] - x[j-1]) / (y[j] - y[j - delta]); + } + } + + // evaluate Q(targetY) + double x0 = 0; + for (int j = end - 1; j >= start; --j) { + x0 = x[j] + x0 * (targetY - y[j]); + } + + return x0; + + } + + /** {@inheritDoc} */ + public double solve(int maxEval, UnivariateRealFunction f, double min, + double max, AllowedSolutions allowedSolutions) { + this.allowed = allowedSolutions; + return super.solve(maxEval, f, min, max); + } + + /** {@inheritDoc} */ + public double solve(int maxEval, UnivariateRealFunction f, double min, + double max, double startValue, + AllowedSolutions allowedSolutions) { + this.allowed = allowedSolutions; + return super.solve(maxEval, f, min, max, startValue); + } + +} diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml index 61bae9a52..aa63d6979 100644 --- a/src/site/xdoc/changes.xml +++ b/src/site/xdoc/changes.xml @@ -52,6 +52,10 @@ The type attribute can be add,update,fix,remove. If the output is not quite correct, check for invisible trailing spaces! --> + + Added a Brent-like solver that has higher (user specified) order and does + bracket selection on the result: BracketingNthOrderBrentSolver. + Added a few shortcut methods and predicates to Dfp (abs, isZero, negativeOrNull, strictlyNegative, positiveOrNull, strictlyPositive). diff --git a/src/site/xdoc/userguide/analysis.xml b/src/site/xdoc/userguide/analysis.xml index daecf80b0..f8fc84b2f 100644 --- a/src/site/xdoc/userguide/analysis.xml +++ b/src/site/xdoc/userguide/analysis.xml @@ -73,6 +73,12 @@ yes no + + bracketing nth order Brent + variable order, guaranteed + yes + yes + Illinois Method univariate real-valued functions @@ -206,7 +212,8 @@ UnivariateRealFunction function = // some user defined function object final double relativeAccuracy = 1.0e-12; final double absoluteAccuracy = 1.0e-8; -UnivariateRealSolver solver = new PegasusSolver(relativeAccuracy, absoluteAccuracy); +final int maxOrder = 5; +UnivariateRealSolver solver = new BracketingNthOrderBrentSolver(relativeAccuracy, absoluteAccuracy, maxOrder); double c = solver.solve(100, function, 1.0, 5.0, AllowedSolutions.LEFT_SIDE);

Force bracketing, by refining a base solution found by a non-bracketing solver: @@ -221,9 +228,9 @@ double c = UnivariateRealSolverUtils.forceSide(100, function, baseRoot, 1.0, 5.0, AllowedSolutions.LEFT_SIDE);

- The BrentSolve uses the Brent-Dekker algorithm which is - fast and robust. This algorithm is recommended for most users. If there - are multiple roots in the interval, or there is a large domain of indeterminacy, the + The BrentSolver uses the Brent-Dekker algorithm which is + fast and robust. If there are multiple roots in the interval, + or there is a large domain of indeterminacy, the algorithm will converge to a random root in the interval without indication that there are problems. Interestingly, the examined text book implementations all disagree in details of the convergence @@ -232,6 +239,15 @@ double c = UnivariateRealSolverUtils.forceSide(100, function, get exactly the same root values as for other implementations of this algorithm.

+

+ The BracketingNthOrderBrentSolver uses an extension of the + Brent-Dekker algorithm which uses inverse nth order polynomial + interpolation instead of inverse quadratic interpolation, and which allows + selection of the side of the convergence interval for result bracketing. + This is now the recommended algorithm for most users since it has the + largest order, doesn't require derivatives, has guaranteed convergence + and allows result bracket selection. +

The SecantSolver uses a straightforward secant algorithm which does not bracket the search and therefore does not diff --git a/src/test/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolverTest.java b/src/test/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolverTest.java new file mode 100644 index 000000000..1488a289b --- /dev/null +++ b/src/test/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolverTest.java @@ -0,0 +1,181 @@ +/* + * 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.solvers; + +import org.apache.commons.math.analysis.DifferentiableUnivariateRealFunction; +import org.apache.commons.math.analysis.QuinticFunction; +import org.apache.commons.math.analysis.UnivariateRealFunction; +import org.apache.commons.math.exception.NumberIsTooSmallException; +import org.apache.commons.math.exception.TooManyEvaluationsException; +import org.apache.commons.math.util.FastMath; +import org.junit.Assert; +import org.junit.Test; + +/** + * Test case for {@link BracketingNthOrderBrentSolver bracketing nth order Brent} solver. + * + * @version $Id$ + */ +public final class BracketingNthOrderBrentSolverTest extends BaseSecantSolverAbstractTest { + /** {@inheritDoc} */ + protected UnivariateRealSolver getSolver() { + return new BracketingNthOrderBrentSolver(); + } + + /** {@inheritDoc} */ + protected int[] getQuinticEvalCounts() { + return new int[] {1, 3, 8, 1, 9, 4, 8, 1, 12, 1, 14}; + } + + @Test(expected=NumberIsTooSmallException.class) + public void testInsufficientOrder1() { + new BracketingNthOrderBrentSolver(1.0e-10, 1); + } + + @Test(expected=NumberIsTooSmallException.class) + public void testInsufficientOrder2() { + new BracketingNthOrderBrentSolver(1.0e-10, 1.0e-10, 1); + } + + @Test(expected=NumberIsTooSmallException.class) + public void testInsufficientOrder3() { + new BracketingNthOrderBrentSolver(1.0e-10, 1.0e-10, 1.0e-10, 1); + } + + @Test + public void testConstructorsOK() { + Assert.assertEquals(2, new BracketingNthOrderBrentSolver(1.0e-10, 2).getMaximalOrder()); + Assert.assertEquals(2, new BracketingNthOrderBrentSolver(1.0e-10, 1.0e-10, 2).getMaximalOrder()); + Assert.assertEquals(2, new BracketingNthOrderBrentSolver(1.0e-10, 1.0e-10, 1.0e-10, 2).getMaximalOrder()); + } + + @Test + public void testConvergenceOnFunctionAccuracy() { + BracketingNthOrderBrentSolver solver = + new BracketingNthOrderBrentSolver(1.0e-12, 1.0e-10, 0.001, 3); + QuinticFunction f = new QuinticFunction(); + double result = solver.solve(20, f, 0.2, 0.9, 0.4, AllowedSolutions.BELOW_SIDE); + Assert.assertEquals(0, f.value(result), solver.getFunctionValueAccuracy()); + Assert.assertTrue(f.value(result) <= 0); + Assert.assertTrue(result - 0.5 > solver.getAbsoluteAccuracy()); + result = solver.solve(20, f, -0.9, -0.2, -0.4, AllowedSolutions.ABOVE_SIDE); + Assert.assertEquals(0, f.value(result), solver.getFunctionValueAccuracy()); + Assert.assertTrue(f.value(result) >= 0); + Assert.assertTrue(result + 0.5 < -solver.getAbsoluteAccuracy()); + } + + @Test + public void testFasterThanNewton() { + // the following test functions come from Beny Neta's paper: + // "Several New Methods for solving Equations" + // intern J. Computer Math Vol 23 pp 265-282 + // available here: http://www.math.nps.navy.mil/~bneta/SeveralNewMethods.PDF + // the reference roots have been computed by the Dfp solver to more than + // 80 digits and checked with emacs (only the first 20 digits are reproduced here) + compare(new TestFunction(0.0, -2, 2) { + public double value(double x) { return FastMath.sin(x) - 0.5 * x; } + public double derivative(double x) { return FastMath.cos(x) - 0.5; } + }); + compare(new TestFunction(6.3087771299726890947, -5, 10) { + public double value(double x) { return FastMath.pow(x, 5) + x - 10000; } + public double derivative(double x) { return 5 * FastMath.pow(x, 4) + 1; } + }); + compare(new TestFunction(9.6335955628326951924, 0.001, 10) { + public double value(double x) { return FastMath.sqrt(x) - 1 / x - 3; } + public double derivative(double x) { return 0.5 / FastMath.sqrt(x) + 1 / (x * x); } + }); + compare(new TestFunction(2.8424389537844470678, -5, 5) { + public double value(double x) { return FastMath.exp(x) + x - 20; } + public double derivative(double x) { return FastMath.exp(x) + 1; } + }); + compare(new TestFunction(8.3094326942315717953, 0.001, 10) { + public double value(double x) { return FastMath.log(x) + FastMath.sqrt(x) - 5; } + public double derivative(double x) { return 1 / x + 0.5 / FastMath.sqrt(x); } + }); + compare(new TestFunction(1.4655712318767680266, -0.5, 1.5) { + public double value(double x) { return (x - 1) * x * x - 1; } + public double derivative(double x) { return (3 * x - 2) * x; } + }); + + } + + private void compare(TestFunction f) { + compare(f, f.getRoot(), f.getMin(), f.getMax()); + } + + private void compare(DifferentiableUnivariateRealFunction f, + double root, double min, double max) { + NewtonSolver newton = new NewtonSolver(1.0e-12); + BracketingNthOrderBrentSolver bracketing = + new BracketingNthOrderBrentSolver(1.0e-12, 1.0e-12, 1.0e-18, 5); + double resultN; + try { + resultN = newton.solve(100, f, min, max); + } catch (TooManyEvaluationsException tmee) { + resultN = Double.NaN; + } + double resultB; + try { + resultB = bracketing.solve(100, f, min, max); + } catch (TooManyEvaluationsException tmee) { + resultB = Double.NaN; + } + Assert.assertEquals(root, resultN, newton.getAbsoluteAccuracy()); + Assert.assertEquals(root, resultB, bracketing.getAbsoluteAccuracy()); + Assert.assertTrue(bracketing.getEvaluations() < newton.getEvaluations()); + } + + private static abstract class TestFunction implements DifferentiableUnivariateRealFunction { + + private final double root; + private final double min; + private final double max; + + protected TestFunction(final double root, final double min, final double max) { + this.root = root; + this.min = min; + this.max = max; + } + + public double getRoot() { + return root; + } + + public double getMin() { + return min; + } + + public double getMax() { + return max; + } + + public abstract double value(double x); + + public abstract double derivative(double x); + + public UnivariateRealFunction derivative() { + return new UnivariateRealFunction() { + public double value(double x) { + return derivative(x); + } + }; + } + + } + +}