Bisection algorithm for concave function

Question

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.

Answer 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.

Answer 2

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).