Root finding via local measurement

TL;DR

This paper introduces two novel numerical methods for root finding that rely solely on derivative measurements at a point, avoiding direct function evaluations, and analyzes their convergence properties.

Contribution

It proposes and characterizes two new derivative-only root-finding algorithms, expanding the toolkit for scenarios with limited function access.

Findings

Both methods converge algebraically with respect to N.

Convergence rate improves with the number of derivatives used.

Methods are applicable under certain conditions for derivative measurements.

Abstract

We consider the problem of numerically identifying roots of a target function - under the constraint that we can only measure the derivatives of the function at a given point, not the function itself. We describe and characterize two methods for doing this: (1) a local-inversion "inching process", where we use local measurements to repeatedly identify approximately how far we need to move to drop the target function by the initial value over N, an input parameter, and (2) an approximate Newton's method, where we estimate the current function value at a given iteration via estimation of the integral of the function's derivative, using N samples. When applicable, both methods converge algebraically with N, with the power of convergence increasing with the number of derivatives applied in the analysis.

Equations

$$dy\equiv \sum _{j=1}^m{a}_{j}d{x}^j+O\left(d{x}^{m+1}\right)$$

$$dx=\sum _{j=1}^m{A}_{j}d{y}^j+O\left(d{y}^{m+1}\right)$$

$$A_n=\frac{1}{n{a}_1^n}\sum _{s_1,s_2,s_3,\cdots }(-1{)}^{s_1+s_2+s_3+\cdots }\frac{n(n+1)\cdots (n-1+s_1+s_2+s_3+\cdots )}{s_1!s_2!s_3!\cdots }{\left(\frac{a_2}{a_1}\right)}^{s_1}{\left(\frac{a_3}{a_1}\right)}^{s_2}\cdots .$$

$$s_1+2s_2+3s_3+\cdots =n-1.$$

$$A_1$$

$$A_2$$

$$A_3$$

$$x_k=x_{k-1}+d{x}_{k-1},\text{ for k∈{1,…,N}}$$

$$r=\frac{1}{4}\frac{(sR{)}^2}{M}.$$

$$d{x}_k=\sum _{j=1}^m{A}_j{\left(-\frac{y_0}{N}\right)}^j+O\left(N^{-(m+1)}\right)$$

$$x_N\sim x_{\ast }+O(N^{-m}).$$

$$x_{\ast }=x_0-{\int }_0^{y_0}\frac{dy}{y^{\prime }(x)},$$

$$dx\sim \frac{-a_1+\sqrt{a_1^2+4a_{2}dy}}{2a_2}+O(d{x}^3).$$

$$x_{k+1}=x_k-\frac{y(x_k)}{y^{\prime }(x_k)}$$

$$y(x_k)=y(x_0)+{\int }_{x_0}^{x_k}{y}^{\prime }(x)dx$$

$${\tilde{x}}_i=x_0+\frac{i}{N-1}(x_k-x_0),\text{ i∈0,…,N−1}.$$

$$x_N\sim x_{\ast }+O(N^{-2\lfloor \frac{m}{2}\rfloor -2})$$

$${\int }_{{\tilde{x}}_0}^{{\tilde{x}}_N}{y}^{\prime }(x)dx={\int }_0^N{y}^{\prime }(\tilde{x}(g))\frac{d\tilde{x}}{dg}dg$$

$$\tilde{x}(g)={\tilde{x}}_g,\text{ for g∈{0,1,…N}}$$

$${\int }_{{\tilde{x}}_0}^{{\tilde{x}}_N}{y}^{\prime }(x)dx\sim \sum _{g=0}^N{y}^{\prime }({\tilde{x}}_g){\frac{d\tilde{x}}{dg}|}_g-\frac{y^{\prime }({\tilde{x}}_0){\frac{d\tilde{x}}{dg}|}_{g=0}+y^{\prime }(x_N){\frac{d\tilde{x}}{dg}|}_{g=N}}{2}+\sum _{k=1}^{\lfloor \frac{m}{2}\rfloor }\frac{B_{2k}}{2k!}\left(G^{(2k-1)}({\tilde{x}}_N)-G^{(2k-1)}({\tilde{x}}_0)\right)$$

$$G(g)=y^{\prime }(\tilde{x}(g))\frac{d\tilde{x}}{dg}$$

$$y=|x{|}^{\gamma }$$

$$x_1=x_0\left(1-\frac{1}{\gamma }\right)$$

