Complex Newton's method

doc/roots/roots.qbk

@@ -29,6 +29,9 @@
    template <class F, class T>
    T schroder_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter);
 
+   template <class F, class Complex>
+   Complex complex_newton(F f, Complex guess, int max_iterations = std::numeric_limits<typename Complex::value_type>::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,17 @@ 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.
+
+Finally, the derivative of /f/ must be continuous at the root or else non-roots can be found; see [@https://math.stackexchange.com/questions/3017766/constructing-newton-iteration-converging-to-non-root here] for an example.
+
+An example usage of `complex_newton` is given in `examples/daubechies_coefficients.cpp`.
+
 [h4 Examples]
 
 See __root_finding_examples.

example/daubechies_coefficients.cpp

@@ -16,6 +16,7 @@
 #include <boost/math/tools/roots.hpp>
 #include <boost/math/special_functions/binomial.hpp>
 #include <boost/multiprecision/cpp_complex.hpp>
+#include <boost/multiprecision/complex128.hpp>
 #include <boost/math/quadrature/gauss_kronrod.hpp>
 
 using std::string;
@@ -24,44 +25,14 @@ 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;
 using boost::math::quadrature::gauss_kronrod;
 using boost::multiprecision::cpp_bin_float_100;
 using boost::multiprecision::cpp_complex_100;
 
-template<class Complex, class F>
-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<Real>::epsilon()))
-       {
-          return xn;
-       }
-
-    } while(iterations);
-
-    if (abs(f) > std::numeric_limits<typename Complex::value_type>::epsilon())
-    {
-        return {std::numeric_limits<typename Complex::value_type>::quiet_NaN(), std::numeric_limits<typename Complex::value_type>::quiet_NaN()};
-    }
-    return xn;
-}
-
 template<class Complex>
 std::vector<std::pair<Complex, Complex>> find_roots(size_t p)
 {
@@ -76,6 +47,8 @@ std::vector<std::pair<Complex, Complex>> find_roots(size_t p)
 
     polynomial<Complex> P(std::move(coeffs));
     polynomial<Complex> Pcopy = P;
+    polynomial<Complex> Pcopy_prime = P.prime();
+    auto orig = [&](Complex z) { return std::make_pair<Complex, Complex>(Pcopy(z), Pcopy_prime(z)); };
 
     polynomial<Complex> P_prime = P.prime();
 
@@ -94,7 +67,33 @@ std::vector<std::pair<Complex, Complex>> find_roots(size_t p)
             return std::make_pair<Complex, Complex>(P(x), P_prime(x));
         };
 
-        Complex r = complex_newton(f, guess, search_radius, 100);
+        Complex r = complex_newton(f, guess);
+        using std::isnan;
+        if(isnan(r.real()))
+        {
+            int i = 50;
+            do {
+                // Try a different guess
+                guess *= Complex(1.0,-1.0);
+                r = complex_newton(f, guess);
+                std::cout << "New guess: " << guess << ", result? " << r << std::endl;
+
+            } while (isnan(r.real()) && i-- > 0);
+
+            if (isnan(r.real()))
+            {
+                std::cout << "Polynomial that killed the process: " << P << std::endl;
+                throw std::logic_error("Newton iteration did not converge");
+            }
+        }
+        // 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);
+        if (isnan(r.real()))
+        {
+            throw std::logic_error("Found a root for the deflated polynomial which is not a root for the original. Indicative of catastrophic numerical error.");
+        }
         // Test the root:
         using std::sqrt;
         Real tol = sqrt(sqrt(std::numeric_limits<Real>::epsilon()));
