2
$\begingroup$

I need to find the (unique) root of a nonlinear function $f(x)$, $x \in \mathbb{R}$. For the record, $f(x)$ is the CDF of a probability density minus a constant $0 < p < 1$ (I am inverting the CDF; see this question). The pdf is defined as a 1-D mixture of Gaussians times a polynomial. $f$ is analytical but relatively expensive to compute, so I would like to speed up the root-finding algorithm.

So far, I have been using fzero in MATLAB, which implements Brent's method. In several places in this forum (see also the question above) I have seen recommended Alefeld-Potra-Shi's method instead, which makes me think it is state-of-the art (I am not particularly familiar with the root-finding literature).

What are the performance differences between Brent's and Alefeld-Potra-Shi's methods, and when is one method expected to be better than the other?

If Alefeld-Potra-Shi is potentially better for my problem, I might work on a MATLAB port of Alefeld-Potra-Shi, which I could not find so far (possibly there is an Octave implementation, and a Julia one here).

Improve this question
$\endgroup$
6
  • 1
    $\begingroup$ If the probabilty function decays sufficiently fast to zero. You could get away in evaluating your function at quadrature points by using domain truncation and chebyshev expansion. Finding roots for the approximating polynomial is relatively cheap and either a very good starting point for a finer iteration or already good enough. E.g. see here. As a bonus, you get all roots in the interval if done correctly, not just the one closest to the starting point for Brent's search. $\endgroup$
    – Bort
    May 31, 2016 at 16:07
  • $\begingroup$ (+1) Thanks @Bort, it seems like an interesting approach. I am concerned though that the cost of computing the approximating polynomial would be pretty much the same as computing the CDF (the CDF is only mildly expensive, the cost stems from a bunch of series expansions and evaluations of erfc). Also, $f$ has only one root, so I don't need a method that returns multiple roots. As a starting point, for the moment I am using a normal approximation to the pdf which seems to work okay. $\endgroup$
    – lacerbi
    May 31, 2016 at 17:40
  • 1
    $\begingroup$ Alefeld, Potra, Shi (M.Com 1993) compares the two methods (see Table 5.2). To compare them directly yourself, you could try using C++ with Boost (which implements Alefeld-Potra-Shi) and GSL (which implements Brent's method). That should give you a very direct comparison for your exact function. $\endgroup$
    – Kirill
    May 31, 2016 at 19:36
  • $\begingroup$ Thanks @Kirill. According to Table 5.2, it seems there is not much difference between Brent's and APS (on average) for the test functions they use; so I probably need to test it on my own problem as you suggest. $\endgroup$
    – lacerbi
    May 31, 2016 at 20:02
  • $\begingroup$ @lacerbi I see. Reading your first question, I thought you might have also multiple roots because of the missing monotonicity. For which purpose do you need the inverse cdf? In case of drawing samples, the exact form might not be needed to be efficient. $\endgroup$
    – Bort
    Jun 2, 2016 at 12:03

1 Answer 1

Reset to default
2
$\begingroup$

The asymptotic order of convergence of Brent's method tends to be either 1.618 and 1.689, of which the Alefeld-Potra-Shi's methods lie directly in between. The main difference of the Alefeld-Potra-Shi's methods are the tightness of the bounds, which Brent's method may fail to give during intermediate iterations.

Overall, I do not recommend Alefeld-Potra-Shi's methods for your purposes. The method is intended to be a black-box sledgehammer for bracketing the root. Brent's method is in most cases 'about as good'.

For a much simpler method which may meet your purposes, Chandrupatla's method may be of interest. It has the same asymptotic nature as Brent's method for simple roots, but has a tendency to converge faster initially due to a more 'intelligent' algorithm for using either bisection or interpolation.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.