Bracketed Halley?

[longemen3000]

Hi, i'm in need of a "bracketed halley" method. basically something like:

fx,dfx,d2fx2 = f(x)
dβ = - (2*fx*dfx)/(2*dfx^2-fx*d2fx2)

# restricted β space
if fx < 0.
  βmax = β
elseif fx > 0.
  βmin = β
end

#updating β
βnew =  β + dβ
if βmin < βnew && βnew < βmax
  β = βnew
else
  dβ = (βmin + βmax) / 2 - β
  β = dβ + β
end

is any method available in Roots.jl to do this?

[longemen3000]

of course, you could generalize and say that any method that does not use brackets could be generalized into a bracketed method, something like Bracketed(inner::AbstractNonBracketingMethod) that selects or rejects an step based on the results of inner

[jverzani]

Does this do what you need: https://github.com/JuliaMath/Roots.jl/blob/master/src/Derivative/lith.jl#L291

I would think that algorithm could be tweaked to include the second derivative, if needed.

[longemen3000]

That seems like the generalization of what I need hahah. It would be great if it available

[jverzani]

This is Chebyshev, not Halley, but if you really need Halley it could be adjusted in a few lines. I'm not sure it is worth including, but perhaps if you like this you could convince me otherwise.

# Order 3 bracketing method
#
using Roots
import Roots: @set!, bracket, adjust_bracket, NullTracks, isissue, log_message
struct ChebyshevBracket <: Roots.AbstractBracketingMethod end
Roots.fn_argout(::ChebyshevBracket) = 3
struct ChebyshevBracketState{T,S} <: Roots.AbstractUnivariateZeroState{T,S}
    xn1::T
    xn0::T
    fxn1::S
    fxn0::S
end

# use xatol, xrtol only, but give some breathing room over the strict ones and cap number of steps
function Roots.default_tolerances(::ChebyshevBracket, ::Type{T}, ::Type{S}) where {T,S}
    xatol = eps(real(T)) * oneunit(real(T))
    xrtol = 2eps(real(T))  # unitless
    atol = 4 * eps(real(float(S))) * oneunit(real(S))
    rtol = 4 * eps(real(float(S))) * one(real(S))
    maxevals = 40
    strict = false
    (xatol, xrtol, atol, rtol, maxevals, strict)
end


function Roots.init_state(M::ChebyshevBracket, F::Roots.Callable_Function, x)
    x₀, x₁ = adjust_bracket(x)
    fx₀, fx₁ = first(F(x₀)), first(F(x₁))
    state = Roots.init_state(M, F, x₀, x₁, fx₀, fx₁)
end

function Roots.init_state(L::ChebyshevBracket, F, x₀, x₁, fx₀, fx₁)
    # keep xn1 the smallest of fxn1, fxn0
    if abs(fx₀) > abs(fx₁)
        state = ChebyshevBracketState(x₁, x₀, fx₁, fx₀)
    else
        state = ChebyshevBracketState(x₀, x₁, fx₀, fx₁)
    end
    state
end

# bracketed Chebyshev
function Roots.update_state(
    M::ChebyshevBracket,
    F,
    o::ChebyshevBracketState{T,S},
    options,
    l=NullTracks(),
) where {T,S}

    c,b,fc,fb = o.xn1, o.xn0, o.fxn1, o.fxn0
    if c < b
        a, fa = c, fc
    else
        a,b,fa,fc = b,c,fb,fc
    end

    fc::S, (Δ, Δ₂) = F(c)

    if isissue(Δ)
        log_message(l, "Issue with computing `Δ`")
        return (o, true)
    end

    # try chebyshev
    # if issue, then newton
    # if issue, then secant
    # if issue, then bisection
    x = c - Δ - 1/2 * Δ^2/Δ₂
    if !(a < x < b)
        x = c - Δ
    end

    if a < x < b
        fx,_ = F(x)
        ā,b̄,c,fā,fb̄,fc = bracket(a,b,x,fa,fb,fx)
    else
        # secant step
        x = a - fa * (b-a)/(fb-fa)
        if x - a < 0.0001 * (b-a) || b-x < 0.0001 * (b-a)
            x = a + (b-a)*0.5
        end

        fx,_ = F(x)
        ā,b̄,c,fā,fb̄,fc = bracket(a,b,x,fa,fb,fx)
    end
    c, fc = abs(fā) < abs(fb̄) ? (ā,fā) : (b̄, fb̄)

    a,b,fa,fb = ā, b̄, fā, fb̄

    @set! o.xn0 = a == c ? b : a
    @set! o.xn1 = c
    @set! o.fxn0 = a == c ? fb : fa
    @set! o.fxn1 = fc

    return o, false
end

[longemen3000]

Seems excelent!, to be fair, the bracketed halley is used because a) we know the second derivatives of the problem, b) the solution has a clear bracket. Apart from that, there isn't any real preference for the second order method other than "Halley uses second derivatives". I'm gonna try that code posted here.

[longemen3000]

the method seems to pass the tests given by the package. in what concerns this issue, the code posted above would fill a hole (niche hole, but a hole either way), but i can use that code directly if this does not appear in a release

[jverzani]

I think I can do something to add this in. I'll keep the issue open for awhile.