1
$\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).

|cite|improve this question
$\endgroup$
  • 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 '16 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 '16 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 '16 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 '16 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 '16 at 12:03

Your Answer

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

Browse other questions tagged or ask your own question.