$$s_1(x)=3x^2-2x^3,x\in (0,1)$$

$$r_n(x)=\frac{2\Gamma \left(n+\frac{3}{2}\right)}{\sqrt{\pi }\Gamma (n+1)}{\int }_0^x(1-u^2{)}^{n}dx,x\in (-1,1)$$

$$s_n(x)=\frac{1}{2}{r}_n(2x-1)-\frac{1}{2}.$$

$$r_n^{(1)}(x)$$

$$f(x_1,x_2)=x_1^4+x_1^2+3x_1{x}_2+x_2^2+7x_1+9x_2.$$

$$\nabla f$$

$$d\nabla f=H\cdot (d{x}_1,d{x}_2{)}^T+O(d{x}^2),$$

$$\displaystyle H=\left(\begin{array}[]{cc}12x_{1}^{2}+2&3\\ 3&2\end{array}\right)$$

Full text

I Introduction

Here, we consider the general, constrained problem of finding a root of a function $y$ that we cannot measure directly, but for which we can measure its derivatives. Specifically, we imagine that we are given a particular point $(x_{0},y_{0})$ on the curve – i.e., an initialization point, where the function value is provided – and we want to use this and available derivatives to find a nearby root $x_{*}$ of $y$. Standard methods for root finding – e.g., Newton’s method and the bisection method – cannot be applied directly in this case because each requires a function evaluation with each iterative search step. We suggest two methods for root finding that avoid this, both of which involve making use of derivatives to estimate how much $y$ has changed as we adjust $x$ to new points, hopefully closer to a root of $y$. The first of these, we call the local inversion method – this proceeds by “inching” towards the root in $N$ steps, at each step adjusting $x$ by the amount needed to approximately decrease $y$ by $y_{0}/N$. The second method, an approximate Newton’s method, makes use of approximate integrals of $y^{\prime}$ to estimate the function value at subsequent root estimates.

In general, we expect the methods we discuss here to provide competitive strategies whenever the target function $y$ – whose root we seek – is hard to evaluate, relative to its derivatives. This can of course occur in various one-off situations, but there are also certain broad classes of problem where this situation can arise. For example, many physical systems are well-modeled via differential equations that relate one or more derivatives of a target function to a specified driving term – Newton’s laws of motion take this form. Given an equation like this, we can often more readily solve for the derivatives of the target than for the target itself – whose evaluation may require numerical integration. One simple example: Consider a flying, EV drone whose battery drain rate is a given function of its instantaneous height $h$, horizontal speed $\partial_{t}r$, and rate of climb $\partial_{t}h$. If the drone moves along a planned trajectory beginning with a full battery, root finding techniques can be applied to the total charge function – the integral of the instantaneous rate of discharge – to forecast when in the future the battery will reach a particular level. With the methods we describe here, this can be done without having to repeatedly evaluate the precise charge at candidate times.

Whereas standard methods for root finding – e.g., Newton’s method and the bisection method – generally converge exponentially quickly in the number of function evaluations applied, the error in the methods we discuss here converge to zero only algebraically. However, the power at which we get convergence increases with the number of derivatives supplied, and so we can often obtain accurate estimates quite quickly. Another virtue of the methods we discuss here is that the simplest forms of the equations we consider are quite easy to implement in code. Given these points, these methods may be of general interest for those who work with applied mathematics.

The remainder of our paper proceeds as follows: In the next two sections, we provide overviews of the two principal methods we discuss as well as a hybrid strategy. In the following section, we cover some example applications. We then conclude with a quick summary discussion.

Note: We have written up the main strategies we discuss here in an open source python package, inchwormrf. This is available for download on pypi, and on github at github.com/efavdb/inchwormrf. The code can be applied directly or used as demo code from which to base other implementations.

II Local inversion

II.1 Algorithm strategy

The first strategy we will discuss is the local inversion method. We suppose we measure $(x_{0},y_{0})$ and seek a root of $y$. Our strategy will be to iteratively work towards the root in $N$ steps, with $x_{k}$ chosen so that $y(x_{k})\approx y_{0}(1-k/N)$, so that $x_{N}$ is an approximate root. At a given step, then, the goal of the algorithm is to inch its way “downhill” by an amount $y_{0}/N$, towards $y=0$. This process will terminate near a root, provided one does indeed sit downhill of the initial $(x_{0},y_{0})$ point – i.e., provided one can get to a root without having to cross any zero derivative points of $y$. Fig. 1 illustrates the idea.

