Root Finding With Derivatives: Newton-Raphson, Halley & Schroeder

Synopsis

#include <boost/math/tools/roots.hpp>

namespace boost{ namespace math{
namespace tools{

template <class F, class T>
T newton_raphson_iterate(F f, T guess, T min, T max, int digits);

template <class F, class T>
T newton_raphson_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter);

template <class F, class T>
T halley_iterate(F f, T guess, T min, T max, int digits);

template <class F, class T>
T halley_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter);

template <class F, class T>
T schroeder_iterate(F f, T guess, T min, T max, int digits);

template <class F, class T>
T schroeder_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter);

}}} // namespaces

Description

These functions newton_raphson_iterate, halley_iterate, and schroeder_iterate all perform iterative root finding using derivatives.

(These functions should be used whenever derivative(s) are available, otherwise use root-finding without derivatives).

newton_raphson_iterate performs second order Newton-Raphson iteration,
halley_iterate andschroeder_iterate perform third order Halley and Schroeder iteration.

The functions all take the same parameters:

Parameters of the root finding functions

F f

Type F must be a callable function object that accepts one parameter and returns a tuple: the tuple type can be any suitable implementation which adhers to the tuple concept in C++11 (for example std::tuple, std::pair, boost::tuple, boost::fusion::tuple).

For the second order iterative methods (Newton Raphson) the tuple should have two elements containing the evaluation of the function and its first derivative respectively.

For the third order methods (Halley and Schroeder) the tuple should have three elements containing the evaluation of the function and its first and second derivatives respectively.

T guess

The initial starting value. A good guess is crucial to quick convergence!

T min

The minimum possible value for the result, this is used as an initial lower bracket.

T max

The maximum possible value for the result, this is used as an initial upper bracket.

int digits

The desired number of binary digits of precision.

uintmax_t& max_iter

An optional maximum number of iterations to perform. Default max_iter is (std::numeric_limits<boost::uintmax_t>::max)(), perhaps 18446744073709551615 which may be more than you want for wait for!
On exit this is set to the actual number of iterations performed: consequently you should always check this value on exit to determine whether iteration terminated because we had converged on the root, or because the maximum number of iterations was exceded.

When using these functions you should note that:

Default max_iter = (std::numeric_limits<boost::uintmax_t>::max)() is effectively 'iterate for ever'!.
The functions may be very sensitive to the initial guess. Typically they converge very rapidly if the initial guess has two or three decimal digits correct. However, convergence can be no better than bisection, or in some rare cases, even worse than bisection if the initial guess is a long way from the correct value and the derivatives are close to zero.
These functions include special cases to handle zero first (and second where appropriate) derivatives, and fall back to bisection in this case. However, it is helpful if functor F is defined to return an arbitrarily small value of the correct sign rather than zero.
If the derivative at the current best guess for the result is infinite (or very close to being infinite) then these functions may terminate prematurely. A large first derivative leads to a very small next step, triggering the termination condition. Derivative-based iteration may not be appropriate in such cases.
If the function is 'Really Well Behaved' (monotonic and has only one root), and only if this is the case then the bracket bounds min and max may as well be set to the widest limits like zero and numeric_limits<T>::max().
If the function more complex and may have more than one root or a pole, then the choice of bounds provides protection against jumping out to seek the 'wrong' root.
If the function is badly behaved then these functions fall back to bisection if the next computed step would take the next value out of bounds. The bounds are updated after each step to ensure this leads to convergence. In this situation, a good initial guess backed up by asymptotically-tight bounds will improve performance no end - rather than relying on bisecting [0, LARGE] for example.
The value of digits is crucial to good performance of these functions, if it is set too high then, at best, you will get one extra (unnecessary) iteration, and, at worst, the last few steps will proceed by bisection. Remember that the returned value can never be more accurate than f(x) can be evaluated, and that if f(x) suffers from cancellation errors as it tends to zero, then the computed steps will be effectively random. The value of digits should be set so that iteration terminates before this point: remember that for second and third order methods the number of correct digits in the result is increasing quite substantially with each iteration, digits should be set by experiment so that the final iteration just takes the next value into the zone where f(x) becomes inaccurate.
A good choice for the maximum number of iterations would be boost::math::policies::get_max_root_iterations<Policy>()).
If you need some diagnostic output to see what is going on, you can #define BOOST_MATH_INSTRUMENT before the #include <boost/math/tools/roots.hpp>, and also ensure that display of all the possibly significant digits with cout.precision(std::numeric_limits<double>::max_digits10): but be warned, this may produce copious output!
Finally: you may well be able to do better than these functions by hand-coding the heuristics used so that they are tailored to a specific function. You may also be able to compute the ratio of derivatives used by these methods more efficiently than computing the derivatives themselves. As ever, algebraic simplification can be a big win.

