Roots.jl

Roots.jl copied to clipboard

It would be good to have a clear idea of the intended use cases for this package.

The solvers are great and fast for large problems run from the REPL. But running find_zero() a million times a second doesn't seem to be a good idea at the moment. Unfortunately that's the most common use case in my field. The setup costs can be larger than the total runtime cost, and can't always be cached, like when running inside a solver that is swapping between Dual numbers and normal numbers, for example.

A tiny, basic secant or zbrent solver with normal args is much faster for these situations, and I keep ending up back there - with my own custom roots package to maintain! It may be that such a separate SimpleRoots package is actually a good idea to divide up this space, and would have a small maintenance footprint anyway.

Any thoughts on that?

Thanks for the question. In working on #116 I added a speedier alternative to the Order0 method, bypassing the machinery of the latter. I can see the merit of some speedier alternatives. This could consist, at first, of bisection64, a42 (a Brent alternative), this Order0 alternative, and a secant method. This only requires a secant implementation, which I'll add to that PR.

Having a separate package seems like it will cause more confusion. Once the ForwardDiffdependency is dropped (I was planning on doing that after merging #116) this package should be fairly lightweight.

Great to hear your'e working on that! Sorry for the constant feature requests ;)

I found that normal args are faster than keyword args, so having a wrapper with keywords and an internal function where all atol etc are normal args is a good solution.

The setup cost on this secant is basically non-existent, maybe a good target:

`
find_zero(f, bracket, ::Secant; atol=1e-10, max_iter=100) = find_zero(f, bracket, atol, max_iter)
find_zero(f, bracket, ::Secant, atol, max_iter) = begin
    local x
    x0, x1 = bracket
    y0, y1 = f(x0), f(x1) 

    for _ in 1:max_iter
        x = x1 - y1 * (x1 - x0) / (y1 - y0)
        if abs(x - x1) < atol return x end
        x0 = x1
        y0 = y1
        x1 = x
        y1 = f(x1)
    end
    error("Root not found")
end

I found that normal args are faster than keyword args

Worth to test on 0.7 where the implementation of keyword args have been overhauled.

Anyone did that testing at any point?