@@ -160,13 +159,13 @@ std::vector<typename Complex::value_type> daubechies_coefficients(std::vector<st
     std::vector<Complex> chosen_roots(p-1);
     for (size_t i = 0; i < p - 1; ++i)
     {
-        if(abs(Qroots[i].first) <= 1)
+        if(norm(Qroots[i].first) <= 1)
         {
             chosen_roots[i] = Qroots[i].first;
         }
         else
         {
-            BOOST_ASSERT(abs(Qroots[i].second) <= 1);
+            BOOST_ASSERT(norm(Qroots[i].second) <= 1);
             chosen_roots[i] = Qroots[i].second;
         }
     }
@@ -209,16 +208,15 @@ std::vector<typename Complex::value_type> daubechies_coefficients(std::vector<st
 
 int main()
 {
-    typedef cpp_complex_100 Complex;
-    // This is probably the most ill-conditioned algorithm I've ever written.
-    for(size_t p = 1; p < 57; ++p)
+    typedef boost::multiprecision::cpp_complex<100> Complex;
+    for(size_t p = 1; p < 200; ++p)
     {
         auto roots = find_roots<Complex>(p);
         auto h = daubechies_coefficients(roots);
         std::cout << "h_" << p << "[] = {";
         for (auto& x : h) {
             std::cout << x << ", ";
         }
-        std::cout << "}\n";
+        std::cout << "} // = h_" << p << "\n\n\n\n";
     }
 }

include/boost/math/tools/roots.hpp

@@ -9,7 +9,7 @@
 #ifdef _MSC_VER
 #pragma once
 #endif
-
+#include <iostream>
 #include <utility>
 #include <boost/config/no_tr1/cmath.hpp>
 #include <stdexcept>
@@ -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 <class T, class U, class V>
@@ -436,9 +436,9 @@ namespace detail{
          // check for out of bounds step:
          if(result < min)
          {
-            T diff = ((fabs(min) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(min))) 
-               ? T(1000) 
-               : (fabs(min) < 1) && (fabs(tools::max_value<T>() * min) < fabs(result)) 
+            T diff = ((fabs(min) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(min)))
+               ? T(1000)
+               : (fabs(min) < 1) && (fabs(tools::max_value<T>() * min) < fabs(result))
                ? ((min < 0) != (result < 0)) ? -tools::max_value<T>() : tools::max_value<T>() : T(result / min);
             if(fabs(diff) < 1)
                diff = 1 / diff;
@@ -563,10 +563,94 @@ 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 even in this somewhat pathological case.
+ * For isolated roots, the problem is so rapidly convergent that this doesn't matter at all.
+ */
+template<class Complex, class F>
+Complex complex_newton(F g, Complex guess, int max_iterations=std::numeric_limits<typename Complex::value_type>::digits)
+{
+    typedef typename Complex::value_type Real;
+    using std::norm;
+    using std::abs;
+    using std::max;
+    // z0, z1, and z2 cannot be the same, in case we immediately need to resort to Muller's Method:
+    Complex z0 = guess + Complex(1,0);
+    Complex z1 = guess + Complex(0,1);
+    Complex z2 = guess;
+
+    do {
+       auto pair = g(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<Complex>(1)+q;
+           Complex A = q*(pair.first - qp1*P1.first + q*P0.first);
+
+           Complex B = (static_cast<Complex>(2)*q+static_cast<Complex>(1))*pair.first - qp1*qp1*P1.first +q*q*P0.first;
+           Complex C = qp1*pair.first;
+           Complex rad = sqrt(B*B - static_cast<Complex>(4)*A*C);
+           Complex denom1 = B + rad;
+           Complex denom2 = B - rad;
+           Complex correction = (z1-z2)*static_cast<Complex>(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);
+       }
+
+       // See: https://math.stackexchange.com/questions/3017766/constructing-newton-iteration-converging-to-non-root
+       // If f' is continuous, then convergence of x_n -> x* implies f(x*) = 0.
+       // This condition approximates this convergence condition by requiring three consecutive iterates to be clustered.
+       Real tol = max(abs(z2)*std::numeric_limits<Real>::epsilon(), std::numeric_limits<Real>::epsilon());
+       bool real_close = abs(z0.real() - z1.real()) < tol && abs(z0.real() - z2.real()) < tol && abs(z1.real() - z2.real()) < tol;
+       bool imag_close = abs(z0.imag() - z1.imag()) < tol && abs(z0.imag() - z2.imag()) < tol && abs(z1.imag() - z2.imag()) < tol;
+       if (real_close && imag_close)
+       {
+           return z2;
+       }
+
+   } while(max_iterations--);
+
+    // 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 < eps
+    // This is somewhat awkward as it isn't scale invariant, but using the Daubechies coefficient example code,
+    // I found this condition generates correct roots, whereas the scale invariant condition discussed here:
+    // https://scicomp.stackexchange.com/questions/30597/defining-a-condition-number-and-termination-criteria-for-newtons-method
+    // allows nonroots to be passed off as roots.
+    auto pair = g(z2);
+    if (abs(pair.first) < sqrt(std::numeric_limits<Real>::epsilon()))
+    {
+        return z2;
+    }
+
+    return {std::numeric_limits<Real>::quiet_NaN(),
+            std::numeric_limits<Real>::quiet_NaN()};
+}
+#endif
 
 } // namespace tools
 } // namespace math
 } // namespace boost
 
 #endif // BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP
-