In order to accomplish the above, we need an approximate local inverse of $y$ that we can employ to identify how far to move in $x$ to drop $y$ by $y_{0}/N$ at each step. Here, we will focus on use of a Taylor series to approximate the inverse, but other methods can also be applied. To begin, we posit an available and locally convergent Taylor series for $y$ at each $x_{k}$,

[TABLE]

We can invert this series to obtain an estimate for the change $dx$ needed to have $y$ adjust by $dy$, as

[TABLE]

where the coefficient $A_{J}$ is given by [2],

[TABLE]

Here, the sum over $s$ values is restricted to partitions of $n-1$,

[TABLE]

The first few $A_{i}$ given by (3) are,

[TABLE]

In a first order approximation, we would use only the expression for $A_{1}$ above in (2). This would correspond to approximating the function $y$ as linear about each $x_{k}$, using just the local slope to estimate how far to move to get to $x_{k+1}$. We can obtain a more refined estimate by using the first two terms, which would correspond to approximating $y$ using a local quadratic, and so on. In practice, it is convenient to “hard-code” the expressions for the first few coefficients, but the number of terms $p(k)$ present in $A_{k}$ increases quickly with $k$, scaling as $\log p(k)\sim\sqrt{n}$ [3]. For general code, then, it’s useful to also develop a method that can directly evaluate the sums (3) for large $n$.

This gives our first strategy: Begin at $(x_{0},y_{0})$, then write

[TABLE]

with $x_{N}$ being the final root estimate and $dx_{k-1}$ being the value returned by (2), in an expansion about $x_{k-1}$, with the shift in $y$ set to $dy=-y_{0}/N$. Example code illustrating this approach with only one derivative is shown in the boxed Algorithm 1.

II.2 Asymptotic error analysis

Upon termination, the local inversion algorithm produces one root estimate, $x_{N}\approx x_{*}$. As noted above, one requirement for this to be a valid root estimate is that there be a root downhill from $x_{0}$, not separated from it by any zero derivative points of $y$. A second condition is that the inverse series (2) provide reasonable estimates for the local inverse at each $x_{k}$. This should hold provided the inverse series converges to the local inverse at each point, within some finite radius sufficiently large to cover the distance between adjacent sample points.

We can identify a sufficient set of conditions for the above to hold as follows: If at the point $x$, $y$ is analytic, with (1) converging within a radius $R$ about $x$ and $|y(\tilde{x})-y(x)|\leq M$ for $\tilde{x}$ within this radius, and $y^{\prime}(x)=s$ is non-zero at $x$, then the inverse of $y$ will also be analytic at $x$, and its Taylor series (2) will converge within a radius of at least [4]

[TABLE]

If the minimum of (9) between $x_{0}$ and $x_{*}$ is some positive value $\tilde{r}>0$, it follows that each of inverse expansions applied in our algorithm will converge between adjacent sample points, provided $N\geq|y_{0}|/\tilde{r}$.

Under the conditions noted above, we can readily characterize the $N$ dependence in the error in the root estimate $x_{N}$. Setting $dy=-\frac{y_{0}}{N}$ and plugging into (2), Taylor’s theorem tells us that the change in $x$ as $y$ moves by $dy$ will be given by

[TABLE]

That is, when using the truncated series expansion to approximate how far we must move to drop $y$ by $y_{0}/N$, we’ll be off by amount that scales like $N^{-(m+1)}$. After taking $N$ steps like this, the errors at each step will accumulate to give a net error in the root estimate that scales like

[TABLE]

That is, the local inversion root estimate’s error converges to zero algebraically with $N$, with a power that increases with $m$, the number of derivatives used in the expansions.

The left plot of Fig. 2 illustrates this scaling in a simple numerical experiment. Here, we’ve applied the local inversion algorithm to find the single real root of $y=x^{5}-3$. Plots of the error versus $N$ are consistent with (11) for each $m$.

II.3 Formal expression for root

If we take the limit of $N\to\infty$ of the $m=1$ version of the local inverse algorithm, we obtain

[TABLE]

a formal, implicit expression for the root we’re estimating. Note that this depends only on the initial point $(x_{0},y_{0})$ and the derivative function, $y^{\prime}$.

