`Muller` algorithm benchmark

[fgittins]

Motivated by #501, I was curious to see how the newly added Muller algorithm (#568) fares in the SciML benchmarks, particularly against the implementation in Roots.jl.

Below are performance results on my local machine (MacBook Pro M3, 16 GB RAM), using BenchmarkTools. Each test solves the same root-finding problem N = 100_000 times (see code below):

Method	Time	Allocations
Roots.jl default	155.340 ms	100,000
ITP	22.850 ms	0
Roots.muller	26.713 ms	0
Muller	11.839 ms	0

Highlights:

• I was able to recreate the expected speed-up of ITP vs. the default solver in Roots.jl, consistent with the SciML benchmarks.
• I found that Roots.muller is significantly faster than the default method (likely because it is non-allocating), but both NonlinearSolve.jl algorithms outperform it.
• Finally, it was interesting that Muller is almost two times faster than ITP for this function. Having run this several times on my own machine, this factor-of-two difference seems consistent. However, this may only reflect the suitability of Muller for this specific problem rather than general performance.

Benchmark script

This is the script that I ran to generate the benchmarks:

using BenchmarkTools
using NonlinearSolve
using Roots

const N = 100_000

myfun(u, level) = u * sin(u) - level

function froots!(out, levels, uspan)
    for i in 1:N
        out[i] = find_zero(x->myfun(x, levels[i]), uspan)
    end
    out
end

function f!(out, levels, uspan)
    for i in 1:N
        prob = IntervalNonlinearProblem{false}(IntervalNonlinearFunction{false}(myfun), uspan, levels[i])
        sol = solve(prob, ITP())
        out[i] = sol.u
    end
    out
end

function groots!(out, levels, uspan)
    xᵢ₋₂, xᵢ = uspan
    xᵢ₋₁ = (xᵢ₋₂ + xᵢ) / 2  # choose midpoint for consistency with Muller()
    for i in 1:N
        out[i] = Roots.muller(x->myfun(x, levels[i]), xᵢ₋₂, xᵢ₋₁, xᵢ)
    end
    out
end

function g!(out, levels, uspan)
    for i in 1:N
        prob = IntervalNonlinearProblem{false}(IntervalNonlinearFunction{false}(myfun), uspan, levels[i])
        sol = solve(prob, Muller())
        out[i] = sol.u
    end
    out
end

function main()
    out = zeros(N)
    levels = 1.5 .* rand(N)
    uspan = (0.0, 2.0)

    @btime froots!($out, $levels, $uspan)
    @btime f!($out, $levels, $uspan)
    @btime groots!($out, $levels, $uspan)
    @btime g!($out, $levels, $uspan)

    nothing
end

main()

If this speed-up is robust, it may be worth adding to the benchmarks.

[ChrisRackauckas]

It would be good to turn that into a more comprehensive benchmark of the bracketing solvers, loop through them and make a table of results.

As for defaulting and such, ITP is generally chosen as a default in things like ODE solver benchmarking because it gets decent speed but it's also using a bisection part so it's very robust to getting floating point accurate. Things like Muller and Newton which don't have the bisection part don't tend to get down to the exact floating point number in ill-conditioned cases. So ITP takes a bit of a hit because of this while balancing having some faster convergence steps. That's why it's the default. But I would expect Muller, Newton, etc. to be faster in almost any case where they are sufficiently robust. It would be nice to setup some interval rootfinding benchmarks that really stress this but alas most of the good cases I know of are buried in ODE solver continuous callbacks.

[fgittins]

It would be good to turn that into a more comprehensive benchmark of the bracketing solvers, loop through them and make a table of results.

Agreed. This would be nice to see.

As for defaulting and such, ITP is generally chosen as a default in things like ODE solver benchmarking because it gets decent speed but it's also using a bisection part so it's very robust to getting floating point accurate.

This makes sense to me. Inevitably there is a balance to be made between robustness and speed. For the ODE solving, it makes sense to err on the side of caution.

[oscardssmith]

specifically a good benchmark function is something like exp(x)-1e-15 since the curvature there is very unkind to lots of methods

[fgittins]

I ran the benchmark for the other BracketingNonlinearSolve.jl algorithms and these are the results:

Method	Time for u sin(u) - level	Time for exp(u) - 1e-15
Alefeld	72.875 ms	127.246 ms
Bisection	70.463 ms	87.490 ms
Brent	36.273 ms	41.182 ms
Falsi	88.493 ms	80.497 ms
ITP	23.465 ms	66.173 ms
Muller	12.215 ms	2.364 ms
Ridder	39.761 ms	39.983 ms

For my tests, Muller was the fastest. It's notable that Muller is an order of magnitude faster than the closest algorithm for the exponential function. This might be explained by the fact that the algorithm approximates the function as quadratic, but I haven't delved into this.

However, as is expected, the performance of many of the algorithms is very sensitive to uspan. When I changed it from (-100.0, 0.0) to (-100.0, 10.0), the results became:

Method	Time for exp(u) - 1e-15
Alefeld	1.635 ms
Bisection	87.534 ms
Brent	63.610 ms
Falsi	81.113 ms
ITP	61.463 ms
Muller	2.364 ms
Ridder	28.026 ms

where Alefeld dramatically improves and takes the lead. But it's notable that Muller stays consistent.

Benchmark script

using BenchmarkTools
using BracketingNonlinearSolve

myfun(u, level) = u * sin(u) - level
myotherfun(u, level) = exp(u) - 1e-15

function f!(out, levels, uspan, alg)
    N = length(out)
    for i in 1:N
        prob = IntervalNonlinearProblem{false}(IntervalNonlinearFunction{false}(myfun), uspan, levels[i])
        sol = solve(prob, alg)
        out[i] = sol.u
    end
    out
end

function g!(out, uspan, alg)
    N = length(out)
    for i in 1:N
        prob = IntervalNonlinearProblem{false}(IntervalNonlinearFunction{false}(myotherfun), uspan)
        sol = solve(prob, alg)
        out[i] = sol.u
    end
    out
end

function main()
    N = 100_000
    out = zeros(N)
    levels = 1.5 .* rand(N)
    uspan = (0.0, 2.0)

    for alg in (Alefeld(), Bisection(), Brent(), Falsi(),
                ITP(), Muller(), Ridder())
        println("Benchmark of $alg for myfun:")
        @btime f!($out, $levels, $uspan, $alg)
    end

    uspan = (-100.0, 0.0)

    for alg in (Alefeld(), Bisection(), Brent(), Falsi(),
                ITP(), Muller(), Ridder())
        println("Benchmark of $alg for myotherfun:")
        @btime g!($out, $uspan, $alg)
    end

    nothing
end

main()

[ChrisRackauckas]

Should also show the final error in the solution too.

[fgittins]

Good suggestion. I checked this for the exponential and Muller is missing the root. This probably isn't too surprising, since it isn't a "bracketing" method.

Method	Error
Alefeld	1.1832913578315177e-30
Bisection	-3.4512664603419266e-29
Brent	1.1832913578315177e-30
Falsi	0.9999999999998214
ITP	1.1832913578315177e-30
Muller	-9.999998071250153e-16
Ridder	1.1832913578315177e-30

[ChrisRackauckas]

This would be great in the SciMLBenchmarks