Make `find_zero` more robust

General Usage
,
1

In my PowerAnalyses.jl package, I’m having trouble in finding a robust way to solve for zero. It works, but it’s not very robust. For example: ArgumentError when effect size is small (<=0.10) · Issue #21 · rikhuijzer/PowerAnalyses.jl · GitHub.

This is the function (from Roots.jl) that crashes (source):

function get_n(T::StatisticalTest; alpha::Real, power::Real, es::Real)
    f(n) = get_alpha(T; es, power, n) - alpha
    initial_value = (2, 1000)
    return find_zero(f, initial_value)
end

Is there anyone here who has tips on making find_zero more robust?

2

Use NonlinearSolve.jl instead. ITP will be a more stable algorithm (and you should get a nice performance boost).

1 Like
3

This is weird because you suggest switching interface with an argument about a better algorithm. Roots.jl has ITP : Reference/API · Roots, so it’ll be IMHO largely easier to simply switch algorithm inside Roots.jl

If the difference is so impressive, the fruit to make Roots.jl faster is probably low-hanging and we should open an issue on Roots.jl

4

I think the issue is the fixed choice of bracketing interval in the package. With 'Bisection(), these can be infinite, but I think the default is an Alefeld, Potra, Shi algorithm, which is much more efficient, but can have issues if the function evaluates to Inf. For robustness, you could use (0,Inf) with Bisection() at the cost of more function calls. If you have a reasonable starting guess, the default will be a hybrid of a secant method and a bracketing method. That might be more robust as well, but I’m not sure how you would get a working initial guess over all the cases considered in the package.

The ITP algorithm is also available in Roots. I wouldn’t expect much more, but you could experiment.

As for performance, Roots is a bit less performant than NonlinearSolve, but not markedly so. Any thoughts on closing that gap would be most welcome.

2 Likes