RootFinding: more robust hybrid newton-raphson

src/Numerics/Numerics.csproj

@@ -110,7 +110,8 @@
     <Compile Include="LinearAlgebra\Generic\Matrix.BCL.cs" />
     <Compile Include="LinearAlgebra\Generic\Vector.BCL.cs" />
     <Compile Include="NonConvergenceException.cs" />
-    <Compile Include="RootFinding\Algorithms\NewtonRaphson.cs" />
+    <Compile Include="RootFinding\Algorithms\HybridNewtonRaphson.cs" />
+    <Compile Include="RootFinding\Algorithms\ZeroCrossingBracketing.cs" />
     <Compile Include="RootFinding\Bracketing.cs" />
     <Compile Include="RootFinding\Algorithms\Brent.cs" />
     <Compile Include="RootFinding\FloatingPointRoots.cs" />

src/Numerics/RootFinding/Algorithms/NewtonRaphson.cs → src/Numerics/RootFinding/Algorithms/HybridNewtonRaphson.cs

@@ -32,20 +32,40 @@
 
 namespace MathNet.Numerics.RootFinding.Algorithms
 {
-    public static class NewtonRaphson
+    public static class HybridNewtonRaphson
     {
         /// <summary>Find a solution of the equation f(x)=0.</summary>
-        /// <remarks>Hybrid Newton-Raphson that when failing falls back to bisection.</remarks>
+        /// <remarks>Hybrid Newton-Raphson that falls back to bisection when overshooting or converging too slow, or to subdivision on lacking bracketing.</remarks>
         /// <exception cref="NonConvergenceException"></exception>
-        public static double FindRoot(Func<double, double> f, Func<double, double> df, double lowerBound, double upperBound, double accuracy = 1e-8, int maxIterations = 100)
+        public static double FindSingleRoot(Func<double, double> f, Func<double, double> df, double lowerBound, double upperBound, double accuracy, int maxIterations, int subdivision)
+        {
+            double root;
+            if (TryFindSingleRoot(f, df, lowerBound, upperBound, accuracy, maxIterations, subdivision, out root))
+            {
+                return root;
+            }
+            throw new NonConvergenceException("The algorithm has exceeded the number of iterations allowed");
+        }
+
+        /// <summary>Find a solution of the equation f(x)=0.</summary>
+        /// <remarks>Hybrid Newton-Raphson that falls back to bisection when overshooting or converging too slow, or to subdivision on lacking bracketing.</remarks>
+        public static bool TryFindSingleRoot(Func<double, double> f, Func<double, double> df, double lowerBound, double upperBound, double accuracy, int maxIterations, int subdivision, out double root)
         {
             double fmin = f(lowerBound);
             double fmax = f(upperBound);
 
-            if (fmin == 0.0) return lowerBound;
-            if (fmax == 0.0) return upperBound;
+            if (fmin == 0.0)
+            {
+                root = lowerBound;
+                return true;
+            }
+            if (fmax == 0.0)
+            {
+                root = upperBound;
+                return true;
+            }
 
-            double root = 0.5*(lowerBound + upperBound);
+            root = 0.5*(lowerBound + upperBound);
             double fx = f(root);
             double lastStep = Math.Abs(upperBound - lowerBound);
             for (int i = 0; i < maxIterations; i++)
@@ -56,35 +76,54 @@ public static double FindRoot(Func<double, double> f, Func<double, double> df, d
                 double step = fx/dfx;
                 root -= step;
 
-                if (root < lowerBound || root > upperBound || Math.Abs(2*fx) > Math.Abs(lastStep*dfx))
+                bool overshoot = root > upperBound, undershoot = root < lowerBound;
+                if (overshoot || undershoot || Math.Abs(2*fx) > Math.Abs(lastStep*dfx))
                 {
-                    // Newton-Raphson step failed ->  bisect instead
+                    // Newton-Raphson step failed
+
+                    // If same signs, try subdivision to scan for zero crossing intervals
+                    if (Math.Sign(fmin) == Math.Sign(fmax) && TryScanForCrossingsWithRoots(f, df, lowerBound, upperBound, accuracy, maxIterations - i - 1, subdivision, out root))
+                    {
+                        return true;
+                    }
+
+                    // Bisection
                     root = 0.5*(upperBound + lowerBound);
                     fx = f(root);
                     lastStep = 0.5*Math.Abs(upperBound - lowerBound);
                     if (Math.Sign(fx) == Math.Sign(fmin))
                     {
                         lowerBound = root;
                         fmin = fx;
+                        if (overshoot)
+                        {
+                            root = upperBound;
+                            fx = fmax;
+                        }
                     }
                     else
                     {
                         upperBound = root;
                         fmax = fx;
+                        if (undershoot)
+                        {
+                            root = lowerBound;
+                            fx = fmin;
+                        }
                     }
                     continue;
                 }
-                
+
                 if (Math.Abs(step) < accuracy)
                 {
-                    return root;
+                    return true;
                 }
 
                 // Evaluation
                 fx = f(root);
                 lastStep = step;
 
-                // Update Bounds
+                // Update bounds
                 if (Math.Sign(fx) != Math.Sign(fmin))
                 {
                     upperBound = root;
@@ -97,11 +136,26 @@ public static double FindRoot(Func<double, double> f, Func<double, double> df, d
                 }
                 else if (Math.Sign(fmin) != Math.Sign(fmax))
                 {
-                    return root;
+                    return true;
                 }
             }
 
-            throw new NonConvergenceException("The algorithm has exceeded the number of iterations allowed");
+            return false;
+        }
+
+        static bool TryScanForCrossingsWithRoots(Func<double, double> f, Func<double, double> df, double lowerBound, double upperBound, double accuracy, int maxIterations, int subdivision, out double root)
+        {
+            var zeroCrossings = ZeroCrossingBracketing.FindIntervalsWithin(f, lowerBound, upperBound, subdivision);
+            foreach (Tuple<double, double> bounds in zeroCrossings)
+            {
+                if (TryFindSingleRoot(f, df, bounds.Item1, bounds.Item2, accuracy, maxIterations, subdivision, out root))
+                {
+                    return true;
+                }
+            }
+
+            root = double.NaN;
+            return false;
         }
     }
 }

src/Numerics/RootFinding/Algorithms/ZeroCrossingBracketing.cs

@@ -0,0 +1,44 @@
+using System;
+using System.Collections.Generic;
+
+namespace MathNet.Numerics.RootFinding.Algorithms
+{
+    public static class ZeroCrossingBracketing
+    {
+        public static IEnumerable<Tuple<double, double>> FindIntervalsWithin(Func<double, double> f, double lowerBound, double upperBound, int parts)
+        {
+            // TODO: Consider binary-style search instead of linear scan
+
+            var fmin = f(lowerBound);
+            var fmax = f(upperBound);
+
+            if (Math.Sign(fmin) != Math.Sign(fmax))
+            {
+                yield return new Tuple<double, double>(lowerBound, upperBound);
+                yield break;
+            }
+
+            double subdiv = (upperBound - lowerBound)/parts;
+            var smin = lowerBound;
+            int sign = Math.Sign(fmin);
+
+            for (int k = 0; k < parts; k++)
+            {
+                var smax = smin + subdiv;
+                var sfmax = f(smax);
+                if (double.IsInfinity(sfmax))
+                {
+                    // expand interval to include pole
+                    smin = smax;
+                    continue;
+                }
+                if (Math.Sign(sfmax) != sign)
+                {
+                    yield return new Tuple<double, double>(smin, smax);
+                    sign = Math.Sign(sfmax);
+                }
+                smin = smax;
+            }
+        }
+    }
+}

src/Numerics/RootFinding/FloatingPointRoots.cs

@@ -35,21 +35,22 @@ namespace MathNet.Numerics.RootFinding
 {
     public static class FloatingPointRoots
     {
-        public static double OfFunction(Func<double, double> f, double lowerBound, double upperBound)
+        public static double OfFunction(Func<double, double> f, double lowerBound, double upperBound, double accuracy = 1e-8)
         {
-            return Brent.FindRoot(f, lowerBound, upperBound, 1e-8, 100);
+            return Brent.FindRoot(f, lowerBound, upperBound, accuracy, 100);
         }
 
-        public static double OfFunctionAndDerivative(Func<double, double> f, Func<double, double> df, double lowerBound, double upperBound)
+        public static double OfFunctionAndDerivative(Func<double, double> f, Func<double, double> df, double lowerBound, double upperBound, double accuracy = 1e-8)
         {
-            try
+            double root;
+            if (HybridNewtonRaphson.TryFindSingleRoot(f, df, lowerBound, upperBound, accuracy, 100, 20, out root))
             {
-                return NewtonRaphson.FindRoot(f, df, lowerBound, upperBound, 1e-8, 100);
-            }
-            catch (NonConvergenceException)
-            {
-                return Brent.FindRoot(f, lowerBound, upperBound, 1e-8, 100);
+                return root;
             }
+
+            return Brent.FindRoot(f, lowerBound, upperBound, accuracy, 100);
+
+            //throw new NonConvergenceException("The algorithm has exceeded the number of iterations allowed");
         }
     }
 }

src/Portable/Portable.csproj

@@ -990,8 +990,11 @@
     <Compile Include="..\Numerics\RootFinding\Algorithms\Brent.cs">
       <Link>RootFinding\Algorithms\Brent.cs</Link>
     </Compile>
-    <Compile Include="..\Numerics\RootFinding\Algorithms\NewtonRaphson.cs">
-      <Link>RootFinding\Algorithms\NewtonRaphson.cs</Link>
+    <Compile Include="..\Numerics\RootFinding\Algorithms\HybridNewtonRaphson.cs">
+      <Link>RootFinding\Algorithms\HybridNewtonRaphson.cs</Link>
+    </Compile>
+    <Compile Include="..\Numerics\RootFinding\Algorithms\ZeroCrossingBracketing.cs">
+      <Link>RootFinding\Algorithms\ZeroCrossingBracketing.cs</Link>
     </Compile>
     <Compile Include="..\Numerics\RootFinding\Bracketing.cs">
       <Link>RootFinding\Bracketing.cs</Link>

src/UnitTests/RootFindingTests/NewtonRaphsonTest.cs

@@ -29,6 +29,7 @@
 // </copyright>
 
 using System;
+using MathNet.Numerics.RootFinding;
 using MathNet.Numerics.RootFinding.Algorithms;
 using NUnit.Framework;
 
@@ -43,38 +44,61 @@ public void MultipleRoots()
             // Roots at -2, 2
             Func<double, double> f1 = x => x * x - 4;
             Func<double, double> df1 = x => 2 * x;
-            Assert.AreEqual(0, f1(NewtonRaphson.FindRoot(f1, df1, -5, 5, 1e-14, 100)));
-            Assert.AreEqual(-2, NewtonRaphson.FindRoot(f1, df1, -5, -1, 1e-14, 100));
-            Assert.AreEqual(2, NewtonRaphson.FindRoot(f1, df1, 1, 4, 1e-14, 100));
-            Assert.AreEqual(0, f1(NewtonRaphson.FindRoot(x => -f1(x), x => -df1(x), -5, 5, 1e-14, 100)));
-            Assert.AreEqual(-2, NewtonRaphson.FindRoot(x => -f1(x), x => -df1(x), -5, -1, 1e-14, 100));
-            Assert.AreEqual(2, NewtonRaphson.FindRoot(x => -f1(x), x => -df1(x), 1, 4, 1e-14, 100));
+            Assert.AreEqual(0, f1(HybridNewtonRaphson.FindSingleRoot(f1, df1, -5, 5, 1e-14, 100, 20)));
+            Assert.AreEqual(-2, HybridNewtonRaphson.FindSingleRoot(f1, df1, -5, -1, 1e-14, 100, 20));
+            Assert.AreEqual(2, HybridNewtonRaphson.FindSingleRoot(f1, df1, 1, 4, 1e-14, 100, 20));
+            Assert.AreEqual(0, f1(HybridNewtonRaphson.FindSingleRoot(x => -f1(x), x => -df1(x), -5, 5, 1e-14, 100, 20)));
+            Assert.AreEqual(-2, HybridNewtonRaphson.FindSingleRoot(x => -f1(x), x => -df1(x), -5, -1, 1e-14, 100, 20));
+            Assert.AreEqual(2, HybridNewtonRaphson.FindSingleRoot(x => -f1(x), x => -df1(x), 1, 4, 1e-14, 100, 20));
 
             // Roots at 3, 4
             Func<double, double> f2 = x => (x - 3) * (x - 4);
             Func<double, double> df2 = x => 2 * x - 7;
-            Assert.AreEqual(0, f2(NewtonRaphson.FindRoot(f2, df2, -5, 5, 1e-14, 100)));
-            Assert.AreEqual(3, NewtonRaphson.FindRoot(f2, df2, -5, 3.5, 1e-14, 100));
-            Assert.AreEqual(4, NewtonRaphson.FindRoot(f2, df2, 3.2, 5, 1e-14, 100));
-            Assert.AreEqual(3, NewtonRaphson.FindRoot(f2, df2, 2.1, 3.9, 0.001, 50), 0.001);
-            Assert.AreEqual(3, NewtonRaphson.FindRoot(f2, df2, 2.1, 3.4, 0.001, 50), 0.001);
+            Assert.AreEqual(0, f2(HybridNewtonRaphson.FindSingleRoot(f2, df2, -5, 5, 1e-14, 100, 20)));
+            Assert.AreEqual(3, HybridNewtonRaphson.FindSingleRoot(f2, df2, -5, 3.5, 1e-14, 100, 20));
+            Assert.AreEqual(4, HybridNewtonRaphson.FindSingleRoot(f2, df2, 3.2, 5, 1e-14, 100, 20));
+            Assert.AreEqual(3, HybridNewtonRaphson.FindSingleRoot(f2, df2, 2.1, 3.9, 0.001, 50, 20), 0.001);
+            Assert.AreEqual(3, HybridNewtonRaphson.FindSingleRoot(f2, df2, 2.1, 3.4, 0.001, 50, 20), 0.001);
         }
 
         [Test]
         public void LocalMinima()
         {
             Func<double, double> f1 = x => x * x * x - 2 * x + 2;
             Func<double, double> df1 = x => 3 * x * x - 2;
-            Assert.AreEqual(0, f1(NewtonRaphson.FindRoot(f1, df1, -5, 5, 1e-14, 100)));
-            Assert.AreEqual(0, f1(NewtonRaphson.FindRoot(f1, df1, -2, 4, 1e-14, 100)));
+            Assert.AreEqual(0, f1(HybridNewtonRaphson.FindSingleRoot(f1, df1, -5, 5, 1e-14, 100, 20)));
+            Assert.AreEqual(0, f1(HybridNewtonRaphson.FindSingleRoot(f1, df1, -2, 4, 1e-14, 100, 20)));
+        }
+
+        [Test]
+        public void Pole()
+        {
+            Func<double, double> f1 = x => 1/(x - 2) + 2;
+            Func<double, double> df1 = x => -1/(x*x - 4*x + 4);
+            Assert.AreEqual(1.5, HybridNewtonRaphson.FindSingleRoot(f1, df1, 1, 2, 1e-14, 100, 20));
+            Assert.AreEqual(1.5, HybridNewtonRaphson.FindSingleRoot(f1, df1, 1, 6, 1e-14, 100, 20));
+            Assert.AreEqual(1.5, FloatingPointRoots.OfFunctionAndDerivative(f1, df1, 1, 6));
+
+            Func<double, double> f2 = x => -1/(x - 2) + 2;
+            Func<double, double> df2 = x => 1/(x*x - 4*x + 4);
+            Assert.AreEqual(2.5, HybridNewtonRaphson.FindSingleRoot(f2, df2, 2, 3, 1e-14, 100, 20));
+            Assert.AreEqual(2.5, HybridNewtonRaphson.FindSingleRoot(f2, df2, -2, 3, 1e-14, 100, 20));
+            Assert.AreEqual(2.5, FloatingPointRoots.OfFunctionAndDerivative(f2, df2, -2, 3));
+
+            Func<double, double> f3 = x => 1/(x - 2) + x + 2;
+            Func<double, double> df3 = x => -1/(x*x - 4*x + 4) + 1;
+            Assert.AreEqual(-Math.Sqrt(3), HybridNewtonRaphson.FindSingleRoot(f3, df3, -2, -1, 1e-14, 100, 20), 1e-14);
+            Assert.AreEqual(Math.Sqrt(3), HybridNewtonRaphson.FindSingleRoot(f3, df3, 1, 1.99, 1e-14, 100, 20));
+            Assert.AreEqual(Math.Sqrt(3), HybridNewtonRaphson.FindSingleRoot(f3, df3, -1.5, 1.99, 1e-14, 100, 20));
+            Assert.AreEqual(Math.Sqrt(3), HybridNewtonRaphson.FindSingleRoot(f3, df3, 1, 6, 1e-14, 100, 20));
         }
 
         [Test]
         public void NoRoot()
         {
             Func<double, double> f1 = x => x * x + 4;
             Func<double, double> df1 = x => 2 * x;
-            Assert.Throws<NonConvergenceException>(() => NewtonRaphson.FindRoot(f1, df1, -5, 5, 1e-14, 50));
+            Assert.Throws<NonConvergenceException>(() => HybridNewtonRaphson.FindSingleRoot(f1, df1, -5, 5, 1e-14, 50, 20));
         }
     }
 }