ITP method

[tecosaur]

From Zulip in early 2022, I developed this ITP implementation that could be turned into a PR.

"""
    ITP(; n₀::Int, κ₁::Real, κ₂::Real)

Use the [ITP method](https://en.wikipedia.org/wiki/ITP_method) to find
a root of a bracketed function, with a convergence rate between 1 and 1.62.

This method was introduced in the paper "An Enhancement of the Bisection Method
Average Performance Preserving Minmax Optimality"
(https://doi.org/10.1145/3423597) by I. F. D. Oliveira and R. H. C. Takahashi.

# Tuning Parameters

The following keyword parameters are accepted.

- `n₀::Int = 1`, the 'slack'. Must not be negative.\n
  When n₀ = 0 the worst-case is identical to that of bisection,
  but increacing n₀ provides greater oppotunity for superlinearity.
- `κ₁::Float64 = 0.1`. Must not be negative.\n
  The recomended value is `0.2/(x₂ - x₁)`.
  Lower values produce tighter asymptotic behaviour, while higher values
  improve the steady-state behaviour when truncation is not helpful.
- `κ₂::Real = 2`. Must lie in [1, 1+ϕ ≈ 2.62).\n
  Higher values allow for a greater convergence rate,
  but also make the method more succeptable to worst-case performance.
  In practice, κ=1,2 seems to work well due to the computational simplicity,
  as κ₂ is used as an exponent in the method.

### Worst Case Performance

n½ + `n₀` iterations, where n½ is the number of iterations using bisection
(n½ = ⌈log2(Δx)/2`tol`⌉).

### Asymptotic Performance

If `f` is twice differentiable and the root is simple,
then with `n₀` > 0 the convergence rate is √`κ₂`.
"""
struct ITP{K1,K2} <: AbstractBracketingAlgorithm where {K1, K2 <: Real}
    n₀::Int
    κ₁::K1
    κ₂::K2
    function ITP(n₀::Int=1, κ₁::K1=0.2, κ₂::K2=2) where {K1, K2 <: Real}
        if κ₁ < 0
            throw(ArgumentError("κ₁ cannot be netagive ($κ₁ ≱ 0)"))
        end
        if !(1 <= κ₂ <= (3+√5)/2)
            throw(ArgumentError(
                string("κ₂ must lie in [1, 1+ϕ) ($κ₂ ",
                       if κ₂ < 1; "≱ 1)" else "≮ 1+ϕ ≈ 2.618)" end)))
        end
        if n₀ < 0
            throw(ArgumentError("n₀ cannot be negative ($n₀ ≱ 0)"))
        end
        ITP{float(K1), K2}(n₀, float(κ₁), κ₂)
    end
end

ITP(; n₀::Int=1, κ₁::K1=0.2, κ₂::K2=2) where {K1, K2 <: Real} =
    ITP{K1,K2}(n₀, κ₁, κ₂)

struct ITPCache{F <: AbstractFloat}
    x̂::F
    ϵₛ::F
end

function alg_cache(::ITP, x₁, x₂, _, _)
    nₘₐₓ = ceil(log2((x₂ - x₁)/2ϵ))
    ϵₛ = ϵ * 2^(nₘₐₓ + n₀)
    x̂ = (x₂ + x₁)/2
    ITPCache(x₁, x₂, x̂, ϵₛ)
end

function perform_step(solver::BracketingImmutableSolver, alg::ITP, _)
    @unpack f, left, right, fl, fr, cache = solver
    x₁, x₂, y₁, y₂ = left, right, fl, fr

    @unpack x̂, ϵₛ = solver.cache

    r = ϵₛ - (x₂ - x₁)/2 # Minimax radius
    δ = κ₁ * (x₂ - x₁)^κ₂ # Truncation error
    # Interpolation, Regula Falsi
    xᵢ = (y₂ * x₁ - y₁ * x₂) / (y₂ - y₁)
    σ = sign(x̂ - xᵢ)
    # Truncation, perturbing xᵢ towards x̂
    xₜ = if δ <= abs(x̂ - xᵢ); xᵢ + σ*δ else x̂ end
    # Projection of xₜ onto the minimax interval [x̂ - r, x̂ + r]
    xₚ = if abs(xₜ - x̂) <= r; xₜ else x̂ - σ*r end
    yₚ = f(xₚ)
    if yₚ > 0
        x₂, y₂ = xₚ, yₚ
    else
        x₁, y₁ = xₚ, yₚ
    end
    x̂ = (x₂ + x₁)/2
    ϵₛ /= 2

    @set! solver.cache.x̂ = x̂
    @set! solver.cache.ϵₛ = ϵₛ

    @set! solver.left = x₁
    @set! solver.right = x₂
    @set! solver.fl = y₁
    @set! solver.fr = y₂

    solver
end

[tecosaur]

Hmmm, it looks like the package has changed a fair bit since this was written. Oh well, it's still probably better to have this issue than not.

[ChrisRackauckas]

This would fit better into SimpleNonlinearSolve.jl. @yash2798 this might be a good one to handle early.

[ChrisRackauckas]

Thanks! It's now in SimpleNonlinearSolve, heavily based on your code.