I have been unable to find any real code implementation of ITP method - Wikipedia. However, it claims to be the bee’s knees - just read the first paragraph of the Wikipedia article:
“In numerical analysis, the ITP method , short for Interpolate Truncate and Project , is the first root-finding algorithm that achieves the superlinear convergence of the secant method[1] while retaining the optimal[2] worst-case performance of the bisection method.[3] It is also the first method with guaranteed average performance strictly better than the bisection method under any continuous distribution.[3] In practice it performs better than traditional interpolation and hybrid based strategies (Brent’s Method, Ridders, Illinois), since it not only converges super-linearly over well behaved functions but also guarantees fast performance under ill-behaved functions where interpolations fail.[3]”
I implemented the pseudocode in the Wikipedia article, and thought I should share it. If anyone wants to put this in a package like Roots.jl, NLSolve.jl or NonlinearSolve.jl, feel free to - I’d be thrilled.
EDIT: provided defaults, and changed argument order to make having defaults more sensible.
EDIT 2: Put the function argument first, diverging further from the Wikipedia pseudocode, but converging towards the conventional (and sensible) ordering. Also added reference to the source of hyperparameter defaults.
"""
ITP(f, a, b, ϵ=eps(), κ₁=0.2/(b-a), κ₂=2, n₀=1)
Find the root of the funtion `f` on the interval [a, b], given a
tolerance ϵ and hyperparameters κ₁, κ₂, n₀. The hyperparameters are
subject to the following constraints:
0 < κ₁ < Inf
1 ≤ κ₂ < 1 + (1 + √5)/2
0 < n₀ < Inf
A Julia implementation of the Interpolate, Truncate Project algorithm,
as detailed in `https://en.wikipedia.org/wiki/ITP_method`.
Hyperparameter defaults are set on the basis of the reasoning provided
in `https://docs.rs/kurbo/0.8.1/kurbo/common/fn.solve_itp.html`
"""
function ITP(f, a, b, ϵ=eps(), κ₁=0.2/(b-a), κ₂=2, n₀=1)
0 < κ₁ < Inf || error("κ₁ must be between 0 and ∞")
1 ≤ κ₂ < 1 + (1 + √5)/2 || error("κ₂ must be between 1 and 1+ (1 + √5)/2 (1 + the golden ratio)")
0 < n₀ < Inf || error("n₀ must be between 0 and ∞")
n_1div2 = ceil(Int, log2((b-a)/2ϵ))
nₘₐₓ = n_1div2 + n₀
y_a = f(a)
y_b = f(b)
sign(y_a) == sign(y_b) && error("sign(f(a)) = sign(f(b)). There is no guaranteed root in the given interval.")
j = 0
while b-a > 2ϵ
# Calculating parameters:
x_1div2 = (a+b)/2
r = ϵ*2^(nₘₐₓ - j) - (b-a)/2
δ = κ₁*(b-a)^κ₂
# Interpolation:
x_f = /(y_b*a - y_a*b, y_b-y_a)
# Truncation:
σ = sign(x_1div2 - x_f)
δ ≤ abs(x_1div2 - x_f) ? (x_t=x_f+σ*δ) : (x_t = x_1div2)
# Projection:
abs(x_t - x_1div2) ≤ r ? (x_ITP = x_t) : (x_ITP = x_1div2 - σ*r)
# Updating Interval:
y_ITP = f(x_ITP)
if y_ITP > 0
b = x_ITP
y_b = y_ITP
elseif y_ITP < 0
a = x_ITP
y_a = y_ITP
else
a = b = x_ITP
end
j += 1
end
return (a+b)/2
end
One could concider adding the following constraint:
ϵ ≤ big(1e-78) && error("A tolerance this small has lead to indefinite iterations in the past. Aborting.")
I don’t have anything but the code in the OP, so it should be easy to test it for anyone who wants to.