Re: Issues with BrentRootFinder, alternative code, performance studies.

Dear Amnon,

 Thanks a lot for your detailed and interesting study in the Brent Method. The algorithms implemented in the BrentMethods.cxx comes from the Numerical Receipt book, it is possible that differs from the original and is not optimal. I had already noticed that the Brent implementation in MathMore, available in the GSLRootFinder class performed better. have you tried also to use that algorithm ? Also, as you said, since it implements TF1::GetX, it is more target for robustness vs performances. It is true basically the same algorithm for finding the minimum (TF1::GetMinimumX ) and the root (TF1::GetX) is used and maybe this is not optimal and we should do as you suggested.

I will look at your code in detail and tested it, and I am willing to include in the next ROOT version. If you also find further improvements or fixes, please let me know,

 Thanks again for your contribution,
Best,
 Lorenzo
On Mar 29, 2012, at 1:22 PM, Amnon Harel wrote:

> Dear ROOT maintainers,
>
> After a bit of study and experiment, I realized that Brent's
> Method is basically the only root finding algorithm worth
> working with.
> Looking into ROOT's BrentRootFinder, several problems are apparent.
> They are listed from the profound to the technical.
>
> 0) First and foremost, despite misleading names, comments, and
> documentation, this code does not implement Brent's Root
> Finding Algorithm.
>
> Brent's algorithm is a bracketing algorithm. Thus, by definition
> it relies heavily on the sign of the user's function [1,2].
> This code calls the functions in BrentMethods with type=4
> and fy=0. Thus these functions use only the absolute value
> of the user's function and are completely unaware of its sign.
>
> This algorithm may have its merits. In particular, since this
> is the default algorithm in TF1::GetX(), it is arguably nice to
> return something useful even if there is no root in the range.
> Presumably this led to the code mutating away from being an
> implementation of Brent's algorithm (a search for f=0) to being
> a search for minimal |f|. I doubt that this derivative
> of Brent's algorithm is a good way to do that, but perhaps
> its merits escape me. Two arguments against it are below.
>
> For what it's worth:
> A) I suggest separating out the concepts of finding f=f0 ("GetX"?)
> and finding a local minimum of |f-f0| ("GetCloseX"?).
> In particular, the documentation of TF1::GetX claims the former,
> but the latter is implemented.
> It's probably better to have two functions, with the f=f0
> function being allowed to report "no f=f0 found".
> B) If searching for min|f-f0|, I suggest to first search for a
> root, which can be done very efficiently, before settling
> for a local minimum of |f|.
> C) My proposed solution also gives some estimate of minimum
> |f-f0| when there is no root. But it's probably less accurate
> than the current algorithm, if I correctly diagnosed all its
> problems.
>
>
> 1) The early grid search is wrong, as a direct results of (0).
> Instead of attempting to bracket the root (as documented correctly)
> it ignores the sign of the function.
> Thus the results are rather random, and only a very fine
> grid tends to, accidently, find a correct bracket.
>
> Here's an example of wrong results from the grid search, the
> ranges in square brackets are from debug prints in MinimStep
> (on entry and exit), and the "cf" numbers are from the user-function
> itself. Npx = 6.
>
> MinimStep type: 4, [ 0 - 39.1663804 ] cf(0)=-1738724.88 cf(7.83327607)=-2379204.37
> cf(15.6665521)=-2739932.05 cf(23.4998282)=-2745476.65 cf(31.3331043)=-2132798.29
> cf(39.1663804)=2261792.81 -> [ 0 - 7.83327607 ]
>
> Whereas the correct bracket would be [ 31.3331043 - 39.1663804 ].
>
>
> 2) Furthermore, in a real Brent's algorithm the time-consuming grid
> search is not needed whenever sign(f(xlow)) != sign(f(xhigh)).
> In the example above:
> [ 0 - 39.1663804 ] is a fine starting bracket.
>
>
> 3) The handling of logStep isn't consistent. With this option:
> - the early exit condition npx<2 leaves the logs in xmin and xmax.
> - the "dx" used to update xmin and xmax before exiting does so
> on a linear scale (while the grid tests used the log scale).
>
>
> 4) The "midpoint" in the early exit clause is not a midpoint.
> "-" --> "+".
>
>
> Alternative code
> ================
>
> An implementation that resolves these issues and vastly out-performs
> the current BrentRootFinder is available in:
>
> http://www-d0.fnal.gov/~aharel/ModBrentMethods.c
> http://www-d0.fnal.gov/~aharel/ModBrentMethods.h
> http://www-d0.fnal.gov/~aharel/ModBrentRootFinder.c
> http://www-d0.fnal.gov/~aharel/ModBrentRootFinder.h
>
> Notes:
> 1) Somewhat recursive grid search when no bracket found
> 2) I could not disentangle the concepts of |g| and g in the
> current implementation, so I based my implementation on [1].
> 3) I don't know what bug was fixed in 11/3/2010, but the condition
> in question isn't the one Brent used ([1-3] all use the same
> condition).
> 4) The code can be further improved. If there is any interest in
> making it part of ROOT, I'm willing to do so.
> - streamlining, at the cost of moving further away from the
> current interface (e.g. have grid function return a bool
> indicating whether a bracket was found).
> - redesign to reduce multiple calls of the user function at
> the same location.
> - the logStep option isn't debugged.
>
>
> Performance and Testing results
> ===============================
>
> Though the functions I use in tests have a particular form, their
> parameter space (7-9 parameters, depending on how you count) is huge
> and they seems quite difficult for these algorithms.
> I run my algorithm with Npx=20. This value was not optimized.
> This is a fair comparison as reducing Npx below ~50 in the old
> algorithm causes it to fail regularly in my tests.
>
> In 10K test cases where at least one root exists, the existing ROOT code
> was better than the proposed code (that is, found a lower |f(x)|) 912
> times, and worst 9088 times. Some of "worst" cases are when it completely
> misses the root.
> The average number of calls of the user function was 127 for the old
> code, and 16 for the new code.
>
> That is, the new code is eight times faster and the probability of
> getting a wrong answer is drastically reduced (due to the nature of
> the problem, it can not be eliminated).
>
> cheers,
> Amnon Harel,
> University of Rochester
>
>
> [1] I do not have Brent's book, but his original code is available in:
> http://www.netlib.org/go/zeroin.f
>
> [2] http://en.wikipedia.org/wiki/Brent's_method
>
> [3] http://math.fullerton.edu/mathews/n2003/BrentMethodMod.html
>
Received on Thu Mar 29 2012 - 14:57:36 CEST