boost/libs/math/doc/internals/roots.qbk

[section:roots Root Finding With Derivatives: Newton-Raphson, Halley & Schroeder]

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

[h4 Description]

These functions all perform iterative root finding using derivatives:

* `newton_raphson_iterate` performs second order
[link math_toolkit.internals1.roots.newton Newton-Raphson iteration],

* `halley_iterate` and`schroeder_iterate` perform third order
[link math_toolkit.internals1.roots.halley Halley] and [link math_toolkit.internals1.roots.schroeder Schroeder] iteration.

The functions all take the same parameters:

[variablelist Parameters of the root finding functions
[[F f] [Type F must be a callable function object that accepts one parameter and
        returns a __tuple:

For the second order iterative methods ([@http://en.wikipedia.org/wiki/Newton_Raphson Newton Raphson])
        the __tuple should have *two* elements containing the evaluation
        of the function and its first derivative.

For the third order methods
([@http://en.wikipedia.org/wiki/Halley%27s_method Halley] and
Schroeder)
        the __tuple should have *three* elements containing the evaluation of
        the function and its first and second derivatives.]]
[[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.]]
[[uintmax_t& max_iter] [An optional maximum number of iterations to perform.  On exit this is set to the actual number of iterations performed.]]
]

When using these functions you should note that:

* Default max_iter = `(std::numeric_limits<boost::uintmax_t>::max)()` is effectively 'iterate for ever'!.
* They 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)
the bracket bounds min and max may as well be set to the widest limits
like zero and `numeric_limits<T>::max()`.
*But if the function more complex and may have more than one root or a pole,
the choice of bounds is protection against jumping out to seek the 'wrong' root.
* 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.  However, a good initial guess backed up by asymptotically-tight
bounds will improve performance no end - rather than relying on bisection.
* 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.
* To get the binary digits of accuracy, use 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.

[h4:newton Newton Raphson Method]
Given an initial guess x0 the subsequent values are computed using:

[equation roots1]

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

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

[h4:halley Halley's Method]

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

[equation roots2]

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.

[h4:schroeder Schroeder's Method]

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

[equation roots3]

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.

[h4 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:

[equation roots4]

To begin with lets solve the problem using Newton-Raphson iterations, we'll
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
   {
      cbrt_functor(T const& target) : a(target)
      { // Constructor stores value to be 'cube-rooted'.
      }
      ``__tuple``<T, T> operator()(T const& z)
      { // z is estimate so far.
         return boost::math::make_tuple(
         z*z*z - a, // return both f(x)
         3 * z*z);  // and f'(x)
      }
   private:
      T a; // to be 'cube-rooted'.
   };

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:

   template <class T>
   T cbrt(T z)
   {
      using namespace std; // for frexp, ldexp, numeric_limits.
      using namespace boost::math::tools;

      int exp;
      frexp(z, &exp); // Get exponent of z (ignore mantissa).
      T min = ldexp(0.5, exp/3);
      T max = ldexp(2.0, exp/3);
      T guess = ldexp(1.0, exp/3); // Rough guess is to divide the exponent by three.
      int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T.
      return newton_raphson_iterate(detail::cbrt_functor<T>(z), guess, min, max, digits);
   }

Using the test data in `libs/math/test/cbrt_test.cpp` this 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 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.

Note also that the above code omits a probably optimization by computing z[sup2],
and reusing it, omits error handling, and does not handle
negative values of z correctly. (These are left as an exercise for the reader!)

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

Now let's adapt the functor slightly to return the second derivative as well:

   template <class T>
   struct cbrt_functor
   {
      cbrt_functor(T const& target) : a(target){}
      ``__tuple``<T, T, T> operator()(T const& z)
      {
         return boost::math::make_tuple(
         z*z*z - a,
         3 * z*z,
         6 * z);
      }
   private:
      T a;
   };

And then adapt the `cbrt` function to use Halley iterations:

   template <class T>
   T cbrt(T z)
   {
      using namespace std;
      using namespace boost::math::tools;

      int exp;
      frexp(z, &exp);
      T min = ldexp(0.5, exp/3);
      T max = ldexp(2.0, exp/3);
      T guess = ldexp(1.0, exp/3);
      int digits = std::numeric_limits<T>::digits / 2;
      return halley_iterate(detail::cbrt_functor<T>(z), guess, min, max, digits);
   }

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

Finally, had we called `cbrt` with [@http://shoup.net/ntl/doc/RR.txt NTL::RR]
set to 1000 bit precision, then full precision 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!]

[endsect] [/section:roots Root Finding With Derivatives]

[/
  Copyright 2006, 2010, 2012 John Maddock and Paul A. Bristow.
  Distributed under the Boost Software License, Version 1.0.
  (See accompanying file LICENSE_1_0.txt or copy at
  http://www.boost.org/LICENSE_1_0.txt).
]