RootFinding: add hybrid newton-raphson/bisection algorithm

src/Numerics/Numerics.csproj

@@ -110,6 +110,7 @@
     <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\Bracketing.cs" />
     <Compile Include="RootFinding\Algorithms\Brent.cs" />
     <Compile Include="RootFinding\FloatingPointRoots.cs" />

src/Numerics/Properties/Resources.resx

@@ -381,4 +381,7 @@
   <data name="RootNotFound" xml:space="preserve">
     <value>The algorithm ended without root in the range.</value>
   </data>
+  <data name="RootMustBeBracketedByBounds" xml:space="preserve">
+    <value>The lower and upper bounds must bracket a single root.</value>
+  </data>
 </root>

src/Numerics/RootFinding/Algorithms/Bisection.cs

@@ -29,41 +29,39 @@
 // </copyright>
 
 using System;
+using MathNet.Numerics.Properties;
 
 namespace MathNet.Numerics.RootFinding.Algorithms
 {
     public static class Bisection
     {
-        public static double FindRootExpand(Func<double, double> f, double guessLowerBound, double guessUpperBound, double accuracy = 1e-5, double expandFactor = 1.6, int maxExpandIteratons = 100)
+        /// <summary>Find a solution of the equation f(x)=0.</summary>
+        /// <exception cref="NonConvergenceException"></exception>
+        public static double FindRootExpand(Func<double, double> f, double guessLowerBound, double guessUpperBound, double accuracy = 1e-8, double expandFactor = 1.6, int maxExpandIteratons = 100)
         {
             Bracketing.Expand(f, ref guessLowerBound, ref guessUpperBound, expandFactor, maxExpandIteratons);
             return FindRoot(f, guessLowerBound, guessUpperBound, accuracy);
         }
 
         /// <summary>Find a solution of the equation f(x)=0.</summary>
-        public static double FindRoot(Func<double, double> f, double lowerBound, double upperBound, double accuracy = 1e-5)
+        /// <exception cref="NonConvergenceException"></exception>
+        public static double FindRoot(Func<double, double> f, double lowerBound, double upperBound, double accuracy = 1e-8)
         {
             double fmin = f(lowerBound);
             double fmax = f(upperBound);
 
-            if (fmin == 0.0)
-                return lowerBound;
-            if (fmax == 0.0)
-                return upperBound;
-
-            ValidateEvaluation(fmin, lowerBound);
-            ValidateEvaluation(fmax, upperBound);
+            if (fmin == 0.0) return lowerBound;
+            if (fmax == 0.0) return upperBound;
 
             if (Math.Sign(fmin) == Math.Sign(fmax))
             {
-                throw new NonConvergenceException("Bounds do not necessarily span a root.");
+                throw new NonConvergenceException(Resources.RootMustBeBracketedByBounds);
             }
 
             while (Math.Abs(fmax - fmin) > 0.5 * accuracy || Math.Abs(upperBound - lowerBound) > 0.5 * Precision.DoubleMachinePrecision)
             {
                 double midpoint = 0.5*(upperBound + lowerBound);
                 double midval = f(midpoint);
-                ValidateEvaluation(midval, midpoint);
 
                 if (Math.Sign(midval) == Math.Sign(fmin))
                 {
@@ -83,13 +81,5 @@ public static double FindRoot(Func<double, double> f, double lowerBound, double
 
             return 0.5*(lowerBound + upperBound);
         }
-
-        static void ValidateEvaluation(double output, double input)
-        {
-            if (Double.IsInfinity(output) || Double.IsInfinity(output))
-            {
-                throw new Exception(String.Format("Objective function returned non-finite result: f({0}) = {1}", input, output));
-            }
-        }
     }
 }

src/Numerics/RootFinding/Algorithms/Brent.cs

@@ -139,8 +139,7 @@ public static double FindRoot(Func<double, double> f, double lowerBound, double
                 froot = f(root);
             }
 
-            // The algorithm has exceeded the number of iterations allowed
-            throw new NonConvergenceException();
+            throw new NonConvergenceException("The algorithm has exceeded the number of iterations allowed");
         }
 
         /// <summary>Helper method useful for preventing rounding errors.</summary>

src/Numerics/RootFinding/Algorithms/NewtonRaphson.cs

@@ -0,0 +1,107 @@
+// <copyright file="NewtonRaphson.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+// 
+// Copyright (c) 2009-2013 Math.NET
+// 
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+// 
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+// 
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using System;
+
+namespace MathNet.Numerics.RootFinding.Algorithms
+{
+    public static class NewtonRaphson
+    {
+        /// <summary>Find a solution of the equation f(x)=0.</summary>
+        /// <remarks>Hybrid Newton-Raphson that when failing falls back to bisection.</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)
+        {
+            double fmin = f(lowerBound);
+            double fmax = f(upperBound);
+
+            if (fmin == 0.0) return lowerBound;
+            if (fmax == 0.0) return upperBound;
+
+            double root = 0.5*(lowerBound + upperBound);
+            double fx = f(root);
+            double lastStep = Math.Abs(upperBound - lowerBound);
+            for (int i = 0; i < maxIterations; i++)
+            {
+                double dfx = df(root);
+
+                // Netwon-Raphson step
+                double step = fx/dfx;
+                root -= step;
+
+                if (root < lowerBound || root > upperBound || Math.Abs(2*fx) > Math.Abs(lastStep*dfx))
+                {
+                    // Newton-Raphson step failed ->  bisect instead
+                    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;
+                    }
+                    else
+                    {
+                        upperBound = root;
+                        fmax = fx;
+                    }
+                    continue;
+                }
+
+                if (Math.Abs(step) < accuracy)
+                {
+                    return root;
+                }
+
+                // Evaluation
+                fx = f(root);
+                lastStep = step;
+
+                // Update Bounds
+                if (Math.Sign(fx) != Math.Sign(fmin))
+                {
+                    upperBound = root;
+                    fmax = fx;
+                }
+                else if (Math.Sign(fx) != Math.Sign(fmax))
+                {
+                    lowerBound = root;
+                    fmin = fx;
+                }
+                else if (Math.Sign(fmin) != Math.Sign(fmax))
+                {
+                    return root;
+                }
+            }
+
+            throw new NonConvergenceException("The algorithm has exceeded the number of iterations allowed");
+        }
+    }
+}

src/Numerics/RootFinding/FloatingPointRoots.cs

@@ -39,5 +39,17 @@ public static double OfFunction(Func<double, double> f, double lowerBound, doubl
         {
             return Brent.FindRoot(f, lowerBound, upperBound, 1e-8, 100);
         }
+
+        public static double OfFunctionAndDerivative(Func<double, double> f, Func<double, double> df, double lowerBound, double upperBound)
+        {
+            try
+            {
+                return NewtonRaphson.FindRoot(f, df, lowerBound, upperBound, 1e-8, 100);
+            }
+            catch (NonConvergenceException)
+            {
+                return Brent.FindRoot(f, lowerBound, upperBound, 1e-8, 100);
+            }
+        }
     }
 }

