From 3ab69d00eefdf8afb1db13a4aecea3a9e1f8f97f Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Sun, 25 Nov 2018 00:19:11 -0700 Subject: [PATCH] Complex Newton's method. Zero derivatives handled by Muller's Method. --- doc/roots/roots.qbk | 17 +++++- example/daubechies_coefficients.cpp | 45 ++++------------ include/boost/math/tools/roots.hpp | 81 +++++++++++++++++++++++++++-- test/test_roots.cpp | 78 +++++++++++++++++++++++++-- 4 files changed, 177 insertions(+), 44 deletions(-) diff --git a/doc/roots/roots.qbk b/doc/roots/roots.qbk index 175a7d888..43ac9ef65 100644 --- a/doc/roots/roots.qbk +++ b/doc/roots/roots.qbk @@ -29,6 +29,9 @@ template T schroder_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter); + template + Complex complex_newton(F f, Complex guess, int max_iterations = std::numeric_limits::digits); + }}} // namespaces boost::math::tools. [h4 Description] @@ -40,9 +43,10 @@ These functions all perform iterative root-finding [*using derivatives]: * `halley_iterate` and `schroder_iterate` perform third-order __halley and __schroder iteration. -The functions all take the same parameters: +* `complex_newton` performs Newton's method on complex-analytic functions. -[variablelist Parameters of the root finding functions + +[variablelist Parameters of the real-valued root finding functions [[F f] [Type F must be a callable function object (or C++ lambda) that accepts one parameter and returns a __tuple_type: @@ -158,6 +162,15 @@ College Park, MD: University of Maryland, Institute for Advanced Computer Studie This method guarantees at least quadratic convergence (the same as Newton's method), and is known to work well in the presence of multiple roots: something that neither Newton nor Halley can do. +The complex Newton method works slightly differently than the rest of the methods: +Since there is no way to bracket roots in the complex plane, +the `min` and `max` arguments are not accepted. +Failure to reach a root is communicated by returning `nan`s. +Remember that if a function has many roots, +then which root the complex Newton's method converges to is essentially impossible to predict a priori; see the Newton's fractal for more information. + +An example usage of `complex_newton` is given in `examples/daubechies_coefficients.cpp`. + [h4 Examples] See __root_finding_examples. diff --git a/example/daubechies_coefficients.cpp b/example/daubechies_coefficients.cpp index 09915dd04..452687c6f 100644 --- a/example/daubechies_coefficients.cpp +++ b/example/daubechies_coefficients.cpp @@ -24,6 +24,7 @@ using boost::math::binomial_coefficient; using boost::math::tools::schroder_iterate; using boost::math::tools::halley_iterate; using boost::math::tools::newton_raphson_iterate; +using boost::math::tools::complex_newton; using boost::math::constants::half; using boost::math::constants::root_two; using boost::math::constants::pi; @@ -31,37 +32,6 @@ using boost::math::quadrature::gauss_kronrod; using boost::multiprecision::cpp_bin_float_100; using boost::multiprecision::cpp_complex_100; -template -Complex complex_newton(F g, Complex guess, typename Complex::value_type search_radius=100, int iterations=30) -{ - typedef typename Complex::value_type Real; - Complex xn = guess; - - Complex f; - Complex df; - do { - - auto pair = g(xn); - f = pair.first; - df = pair.second; - - --iterations; - - xn = xn - (f/df); - if (abs(f) < sqrt(std::numeric_limits::epsilon())) - { - return xn; - } - - } while(iterations); - - if (abs(f) > std::numeric_limits::epsilon()) - { - return {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}; - } - return xn; -} - template std::vector> find_roots(size_t p) { @@ -76,6 +46,8 @@ std::vector> find_roots(size_t p) polynomial P(std::move(coeffs)); polynomial Pcopy = P; + polynomial Pcopy_prime = P.prime(); + auto orig = [&](Complex z) { return std::make_pair(Pcopy(z), Pcopy_prime(z)); }; polynomial P_prime = P.prime(); @@ -94,7 +66,11 @@ std::vector> find_roots(size_t p) return std::make_pair(P(x), P_prime(x)); }; - Complex r = complex_newton(f, guess, search_radius, 100); + Complex r = complex_newton(f, guess); + // Refine r with the original function. + // We only use the polynomial division to ensure we don't get the same root over and over. + // However, the division induces error which can grow quickly-or slowly! See Numerical Recipes, section 9.5.1. + r = complex_newton(orig, r); // Test the root: using std::sqrt; Real tol = sqrt(sqrt(std::numeric_limits::epsilon())); @@ -210,8 +186,7 @@ std::vector daubechies_coefficients(std::vector(p); auto h = daubechies_coefficients(roots); @@ -219,6 +194,6 @@ int main() for (auto& x : h) { std::cout << x << ", "; } - std::cout << "}\n"; + std::cout << "}\n\n\n\n"; } } diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp index 81e846613..c64587b4a 100644 --- a/include/boost/math/tools/roots.hpp +++ b/include/boost/math/tools/roots.hpp @@ -9,7 +9,7 @@ #ifdef _MSC_VER #pragma once #endif - +#include #include #include #include @@ -65,7 +65,7 @@ inline void unpack_0(const Tuple& t, T& val) BOOST_MATH_NOEXCEPT(T) { using dummy::get; // Rely on ADL to find the correct overload of get: - val = get<0>(t); + val = get<0>(t); } template @@ -559,10 +559,85 @@ inline T schroeder_iterate(F f, T guess, T min, T max, int digits) BOOST_NOEXCEP return schroder_iterate(f, guess, min, max, digits, m); } +#ifndef BOOST_NO_CXX11_AUTO_DECLARATIONS +/* Why do we set the default maximum number of iterations to the number of digits in the type? + * Because for double roots, the number of digits increases linearly with the number of iterations, so this default should recover full precision. + * For isolated roots, the problem is so rapidly convergent that this doesn't matter at all. + * For nonexistence, the problem is harder. + */ +template +Complex complex_newton(F g, Complex guess, int max_iterations=std::numeric_limits::digits) +{ + typedef typename Complex::value_type Real; + using std::norm; + Complex z0 = guess/static_cast(4); + Complex z1 = guess/static_cast(2); + Complex z2 = guess; + do { + auto pair = g(z2); + if (abs(pair.first) < 10*std::numeric_limits::epsilon()) + { + // Computed it, might as well use it: + if (abs(pair.second) > sqrt(std::numeric_limits::epsilon())) + { + return z2 - pair.first/pair.second; + } + return z2; + } + if (norm(pair.second) == 0) + { + // Muller's method. Notation follows Numerical Recipes, 9.5.2: + Complex q = (z2 - z1)/(z1 - z0); + auto P0 = g(z0); + auto P1 = g(z1); + Complex qp1 = static_cast(1)+q; + Complex A = q*(pair.first - qp1*P1.first + q*P0.first); + + Complex B = (static_cast(2)*q+static_cast(1))*pair.first - qp1*qp1*P1.first +q*q*P0.first; + Complex C = qp1*pair.first; + Complex rad = sqrt(B*B - static_cast(4)*A*C); + Complex denom1 = B + rad; + Complex denom2 = B - rad; + Complex correction = (z1-z2)*static_cast(2)*C; + if (norm(denom1) > norm(denom2)) + { + correction /= denom1; + } else { + correction /= denom2; + } + + z0 = z1; + z1 = z2; + z2 = z2 + correction; + } + else + { + z0 = z1; + z1 = z2; + z2 = z2 - (pair.first/pair.second); + } + + } while(max_iterations--); + + // At this point, maybe we should use complex trapezoidal quadrature to see if there are any roots? + // https://www.ams.org/journals/mcom/1967-21-100/S0025-5718-1967-0228165-4/S0025-5718-1967-0228165-4.pdf + // Exposition of the above: http://www.chebfun.org/examples/roots/ComplexRoots.html + + // The idea is that if we can get abs(f) < eps, we should, but if we go through all these iterations + // and abs(f) < sqrt(eps), then roundoff error simply does not allow that we can evaluate f to < esp + auto pair = g(z2); + if (abs(pair.first) < sqrt(std::numeric_limits::epsilon())) + { + return z2; + } + + return {std::numeric_limits::quiet_NaN(), + std::numeric_limits::quiet_NaN()}; +} +#endif } // namespace tools } // namespace math } // namespace boost #endif // BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP - diff --git a/test/test_roots.cpp b/test/test_roots.cpp index c05392377..403055634 100644 --- a/test/test_roots.cpp +++ b/test/test_roots.cpp @@ -11,14 +11,23 @@ #include #include #include +#include #include +#include #include #include #include +#include #include "table_type.hpp" #include #include +#if BOOST_HAS_FLOAT128 +#include +#endif + +#include + #define BOOST_CHECK_CLOSE_EX(a, b, prec, i) \ {\ unsigned int failures = boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed;\ @@ -317,13 +326,74 @@ void test_beta(T, const char* /* name */) test_inverses(ibeta_large_data); } +#ifndef BOOST_NO_CXX11_AUTO_DECLARATIONS +template +void test_complex_newton() +{ + typedef typename Complex::value_type Real; + std::cout << "Testing complex Newton's Method on type " << boost::typeindex::type_id().pretty_name() << "\n"; + using boost::math::tools::complex_newton; + using boost::math::tools::polynomial; + using boost::math::constants::half; + + Real tol = std::numeric_limits::epsilon(); + // p(z) = z^2 + 1, roots: \pm i. + polynomial p{{1,0}, {0, 0}, {1,0}}; + Complex guess{1,1}; + polynomial p_prime = p.prime(); + auto f = [&](Complex z) { return std::make_pair(p(z), p_prime(z)); }; + Complex root = complex_newton(f, guess); + + BOOST_CHECK(abs(root.real()) <= tol); + BOOST_CHECK_CLOSE(root.imag(), 1, tol); + + guess = -guess; + root = complex_newton(f, guess); + BOOST_CHECK(abs(root.real()) <= tol); + BOOST_CHECK_CLOSE(root.imag(), -1, tol); + + // Test that double roots are handled correctly-as correctly as possible. + // Convergence at a double root is not quadratic. + // This sets p = (z-i)^2: + p = polynomial({{-1,0}, {0,-2}, {1,0}}); + p_prime = p.prime(); + guess = -guess; + auto g = [&](Complex z) { return std::make_pair(p(z), p_prime(z)); }; + root = complex_newton(g, guess); + BOOST_CHECK(abs(root.real()) < 10*sqrt(tol)); + BOOST_CHECK_CLOSE(root.imag(), 1, tol); + + // Test that zero derivatives are handled. + // p(z) = z^2 + iz + 1 + p = polynomial({{1,0}, {0,1}, {1,0}}); + // p'(z) = 2z + i + p_prime = p.prime(); + guess = Complex(0,-boost::math::constants::half()); + auto g2 = [&](Complex z) { return std::make_pair(p(z), p_prime(z)); }; + root = complex_newton(g2, guess); + + // Here's the other root, in case code changes cause it to be found: + //Complex expected_root1{0, half()*(sqrt(static_cast(5)) - static_cast(1))}; + Complex expected_root2{0, -half()*(sqrt(static_cast(5)) + static_cast(1))}; + + BOOST_CHECK_CLOSE(expected_root2.imag(),root.imag(), tol); + BOOST_CHECK(abs(root.real()) < tol); +} +#endif + BOOST_AUTO_TEST_CASE( test_main ) { + test_beta(0.1, "double"); - // bug reports: - boost::math::skew_normal_distribution<> dist(2.0, 1.0, -2.5); - BOOST_CHECK(boost::math::isfinite(quantile(dist, 0.075))); - +#ifndef BOOST_NO_CXX11_AUTO_DECLARATIONS + test_complex_newton>(); + test_complex_newton>(); + test_complex_newton>(); + test_complex_newton(); +#if BOOST_HAS_FLOAT128 + test_complex_newton(); +#endif +#endif }