Newton Raphson Method

Given an initial guess x0 the subsequent values are computed using:

Out of bounds steps revert to bisection of the current bounds.

Under ideal conditions, the number of correct digits doubles with each iteration.

Halley's Method

Given an initial guess x0 the subsequent values are computed using:

Over-compensation by the second derivative (one which would proceed in the wrong direction) causes the method to revert to a Newton-Raphson step.

Out-of-bounds steps revert to bisection of the current bounds.

Under ideal conditions, the number of correct digits trebles with each iteration.

Schroeder's Method

Given an initial guess x0 the subsequent values are computed using:

Over-compensation by the second derivative (one which would proceed in the wrong direction) causes the method to revert to a Newton-Raphson step. Likewise a Newton step is used whenever that Newton step would change the next value by more than 10%.

Out of bounds steps revert to bisection of the current bounds.

Under ideal conditions, the number of correct digits trebles with each iteration.

Example

Let's suppose we want to find the cube root of a number: the equation we want to solve along with its derivatives are:

To begin with let us solve the problem using Newton-Raphson iterations. We will begin by defining a function object (functor) that returns the evaluation of the function to solve, along with its first derivative f'(x):

template <class T>
struct cbrt_functor_1stderiv
{ // Functor returning function and 1st derivative.

  cbrt_functor_1stderiv(T const& target) : value(target)
  { // Constructor stores the value to be 'cube_rooted'.
  }

  std::pair<T, T> operator()(T const& z) // z is best estimate so far.
  { // Return both f(x) and derivative f'(x).
    T fx = z*z*z - value; // Difference estimate fx = x^3 - value.
    T d1x =  3 * z*z; // 1st derivative d1x = 3x^2.
    return std::make_pair(fx, d1x); // 'return' both fx and d1x.
  }
private:
  T value; // to be 'cube_rooted'.
}; // cbrt_functor_1stderiv

Implementing the cube root is fairly trivial now, the hardest part is finding a good approximation to begin with: in this case, we'll just divide the exponent by three. (There are better but more complex 'guess' algorithms used in 'real-life'.)

template <class T>
T cbrt_1deriv(T x)
{ //! \return cube root of x using 1st derivative and Newton_Raphson.
  using namespace std; // For frexp, ldexp, numeric_limits.
  using namespace boost::math::tools; // For newton_raphson_iterate.

  int exponent;
  frexp(x, &exponent); // Get exponent of x (ignore mantissa).
  T guess = ldexp(1., exponent/3); // Rough guess is to divide the exponent by three.
  // Set an initial bracket interval.
  T min = ldexp(0.5, exponent/3); // Minimum possible value is half our guess.
  T max = ldexp(2., exponent/3);// Maximum possible value is twice our guess.
  // digits used to control how accurate to try to make the result.
  int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T.

  boost::uintmax_t maxit = 20; // Optionally limit the number of iterations.
  //cout << "Max Iterations " << maxit << endl; //
  T result = newton_raphson_iterate(cbrt_functor_1stderiv<T>(x), guess, min, max, digits, maxit);
  // Can check and show how many iterations (updated by newton_raphson_iterate).
  // cout << "Iterations " << maxit << endl;
  return result;
} // cbrt_1deriv

