Optimization: Minimally working Conjugate Gradient, BFGS, and Newton minimizers

src/Numerics/Numerics.csproj

@@ -262,7 +262,7 @@
     <Compile Include="RootFinding\Brent.cs" />
     <Compile Include="FindRoots.cs" />
     <Compile Include="RootFinding\Bisection.cs" />
-    <Compile Include="Optimization\BFGS.cs" />
+    <Compile Include="Optimization\BfgsMinimizer.cs" />
     <Compile Include="Optimization\ConjugateGradientMinimizer.cs" />
     <Compile Include="Optimization\GoldenSectionMinimizer.cs" />
     <Compile Include="Optimization\ObjectiveFunction.cs" />

src/Numerics/Optimization/BFGS.cs

@@ -1,11 +0,0 @@
-using System;
-using System.Collections.Generic;
-using System.Linq;
-using System.Text;
-
-namespace MathNet.Numerics.Optimization
-{
-    class BFGS
-    {
-    }
-}

src/Numerics/Optimization/BfgsMinimizer.cs

@@ -0,0 +1,124 @@
+using System;
+using MathNet.Numerics.LinearAlgebra;
+using MathNet.Numerics.Optimization.LineSearch;
+
+namespace MathNet.Numerics.Optimization
+{
+    public class BfgsMinimizer
+    {
+        public double GradientTolerance { get; set; }
+        public int MaximumIterations { get; set; }
+
+        public BfgsMinimizer(double gradientTolerance, int maximumIterations)
+        {
+            GradientTolerance = gradientTolerance;
+            MaximumIterations = maximumIterations;
+        }
+
+        public MinimizationResult FindMinimum(IObjectiveFunction objective, Vector<double> initialGuess)
+        {
+            if (!objective.IsGradientSupported)
+                throw new IncompatibleObjectiveException("Gradient not supported in objective function, but required for BFGS minimization.");
+
+            objective.EvaluateAt(initialGuess);
+
+            ValidateGradient(objective);
+
+            // Check that we're not already done
+            if (ExitCriteriaSatisfied(objective.Point, objective.Gradient))
+                return new MinimizationResult(objective, 0, MinimizationResult.ExitCondition.AbsoluteGradient);
+
+            // Set up line search algorithm
+            var lineSearcher = new WeakWolfeLineSearch(1e-4, 0.9, 1000);
+
+            // First step
+            var inversePseudoHessian = CreateMatrix.DenseIdentity<double>(initialGuess.Count);
+            var searchDirection = -objective.Gradient;
+            var stepSize = 100 * GradientTolerance / (searchDirection * searchDirection);
+
+            var previousGradient = objective.Gradient;
+
+            LineSearchResult result;
+            try
+            {
+                result = lineSearcher.FindConformingStep(objective, searchDirection, stepSize);
+            }
+            catch (Exception e)
+            {
+                throw new InnerOptimizationException("Line search failed.", e);
+            }
+
+            objective = result.FunctionInfoAtMinimum;
+            ValidateGradient(objective);
+
+            var gradient = objective.Gradient;
+            var step = objective.Point - initialGuess;
+            stepSize = result.FinalStep;
+
+            // Subsequent steps
+            int iterations = 1;
+            int totalLineSearchSteps = result.Iterations;
+            int iterationsWithNontrivialLineSearch = result.Iterations > 0 ? 0 : 1;
+            while (!ExitCriteriaSatisfied(objective.Point, objective.Gradient) && iterations < MaximumIterations)
+            {
+                var y = objective.Gradient - previousGradient;
+
+                double sy = step * y;
+                inversePseudoHessian = inversePseudoHessian + ((sy + y * inversePseudoHessian * y) / Math.Pow(sy, 2.0)) * step.OuterProduct(step) - ( (inversePseudoHessian * y.ToColumnMatrix())*step.ToRowMatrix() + step.ToColumnMatrix()*(y.ToRowMatrix() * inversePseudoHessian)) * (1.0 / sy);
+
+                searchDirection = -inversePseudoHessian * objective.Gradient;
+
+                if (searchDirection * objective.Gradient >= 0)
+                {
+                    searchDirection = -objective.Gradient;
+                    inversePseudoHessian = CreateMatrix.DenseIdentity<double>(initialGuess.Count);
+                }
+
+                previousGradient = objective.Gradient;
+                var previousPoint = objective.Point;
+
+                try
+                {
+                    result = lineSearcher.FindConformingStep(objective, searchDirection, 1.0);
+                }
+                catch (Exception e)
+                {
+                    throw new InnerOptimizationException("Line search failed.", e);
+                }
+
+                iterationsWithNontrivialLineSearch += result.Iterations > 0 ? 1 : 0;
+                totalLineSearchSteps += result.Iterations;
+                stepSize = result.FinalStep;
+                step = result.FunctionInfoAtMinimum.Point - previousPoint;
+                objective = result.FunctionInfoAtMinimum;
+
+                iterations += 1;
+            }
+
+            if (iterations == this.MaximumIterations)
+                throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", MaximumIterations));
+
+            return new MinimizationWithLineSearchResult(objective, iterations, MinimizationResult.ExitCondition.AbsoluteGradient, totalLineSearchSteps, iterationsWithNontrivialLineSearch);
+        }
+
+        private bool ExitCriteriaSatisfied(Vector<double> candidatePoint, Vector<double> gradient)
+        {
+            return gradient.Norm(2.0) < this.GradientTolerance;
+        }
+
+        private void ValidateGradient(IObjectiveFunction objective)
+        {
+            foreach (var x in objective.Gradient)
+            {
+                if (Double.IsNaN(x) || Double.IsInfinity(x))
+                    throw new EvaluationException("Non-finite gradient returned.", objective);
+            }
+        }
+
+        private void ValidateObjective(IObjectiveFunction objective)
+        {
+            if (Double.IsNaN(objective.Value) || Double.IsInfinity(objective.Value))
+                throw new EvaluationException("Non-finite objective function returned.", objective);
+        }
+    }
+}

