Issues with BrentRootFinder, alternative code, performance studies.

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:

1. 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".
2. 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|.
3. 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.
4. 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 - 13:22:53 CEST