Using the test data in ../../test/cbrt_data.ipp, the test ../../test/test_cbrt.cpp found the cube root exact to the last digit in every case, and in no more than 6 iterations at double precision.

However, you will note that a high (maximum) precision was used in this example, exactly what was warned against earlier on in these docs! In this particular case it is possible to compute f(x) exactly and without undue cancellation error, so a high limit is not too much of an issue. However, reducing the limit to std::numeric_limits<T>::digits * 2 / 3 gave full precision in all but one of the test cases (and that one was out by just one bit). The maximum number of iterations remained 6, but in most cases was reduced by one.

So we might decide to trade the accuracy and limit iterations for speed.

int get_digits = (digits * 2) /3; // Two thirds of maximum possible accuracy.
boost::uintmax_t maxit = 7; // Limit the number of iterations.

Note also that the above code omits an optimization by computing z², and reusing it, omits error handling, and does not handle negative values of z correctly. (As is traditional, these are left as an exercise for the reader!)

The boost::math::cbrt function cbrt also includes these and other improvements.

Now let's adapt the functor slightly to return the second derivative f''(x) as well:

template <class T>
struct cbrt_functor_2deriv
{ // Functor returning both 1st and 2nd derivatives.
  cbrt_functor_2deriv(T const& to_find_root_of) : value(to_find_root_of)
  { // Constructor stores value to find root of, for example:
  }

  // using boost::math::tuple; // to return three values.
  std::tuple<T, T, T> operator()(T const& x)
  { // Return both f(x) and f'(x) and f''(x).
    T fx = x*x*x - value; // Difference (estimate x^3 - value).
    T dx = 3 * x*x; // 1st derivative = 3x^2.
    T d2x = 6 * x; // 2nd derivative = 6x.
    return std:: make_tuple(fx, dx, d2x); // 'return' fx, dx and d2x.
  }
private:
  T value; // to be 'cube_rooted'.
}; // struct cbrt_functor_2deriv

And then adapt the cbrt function to use Halley iterations:

template <class T>
T cbrt_2deriv(T x)
{ // return cube root of x using 1st and 2nd derivatives and Halley.
  using namespace std;  // Help ADL of std functions.
  using namespace boost::math; // halley_iterate

  int exponent;
  frexp(x, &exponent); // Get exponent of z (ignore mantissa).
  T guess = ldexp(1., exponent/3); // Rough guess is to divide the exponent by three.
  T min = ldexp(0.5, exponent/3); // Minimum possible value is half our guess.
  T max = ldexp(2., exponent/3); // Maximum possible value is twice our guess.

  int digits = std::numeric_limits<T>::digits /2 ; // Half maximum possible binary digits accuracy for type T.
  boost::uintmax_t maxit = 10;
  T result = halley_iterate(cbrt_functor_2deriv<T>(x), guess, min, max, digits, maxit);
  // Can show how many iterations (updated by halley_iterate).
  // cout << "Iterations " << maxit << endl;
  return result;
} // cbrt_2deriv(x)

Note that the iterations are set to stop at just one-half of full precision, and yet, even so, not one of the test cases had a single bit wrong. What's more, the maximum number of iterations was now just 4.

Just to complete the picture, we could have called schroeder_iterate in the last example: and in fact it makes no difference to the accuracy or number of iterations in this particular case. However, the relative performance of these two methods may vary depending upon the nature of f(x), and the accuracy to which the initial guess can be computed. There appear to be no generalisations that can be made except "try them and see".

We could also have used a C++ lambda expression in place of the function object.

Finally, had we called cbrt with Boost.Multiprecision or NTL::RR set to 1000 bit precision, then full precision still can be obtained with just 7 iterations. To put that in perspective, an increase in precision by a factor of 20 has less than doubled the number of iterations. That just goes to emphasise that most of the iterations are used up getting the first few digits correct: after that these methods can churn out further digits with remarkable efficiency.

Or to put it another way: nothing beats a really good initial guess!

See root_finding_example.cpp for full source code.