BrentOptimizer.java
/*
* 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.math4.legacy.optim.univariate;
import org.apache.commons.math4.legacy.analysis.UnivariateFunction;
import org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException;
import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
import org.apache.commons.math4.legacy.optim.ConvergenceChecker;
import org.apache.commons.math4.legacy.optim.nonlinear.scalar.GoalType;
import org.apache.commons.math4.core.jdkmath.JdkMath;
import org.apache.commons.numbers.core.Precision;
/**
* 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.
* <br>
* This code is an adaptation, partly based on the Python code from SciPy
* (module "optimize.py" v0.5); the original algorithm is also modified
* <ul>
* <li>to use an initial guess provided by the user,</li>
* <li>to ensure that the best point encountered is the one returned.</li>
* </ul>
*
* @since 2.0
*/
public class BrentOptimizer extends UnivariateOptimizer {
/**
* Golden section.
*/
private static final double GOLDEN_SECTION = 0.5 * (3 - JdkMath.sqrt(5));
/**
* Minimum relative tolerance.
*/
private static final double MIN_RELATIVE_TOLERANCE = 2 * JdkMath.ulp(1d);
/**
* Relative threshold.
*/
private final double relativeThreshold;
/**
* Absolute threshold.
*/
private final double absoluteThreshold;
/**
* The arguments are used implement the original stopping criterion
* of Brent's algorithm.
* {@code abs} and {@code rel} define a tolerance
* {@code tol = rel |x| + abs}. {@code rel} should be no smaller than
* <em>2 macheps</em> and preferably not much less than <em>sqrt(macheps)</em>,
* where <em>macheps</em> is the relative machine precision. {@code abs} must
* be positive.
*
* @param rel Relative threshold.
* @param abs Absolute threshold.
* @param checker Additional, user-defined, convergence checking
* procedure.
* @throws NotStrictlyPositiveException if {@code abs <= 0}.
* @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
*/
public BrentOptimizer(double rel,
double abs,
ConvergenceChecker<UnivariatePointValuePair> checker) {
super(checker);
if (rel < MIN_RELATIVE_TOLERANCE) {
throw new NumberIsTooSmallException(rel, MIN_RELATIVE_TOLERANCE, true);
}
if (abs <= 0) {
throw new NotStrictlyPositiveException(abs);
}
relativeThreshold = rel;
absoluteThreshold = abs;
}
/**
* The arguments are used for implementing the original stopping criterion
* of Brent's algorithm.
* {@code abs} and {@code rel} define a tolerance
* {@code tol = rel |x| + abs}. {@code rel} should be no smaller than
* <em>2 macheps</em> and preferably not much less than <em>sqrt(macheps)</em>,
* where <em>macheps</em> is the relative machine precision. {@code abs} must
* be positive.
*
* @param rel Relative threshold.
* @param abs Absolute threshold.
* @throws NotStrictlyPositiveException if {@code abs <= 0}.
* @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
*/
public BrentOptimizer(double rel,
double abs) {
this(rel, abs, null);
}
/** {@inheritDoc} */
@Override
protected UnivariatePointValuePair doOptimize() {
final boolean isMinim = getGoalType() == GoalType.MINIMIZE;
final double lo = getMin();
final double mid = getStartValue();
final double hi = getMax();
final UnivariateFunction func = getObjectiveFunction();
// Optional additional convergence criteria.
final ConvergenceChecker<UnivariatePointValuePair> checker
= getConvergenceChecker();
double a;
double b;
if (lo < hi) {
a = lo;
b = hi;
} else {
a = hi;
b = lo;
}
double x = mid;
double v = x;
double w = x;
double d = 0;
double e = 0;
double fx = func.value(x);
if (!isMinim) {
fx = -fx;
}
double fv = fx;
double fw = fx;
UnivariatePointValuePair previous = null;
UnivariatePointValuePair current
= new UnivariatePointValuePair(x, isMinim ? fx : -fx);
// Best point encountered so far (which is the initial guess).
UnivariatePointValuePair best = current;
while (true) {
final double m = 0.5 * (a + b);
final double tol1 = relativeThreshold * JdkMath.abs(x) + absoluteThreshold;
final double tol2 = 2 * tol1;
// Default stopping criterion.
final boolean stop = JdkMath.abs(x - m) <= tol2 - 0.5 * (b - a);
if (!stop) {
double p = 0;
double q = 0;
double r = 0;
double u = 0;
if (JdkMath.abs(e) > tol1) { // 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 (p > q * (a - x) &&
p < q * (b - x) &&
JdkMath.abs(p) < JdkMath.abs(0.5 * q * r)) {
// Parabolic interpolation step.
d = p / q;
u = x + d;
// f must not be evaluated too close to a or b.
if (u - a < tol2 || b - u < tol2) {
if (x <= m) {
d = tol1;
} else {
d = -tol1;
}
}
} else {
// Golden section step.
if (x < m) {
e = b - x;
} else {
e = a - x;
}
d = GOLDEN_SECTION * e;
}
} else {
// Golden section step.
if (x < m) {
e = b - x;
} else {
e = a - x;
}
d = GOLDEN_SECTION * e;
}
// Update by at least "tol1".
if (JdkMath.abs(d) < tol1) {
if (d >= 0) {
u = x + tol1;
} else {
u = x - tol1;
}
} else {
u = x + d;
}
double fu = func.value(u);
if (!isMinim) {
fu = -fu;
}
// User-defined convergence checker.
previous = current;
current = new UnivariatePointValuePair(u, isMinim ? fu : -fu);
best = best(best,
best(previous,
current,
isMinim),
isMinim);
if (checker != null && checker.converged(getIterations(), previous, current)) {
return best;
}
// 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 ||
Precision.equals(w, x)) {
v = w;
fv = fw;
w = u;
fw = fu;
} else if (fu <= fv ||
Precision.equals(v, x) ||
Precision.equals(v, w)) {
v = u;
fv = fu;
}
}
} else { // Default termination (Brent's criterion).
return best(best,
best(previous,
current,
isMinim),
isMinim);
}
incrementIterationCount();
}
}
/**
* Selects the best of two points.
*
* @param a Point and value.
* @param b Point and value.
* @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}. When {@code a} and {@code b} have the same function
* value, {@code a} is returned.
*/
private UnivariatePointValuePair best(UnivariatePointValuePair a,
UnivariatePointValuePair b,
boolean isMinim) {
if (a == null) {
return b;
}
if (b == null) {
return a;
}
if (isMinim) {
return a.getValue() <= b.getValue() ? a : b;
} else {
return a.getValue() >= b.getValue() ? a : b;
}
}
}