2

Objective
I need to find roots of $$f(x)=c$$ in interval $[a,b]$, where

  • $f(a)=0$ and $c<f(b)<1$
  • $f(x)$ is unknown outside of the interval
  • $f(x)$ is monotonuous and has decreasing derivative (i.e., $f''(x)<0$)

(Overall it resembles flipped hockey stick.)

What I tried
Right now I am using bisection method, which is somewhat slow for my purposes, as evaluating the function is expensive.

I tried other bracketing methods:

  • false position (regula falsi) is slower due to the stagnant bound. Illinois improves performance a bit, but still slower than bisection.
  • Ridder's method converges in less steps, but overall slower, as it requires more function evaluations (perhaps also due to the stagnant bound).

Question
Could I use my knowledge of the function properties outline in the beginning, to achieve faster convergence (less function evaluations)?

Possibilities. As an idea, one could try using transformation: $$F(x) = \log (1-f(x)), C=\log (1-c),$$ In a hope that false position converges faster for equation $$F(x)=C.$$

Supporting information: more details about the curve can be found here. Mine is not exactly teh same, but very similar. Calculating it directly is expensive, and one resorts to resampling and averaging.

Roger V.
  • 763
  • Do you have an upper bound on $f''$? – Lazy Sep 15 '21 at 15:13
  • Also is $c$ fixed? Is $f$ fixed? – Lazy Sep 15 '21 at 15:19
  • @Lazy Yes, I could calculate $f''$ at the endpoints. – Roger V. Sep 15 '21 at 15:20
  • @Lazy $c$ is a known number, $f(x)$ is a known function (there is an additional level of complexity there, but it would take us away from the question). – Roger V. Sep 15 '21 at 15:22
  • 1
    Did you use a proper implementation of the Illinois method like in https://math.stackexchange.com/questions/3551775/bracketing-root-finding-methods? The en-wiki variant is somewhat short on details, the de-wiki article is better in this regard. – Lutz Lehmann Sep 15 '21 at 15:34
  • @LutzLehmann I based my implementation on the example given here, but thanks for the link - I will check, if something is missing. – Roger V. Sep 15 '21 at 15:44
  • 1
    That link is not more informative than the en-wiki. The key point that is missing in both is the idea to continue the halving of the stalling function value until that boundary moves. This is especially important in the initial "global search" phase with hockeystick-like functions. – Lutz Lehmann Sep 15 '21 at 15:52
  • @RogerVadim I mean this like: Do you want to solve this for multiple points or for multiple functions? Are you generally able to calculate f''? – Lazy Sep 15 '21 at 15:55
  • @Lazy I can't calculate the derivatives. If you want more details: in principle, I know the analytical form of the function, but it is unwieldy, so instead I calculate its values via resampling. – Roger V. Sep 15 '21 at 16:26
  • 1
    If you have a function table, then you can also use inverse interpolation with some spline interpolation function. – Lutz Lehmann Sep 15 '21 at 17:01
  • @RogerVadim Ok. My suggestion would be: If you can somewhat estimate $f''$ on an interval $[c,d]\subset[a,b]$ then you could approximate $f$ by a second degree polynom $a_0+a_1x+a_2x^2$ so that $2a_2\approx f''$ on $[c,d]$ and the functions are identical in the edges. Now we can solve $a_0+a_1x+a_2x^2=c$ and use the solution $x$ as new bisection point. This might speed up the bisection process quite a bit, depending on how behaved your $f$ is.

    If $f''$ behaves it might suffice to calculate $f''$ for some points in $[a,b]$.

     – Lazy Sep 15 '21 at 19:56
  • We might for example use this to invert the sine on $[0,\pi/2]$: By evaluating $f''$ on a few points I can solve $\sin(x)=9/10$ with a precision of $0.0000001$ in $10$ iterations, while regular bisection requires $22$. The problem here is: How do we get a good idea about $f''$ efficiently? If we get a sufficiently? – Lazy Sep 15 '21 at 20:08
  • If we can get an upper and a lower bound on $f''$, then we can use this to simultaneously calculate two bounds at the same time, one lower and one upper. Because if $f'' > C$ then if we approximate $f$ by assuming $f''=C$ then we wil get a solution that is too large, while if $f'' < C$ we get a solution that is too small. – Lazy Sep 15 '21 at 20:08
  • Sorry, the $<,>$ should be the other way round. – Lazy Sep 15 '21 at 20:22