src/Portable/Portable.csproj

@@ -990,6 +990,9 @@
     <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>
     <Compile Include="..\Numerics\RootFinding\Bracketing.cs">
       <Link>RootFinding\Bracketing.cs</Link>
     </Compile>

src/UnitTests/RootFindingTests/BisectionTest.cs

@@ -46,6 +46,9 @@ public void MultipleRoots()
             Assert.AreEqual(0, f1(Bisection.FindRootExpand(f1, 3, 4, 1e-14)));
             Assert.AreEqual(-2, Bisection.FindRoot(f1, -5, -1, 1e-14));
             Assert.AreEqual(2, Bisection.FindRoot(f1, 1, 4, 1e-14));
+            Assert.AreEqual(0, f1(Bisection.FindRoot(x => -f1(x), 0, 5, 1e-14)));
+            Assert.AreEqual(-2, Bisection.FindRoot(x => -f1(x), -5, -1, 1e-14));
+            Assert.AreEqual(2, Bisection.FindRoot(x => -f1(x), 1, 4, 1e-14));
 
             // Roots at 3, 4
             Func<double, double> f2 = x => (x - 3) * (x - 4);
@@ -56,6 +59,14 @@ public void MultipleRoots()
             Assert.AreEqual(3, Bisection.FindRoot(f2, 2.1, 3.4, 0.001), 0.001);
         }
 
+        [Test]
+        public void LocalMinima()
+        {
+            Func<double, double> f1 = x => x * x * x - 2 * x + 2;
+            Assert.AreEqual(0, f1(Bisection.FindRoot(f1, -5, 5, 1e-14)), 1e-14);
+            Assert.AreEqual(0, f1(Bisection.FindRoot(f1, -2, 4, 1e-14)), 1e-14);
+        }
+
         [Test]
         public void NoRoot()
         {

src/UnitTests/RootFindingTests/BrentTest.cs

@@ -45,6 +45,9 @@ public void MultipleRoots()
             Assert.AreEqual(0, f1(Brent.FindRoot(f1, -5, 5, 1e-14, 100)));
             Assert.AreEqual(-2, Brent.FindRoot(f1, -5, -1, 1e-14, 100));
             Assert.AreEqual(2, Brent.FindRoot(f1, 1, 4, 1e-14, 100));
+            Assert.AreEqual(0, f1(Brent.FindRoot(x => -f1(x), -5, 5, 1e-14, 100)));
+            Assert.AreEqual(-2, Brent.FindRoot(x => -f1(x), -5, -1, 1e-14, 100));
+            Assert.AreEqual(2, Brent.FindRoot(x => -f1(x), 1, 4, 1e-14, 100));
 
             // Roots at 3, 4
             Func<double, double> f2 = x => (x - 3)*(x - 4);
@@ -55,6 +58,14 @@ public void MultipleRoots()
             Assert.AreEqual(3, Brent.FindRoot(f2, 2.1, 3.4, 0.001, 50), 0.001);
         }
 
+        [Test]
+        public void LocalMinima()
+        {
+            Func<double, double> f1 = x => x * x * x - 2 * x + 2;
+            Assert.AreEqual(0, f1(Brent.FindRoot(f1, -5, 5, 1e-14, 100)), 1e-14);
+            Assert.AreEqual(0, f1(Brent.FindRoot(f1, -2, 4, 1e-14, 100)), 1e-14);
+        }
+
         [Test]
         public void NoRoot()
         {