Although we won’t pursue this here in detail, we note that (12) provides a convenient starting point for identifying many methods for estimating roots. For example, we can recover Newton’s method by using a linear approximation to the integrand, expanding about $(x_{0},y_{0})$, integrating this, then iterating. Using a discrete approximation to the integral, we can recover the local inversion algorithm above. We can also combine these ideas to obtain an iterated version of the local inversion method.

A special class of functions for which (12) can be quite useful are those for which we can obtain an explicit expression for $y^{\prime}$ in terms of $y$. In this case, (12) becomes an explicit integral which can be approximated in various ways, sometimes giving results that converge much more quickly than does the general local inversion method. For example, with $y=\cos(x)$ and $y^{\prime}=-\sqrt{1-y^{2}}$, (12) returns a familiar integral for $\arccos$. In a numerical experiment, we applied Gaussian quadrature to obtain an estimate for the integral from $y=1/\sqrt{2}$ to [math] – giving estimates for $\pi/4$. The errors in these estimates converged exponentially to zero with the number of samples $N$ used to estimate the integral – beating the algebraic convergence rate of local inversion.

II.4 Behavior near a zero derivative

When a zero derivative point separates $x_{0}$ and $x_{*}$, the local inversion method will not return a valid root estimate. However, we can still apply the local inversion strategy if either $x_{0}$ or $x_{*}$ is itself a zero derivative point. The challenge in this case is that the inverse Taylor series (2) will not converge at the zero derivative point, and some other method of approximating the inverse will be required near this point.

To sketch how this can sometimes be done, we consider here the case where the first derivative of $y$ is zero at $x_{0}$, but the second derivative is not. This situation occurs in some physical applications. E.g., it can arise when seeking turning points of an oscillator, moving about a local minimum at $x_{0}$. To estimate the change in $x$ needed for $y$ to move by $dy$ in this case, we can apply the quadratic formula, which gives

[TABLE]

If we then run the local inversion algorithm with (13) used in place of the inverse Taylor expansion, we can again obtain good root estimates. In this case, applying an error analysis like that we applied above, one can see that the error here will be dominated by the region near the zero derivative, giving a net root estimate error that scales like $N^{-3/2}$. If there is no zero derivative point, we can still use (13) as an alternative to (2). In this case, the root estimates that result again have errors going to zero as $N^{-2}$ – the rate of convergence that we get using series inversion with two derivatives.

To improve on this approach using higher derivative information, one can simply carry out a full asymptotic expansion of (1) near the zero derivative point. This can be applied to obtain more refined estimates for the $dx$ needed at each step in this region. Away from the zero derivative point this local asymptotic expansion will not converge, and one must switch over to the inverse Taylor expansion (2). The resulting algorithm is a little more complicated, but does allow one to obtain root estimates whose errors go to zero more quickly with $N$.

This completes our overview of the local inversion root finding strategy. We now turn to our second approach, an approximate Newton hop strategy.

III approximate Newton hops

III.1 Basic strategy

Here, we discuss our second method for root finding – approximate Newton hops. Recall that in the standard Newton’s method, one iteratively improves upon a root estimate by fitting a line to the function $y$ about the current position, then moving to the root of that line. That is, we take

[TABLE]

Iteratively applying (14), the $x_{k}$ values often quickly approach a root of the function $y$, with the error in the root estimate converging exponentially quickly to zero.

We can’t apply Newton’s method directly under the condition that we consider here, that $y$ is difficult to measure. However, we assume that $y^{\prime}$ can be evaluated relatively quickly, and so we can apply Newton’s method if we replace the $y$ value in (14) above by an approximate integral of $y^{\prime}$, writing

[TABLE]

A simple strategy we can apply to evaluate the above is to sample $y^{\prime}$ at the $N$ equally-spaced points

[TABLE]

Given samples of $y^{\prime}$ at these points, we can then apply the trapezoid rule to estimate the integral in (15). This gives an estimate for $y(x_{k})$ that is accurate to $O(N^{-2})$. Plugging this into (14), we can then obtain an approximate, Newton hop root update estimate. If we iterate, we then obtain an approximate Newton’s method. We find that this approach often brings the system to within the $O(N^{-2})$ “noise floor” of the root within $5$ to $10$ iterations. If $y^{\prime}\not=0$ at the root, this then gives an estimate for the root that is also accurate to $O(N^{-2})$. A simple python implementation of this strategy is given in the boxed Algorithm 2.

III.2 Refinements

