Roots.jl
Roots.jl copied to clipboard
Repetitive solving
I have a non-increasing function $f$ on $\mathbb R_+$ that starts at $f(0) =1$ and ends up at $f(Inf) = 0$. I want to compute the inverse of the function, and I am currently using the folowing code:
f_inv(y) = find_zero(x -> f(x) - y,(0,Inf))
Then i am repetitively calling this inversion as, e.g, :
samples = rand(1e6)
f_inv.(samples)
Would it be possible to reduce the cost of subsequent evaluations ? I guess that, since this is using a bisection method, we could re-start each iteration using the previous bisection. Is there something like that implemented in Roots.jl ?
is your function continuous? it is differentiable?, if a) is true, maybe you could benefit from a secant method, if it is 1) and 2), you could use https://juliamath.github.io/Roots.jl/stable/reference/#Roots.LithBoonkkampIJzermanBracket
Yes, either of those should work, though a bracketing method does offer guarantees.
For a bracketing method, you might gain a few function evaluations by using the interval (f(x_last), oo) as the function is increasing, so if you are looking at finding the function at an increasing selection of points you need not start at 0.
I guess, this is most efficient if you were to find the inverse at select points and then interpolate between those to fill out the description of the curve.
No it is not guaranteed to be continuous, nor differentiable. But it is guaranteed to be convex and non-increasing in the loose sense, that is, $\forall x, y$ such that $y > x$, $\forall \alpha \in [0,1]$,
- $f(y) \le f(x)$
- $\alpha * f(y) + (1-\alpha) * f(x) \le f(\alpha t + (1 - \alpha) x)$.
Those properties are strong, but do not ensure continuity : there might be a point $x$ so that $f(x-) \neq f(x+)$. in which case, I want f_inv to systematically return the same side (which one I need I do not remember RN).
Due to this non-continuity, I thought that the bisection was my best bet, but maybe i'm wrong ?
I do want to store the previously-found inverses so that the bisection can restart from a smaller interval than (0,Inf)
each time, which was the core of my question