src/Numerics/Optimization/ConjugateGradientMinimizer.cs

@@ -49,12 +49,12 @@ namespace MathNet.Numerics.Optimization
             objective = result.FunctionInfoAtMinimum;
             ValidateGradient(objective);
 
-            double stepSize = (objective.Point - initialGuess).Norm(2.0);
+            double stepSize = result.FinalStep;
 
             // Subsequent steps
             int iterations = 1;
             int totalLineSearchSteps = result.Iterations;
-            int noLineSearchIterations = result.Iterations > 0 ? 0 : 1;
+            int iterationsWithNontrivialLineSearch = result.Iterations > 0 ? 0 : 1;
             int steepestDescentResets = 0;
             while (!ExitCriteriaSatisfied(objective.Point, objective.Gradient) && iterations < MaximumIterations)
             {
@@ -77,10 +77,9 @@ namespace MathNet.Numerics.Optimization
                     throw new InnerOptimizationException("Line search failed.", e);
                 }
 
-                noLineSearchIterations += result.Iterations == 0 ? 1 : 0;
+                iterationsWithNontrivialLineSearch += result.Iterations == 0 ? 1 : 0;
                 totalLineSearchSteps += result.Iterations;
-				stepSize = (result.FunctionInfoAtMinimum.Point - objective.Point).Norm(2.0);
-
+				stepSize = result.FinalStep;
                 objective = result.FunctionInfoAtMinimum;
                 iterations += 1;
             }
@@ -90,7 +89,7 @@ namespace MathNet.Numerics.Optimization
                 throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", MaximumIterations));
             }
 
-            return new MinimizationWithLineSearchResult(objective, iterations, MinimizationResult.ExitCondition.AbsoluteGradient, totalLineSearchSteps, noLineSearchIterations);
+            return new MinimizationWithLineSearchResult(objective, iterations, MinimizationResult.ExitCondition.AbsoluteGradient, totalLineSearchSteps, iterationsWithNontrivialLineSearch);
         }
 
         bool ExitCriteriaSatisfied(Vector<double> candidatePoint, Vector<double> gradient)

src/UnitTests/OptimizationTests/TestBfgsMinimizer.cs

@@ -0,0 +1,44 @@
+using System;
+using MathNet.Numerics.LinearAlgebra.Double;
+using MathNet.Numerics.Optimization;
+using NUnit.Framework;
+
+namespace MathNet.Numerics.UnitTests.OptimizationTests
+{
+    [TestFixture]
+    public class TestBfgsMinimizer
+    {
+        [Test]
+        public void FindMinimum_Rosenbrock_Easy()
+        {
+            var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
+            var solver = new BfgsMinimizer(1e-5, 1000);
+            var result = solver.FindMinimum(obj, new DenseVector(new[] { 1.2, 1.2 }));
+
+            Assert.That(Math.Abs(result.MinimizingPoint[0] - 1.0), Is.LessThan(1e-3));
+            Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3));
+        }
+
+        [Test]
+        public void FindMinimum_Rosenbrock_Hard()
+        {
+            var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
+            var solver = new BfgsMinimizer(1e-5, 1000);
+            var result = solver.FindMinimum(obj, new DenseVector(new[] { -1.2, 1.0 }));
+
+            Assert.That(Math.Abs(result.MinimizingPoint[0] - 1.0), Is.LessThan(1e-3));
+            Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3));
+        }
+
+        [Test]
+        public void FindMinimum_Rosenbrock_Overton()
+        {
+            var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
+            var solver = new BfgsMinimizer(1e-5, 1000);
+            var result = solver.FindMinimum(obj, new DenseVector(new[] { -0.9, -0.5 }));
+
+            Assert.That(Math.Abs(result.MinimizingPoint[0] - 1.0), Is.LessThan(1e-3));
+            Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3));
+        }
+    }
+}

