WSJT-X/boost/libs/math/doc/roots/root_finding_examples.qbk

[section:root_finding_examples Examples of Root-Finding (with and without derivatives)]

[import ../../example/root_finding_example.cpp]
[import ../../example/root_finding_n_example.cpp]
[import ../../example/root_finding_multiprecision_example.cpp]

The examples demonstrate how to use the various tools for
[@http://en.wikipedia.org/wiki/Root-finding_algorithm root finding].

We start with the simple cube root function `cbrt` ( C++ standard function name
[@http://en.cppreference.com/w/cpp/numeric/math/cbrt cbrt])
showing root finding __cbrt_no_derivatives.

We then show how use of derivatives can improve the speed of convergence.

(But these examples are only a demonstration and do not try to make
the ultimate improvements of an 'industrial-strength'
implementation, for example, of `boost::math::cbrt`, mainly by using a better computed initial 'guess'
at [@boost:/libs/math/include/boost/math/special_functions/cbrt.hpp  cbrt.hpp]).

Then we show how a higher root (__fifth_root) [super 5][radic] can be computed,
and in
[@../../example/root_finding_n_example.cpp root_finding_n_example.cpp]
a generic method for the __nth_root that constructs the derivatives at compile-time.

These methods should be applicable to other functions that can be differentiated easily.

[section:cbrt_eg Finding the Cubed Root With and Without Derivatives]

First some `#includes` that will be needed.

[root_finding_include_1]

[tip For clarity, `using` statements are provided to list what functions are being used in this example:
you can, of course, partly or fully qualify the names in other ways.
(For your application, you may wish to extract some parts into header files,
but you should never use `using` statements globally in header files).]

Let's suppose we want to find the root of a number ['a], and to start, compute the cube root.

So the equation we want to solve is:

__spaces ['f(x) = x[cubed] -a]

We will first solve this without using any information
about the slope or curvature of the cube root function.

Fortunately, the cube-root function is 'Really Well Behaved' in that it is monotonic
and has only one root (we leave negative values 'as an exercise for the student').

We then show how adding what we can know about this function, first just the slope
or 1st derivative ['f'(x)], will speed homing in on the solution.

Lastly, we show how adding the curvature ['f''(x)] too will speed convergence even more.

[h3:cbrt_no_derivatives Cube root function without derivatives]

First we define a function object (functor):

[root_finding_noderiv_1]

Implementing the cube-root function itself 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'.)

[root_finding_noderiv_2]

This snippet from `main()` in [@../../example/root_finding_example.cpp root_finding_example.cpp]
shows how it can be used.

[root_finding_main_1]

[pre
  cbrt_noderiv(27) = 3
  cbrt_noderiv(28) = 3.0365889718756618
]

The result of `bracket_and_solve_root` is a [@http://www.cplusplus.com/reference/utility/pair/ pair]
of values that could be displayed.

The number of bits separating them can be found using `float_distance(r.first, r.second)`.
The distance is zero (closest representable) for 3[super 3] = 27
but `float_distance(r.first, r.second) = 3` for cube root of 28 with this function.
The result (avoiding overflow) is midway between these two values.

[h3:cbrt_1st_derivative Cube root function with 1st derivative (slope)]

We now solve the same problem, but using more information about the function,
to show how this can speed up finding the best estimate of the root.

For the root function, the 1st differential (the slope of the tangent to a curve at any point) is known.

This algorithm is similar to this [@http://en.wikipedia.org/wiki/Nth_root_algorithm nth root algorithm].

If you need some reminders, then
[@http://en.wikipedia.org/wiki/Derivative#Derivatives_of_elementary_functions derivatives of elementary functions]
may help.

Using the rule that the derivative of ['x[super n]] for positive n (actually all nonzero n) is ['n x[super n-1]],
allows us to get the 1st differential as ['3x[super 2]].

To see how this extra information is used to find a root, view
[@http://en.wikipedia.org/wiki/Newton%27s_method Newton-Raphson iterations]
and the [@http://en.wikipedia.org/wiki/Newton%27s_method#mediaviewer/File:NewtonIteration_Ani.gif animation].

We define a better functor `cbrt_functor_deriv` that returns
both the evaluation of the function to solve, along with its first derivative:

To '['return]' two values, we use a [@http://en.cppreference.com/w/cpp/utility/pair std::pair]
of floating-point values.

[root_finding_1_deriv_1]

The result of [@boost:/libs/math/include/boost/math/tools/roots.hpp `newton_raphson_iterate`]
function is a single value.

[tip There is a compromise between accuracy and speed when chosing the value of `digits`.
It is tempting to simply chose `std::numeric_limits<T>::digits`,
but this may mean some inefficient and unnecessary iterations as the function thrashes around
trying to locate the last bit.  In theory, since the precision doubles with each step
it is sufficient to stop when half the bits are correct: as the last step will have doubled
that to full precision.  Of course the function has no way to tell if that is actually the case
unless it does one more step to be sure.  In practice setting the precision to slightly more
than `std::numeric_limits<T>::digits / 2` is a good choice.]

Note that it is up to the caller of the function to check the iteration count
after the call to see if iteration stoped as a result of running out of iterations
rather than meeting the required precision.

Using the test data in [@../../test/test_cbrt.cpp  /test/test_cbrt.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 probable optimization by computing z[sup2]
and reusing it, omits error handling, and does not handle
negative values of z correctly. (These are left as the customary exercise for the reader!)

The `boost::math::cbrt` function also includes these and other improvements:
most importantly it uses a much better initial guess which reduces the iteration count to
just 1 in almost all cases.

[h3:cbrt_2_derivatives Cube root with 1st & 2nd derivative (slope & curvature)]

Next we define yet another even better functor `cbrt_functor_2deriv` that returns
both the evaluation of the function to solve,
along with its first [*and second] derivative:

__spaces['f''(x) = 6x]

using information about both slope and curvature to speed convergence.

To [''return'] three values, we use a `tuple` of three floating-point values:
[root_finding_2deriv_1]

The function `halley_iterate` also returns a single value,
and the number of iterations will reveal if it met the convergence criterion set by `get_digits`.

The no-derivative method gives a result of

  cbrt_noderiv(28) = 3.0365889718756618

with a 3 bits distance between the bracketed values, whereas the derivative methods both converge to a single value

  cbrt_2deriv(28) = 3.0365889718756627

which we can compare with the [@boost:/libs/math/doc/html/math_toolkit/powers/cbrt.html boost::math::cbrt]

  cbrt(28)              = 3.0365889718756627

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
[link math_toolkit.roots.roots_deriv.schroder `schroder_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 (about 300 decimal digits),
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!]

Full code of this example is at
[@../../example/root_finding_example.cpp  root_finding_example.cpp],

[endsect]

[section:lambda Using C++11 Lambda's]

Since all the root finding functions accept a function-object, they can be made to
work (often in a lot less code) with C++11 lambda's.  Here's the much reduced code for our "toy" cube root function:

[root_finding_2deriv_lambda]

Full code of this example is at
[@../../example/root_finding_example.cpp  root_finding_example.cpp],

[endsect]

[section:5th_root_eg Computing the Fifth Root]

Let's now suppose we want to find the [*fifth root] of a number ['a].

The equation we want to solve is :

__spaces['f](x) = ['x[super 5] -a]

If your differentiation is a little rusty
(or you are faced with an function whose complexity makes differentiation daunting),
then you can get help, for example, from the invaluable
[@http://www.wolframalpha.com/ WolframAlpha site.]

For example, entering the commmand: `differentiate x ^ 5`

or the Wolfram Language command: ` D[x ^ 5, x]`

gives the output: `d/dx(x ^ 5) = 5 x ^ 4`

and to get the second differential, enter: `second differentiate x ^ 5`

or the Wolfram Language command: `D[x ^ 5, { x, 2 }]`

to get the output: `d ^ 2 / dx ^ 2(x ^ 5) = 20 x ^ 3`

To get a reference value, we can enter: [^fifth root 3126]

or: `N[3126 ^ (1 / 5), 50]`

to get a result with a precision of 50 decimal digits:

5.0003199590478625588206333405631053401128722314376

(We could also get a reference value using __multiprecision_root).

The 1st and 2nd derivatives of x[super 5] are:

__spaces['f]\'(x) = 5x[super 4]

__spaces['f]\'\'(x) = 20x[super 3]

[root_finding_fifth_functor_2deriv]
[root_finding_fifth_2deriv]

Full code of this example is at
[@../../example/root_finding_example.cpp  root_finding_example.cpp] and
[@../../example/root_finding_n_example.cpp  root_finding_n_example.cpp].

[endsect]

[section:multiprecision_root  Root-finding using Boost.Multiprecision]

The apocryphally astute reader might, by now, be asking "How do we know if this computes the 'right' answer?".

For most values, there is, sadly, no 'right' answer.
This is because values can only rarely be ['exactly represented] by C++ floating-point types.
What we do want is the 'best' representation - one that is the nearest __representable value.
(For more about how numbers are represented see __floating_point).

Of course, we might start with finding an external reference source like
__WolframAlpha, as above, but this is not always possible.

Another way to reassure is to compute 'reference' values at higher precision
with which to compare the results of our iterative computations using built-in like `double`.
They should agree within the tolerance that was set.

The result  of `static_cast`ing to `double` from a higher-precision type like `cpp_bin_float_50` is guaranteed
to be the [*nearest representable] `double` value.

For example, the cube root functions in our example for `cbrt(28.)` compute

`std::cbrt<double>(28.)  = 3.0365889718756627`

WolframAlpha says  `3.036588971875662519420809578505669635581453977248111123242141...`

`static_cast<double>(3.03658897187566251942080957850) = 3.0365889718756627`

This example `cbrt(28.)   =  3.0365889718756627`

[tip To ensure that all potentially significant decimal digits are displayed use `std::numeric_limits<T>::max_digits10`
(or if not available on older platforms or compilers use `2+std::numeric_limits<double>::digits*3010/10000`).[br]

Ideally, values should agree to `std::numeric-limits<T>::digits10` decimal digits.

This also means that a 'reference' value to be [*input] or `static_cast` should have
at least `max_digits10` decimal digits (17 for 64-bit `double`).
]

If we wish to compute [*higher-precision values] then, on some platforms, we may be able to use `long double`
with a higher precision than `double` to compare with the very common `double`
and/or a more efficient built-in quad floating-point type like `__float128`.

Almost all platforms can easily use __multiprecision,
for example, __cpp_dec_float or a binary type __cpp_bin_float types,
to compute values at very much higher precision.

[note With multiprecision types, it is debatable whether to use the type `T` for computing the initial guesses.
Type `double` is like to be accurate enough for the method used in these examples.
This would limit the exponent range of possible values to that of `double`.
There is also the cost of conversion to and from type `T` to consider.
In these examples, `double` is used via `typedef double guess_type`.]

Since the functors and functions used above are templated on the value type,
we can very simply use them with any of the __multiprecision types.  As a reminder,
here's our toy cube root function using 2 derivatives and C++11 lambda functions to find the root:

[root_finding_2deriv_lambda]

Some examples below are 50 decimal digit decimal and binary types
(and on some platforms a much faster `float128` or `quad_float` type )
that we can use with these includes:

[root_finding_multiprecision_include_1]

Some using statements simplify their use:

[root_finding_multiprecision_example_1]

They can be used thus:

[root_finding_multiprecision_example_2]

A reference value computed by __WolframAlpha is

   N[2^(1/3), 50]  1.2599210498948731647672106072782283505702514647015

which agrees exactly.

To [*show] values to their full precision, it is necessary to adjust the `std::ostream` `precision` to suit the type, for example:

[root_finding_multiprecision_show_1]

[root_finding_multiprecision_example_3]

which outputs:

[pre
cbrt(2) = 1.2599210498948731647672106072782283505702514647015

value = 2, cube root =1.25992104989487
value = 2, cube root =1.25992104989487
value = 2, cube root =1.2599210498948731647672106072782283505702514647015
]

[tip Be [*very careful] about the floating-point type `T` that is passed to the root-finding function.
Carelessly passing a integer by writing
`cpp_dec_float_50  r = cbrt_2deriv(2);` or `show_cube_root(2);`
will provoke many warnings and compile errors.

Even `show_cube_root(2.F);` will produce warnings because `typedef double guess_type` defines the type
used to compute the guess and bracket values as `double`.

Even more treacherous is passing a `double` as in `cpp_dec_float_50  r = cbrt_2deriv(2.);`
which silently gives the 'wrong' result, computing a `double` result and [*then] converting to `cpp_dec_float_50`!
All digits beyond `max_digits10` will be incorrect.
Making the `cbrt` type explicit with `cbrt_2deriv<cpp_dec_float_50>(2.);` will give you the desired 50 decimal digit precision result.
] [/tip]

Full code of this example is at
[@../../example/root_finding_multiprecision_example.cpp  root_finding_multiprecision_example.cpp].

[endsect]

[section:nth_root Generalizing to Compute the nth root]

If desired, we can now further generalize to compute the ['n]th root by computing the derivatives [*at compile-time]
using the rules for differentiation and `boost::math::pow<N>`
where template parameter `N` is an integer and a compile time constant.  Our functor and function now have an additional template parameter `N`,
for the root required.

[note Since the powers and derivatives are fixed at compile time, the resulting code is as efficient as as if hand-coded as the cube and fifth-root examples above.
A good compiler should also optimise any repeated multiplications.]

Our ['n]th root functor is

[root_finding_nth_functor_2deriv]

and our ['n]th root  function is

[root_finding_nth_function_2deriv]

[root_finding_n_example_2]

produces an output similar to this

[root_finding_example_output_1]

[tip Take care with the type passed to the function. It is best to pass a `double` or greater-precision floating-point type.

Passing an integer value, for example, `nth_2deriv<5>(2)` will be rejected, while `nth_2deriv<5, double>(2)` converts the integer to `double`.

Avoid passing a `float` value that will provoke warnings (actually spurious) from the compiler about potential loss of data,
as noted above.]

[warning Asking for unreasonable roots, for example, `show_nth_root<1000000>(2.);` may lead to
[@http://en.wikipedia.org/wiki/Loss_of_significance Loss of significance] like
`Type double value = 2, 1000000th root = 1.00000069314783`.
Use of the the `pow` function is more sensible for this unusual need.
]

Full code of this example is at
[@../../example/root_finding_n_example.cpp  root_finding_n_example.cpp].

[endsect]

[section:elliptic_eg A More complex example - Inverting the Elliptic Integrals]

The arc length of an ellipse with radii ['a] and ['b] is given by:

[pre L(a, b) = 4aE(k)]

with:

[pre k = [sqrt](1 - b[super 2]/a[super 2])]

where ['E(k)] is the complete elliptic integral of the second kind - see __ellint_2.

Let's suppose we know the arc length and one radii, we can then calculate the other
radius by inverting the formula above.  We'll begin by encoding the above formula
into a functor that our root-finding algorithms can call.

Note that while not
completely obvious from the formula above, the function is completely symmetrical
in the two radii - which can be interchanged at will - in this case we need to
make sure that `a >= b` so that we don't accidentally take the square root of a negative number:

[import ../../example/root_elliptic_finding.cpp]

[elliptic_noderv_func]

We'll also need a decent estimate to start searching from, the approximation:

[pre L(a, b) [approx] 4[sqrt](a[super 2] + b[super 2])]

Is easily inverted to give us what we need, which using derivative-free root
finding leads to the algorithm:

[elliptic_root_noderiv]

This function generally finds the root within 8-10 iterations, so given that the runtime
is completely dominated by the cost of calling the ellliptic integral it would be nice to
reduce that count somewhat. We'll try to do that by using a derivative-based method;
the derivatives of this function are rather hard to work out by hand, but fortunately
[@http://www.wolframalpha.com/input/?i=d%2Fda+\[4+*+a+*+EllipticE%281+-+b^2%2Fa^2%29\]
Wolfram Alpha] can do the grunt work for us to give:

[pre d/da L(a, b) = 4(a[super 2]E(k) - b[super 2]K(k)) / (a[super 2] - b[super 2])]

Note that now we have [*two] elliptic integral calls to get the derivative, so our
functor will be at least twice as expensive to call as the derivative-free one above:
we'll have to reduce the iteration count quite substantially to make a difference!

Here's the revised functor:

[elliptic_1deriv_func]

The root-finding code is now almost the same as before, but we'll make use of
Newton-iteration to get the result:

[elliptic_1deriv]

The number of iterations required for `double` precision is now usually around 4 -
so we've slightly more than halved the number of iterations, but made the
functor twice as expensive to call!

Interestingly though, the second derivative requires no more expensive
elliptic integral calls than the first does, in other words it comes
essentially "for free", in which case we might as well make use of it
and use Halley-iteration.  This is quite a typical situation when
inverting special-functions.  Here's the revised functor:

[elliptic_2deriv_func]

The actual root-finding code is almost the same as before, except we can
use Halley, rather than Newton iteration:

[elliptic_2deriv]

While this function uses only slightly fewer iterations (typically around 3)
to find the root, compared to the original derivative-free method, we've moved from
8-10 elliptic integral calls to 6.

Full code of this example is at
[@../../example/root_elliptic_finding.cpp  root_elliptic_finding.cpp].

[endsect]


[endsect] [/section:root_examples Examples of Root Finding (with and without derivatives)]

[section:bad_guess The Effect of a Poor Initial Guess]

It's instructive to take our "toy" example algorithms, and use deliberately bad initial guesses to see how the
various root finding algorithms fair.  We'll start with the cubed root, and using the cube root of 500 as the test case:

[table
[[Initial Guess=][-500% ([approx]1.323)][-100% ([approx]3.97)][-50% ([approx]3.96)][-20% ([approx]6.35)][-10% ([approx]7.14)][-5% ([approx]7.54)][5% ([approx]8.33)][10% ([approx]8.73)][20% ([approx]9.52)][50% ([approx]11.91)][100% ([approx]15.87)][500 ([approx]47.6)]]
[[bracket_and_solve_root][12][8][8][10][11][11][11][11][11][11][7][13]]
[[newton_iterate][12][7][7][5][5][4][4][5][5][6][7][9]]
[[halley_iterate][7][4][4][3][3][3][3][3][3][4][4][6]]
[[schroder_iterate][11][6][6][4][3][3][3][3][4][5][5][8]]
]

As you can see `bracket_and_solve_root` is relatively insensitive to starting location - as long as you don't start many orders of magnitude away from the root it will
take roughly the same number of steps to bracket the root and solve it.  On the other hand the derivative-based methods are slow to start, but once they have some digits
correct they increase precision exceptionally fast: they are therefore quite sensitive to the initial starting location.

The next table shows the number of iterations required to find the second radius of an ellipse with first radius 50 and arc-length 500:

[table
[[Initial Guess=][-500% ([approx]20.6)][-100% ([approx]61.81)][-50% ([approx]61.81)][-20% ([approx]98.9)][-10% ([approx]111.3)][-5% ([approx]117.4)][5% ([approx]129.8)][10% ([approx]136)][20% ([approx]148.3)][50% ([approx]185.4)][100% ([approx]247.2)][500 ([approx]741.7)]]
[[bracket_and_solve_root][11][5][5][8][8][7][7][8][9][8][6][10]]
[[newton_iterate][4][4][4][3][3][3][3][3][3][4][4][4]]
[[halley_iterate][4][3][3][3][3][2][2][3][3][3][3][3]]
[[schroder_iterate][4][3][3][3][3][2][2][3][3][3][3][3]]
]

Interestingly this function is much more resistant to a poor initial guess when using derivatives.

[endsect]

[section:bad_roots Examples Where Root Finding Goes Wrong]

There are many reasons why root root finding can fail, here are just a few of the more common examples:

[h3 Local Minima]

If you start in the wrong place, such as z[sub 0] here:

[$../roots/bad_root_1.svg]

Then almost any root-finding algorithm will descend into a local minima rather than find the root.

[h3 Flatlining]

In this example, we're starting from a location (z[sub 0]) where the first derivative is essentially zero:

[$../roots/bad_root_2.svg]

In this situation the next iteration will shoot off to infinity (assuming we're using derivatives that is).  Our
code guards against this by insisting that the root is always bracketed, and then never stepping outside those bounds.
In a case like this, no root finding algorithm can do better than bisecting until the root is found.

Note that there is no scale on the graph, we have seen examples of this situation occur in practice ['even when
several decimal places of the initial guess z[sub 0] are correct.]

This is really a special case of a more common situation where root finding with derivatives is ['divergent].  Consider
starting at z[sub 0] in this case:

[$../roots/bad_root_4.svg]

An initial Newton step would take you further from the root than you started, as will all subsequent steps.

[h3 Micro-stepping / Non-convergence]

Consider starting at z[sub 0] in this situation:

[$../roots/bad_root_3.svg]

The first derivative is essentially infinite, and the second close to zero (and so offers no correction if we use it),
as a result we take a very small first step.  In the worst case situation, the first step is so small
- perhaps even so small that subtracting from z[sub 0] has no effect at the current working precision - that our algorithm
will assume we are at the root already and terminate.  Otherwise we will take lot's of very small steps which never converge
on the root: our algorithms will protect against that by reverting to bisection.

An example of this situation would be trying to find the root of e[super -1/z[super 2]] - this function has a single
root at ['z = 0], but for ['z[sub 0] < 0] neither Newton nor Halley steps will ever converge on the root, and for ['z[sub 0] > 0]
the steps are actually divergent.

[endsect]

[/
  Copyright 2015 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).
]