Roots.jl icon indicating copy to clipboard operation
Roots.jl copied to clipboard

Repetitive solving

Open lrnv opened this issue 1 year ago • 3 comments

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 ?

lrnv avatar Nov 09 '23 12:11 lrnv

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

longemen3000 avatar Nov 09 '23 22:11 longemen3000

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.

jverzani avatar Nov 09 '23 22:11 jverzani

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

lrnv avatar Nov 09 '23 22:11 lrnv