Brent's method in Julia

Brent's method or Wijngaarden-Brent-Dekker method is a root-finding algorithm which combines the bisection method, the secant method and inverse quadratic interpolation. This method always converges as long as the values of the function are computable within a given region containing a root.

Given a function \(f(x)\) and the bracket \([x_0, x_1]\) two new points, \(x_2\) and \(x_3\), are initialized with the \(x_1\) value. \(x_0\) and \(x_1\) are swapped if \(|f(x_0)| < |f(x_1)|\).

Then, in each iteration if the evaluation of the points \(x_0\), \(x_1\) and \(x_2\) are different (according to a certain tolerance) the inverse quadratic interpolation is used to get the new guess \(x\). Otherwise, the linear interpolation (secant method) is used to obtain the guess. After that, if any of the following conditions are satisfied \(x\) will be redefined using the bisection method:

• If \(x\) does not lie between \((3x_0+x_1)/4\) and \(x_1\).
• If in the previous iteration the bisection method was used or it is the first iteration and \(|x-x_1| \geq {1 \over 2} \,|x_1-x_2|\).
• If in the previous iteration the bisection method was not used and \(|x-x_1| \geq {1 \over 2} \,|x_2-x_3|\).
• If in the previous iteration the bisection method was used or it is the first iteration and \(|x_1-x_2| < |\delta|\), where \(\delta\) is a tolerance factor.
• If in the previous iteration the bisection method was not used and \(|x_2-x_3| < |\delta|\), where \(\delta\) is a tolerance factor.

We define \(\delta\) as \(2 \epsilon x_1\), where \(\epsilon\) is the machine epsilon.

\(x_3\) and \(x_2\) are redefined in each iteration with \(x_2\) and \(x_1\) value, respectively, and the new guess \(x\) will be set as \(x_1\) if \(f(x_0)f(x)<0\) or as \(x_2\) otherwise. Again, \(x_0\) and \(x_1\) are swapped if \(|f(x_0)| < |f(x_1)|\). The algorithm converges when \(f(x)\) or \(|x_1-x_0|\) are small enough, both according to tolerance factors.

In the following implementation, the inverse quadratic interpolation is applied directly. In other cases, like the implementation in Numerical recipes, used for example in Boost, the Lagrange polynomial is reduced defining the variables \(p\), \(q\), \(r\), \(s\) and \(t\) as explained in MathWorld and \(x\) value is not overwritten with the bisection method, but modified. For simplicity of the code, here the inverse quadratic interpolation is applied directly as in the entry Inverse quadratic interpolation in Julia and the new guess is overwritten if needed.

In Julia:

function brent(f::Function, x0::Number, x1::Number, args::Tuple=();
               xtol::AbstractFloat=1e-7, ytol=2eps(Float64),
               maxiter::Integer=50)
    EPS = eps(Float64)
    y0 = f(x0,args...)
    y1 = f(x1,args...)
    if abs(y0) < abs(y1)
        # Swap lower and upper bounds.
        x0, x1 = x1, x0
        y0, y1 = y1, y0
    end
    x2 = x0
    y2 = y0
    x3 = x2
    bisection = true
    for _ in 1:maxiter
        # x-tolerance.
        if abs(x1-x0) < xtol
            return x1
        end

        # Use inverse quadratic interpolation if f(x0)!=f(x1)!=f(x2)
        # and linear interpolation (secant method) otherwise.
        if abs(y0-y2) > ytol && abs(y1-y2) > ytol
            x = x0*y1*y2/((y0-y1)*(y0-y2)) +
                x1*y0*y2/((y1-y0)*(y1-y2)) +
                x2*y0*y1/((y2-y0)*(y2-y1))
        else
            x = x1 - y1 * (x1-x0)/(y1-y0)
        end

        # Use bisection method if satisfies the conditions.
        delta = abs(2EPS*abs(x1))
        min1 = abs(x-x1)
        min2 = abs(x1-x2)
        min3 = abs(x2-x3)
        if (x < (3x0+x1)/4 && x > x1) ||
           (bisection && min1 >= min2/2) ||
           (!bisection && min1 >= min3/2) ||
           (bisection && min2 < delta) ||
           (!bisection && min3 < delta)
            x = (x0+x1)/2
            bisection = true
        else
            bisection = false
        end

        y = f(x,args...)
        # y-tolerance.
        if abs(y) < ytol
            return x
        end
        x3 = x2
        x2 = x1
        if sign(y0) != sign(y)
            x1 = x
            y1 = y
        else
            x0 = x
            y0 = y
        end
        if abs(y0) < abs(y1)
            # Swap lower and upper bounds.
            x0, x1 = x1, x0
            y0, y1 = y1, y0
        end
    end
    error("Max iteration exceeded")
end

As an example, for the function \(f(x)=x^4-2x^2+1/4 \quad (0 \leq x \leq 1)\), the solution is \(\sqrt{1-\sqrt{3}/2}\):

julia> (f(x)=x^4-2x^2+1/4; x0=0; x1=1; brent(f,x0,x1))
0.3660254037844386
julia> sqrt(1-sqrt(3)/2)
0.3660254037844387