src/UnitTests/OptimizationTests/TestRosenbrockFunction.cs

@@ -0,0 +1,55 @@
+using System;
+using MathNet.Numerics.LinearAlgebra.Double;
+using NUnit.Framework;
+
+namespace MathNet.Numerics.UnitTests.OptimizationTests
+{
+    [TestFixture]
+    class TestRosenbrockFunction
+    {
+        [Test]
+        public void TestGradient()
+        {
+            var input = new DenseVector(new[]{ -0.9, -0.5 } );
+
+            var v1 = RosenbrockFunction.Value(input);
+            var g = RosenbrockFunction.Gradient(input);
+
+            var eps = 1e-5;
+            var eps0 = (new DenseVector(new[] { 1.0, 0.0 })) * eps;
+            var eps1 = (new DenseVector(new[] { 0.0, 1.0 })) * eps;
+
+            var g0 = (RosenbrockFunction.Value(input + eps0) - RosenbrockFunction.Value(input - eps0)) / (2 * eps);
+            var g1 = (RosenbrockFunction.Value(input + eps1) - RosenbrockFunction.Value(input - eps1)) / (2 * eps);
+
+            Assert.That(Math.Abs(g0 - g[0]) < 1e-3);
+            Assert.That(Math.Abs(g1 - g[1]) < 1e-3);
+        }
+
+        [Test]
+        public void TestHessian()
+        {
+            var input = new DenseVector(new[] { -0.9, -0.5 });
+
+            var v1 = RosenbrockFunction.Value(input);
+            var h = RosenbrockFunction.Hessian(input);
+
+            var eps = 1e-5;
+
+            var eps0 = (new DenseVector(new[] { 1.0, 0.0 })) * eps;
+            var eps1 = (new DenseVector(new[] { 0.0, 1.0 })) * eps;
+
+            var epsuu = (new DenseVector(new[] { 1.0, 1.0 })) * eps;
+            var epsud = (new DenseVector(new[] { 1.0, -1.0 })) * eps;
+
+            var h00 = (RosenbrockFunction.Value(input + eps0) - 2*RosenbrockFunction.Value(input) + RosenbrockFunction.Value(input - eps0)) / (eps*eps);
+            var h11 = (RosenbrockFunction.Value(input + eps1) - 2 * RosenbrockFunction.Value(input) + RosenbrockFunction.Value(input - eps1)) / (eps * eps);
+            var h01 = (RosenbrockFunction.Value(input + epsuu) - RosenbrockFunction.Value(input + epsud) - RosenbrockFunction.Value(input - epsud) + RosenbrockFunction.Value(input - epsuu)) / (4*eps * eps);
+
+            Assert.That(Math.Abs(h00 - h[0,0]) < 1e-3);
+            Assert.That(Math.Abs(h11 - h[1,1]) < 1e-3);
+            Assert.That(Math.Abs(h01 - h[0, 1]) < 1e-3);
+            Assert.That(Math.Abs(h01 - h[1, 0]) < 1e-3);
+        }
+    }
+}

src/UnitTests/UnitTests.csproj

@@ -372,7 +372,10 @@
     <Compile Include="OptimizationTests\BfgsTest.cs" />
     <Compile Include="RootFindingTests\BisectionTest.cs" />
     <Compile Include="OptimizationTests\RosenbrockFunction.cs" />
+    <Compile Include="OptimizationTests\TestBfgsMinimizer.cs" />
     <Compile Include="OptimizationTests\TestConjugateGradientMinimizer.cs" />
+    <Compile Include="OptimizationTests\TestNewtonMinimizer.cs" />
+    <Compile Include="OptimizationTests\TestRosenbrockFunction.cs" />
     <Compile Include="PermutationTest.cs" />
     <Compile Include="PrecisionTest.cs" />
     <Compile Include="Properties\AssemblyInfo.cs" />