We can obtain more quickly converging versions of the approximate Newton hop strategy if we make use of higher order derivatives. The Euler-Maclaurin formula provides a convenient method for incorporating this information [6]. We quote the formula below in a special limit, but first discuss the convergence rates that result from its application: If $m$ derivatives are supplied and we use $N$ points to estimate the integral at right in (15), this approach will give root estimates that converge to the correct values as

[TABLE]

For $m=1$, the error term goes as $O(N^{-2})$, matching that of the trapezoid rule we discussed just above. For even $m\geq 2$, we get $O(N^{-m-2})$. That is, the error in this case goes to zero with two extra powers of $N$, relative to that of the local inversion method, (11). This can give us a very significant speed up.

In our applications of the Euler-Maclaurin formula, we are often interested in a slightly more general situation than that noted above, where the samples are always equally spaced. To that end, we’ll posit now that we have samples of $y^{\prime}$ at an ordered set of values $\tilde{x}_{i}$, with $i\in(0,\ldots,N)$. Direct application of the formula requires evenly spaced samples. To apply it in this case, then, we consider a change of variables, writing

[TABLE]

Here, $\tilde{x}(g)$ is now some interpolating function that maps the domain $[0,N]$ to $[\tilde{x}_{0},\tilde{x}_{N}]$, satisfying

[TABLE]

With this change of variables, the integral at right in (18) is now sampled at an evenly spaced set of points, and we can use the Euler-Maclaurin formula to estimate its value. If we plan to use $m$ derivatives in our analysis, the interpolating function $\tilde{x}$ should also have at least $m$ continuous derivatives. In our open source implementation, we have used a polynomial spline of degree $2\lfloor m/2\rfloor+1$ for $m>1$ and a cubic spline for $m=1$.

We now quote the Euler-Maclaurin formula, as applied to the right side of (18). This reads,

[TABLE]

Here $B_{k}$ is the $k$-th Bernoulli number and $G$ is the function of $g$ given by

[TABLE]

Notice from (20) that the higher order derivatives of $G$ are only needed at the boundaries of integration, where $g=0,N$. We can evaluate these higher-order derivatives using the Faa di Bruno theorem [5], which generalizes the first derivative chain rule for derivatives of composite functions. Each extra derivative with respect to $g$ gives another factor scaling like $O(N^{-1})$, resulting in the scaling of the error quoted in (17) at fixed $m$.

With the formulas noted above, we can implement more refined versions of the approximate Newton hop method: With each iteration, we apply (20) to approximate the current $y$ value, then plug this into (14) to obtain a new root estimate. The samples used to estimate the integral of $y^{\prime}$ can be equally spaced, or not. Again, the former situation is more convenient to code up, but the latter might be used to allow for more efficient sampling: E.g., one might re-use prior samples of $y^{\prime}$ in each iteration, perhaps adding only a single new sample at the current root estimate. Doing this, the samples will no longer each be evenly spaced, which will force use of an interpolation function. However, this can allow for a significant speed up when $y^{\prime}$ evaluations are also expensive.

A hybrid local-inversion, final Newton hop approach can also be taken: First carrying out the local inversion process with $m$ derivatives, will get us to a root estimate that is within $O(N^{-m})$ of the true root. If we save the $y^{\prime}$ values observed throughout this process, we can then carry out a single, final Newton hop using the estimate (20) for the $y$ value at termination of the local inversion process. This will then bring us to a root estimate accurate to $O(N^{-2\lfloor\frac{m}{2}\rfloor-2})$, as in (17).

In our open source package, we include a higher order Newton’s method, but for simplicity have only coded up the case where $N$ new samples of $y^{\prime}$ are evaluated with each iteration, always evenly spaced. When $y^{\prime}$ is easy to evaluate, this approach is much faster than that where a new interpolation is required with each iteration. However, to support the situation where $y^{\prime}$ is also expensive to evaluate, we have also implemented the hybrid local inversion-Newton method. This requires only one interpolation to be carried out in the final step, makes efficient use of $y^{\prime}$ evaluations, and results in convergence rates consistent with our analysis here. The right plot of Fig. 2 illustrates this rate of convergence in a simple application.

IV Example applications

IV.1 Functions with high curvature

Here, we review the point that Newton’s method will sometimes not converge for functions that are highly curved, but that in cases like this the local inversion method generally will. Consider the function

[TABLE]

where $\gamma>0$ is a fixed parameter. If we initialize Newton’s method at the point $(x_{0},x_{0}^{\gamma})$, with $x_{0}>0$, a bit of algebra shows that Newton’s method will move us to the point

