Author: luc
Date: Mon Mar 17 15:14:07 2014
New Revision: 1578428
URL: http://svn.apache.org/r1578428
Log:
Improved brackting utility for univariate solvers.
Bracketing utility for univariate root solvers now returns a tighter
interval than before. It also allows choosing the search interval
expansion rate, supporting both linear and asymptotically exponential
rates.
Modified:
commons/proper/math/trunk/src/changes/changes.xml
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/solvers/UnivariateSolverUtils.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/solvers/UnivariateSolverUtilsTest.java
Modified: commons/proper/math/trunk/src/changes/changes.xml
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/changes/changes.xml?rev=1578428&r1=1578427&r2=1578428&view=diff
==============================================================================
--- commons/proper/math/trunk/src/changes/changes.xml (original)
+++ commons/proper/math/trunk/src/changes/changes.xml Mon Mar 17 15:14:07 2014
@@ -51,6 +51,11 @@ If the output is not quite correct, chec
</properties>
<body>
<release version="3.3" date="TBD" description="TBD">
+ <action dev="luc" type="update" >
+ Bracketing utility for univariate root solvers returns a tighter
interval than before.
+ It also allows choosing the search interval expansion rate, supporting
both linear
+ and asymptotically exponential rates.
+ </action>
<action dev="luc" type="fix" issue="MATH-1107" due-to="Bruce A Johnson">
Prevent penalties to grow multiplicatively in CMAES for out of bounds
points.
</action>
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/solvers/UnivariateSolverUtils.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/solvers/UnivariateSolverUtils.java?rev=1578428&r1=1578427&r2=1578428&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/solvers/UnivariateSolverUtils.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/solvers/UnivariateSolverUtils.java
Mon Mar 17 15:14:07 2014
@@ -171,31 +171,16 @@ public class UnivariateSolverUtils {
}
/**
- * This method attempts to find two values a and b satisfying <ul>
- * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
- * <li> <code> f(a) * f(b) < 0 </code></li>
- * </ul>
- * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
- * and <code>b</code> bracket a root of f.
- * <p>
- * The algorithm starts by setting
- * <code>a := initial -1; b := initial +1,</code> examines the value of the
- * function at <code>a</code> and <code>b</code> and keeps moving
- * the endpoints out by one unit each time through a loop that terminates
- * when one of the following happens: <ul>
- * <li> <code> f(a) * f(b) < 0 </code> -- success!</li>
- * <li> <code> a = lower </code> and <code> b = upper</code>
- * -- NoBracketingException </li>
- * <li> <code> Integer.MAX_VALUE</code> iterations elapse
- * -- NoBracketingException </li>
- * </ul></p>
- * <p>
+ * This method simply calls {@link #bracket(UnivariateFunction, double,
double, double,
+ * double, double, int) bracket(function, initial, lowerBound, upperBound,
q, r, maximumIterations)}
+ * with {@code q} and {@code r} set to 1.0 and {@code maximumIterations}
set to {@code Integer.MAX_VALUE}.
* <strong>Note: </strong> this method can take
* <code>Integer.MAX_VALUE</code> iterations to throw a
* <code>ConvergenceException.</code> Unless you are confident that there
* is a root between <code>lowerBound</code> and <code>upperBound</code>
* near <code>initial,</code> it is better to use
- * {@link #bracket(UnivariateFunction, double, double, double, int)},
+ * {@link #bracket(UnivariateFunction, double, double, double,
+ * double, int) bracket(function, initial, lowerBound, upperBound, delta,
maximumIterations)},
* explicitly specifying the maximum number of iterations.</p>
*
* @param function Function.
@@ -215,28 +200,13 @@ public class UnivariateSolverUtils {
throws NullArgumentException,
NotStrictlyPositiveException,
NoBracketingException {
- return bracket(function, initial, lowerBound, upperBound,
Integer.MAX_VALUE);
+ return bracket(function, initial, lowerBound, upperBound, 1.0, 1.0,
Integer.MAX_VALUE);
}
/**
- * This method attempts to find two values a and b satisfying <ul>
- * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
- * <li> <code> f(a) * f(b) <= 0 </code> </li>
- * </ul>
- * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
- * and <code>b</code> bracket a root of f.
- * <p>
- * The algorithm starts by setting
- * <code>a := initial -1; b := initial +1,</code> examines the value of the
- * function at <code>a</code> and <code>b</code> and keeps moving
- * the endpoints out by one unit each time through a loop that terminates
- * when one of the following happens: <ul>
- * <li> <code> f(a) * f(b) <= 0 </code> -- success!</li>
- * <li> <code> a = lower </code> and <code> b = upper</code>
- * -- NoBracketingException </li>
- * <li> <code> maximumIterations</code> iterations elapse
- * -- NoBracketingException </li></ul></p>
- *
+ * This method simply calls {@link #bracket(UnivariateFunction, double,
double, double,
+ * double, double, int) bracket(function, initial, lowerBound, upperBound,
q, r, maximumIterations)}
+ * with {@code q} and {@code r} set to 1.0.
* @param function Function.
* @param initial Initial midpoint of interval being expanded to
* bracket a root.
@@ -257,38 +227,132 @@ public class UnivariateSolverUtils {
throws NullArgumentException,
NotStrictlyPositiveException,
NoBracketingException {
+ return bracket(function, initial, lowerBound, upperBound, 1.0, 1.0,
maximumIterations);
+ }
+
+ /**
+ * This method attempts to find two values a and b satisfying <ul>
+ * <li> {@code lowerBound <= a < initial < b <= upperBound} </li>
+ * <li> {@code f(a) * f(b) <= 0} </li>
+ * </ul>
+ * If {@code f} is continuous on {@code [a,b]}, this means that {@code a}
+ * and {@code b} bracket a root of {@code f}.
+ * <p>
+ * The algorithm checks the sign of \( f(l_k) \) and \( f(u_k) \) for
increasing
+ * values of k, where \( l_k = max(lower, initial - \delta_k) \),
+ * \( u_k = min(upper, initial + \delta_k) \), using recurrence
+ * \( \delta_{k+1} = r \delta_k + q, \delta_0 = 0\) and starting search
with \( k=1 \).
+ * The algorithm stops when one of the following happens: <ul>
+ * <li> at least one positive and one negative value have been found --
success!</li>
+ * <li> both endpoints have reached their respective limites --
NoBracketingException </li>
+ * <li> {@code maximumIterations} iterations elapse --
NoBracketingException </li></ul></p>
+ * <p>
+ * If different signs are found at first iteration ({@code k=1}), then the
returned
+ * interval will be \( [a, b] = [l_1, u_1] \). If different signs are
found at a later
+ * iteration ({code k>1}, then the returned interval will be either
+ * \( [a, b] = [l_{k+1}, l_{k}] \) or ( [a, b] = [u_{k}, u_{k+1}] \). A
root solver called
+ * with these parameters will therefore start with the smallest bracketing
interval known
+ * at this step.
+ * </p>
+ * <p>
+ * Interval expansion rate is tuned by changing the recurrence parameters
{@code r} and
+ * {@code q}. When the multiplicative factor {@code r} is set to 1, the
sequence is a
+ * simple arithmetic sequence with linear increase. When the
multiplicative factor {@code r}
+ * is larger than 1, the sequence has an asymtotically exponential rate.
Note than the
+ * additive parameter {@code q} should never be set to zero, otherwise the
interval would
+ * degenerate to the single initial point for all values of {@code k}.
+ * </p>
+ * <p>
+ * As a rule of thumb, when the location of the root is expected to be
approximately known
+ * within some error margin, {@code r} should be set to 1 and {@code q}
should be set to the
+ * order of magnitude of the error margin. When the location of the root
is really a wild guess,
+ * then {@code r} should be set to a value larger than 1 (typically 2 to
double the interval
+ * length at each iteration) and {@code q} should be set according to half
the initial
+ * search interval length.
+ * </p>
+ * <p>
+ * As an example, if we consider the trivial function {@code f(x) = 1 - x}
and use
+ * {@code initial = 4}, {@code r = 1}, {@code q = 2}, the algorithm will
compute
+ * {@code f(4-2) = f(2) = -1} and {@code f(4+2) = f(6) = -5} for {@code k
= 1}, then
+ * {@code f(4-4) = f(0) = +1} and {@code f(4+4) = f(8) = -7} for {@code k
= 2}. Then it will
+ * return the interval {@code [0, 2]} as the smallest one known to be
bracketing the root.
+ * As shown by this example, the initial value (here {@code 4}) may lie
outside of the returned
+ * bracketing interval.
+ * </p>
+ * @param function function to check
+ * @param initial Initial midpoint of interval being expanded to
+ * bracket a root.
+ * @param lowerBound Lower bound (a is never lower than this value).
+ * @param upperBound Upper bound (b never is greater than this
+ * value).
+ * @param q additive offset used to compute bounds sequence (must be
strictly positive)
+ * @param r multiplicative factor used to compute bounds sequence
+ * @param maximumIterations Maximum number of iterations to perform
+ * @return a two element array holding the bracketing values.
+ * @exception NoBracketingException if function cannot be bracketed in the
search interval
+ */
+ public static double[] bracket(final UnivariateFunction function, final
double initial,
+ final double lowerBound, final double
upperBound,
+ final double q, final double r, final int
maximumIterations)
+ throws NoBracketingException {
+
if (function == null) {
throw new NullArgumentException(LocalizedFormats.FUNCTION);
}
+ if (q <= 0) {
+ throw new NotStrictlyPositiveException(q);
+ }
if (maximumIterations <= 0) {
throw new
NotStrictlyPositiveException(LocalizedFormats.INVALID_MAX_ITERATIONS,
maximumIterations);
}
verifySequence(lowerBound, initial, upperBound);
- double a = initial;
- double b = initial;
- double fa;
- double fb;
- int numIterations = 0;
-
- do {
- a = FastMath.max(a - 1.0, lowerBound);
- b = FastMath.min(b + 1.0, upperBound);
- fa = function.value(a);
-
- fb = function.value(b);
- ++numIterations;
- } while ((fa * fb > 0.0) && (numIterations < maximumIterations) &&
- ((a > lowerBound) || (b < upperBound)));
-
- if (fa * fb > 0.0) {
- throw new NoBracketingException(LocalizedFormats.FAILED_BRACKETING,
- a, b, fa, fb,
- numIterations, maximumIterations,
initial,
- lowerBound, upperBound);
+ // initialize the recurrence
+ double a = initial;
+ double b = initial;
+ double fa = Double.NaN;
+ double fb = Double.NaN;
+ double delta = 0;
+
+ for (int numIterations = 0;
+ (numIterations < maximumIterations) && (a > lowerBound || b >
upperBound);
+ ++numIterations) {
+
+ final double previousA = a;
+ final double previousFa = fa;
+ final double previousB = b;
+ final double previousFb = fb;
+
+ delta = r * delta + q;
+ a = FastMath.max(initial - delta, lowerBound);
+ b = FastMath.min(initial + delta, upperBound);
+ fa = function.value(a);
+ fb = function.value(b);
+
+ if (numIterations == 0) {
+ // at first iteration, we don't have a previous interval
+ // we simply compare both sides of the initial interval
+ if (fa * fb <= 0) {
+ // the first interval already brackets a root
+ return new double[] { a, b };
+ }
+ } else {
+ // we have a previous interval with constant sign and expand
it,
+ // we expect sign changes to occur at boundaries
+ if (fa * previousFa <= 0) {
+ // sign change detected at near lower bound
+ return new double[] { a, previousA };
+ } else if (fb * previousFb <= 0) {
+ // sign change detected at near upper bound
+ return new double[] { previousB, b };
+ }
+ }
+
}
- return new double[] {a, b};
+ // no bracketing found
+ throw new NoBracketingException(a, b, fa, fb);
+
}
/**
Modified:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/solvers/UnivariateSolverUtilsTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/solvers/UnivariateSolverUtilsTest.java?rev=1578428&r1=1578427&r2=1578428&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/solvers/UnivariateSolverUtilsTest.java
(original)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/solvers/UnivariateSolverUtilsTest.java
Mon Mar 17 15:14:07 2014
@@ -21,6 +21,7 @@ import org.apache.commons.math3.analysis
import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.analysis.function.Sin;
import org.apache.commons.math3.exception.MathIllegalArgumentException;
+import org.apache.commons.math3.exception.NoBracketingException;
import org.apache.commons.math3.util.FastMath;
import org.junit.Assert;
import org.junit.Test;
@@ -87,6 +88,56 @@ public class UnivariateSolverUtilsTest {
}
@Test
+ public void testBracketCentered() {
+ double initial = 0.1;
+ double[] result = UnivariateSolverUtils.bracket(sin, initial, -2.0,
2.0, 0.2, 1.0, 100);
+ Assert.assertTrue(result[0] < initial);
+ Assert.assertTrue(result[1] > initial);
+ Assert.assertTrue(sin.value(result[0]) < 0);
+ Assert.assertTrue(sin.value(result[1]) > 0);
+ }
+
+ @Test
+ public void testBracketLow() {
+ double initial = 0.5;
+ double[] result = UnivariateSolverUtils.bracket(sin, initial, -2.0,
2.0, 0.2, 1.0, 100);
+ Assert.assertTrue(result[0] < initial);
+ Assert.assertTrue(result[1] < initial);
+ Assert.assertTrue(sin.value(result[0]) < 0);
+ Assert.assertTrue(sin.value(result[1]) > 0);
+ }
+
+ @Test
+ public void testBracketHigh(){
+ double initial = -0.5;
+ double[] result = UnivariateSolverUtils.bracket(sin, initial, -2.0,
2.0, 0.2, 1.0, 100);
+ Assert.assertTrue(result[0] > initial);
+ Assert.assertTrue(result[1] > initial);
+ Assert.assertTrue(sin.value(result[0]) < 0);
+ Assert.assertTrue(sin.value(result[1]) > 0);
+ }
+
+ @Test(expected=NoBracketingException.class)
+ public void testBracketLinear(){
+ UnivariateSolverUtils.bracket(new UnivariateFunction() {
+ public double value(double x) {
+ return 1 - x;
+ }
+ }, 1000, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY, 1.0, 1.0,
100);
+ }
+
+ @Test
+ public void testBracketExponential(){
+ double[] result = UnivariateSolverUtils.bracket(new
UnivariateFunction() {
+ public double value(double x) {
+ return 1 - x;
+ }
+ }, 1000, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY, 1.0, 2.0,
10);
+ Assert.assertTrue(result[0] <= 1);
+ Assert.assertTrue(result[1] >= 1);
+ }
+
+ @Test
public void testBracketEndpointRoot() {
double[] result = UnivariateSolverUtils.bracket(sin, 1.5, 0, 2.0);
Assert.assertEquals(0.0, sin.value(result[0]), 1.0e-15);
@@ -97,18 +148,28 @@ public class UnivariateSolverUtilsTest {
public void testNullFunction() {
UnivariateSolverUtils.bracket(null, 1.5, 0, 2.0);
}
-
+
@Test(expected=MathIllegalArgumentException.class)
public void testBadInitial() {
UnivariateSolverUtils.bracket(sin, 2.5, 0, 2.0);
}
-
+
+ @Test(expected=MathIllegalArgumentException.class)
+ public void testBadAdditive() {
+ UnivariateSolverUtils.bracket(sin, 1.0, -2.0, 3.0, -1.0, 1.0, 100);
+ }
+
+ @Test(expected=NoBracketingException.class)
+ public void testIterationExceeded() {
+ UnivariateSolverUtils.bracket(sin, 1.0, -2.0, 3.0, 1.0e-5, 1.0, 100);
+ }
+
@Test(expected=MathIllegalArgumentException.class)
public void testBadEndpoints() {
// endpoints not valid
UnivariateSolverUtils.bracket(sin, 1.5, 2.0, 1.0);
}
-
+
@Test(expected=MathIllegalArgumentException.class)
public void testBadMaximumIterations() {
// bad maximum iterations
