490 lines
17 KiB
Python
490 lines
17 KiB
Python
"""Classical and state-of-the-art bracketing root finders, all with a common
|
|
signature and a shared function-evaluation counter.
|
|
|
|
Every method returns (x, nfev, converged).
|
|
|
|
Included:
|
|
bisection - the null model
|
|
regula_falsi - the raw "step proportional to dissimilarity" rule
|
|
illinois - Hafner/Dowell 1971 stalling fix (halve the stale end)
|
|
pegasus - Dowell & Jarratt 1972
|
|
anderson_bjorck - Anderson & Bjorck 1973
|
|
ridders - Ridders 1979
|
|
brent - Brent 1973 (own implementation, verified vs scipy)
|
|
toms748 - Alefeld, Potra & Shi 1995 (via scipy)
|
|
itp - Oliveira & Takahashi, ACM TOMS 2021 (own implementation)
|
|
|
|
The ITP method is the current state of the art for bracketed root finding:
|
|
superlinear average performance with a worst case that is provably no worse
|
|
than bisection. Note that it carries THREE tuning parameters (kappa1, kappa2,
|
|
n0); the authors' recommended defaults are used here and are reported as such.
|
|
"""
|
|
from __future__ import annotations
|
|
|
|
import math
|
|
|
|
from scipy.optimize import toms748 as _scipy_toms748
|
|
|
|
|
|
class Counter:
|
|
"""Counts calls to f and caches nothing (evaluations are the cost model)."""
|
|
|
|
__slots__ = ("f", "n", "limit")
|
|
|
|
def __init__(self, f, limit):
|
|
self.f = f
|
|
self.n = 0
|
|
self.limit = limit
|
|
|
|
def __call__(self, x):
|
|
if self.n >= self.limit:
|
|
raise _Budget()
|
|
self.n += 1
|
|
return self.f(x)
|
|
|
|
|
|
class _Budget(Exception):
|
|
pass
|
|
|
|
|
|
def _converged(a, b, xtol):
|
|
return abs(b - a) <= xtol
|
|
|
|
|
|
# --------------------------------------------------------------------------
|
|
# null model
|
|
# --------------------------------------------------------------------------
|
|
def bisection(f, a, b, xtol, maxfev):
|
|
g = Counter(f, maxfev)
|
|
try:
|
|
fa, fb = g(a), g(b)
|
|
if fa == 0.0:
|
|
return a, g.n, True
|
|
if fb == 0.0:
|
|
return b, g.n, True
|
|
x = 0.5 * (a + b)
|
|
while not _converged(a, b, xtol):
|
|
x = 0.5 * (a + b)
|
|
fx = g(x)
|
|
if fx == 0.0:
|
|
return x, g.n, True
|
|
if fa * fx < 0.0:
|
|
b, fb = x, fx
|
|
else:
|
|
a, fa = x, fx
|
|
return 0.5 * (a + b), g.n, True
|
|
except _Budget:
|
|
return 0.5 * (a + b), g.n, False
|
|
|
|
|
|
# --------------------------------------------------------------------------
|
|
# the false-position family (all are "similarity-proportional steps" with
|
|
# different stale-endpoint corrections)
|
|
# --------------------------------------------------------------------------
|
|
def _false_position_family(f, a, b, xtol, maxfev, correction):
|
|
g = Counter(f, maxfev)
|
|
x = 0.5 * (a + b) # defined before the first evaluation: a budget that
|
|
# runs out on f(a) must give an honest failure, not
|
|
# an UnboundLocalError
|
|
try:
|
|
fa, fb = g(a), g(b)
|
|
if fa == 0.0:
|
|
return a, g.n, True
|
|
if fb == 0.0:
|
|
return b, g.n, True
|
|
while not _converged(a, b, xtol):
|
|
denom = fb - fa
|
|
if denom == 0.0:
|
|
x = 0.5 * (a + b)
|
|
else:
|
|
x = (a * fb - b * fa) / denom
|
|
lo, hi = (a, b) if a < b else (b, a)
|
|
if not (lo < x < hi) or not math.isfinite(x):
|
|
x = 0.5 * (a + b)
|
|
fx = g(x)
|
|
if fx == 0.0:
|
|
return x, g.n, True
|
|
if fa * fx < 0.0:
|
|
# b replaced; a became stale -> correct fa
|
|
fa = correction(fa, fb, fx)
|
|
b, fb = x, fx
|
|
else:
|
|
fb = correction(fb, fa, fx)
|
|
a, fa = x, fx
|
|
return 0.5 * (a + b), g.n, True
|
|
except _Budget:
|
|
return x, g.n, False
|
|
|
|
|
|
def regula_falsi(f, a, b, xtol, maxfev):
|
|
return _false_position_family(f, a, b, xtol, maxfev,
|
|
lambda f_stale, f_moved, fx: f_stale)
|
|
|
|
|
|
def illinois(f, a, b, xtol, maxfev):
|
|
return _false_position_family(f, a, b, xtol, maxfev,
|
|
lambda f_stale, f_moved, fx: 0.5 * f_stale)
|
|
|
|
|
|
def pegasus(f, a, b, xtol, maxfev):
|
|
return _false_position_family(
|
|
f, a, b, xtol, maxfev,
|
|
lambda f_stale, f_moved, fx: f_stale * f_moved / (f_moved + fx))
|
|
|
|
|
|
def _ab_corr(f_stale, f_moved, fx):
|
|
m = 1.0 - fx / f_moved
|
|
if m <= 0.0:
|
|
m = 0.5
|
|
return f_stale * m
|
|
|
|
|
|
def anderson_bjorck(f, a, b, xtol, maxfev):
|
|
return _false_position_family(f, a, b, xtol, maxfev, _ab_corr)
|
|
|
|
|
|
# --------------------------------------------------------------------------
|
|
# Ridders 1979
|
|
# --------------------------------------------------------------------------
|
|
def ridders(f, a, b, xtol, maxfev):
|
|
g = Counter(f, maxfev)
|
|
x = 0.5 * (a + b)
|
|
try:
|
|
fa, fb = g(a), g(b)
|
|
if fa == 0.0:
|
|
return a, g.n, True
|
|
if fb == 0.0:
|
|
return b, g.n, True
|
|
while not _converged(a, b, xtol):
|
|
c = 0.5 * (a + b)
|
|
fc = g(c)
|
|
if fc == 0.0:
|
|
return c, g.n, True
|
|
d = fc * fc - fa * fb
|
|
if d <= 0.0:
|
|
x = c
|
|
else:
|
|
s = math.sqrt(d)
|
|
x = c + (c - a) * (math.copysign(1.0, fa - fb) * fc) / s
|
|
fx = g(x)
|
|
if fx == 0.0:
|
|
return x, g.n, True
|
|
# rebuild the tightest bracket among a, c, x, b
|
|
if fc * fx < 0.0:
|
|
a, fa, b, fb = (c, fc, x, fx) if c < x else (x, fx, c, fc)
|
|
elif fa * fx < 0.0:
|
|
b, fb = x, fx
|
|
else:
|
|
a, fa = x, fx
|
|
return 0.5 * (a + b), g.n, True
|
|
except _Budget:
|
|
return x, g.n, False
|
|
|
|
|
|
# --------------------------------------------------------------------------
|
|
# Brent 1973
|
|
# --------------------------------------------------------------------------
|
|
def brent(f, a, b, xtol, maxfev, stats=None):
|
|
"""stats: optional dict; on return it carries 'bisect' = how many of the
|
|
steps taken were Brent's midpoint fallback rather than an interpolation.
|
|
Used to explain, rather than assert, where the evaluations went."""
|
|
g = Counter(f, maxfev)
|
|
n_bisect = 0
|
|
b_ = b
|
|
try:
|
|
fa, fb = g(a), g(b)
|
|
if fa == 0.0:
|
|
return a, g.n, True
|
|
if fb == 0.0:
|
|
return b, g.n, True
|
|
if stats is not None:
|
|
stats['bisect'] = 0
|
|
c, fc = a, fa
|
|
d = e = b - a
|
|
while True:
|
|
if fb * fc > 0.0:
|
|
c, fc = a, fa
|
|
d = e = b - a
|
|
if abs(fc) < abs(fb):
|
|
a, fa = b, fb
|
|
b, fb = c, fc
|
|
c, fc = a, fa
|
|
tol1 = 2.0 * 2.22e-16 * abs(b) + 0.5 * xtol
|
|
xm = 0.5 * (c - b)
|
|
if abs(xm) <= tol1 or fb == 0.0:
|
|
if stats is not None:
|
|
stats['bisect'] = n_bisect
|
|
return b, g.n, True
|
|
if abs(e) >= tol1 and abs(fa) > abs(fb):
|
|
s = fb / fa
|
|
if a == c:
|
|
p = 2.0 * xm * s
|
|
q = 1.0 - s
|
|
else:
|
|
q = fa / fc
|
|
r = fb / fc
|
|
p = s * (2.0 * xm * q * (q - r) - (b - a) * (r - 1.0))
|
|
q = (q - 1.0) * (r - 1.0) * (s - 1.0)
|
|
if p > 0.0:
|
|
q = -q
|
|
p = abs(p)
|
|
if 2.0 * p < min(3.0 * xm * q - abs(tol1 * q), abs(e * q)):
|
|
e = d
|
|
d = p / q
|
|
else:
|
|
d = xm
|
|
e = d
|
|
n_bisect += 1
|
|
else:
|
|
d = xm
|
|
e = d
|
|
n_bisect += 1
|
|
a, fa = b, fb
|
|
b = b + (d if abs(d) > tol1 else math.copysign(tol1, xm))
|
|
b_ = b
|
|
fb = g(b)
|
|
except _Budget:
|
|
if stats is not None:
|
|
stats['bisect'] = n_bisect
|
|
return b_, g.n, False
|
|
|
|
|
|
# --------------------------------------------------------------------------
|
|
# TOMS 748 (Alefeld, Potra & Shi 1995) - via scipy, counted by the wrapper
|
|
# --------------------------------------------------------------------------
|
|
def toms748(f, a, b, xtol, maxfev):
|
|
g = Counter(f, maxfev)
|
|
try:
|
|
x, r = _scipy_toms748(g, a, b, xtol=xtol, rtol=8.9e-16,
|
|
maxiter=maxfev, full_output=True, disp=False)
|
|
return x, g.n, bool(r.converged)
|
|
except _Budget:
|
|
return 0.5 * (a + b), g.n, False
|
|
except Exception:
|
|
return 0.5 * (a + b), g.n, False
|
|
|
|
|
|
# --------------------------------------------------------------------------
|
|
# ITP - Oliveira & Takahashi, ACM TOMS 47(1), 2021
|
|
# Three tuning parameters, authors' recommended defaults.
|
|
# --------------------------------------------------------------------------
|
|
def itp(f, a, b, xtol, maxfev, kappa1=None, kappa2=2.0, n0=1):
|
|
"""ITP is natively written against a half-width epsilon: it terminates on
|
|
(b - a) <= 2*eps. Every other method here terminates on (b - a) <= xtol.
|
|
Feeding xtol straight in as eps would let ITP stop at twice the bracket
|
|
width of its competitors and win the benchmark on the stopping rule rather
|
|
than on the search. eps = xtol/2 is therefore used, which makes the final
|
|
bracket identical across all methods. Requires a < b."""
|
|
if not a < b:
|
|
raise ValueError("itp requires a < b")
|
|
g = Counter(f, maxfev)
|
|
x = 0.5 * (a + b)
|
|
eps = 0.5 * xtol
|
|
try:
|
|
if kappa1 is None:
|
|
kappa1 = 0.2 / (b - a)
|
|
fa, fb = g(a), g(b)
|
|
if fa == 0.0:
|
|
return a, g.n, True
|
|
if fb == 0.0:
|
|
return b, g.n, True
|
|
n_half = math.ceil(math.log2((b - a) / (2.0 * eps))) if b - a > 2 * eps else 0
|
|
n_max = n_half + n0
|
|
j = 0
|
|
while (b - a) > 2.0 * eps:
|
|
# interpolation (regula falsi)
|
|
x_f = (b * fa - a * fb) / (fa - fb)
|
|
# truncation
|
|
x_half = 0.5 * (a + b)
|
|
sigma = math.copysign(1.0, x_half - x_f)
|
|
delta = kappa1 * (b - a) ** kappa2
|
|
x_t = x_f + sigma * delta if delta <= abs(x_half - x_f) else x_half
|
|
# projection
|
|
r = 2.0 ** (n_max - j) * eps - 0.5 * (b - a)
|
|
x = x_t if abs(x_t - x_half) <= r else x_half - sigma * r
|
|
fx = g(x)
|
|
if fx == 0.0:
|
|
return x, g.n, True
|
|
if fx * fa < 0.0:
|
|
b, fb = x, fx
|
|
else:
|
|
a, fa = x, fx
|
|
j += 1
|
|
return 0.5 * (a + b), g.n, True
|
|
except _Budget:
|
|
return x, g.n, False
|
|
|
|
|
|
# --------------------------------------------------------------------------
|
|
# Chandrupatla 1997 - the hybrid quadratic/bisection method. It owns one of
|
|
# the two benchmark sets used here and is the closest relative of the proposed
|
|
# rule, so its absence from the comparison would be the first thing a reader
|
|
# asked about.
|
|
# --------------------------------------------------------------------------
|
|
def chandrupatla(f, a, b, xtol, maxfev):
|
|
g = Counter(f, maxfev)
|
|
x1, x2 = b, a
|
|
# defined before the first evaluation: a budget exhausted on f(b) must be
|
|
# an honest failure, not an UnboundLocalError in the handler
|
|
f1 = f2 = 0.0
|
|
try:
|
|
f1, f2 = g(x1), g(x2)
|
|
if f1 == 0.0:
|
|
return x1, g.n, True
|
|
if f2 == 0.0:
|
|
return x2, g.n, True
|
|
x3, f3 = x2, f2
|
|
t = 0.5
|
|
while True:
|
|
x = x1 + t * (x2 - x1)
|
|
fx = g(x)
|
|
if (fx > 0.0) == (f1 > 0.0):
|
|
x3, f3 = x1, f1
|
|
x1, f1 = x, fx
|
|
else:
|
|
x3, f3 = x2, f2
|
|
x2, f2 = x1, f1
|
|
x1, f1 = x, fx
|
|
xm, fm = (x1, f1) if abs(f1) < abs(f2) else (x2, f2)
|
|
if fm == 0.0:
|
|
return xm, g.n, True
|
|
tol = 2.0 * 2.22e-16 * abs(xm) + xtol
|
|
span = abs(x2 - x1)
|
|
if span == 0.0:
|
|
return xm, g.n, True
|
|
tl = tol / span
|
|
if tl > 0.5:
|
|
return xm, g.n, True
|
|
# accept inverse quadratic interpolation only inside the region
|
|
# where Chandrupatla shows it is safe
|
|
t = 0.5
|
|
if x3 != x2 and f3 != f2:
|
|
xi = (x1 - x2) / (x3 - x2)
|
|
ph = (f1 - f2) / (f3 - f2)
|
|
if 0.0 < xi < 1.0 and 1.0 - math.sqrt(1.0 - xi) < ph < math.sqrt(xi):
|
|
d1, d2, d3 = f1 - f2, f1 - f3, f2 - f3
|
|
if d1 != 0.0 and d2 != 0.0 and d3 != 0.0:
|
|
t = (f1 / d1) * (f3 / (-d3)) + \
|
|
((x3 - x1) / (x2 - x1)) * (f1 / d2) * (f2 / d3)
|
|
if not math.isfinite(t):
|
|
t = 0.5
|
|
t = min(max(t, tl), 1.0 - tl)
|
|
except _Budget:
|
|
return (x1 if abs(f1) < abs(f2) else x2), g.n, False
|
|
|
|
|
|
# --------------------------------------------------------------------------
|
|
# Dekker 1969 - secant with a bisection fallback, the direct ancestor of Brent
|
|
# --------------------------------------------------------------------------
|
|
def dekker(f, a, b, xtol, maxfev):
|
|
g = Counter(f, maxfev)
|
|
bb = b
|
|
try:
|
|
fa, fb = g(a), g(b)
|
|
if fa == 0.0:
|
|
return a, g.n, True
|
|
if fb == 0.0:
|
|
return b, g.n, True
|
|
if abs(fa) < abs(fb):
|
|
a, b, fa, fb = b, a, fb, fa
|
|
c, fc = a, fa
|
|
while abs(b - c) > xtol:
|
|
m = 0.5 * (c + b)
|
|
if fb != fa:
|
|
s = b - fb * (b - a) / (fb - fa)
|
|
else:
|
|
s = m
|
|
lo, hi = (b, m) if b < m else (m, b)
|
|
bnew = s if lo < s < hi else m
|
|
a, fa = b, fb
|
|
b = bnew
|
|
bb = b
|
|
fb = g(b)
|
|
if (fb > 0.0) == (fc > 0.0):
|
|
c, fc = a, fa
|
|
if abs(fc) < abs(fb):
|
|
b, c = c, b
|
|
fb, fc = fc, fb
|
|
if fb == 0.0:
|
|
return b, g.n, True
|
|
return b, g.n, True
|
|
except _Budget:
|
|
return bb, g.n, False
|
|
|
|
|
|
# --------------------------------------------------------------------------
|
|
# Library implementations. Where SciPy ships the method, its implementation
|
|
# is the one a reader can check and the one a user would actually run, so it
|
|
# belongs in the comparison. Ours are kept beside them for two reasons: the
|
|
# timing section needs pure Python throughout, and a competitor is entitled to
|
|
# the better of the two counts (see `best_of` below).
|
|
#
|
|
# Measured difference on the 199 problems: brentq 7004 against our 6991,
|
|
# _chandrupatla 5228 against our 5209 - that is, OUR versions are slightly
|
|
# cheaper, so taking the library number alone would flatter us.
|
|
# --------------------------------------------------------------------------
|
|
def brent_scipy(f, a, b, xtol, maxfev):
|
|
from scipy.optimize import brentq
|
|
g = Counter(f, maxfev)
|
|
try:
|
|
x = brentq(g, a, b, xtol=xtol, rtol=8.9e-16, maxiter=maxfev)
|
|
return x, g.n, True
|
|
except _Budget:
|
|
return 0.5 * (a + b), g.n, False
|
|
except Exception:
|
|
return 0.5 * (a + b), g.n, False
|
|
|
|
|
|
def chandrupatla_scipy(f, a, b, xtol, maxfev):
|
|
import numpy as np
|
|
from scipy.optimize._chandrupatla import _chandrupatla
|
|
g = Counter(f, maxfev)
|
|
|
|
def vec(x):
|
|
return np.asarray([g(float(xi)) for xi in np.atleast_1d(x)]).reshape(
|
|
np.shape(x))
|
|
|
|
try:
|
|
res = _chandrupatla(vec, a, b, xatol=xtol, xrtol=8.9e-16,
|
|
fatol=0.0, frtol=0.0, maxiter=maxfev)
|
|
return (float(np.asarray(res.x).ravel()[0]), g.n,
|
|
bool(np.all(np.asarray(res.success))))
|
|
except _Budget:
|
|
return 0.5 * (a + b), g.n, False
|
|
except Exception:
|
|
return 0.5 * (a + b), g.n, False
|
|
|
|
|
|
def best_of(*methods):
|
|
"""Give the competitor the cheaper of its implementations on every
|
|
problem. The same concession the paper already makes to ITP by reporting
|
|
a per-problem oracle over its constants."""
|
|
def run(f, a, b, xtol, maxfev):
|
|
best = None
|
|
for m in methods:
|
|
try:
|
|
x, n, conv = m(f, a, b, xtol, maxfev)
|
|
except Exception:
|
|
continue
|
|
if conv and (best is None or n < best[1]):
|
|
best = (x, n, True)
|
|
return best if best is not None else (0.5 * (a + b), maxfev, False)
|
|
return run
|
|
|
|
|
|
BASELINES = {
|
|
"chandrupatla": chandrupatla,
|
|
"chandrupatla_scipy": chandrupatla_scipy,
|
|
"brent_scipy": brent_scipy,
|
|
"dekker": dekker,
|
|
"bisection": bisection,
|
|
"regula_falsi": regula_falsi,
|
|
"illinois": illinois,
|
|
"pegasus": pegasus,
|
|
"anderson_bjorck": anderson_bjorck,
|
|
"ridders": ridders,
|
|
"brent": brent,
|
|
"toms748": toms748,
|
|
"itp": itp,
|
|
}
|