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

On Thu, Mar 29, 2012 at 11:22:44AM +0000, Amnon Harel wrote:
>
> After a bit of study and experiment, I realized that Brent's
> Method is basically the only root finding algorithm worth working with.
> ... several problems are ... are listed from the profound to the technical.
>

FWIW, I used to work with algorithms for finding minimums and maximums.

My profound observation can be summed by the question: "how come there is no one universal algorithm that works with any function under every possible conditions?!?". (Think about it).

Several previous lives ago I also studied mathematical proofs for these algorithms, and basically they all involve the proverbial spherical cows in vacuum. For any kind of realistic function, they are euristics at best. Unlike your Grand-father's sqrt() algorithm, they do not have to give the correct answer always every time.

These proofs usually include claims on convergence speed. But the proofs require such unrealistically "nice" functions, that in practice all convergence speed claims do not hold and all theoretically "fast" methods in practice converge much slower than expected (or more usually go straight into the la-la land). Then they have to be kludged in various ways to make them work at all (enter Brent's method).

I also have looked at many implementations of root and minimum finding algorithms, and most do not implement the actual algorithm after which they are named (as described in words, or as implemented in Numerical Recepies or in Netlib). I think this is because the algorithms always have to be kludged to work with any realistic function and every implementer/kludger only does just enough to "make it work" on his particular function at hand.

I never had to find the root (or minimum) of functions like "y=x^2".

The functions I worked with are usually defined on an interpolated grid (computing the value of the function at each grid point involves integration of charged particle trajectories through a magnetic field) and I find that anything too fancy (i.e. with gradients, etc) quickly goes ashtray because the function is "too noisy".

In my experience, only the grid-search methods are fully reliable (i.e. binary section, golden section, etc). If they are too slow, one can try to accelerate them by using derivatives, interpolation and extrapolation. The result would look very much like Brent's method (binary section accelerated by linear or quadratic interpolation).

But in accordance with Murphy's law, all speedups would usually slow you down, and that's what I saw in practice - nothing worked better for me than the good old golden section method.

(This is for 1-dimension functions. For N-dimenstion, the same holds - in the order of robustness: grid search, simplex/complex methods, gradient methods, and forget about using second derivatives and stuff).

K.O.

>
> 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
>

-- 
Konstantin Olchanski
Data Acquisition Systems: The Bytes Must Flow!
Email: olchansk-at-triumf-dot-ca
Snail mail: 4004 Wesbrook Mall, TRIUMF, Vancouver, B.C., V6T 2A3, Canada