2 Answers2

1

I shall suppose that the computation of $f'(x)$ is more expensive that the computation of $f(x)$ itself.

For this kind of problems, I extensively used the so-called RTSAFE subroutine from "Numerical Recipes" (have a look at the documentation here on chapter $9$. Just adapt the FUND subroutine to return the derivative by finite difference.

This corresponds to a combination of bisection and Newton steps.

Claude Leibovici
  • 289,558
0

The largest bottleneck is the function evaluations, but luckily for us the root can be bracketed with something like bisection, ensuring we can avoid interpolation iterations if we know they will likely be unhelpful.

This is exactly where Chandrupatla's method shines. You can check here for an example implementation.

Output of the above example:

\begin{array}{r|r} x & f(x)=x^2-2 \\ \hline 0.0000000000000000 & -2.0000000000000000 \\ 2.0000000000000000 & 2.0000000000000000 \\ 1.0000000000000000 & -1.0000000000000000 \\ 1.5000000000000000 & 0.2500000000000000 \\ 1.4095238095238036 & -0.0132426303855042 \\ 1.4142641817017454 & 0.0001431756445074 \\ 1.4142135574926071 & -0.0000000138041043 \\ 1.4142135623731014 & 0.0000000000000178 \\ 1.4142135623730892 & -0.0000000000000167 \\ 1.4142135623730951 & 0.0000000000000004 \end{array}

It can be seen that Chandrupatla's method intelligently chooses to use bisection during the initial iterations to ensure convergence while it uses interpolation for later iterations where it is confident it should be accurate.

The downside is that in some cases it does not alternate sides of the root each iteration. This leads to a slower order of convergence. When alternating, the order of convergence is $\approx1.839$. When not alternating, the order of convergence is $\approx1.618$. This means the number of digits accurate increases by roughly $83.9\%$ and $61.8\%$ respectively every iteration.

If your function is stochastic i.e. you are estimating $f(x)$ as an average of $f_i(x)$ for many values of $i$, then you may want to resort to stochastic root-finding algorithms, at least for the initial iterations.

You can check here for an example implementation of the Robbins-Monro algorithm with Polyak-Ruppert averaging.

Output of the above example:

\begin{array}{r|r} x & f(x)=x^2-2~~(\pm0.1) \\ \hline 1.0000000000000000 & -1.0000000000000000 \\ 1.3208425751018684 & -0.2553748917982650 \\ 1.2931735764314634 & -0.3277021012194583 \\ 1.3337128053831124 & -0.2212101527571080 \\ 1.3501779010773112 & -0.1770196354424665 \\ 1.3656025721082210 & -0.1351296150514110 \\ 1.3758861134900708 & -0.1069374027051879 \\ 1.3847629867055029 & -0.0824314706504552 \\ 1.3880601300384776 & -0.0732890753975646 \\ 1.3900420572765071 & -0.0677830790024958 \\ 1.3914381380169945 & -0.0638999080717995 \\ 1.3909224968549345 & -0.0653346077428347 \\ 1.3916674836639233 & -0.0632616149125236 \\ 1.3918200228103006 & -0.0628370241043343 \\ 1.3941212662850380 & -0.0564258948918022 \\ 1.3964265279275152 & -0.0499929521003046 \\ 1.3987379168844518 & -0.0435322398697444 \\ 1.4008372481138009 & -0.0376550042969532 \\ 1.4009627374051223 & -0.0373034084023462 \\ 1.4012682932051819 & -0.0364471704578364 \\ 1.4010721993805688 & -0.0369966921228957 \\ 1.4015011630538421 & -0.0357944899587279 \end{array}

Once you feel like you have converged well enough for the stochastic approximation, you can use that estimate to kick-start your bracketing method to reach high-precision results in only a few additional iterations (at the cost of needing to spend more computations on those iterations).

Simply Beautiful Art
  • 76,603