[TABLE]

The distance of $x_{1}$ from the origin is a factor of $|1-\gamma^{-1}|$ times that of $x_{0}$. Repeated applications of Newton’s method will result in $|x_{k}|=x_{0}|1-\gamma^{-1}|^{k}$. If $\gamma>1/2$, the method converges exponentially quickly to the root. However, if $\gamma<1/2$, each iteration drives the estimate further away from the root of $y$, due to “overshooting”.

High curvature is not an issue for the local inversion method: By inching along, it can quickly recalibrate to changes in slope, preventing overshooting. Note, however, that applying a final, approximate Newton hop after the local inversion process will no longer improve convergence in cases of high curvature such as this.

IV.2 Smoothstep

The smoothstep function [1]

[TABLE]

is an “s-shaped” curve that is commonly used in computer graphics to generate smooth transitions, both in shading applications and for specifying motions. This function passes through $(0,0)$ and $(1,1)$ and has zero derivatives at both of these points. More generally, one can consider the degree-$(2n+1)$ polynomial $s_{n}$ that passes through these points and has $n\geq 1$ zero derivatives at both end points. Simple inverse functions are not available for $n>1$, so to determine where $s_{n}(x)$ takes on a specific value, one must often resort to root finding.

The following identity makes our present approach a tidy one for identifying inverse values for general $n$: If we define,

[TABLE]

then

[TABLE]

Whereas evaluation of $r_{n}(x)$ requires us to evaluate the integral above, the first derivative of $r_{n}$ is simply

[TABLE]

Higher order derivatives of $r_{n}$ are easy to obtain from this last line as well. Given these expressions, we can easily invert $r_{n}$ for any value of $n$ using our methods. Fig. 3 provides an illustration.

IV.3 Multiple dimensions

Solutions to coupled sets of equations can also be treated using the local measurement method. We illustrate this here by example. For simplicity, we’ll consider a cost function that can be written out by hand, but understand that the method we describe will be most useful in cases where the derivatives of a cost function are more easy to evaluate than the cost function itself. The function we’ll consider is

[TABLE]

If we wish to find a local extremum of this function, we can set the gradient of the above to zero, which gives

[TABLE]

a coupled set of two equations for $x_{1}$ and $x_{2}$.

To apply the local inversion technique over $N$ steps, we need to expand the gradient about a given position. To first order, this gives

[TABLE]

where

[TABLE]

is the Hessian of (28), and

[TABLE]

is the target gradient shift per step. Inverting (30), we obtain

[TABLE]

There is a single real solution to (29) at $(x_{1},x_{2})\approx(1.35172698,-6.25699505)$. Running the local inversion algorithm (35) finds this solution with error decreasing like $O(N^{-1})$.

To improve upon the above solution, we can either include more terms in the expansion (30), or we can apply a final, approximate Newton hop to refine the solution after the inching process. To carry that out, we need to consider how to estimate the final values of our gradient function. To that end, we note that the integral we are concerned with is now a path integral

[TABLE]

Taking a discrete approximation to this integral, and then applying a final Newton hop, we obtain estimates with errors that converge to zero like $O(N^{-2})$.

V Discussion

Here, we have introduced two methods for evaluating roots of a function that require only the ability to evaluate the derivative(s) of the function in question, not the function itself. This approach might be computationally more convenient than standard methods whenever evaluation of the function itself is relatively inconvenient for some reason. The main cost of applying these methods is that convergence is algebraic in the number of steps taken, rather than exponential – the typical convergence form of more familiar methods for root finding, such as Newton’s method and the bi-section method. However, one virtue of these approaches is that higher order derivative information is very easily incorporated, with each extra derivative provided generally offering an increase to the convergence rate. This property can allow for the relatively slow rate of convergence concern to be mitigated.

The three main approaches that we have detailed here are as follows: (1) the local inversion method, (2) the approximate Newton’s method, and (3) the hybrid local inversion-Newton method. The first of these should be preferred when looking at a function that is highly curved, since Newton’s method is subject to overshooting in this case. The second will converge more quickly for functions that are not too highly curved, and so is preferable in this case. Finally, the third method provides a convenient choice when working with functions where Newton’s method applies, but whose derivatives are also somewhat costly to evaluate: Our implementation here gives the convergence rate of Newton’s method and also minimizes the number of derivative calls required.