Hi everyone,

I have found a bug in the function brent in section 10.2

The line that reads:
if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))

should read:
if (fabs(p) >= fabs(0.5*q*etemp) || p >= q*(a-x) || p >= q*(b-x))

The "less than or equal to" sign should be "greater than or equal to. I confirmed this bug by referring to Brent's original book, which contains the Algol code for this algorithm.

In the current form the algorithm fails in many cases through a divide by zero error:eek:. I am surprised that it has not been reported until now.

Calculon.

Hi,

Brent's algorithm is given correctly in NR. The problem is a misprint in the Algol algorithm in his book, on p. 80. As you point out, it has the opposite inequality for the (a-x) piece. However, on p. 74 where he describes the algorithm, he says that you use golden section if x+p/q is not inside (a,b). You can easily check that this gives the inequalities as coded in NR. (Remember that q>0.)

A Fortran version of the algorithm by Brent is at http://www.netlib.org/go/fmin.f
You can see there that he has corrected his Algol version. Similarly, the book by Forsythe, Macolm and Moler "Computer Methods for Mathematical Computations" has a Fortran version of fmin that agrees with the Netlib version.

So this is definitely Brent's algorithm. I'm not sure why you might be having problems - feel free to submit the simplest example you have that shows the problem.

Best,
Saul Teukolsky

I have been trying to use a version of Brent's algorithm for line minimisation in a variety of multidimensional optimisation algorithms.

The problem with the NR implementation of Brent's algorithm is that there is no explicit checking that the three points used to fit the quadratic are not in fact co-linear. This causes the variable "p" to become zero. When used in multidimensional cases where the minimum is against the search direction, this fails all of the bad fit tests and the program attempts a quadratic step. Of course this leads to the divide by zero error mentioned in my previous post.

This is really for the most part a problem in optimisation algorithms where the minimum along the search line occurs against the search direction. This is most obviously possible in Powell's algorithm, where the direction of the search is arbitrary. It is not as much of a problem in other algorithms which more carefully select the line search direction. However, I think that the given implementation of Brent's algorithm in NR is unstable when used in Powell optimisations.

The needed fix to the algorithm is however very simple. All one needs to do is check to see if the quadratic fit points are co-linear, and if so take the golden section step. The minimal way to do this in the code is to modify the line:

if(fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))

to:

if(fabs(p) < ZERO || fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))

Here ZERO is a number that is practiacally zero on the scale of the problem that you are dealing with. In my code I take the value of 1E-25. This vastly improves the stability of the algorithm at very little computational cost.

In my opinion this should be included in any implementation of the algorithm, unless it can be guaranteed that the algorithm will work when the three points used to fit the qudaratic are co-linear.

Calculon.