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

(updated by cdrnet)

src/Numerics/Numerics.csproj

@@ -195,7 +195,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\BisectionRootFinder.cs" />
     <Compile Include="Optimization\ConjugateGradientMinimizer.cs" />
     <Compile Include="Optimization\Exceptions.cs" />
@@ -203,6 +203,7 @@
     <Compile Include="Optimization\LineSearchingMinimizerOutput.cs" />
     <Compile Include="Optimization\LineSearchOutput.cs" />
     <Compile Include="Optimization\MinimizationOutput.cs" />
+    <Compile Include="Optimization\NewtonMinimizer.cs" />
     <Compile Include="Optimization\ObjectiveChecker.cs" />
     <Compile Include="Optimization\ObjectiveFunction.cs" />
     <Compile Include="Optimization\OptimizationResult.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,131 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using MathNet.Numerics.LinearAlgebra;
+
+namespace MathNet.Numerics.Optimization
+{
+    public class BfgsMinimizer
+    {
+        public double GradientTolerance { get; set; }
+        public int MaximumIterations { get; set; }
+
+        public BfgsMinimizer(double gradient_tolerance, int maximum_iterations)
+        {
+            this.GradientTolerance = gradient_tolerance;
+            this.MaximumIterations = maximum_iterations;
+        }
+
+        public MinimizationOutput FindMinimum(IObjectiveFunction objective, Vector<double> initial_guess)
+        {
+            if (!objective.GradientSupported)
+                throw new IncompatibleObjectiveException("Gradient not supported in objective function, but required for BFGS minimization.");
+
+            if (!(objective is ObjectiveChecker))
+                objective = new ObjectiveChecker(objective, this.ValidateObjective, this.ValidateGradient, null);
+
+            IEvaluation initial_eval = objective.Evaluate(initial_guess);            
+            
+            // Check that we're not already done
+            if (this.ExitCriteriaSatisfied(initial_guess, initial_eval.Gradient))
+                return new MinimizationOutput(initial_eval, 0);
+            
+            // Set up line search algorithm
+            var line_searcher = new WeakWolfeLineSearch(1e-4, 0.9, 1000);
+
+            // Declare state variables
+            IEvaluation candidate_point;
+            double step_size;
+            Vector<double> gradient, previous_gradient, step, search_direction;
+            Matrix<double> inverse_pseudo_hessian;
+
+            // First step
+            inverse_pseudo_hessian = Matrix<double>.Build.DiagonalIdentity(initial_guess.Count);
+            search_direction = -initial_eval.Gradient;
+            step_size = 100 * this.GradientTolerance / (search_direction * search_direction);
+                        
+            LineSearchOutput result;
+            try 
+            {
+                result = line_searcher.FindConformingStep(objective, initial_eval, search_direction, step_size);
+            } 
+            catch (Exception e) 
+            {
+                throw new InnerOptimizationException("Line search failed.", e);
+            }
+
+            candidate_point = result.FunctionInfoAtMinimum;
+            gradient = candidate_point.Gradient;
+            previous_gradient = initial_eval.Gradient;
+            step = candidate_point.Point - initial_guess;
+            step_size = result.FinalStep;
+
+            // Subsequent steps
+            int iterations = 1;
+            int total_line_search_steps = result.Iterations;
+            int iterations_with_nontrivial_line_search = result.Iterations > 0 ? 0 : 1;
+            int steepest_descent_resets = 0;
+            while (!this.ExitCriteriaSatisfied(candidate_point.Point, candidate_point.Gradient) && iterations < this.MaximumIterations)
+            {
+                var y = candidate_point.Gradient - previous_gradient;
+
+                double sy = step * y;
+                inverse_pseudo_hessian = inverse_pseudo_hessian + ((sy + y * inverse_pseudo_hessian * y) / Math.Pow(sy, 2.0)) * step.OuterProduct(step) - ( (inverse_pseudo_hessian * y.ToColumnMatrix())*step.ToRowMatrix() + step.ToColumnMatrix()*(y.ToRowMatrix() * inverse_pseudo_hessian)) * (1.0 / sy);
+
+               search_direction = -inverse_pseudo_hessian * candidate_point.Gradient;
+
+                if (search_direction * candidate_point.Gradient >= 0)
+                {
+                    search_direction = -candidate_point.Gradient;
+                    inverse_pseudo_hessian = Matrix<double>.Build.DiagonalIdentity(initial_guess.Count);
+                    steepest_descent_resets += 1;
+                }
+
+                try
+                {
+                    result = line_searcher.FindConformingStep(objective, candidate_point, search_direction, 1.0);
+                }
+                catch (Exception e)
+                {
+                    throw new InnerOptimizationException("Line search failed.", e);
+                }
+
+                iterations_with_nontrivial_line_search += result.Iterations > 0 ? 1 : 0;
+                total_line_search_steps += result.Iterations;
+
+                step_size = result.FinalStep;
+                step = result.FunctionInfoAtMinimum.Point - candidate_point.Point;
+                previous_gradient = candidate_point.Gradient;
+                candidate_point = result.FunctionInfoAtMinimum;
+
+                iterations += 1;
+            }
+
+            if (iterations == this.MaximumIterations)
+                throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", this.MaximumIterations));
+
+            return new MinimizationWithLineSearchOutput(candidate_point, iterations, total_line_search_steps, iterations_with_nontrivial_line_search);
+        }
+
+        private bool ExitCriteriaSatisfied(Vector<double> candidate_point, Vector<double> gradient)
+        {
+            return gradient.Norm(2.0) < this.GradientTolerance;
+        }
+
+        private void ValidateGradient(Vector<double> gradient, Vector<double> input)
+        {
+            foreach (var x in gradient)
+            {
+                if (Double.IsNaN(x) || Double.IsInfinity(x))
+                    throw new EvaluationException("Non-finite gradient returned.");
+            }
+        }
+
+        private void ValidateObjective(double objective, Vector<double> input)
+        {
+            if (Double.IsNaN(objective) || Double.IsInfinity(objective))
+                throw new EvaluationException("Non-finite objective function returned.");
+        }
+    }
+}

src/Numerics/Optimization/ConjugateGradientMinimizer.cs

@@ -23,7 +23,8 @@ namespace MathNet.Numerics.Optimization
             if (!objective.GradientSupported)
                 throw new IncompatibleObjectiveException("Gradient not supported in objective function, but required for ConjugateGradient minimization.");
 
-            objective = new ObjectiveChecker(objective, this.ValidateObjective, this.ValidateGradient, null);
+            if (!(objective is ObjectiveChecker))
+                objective = new ObjectiveChecker(objective, this.ValidateObjective, this.ValidateGradient, null);
 
             IEvaluation initial_eval = objective.Evaluate(initial_guess);
             var gradient = initial_eval.Gradient;            
@@ -53,14 +54,14 @@ namespace MathNet.Numerics.Optimization
                 throw new InnerOptimizationException("Line search failed.", e);
             }
 
-            candidate_point = result.FunctionInfoAtMinimum;            
-            
-			double step_size = (candidate_point.Point - initial_guess).Norm(2.0);
+            candidate_point = result.FunctionInfoAtMinimum;
+
+            double step_size = result.FinalStep;
             
             // Subsequent steps
             int iterations = 1;
             int total_line_search_steps = result.Iterations;
-            int no_line_search_iterations = result.Iterations > 0 ? 0 : 1;
+            int iterations_with_nontrivial_line_search = result.Iterations > 0 ? 0 : 1;
             int steepest_descent_resets = 0;
             while (!this.ExitCriteriaSatisfied(candidate_point.Point, candidate_point.Gradient) && iterations < this.MaximumIterations)
             {
@@ -93,10 +94,10 @@ namespace MathNet.Numerics.Optimization
                     throw new InnerOptimizationException("Line search failed.", e);
                 }
 
-                no_line_search_iterations += result.Iterations == 0 ? 1 : 0;
+                iterations_with_nontrivial_line_search += result.Iterations > 0 ? 1 : 0;
                 total_line_search_steps += result.Iterations;
 
-				step_size = (result.FunctionInfoAtMinimum.Point - candidate_point.Point).Norm (2.0);
+                step_size = result.FinalStep;
                 candidate_point = result.FunctionInfoAtMinimum;
 
                 iterations += 1;
@@ -105,7 +106,7 @@ namespace MathNet.Numerics.Optimization
             if (iterations == this.MaximumIterations)
                 throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", this.MaximumIterations));
 
-            return new MinimizationWithLineSearchOutput(candidate_point, iterations, total_line_search_steps, no_line_search_iterations);
+            return new MinimizationWithLineSearchOutput(candidate_point, iterations, total_line_search_steps, iterations_with_nontrivial_line_search);
         }
 
         private bool ExitCriteriaSatisfied(Vector<double> candidate_point, Vector<double> gradient)

src/Numerics/Optimization/LineSearchingMinimizerOutput.cs

@@ -8,13 +8,13 @@ namespace MathNet.Numerics.Optimization
     public class MinimizationWithLineSearchOutput : MinimizationOutput
     {
         public int TotalLineSearchIterations { get; private set; }
-        public int IterationsWithNoLineSearch { get; private set; }
+        public int IterationsWithNonTrivialLineSearch { get; private set; }
 
-        public MinimizationWithLineSearchOutput(IEvaluation function_info, int iterations, int total_line_search_iterations, int iterations_with_no_line_search)
+        public MinimizationWithLineSearchOutput(IEvaluation function_info, int iterations, int total_line_search_iterations, int iterations_with_non_trivial_line_search)
             : base(function_info, iterations)
         {
             this.TotalLineSearchIterations = total_line_search_iterations;
-            this.IterationsWithNoLineSearch = iterations_with_no_line_search;
+            this.IterationsWithNonTrivialLineSearch = iterations_with_non_trivial_line_search;
         }
     }
 }

src/Numerics/Optimization/NewtonMinimizer.cs

@@ -0,0 +1,128 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using MathNet.Numerics.LinearAlgebra;
+using LU = MathNet.Numerics.LinearAlgebra.Factorization.LU<double>;
+
+namespace MathNet.Numerics.Optimization
+{
+    public class NewtonMinimizer
+    {
+        public double GradientTolerance { get; set; }
+        public int MaximumIterations { get; set; }
+        public bool UseLineSearch { get; set; }
+
+        public NewtonMinimizer(double gradient_tolerance, int maximum_iterations, bool use_line_search=false)
+        {
+            this.GradientTolerance = gradient_tolerance;
+            this.MaximumIterations = maximum_iterations;
+            this.UseLineSearch = use_line_search;
+        }
+
+        public MinimizationOutput FindMinimum(IObjectiveFunction objective, Vector<double> initial_guess)
+        {
+            if (!objective.GradientSupported)
+                throw new IncompatibleObjectiveException("Gradient not supported in objective function, but required for Newton minimization.");
+
+            if (!objective.HessianSupported)
+                throw new IncompatibleObjectiveException("Hessian not supported in objective function, but required for Newton minimization.");
+
+            if (!(objective is ObjectiveChecker))
+                objective = new ObjectiveChecker(objective, this.ValidateObjective, this.ValidateGradient, this.ValidateHessian);
+
+            IEvaluation initial_eval = objective.Evaluate(initial_guess);            
+            
+            // Check that we're not already done
+            if (this.ExitCriteriaSatisfied(initial_guess, initial_eval.Gradient))
+                return new MinimizationOutput(initial_eval, 0);
+            
+            // Set up line search algorithm            
+            var line_searcher = new WeakWolfeLineSearch(1e-4, 0.9, 1000);
+
+            // Declare state variables
+            IEvaluation candidate_point = initial_eval;            
+            Vector<double> search_direction;            
+            LineSearchOutput result;
+
+            // Subsequent steps
+            int iterations = 0;
+            int total_line_search_steps = 0;
+            int iterations_with_nontrivial_line_search = 0;
+            int steepest_descent_resets = 0;
+            bool tmp_line_search = false;
+            while (!this.ExitCriteriaSatisfied(candidate_point.Point, candidate_point.Gradient) && iterations < this.MaximumIterations)
+            {
+
+                search_direction = candidate_point.Hessian.LU().Solve(-candidate_point.Gradient);
+
+                if (search_direction * candidate_point.Gradient >= 0)
+                {
+                    search_direction = -candidate_point.Gradient;                    
+                    steepest_descent_resets += 1;
+                    tmp_line_search = true;
+                }
+
+                if (this.UseLineSearch || tmp_line_search)
+                {
+                    try
+                    {
+                        result = line_searcher.FindConformingStep(objective, candidate_point, search_direction, 1.0);
+                    }
+                    catch (Exception e)
+                    {
+                        throw new InnerOptimizationException("Line search failed.", e);
+                    }
+                    iterations_with_nontrivial_line_search += result.Iterations > 0 ? 1 : 0;
+                    total_line_search_steps += result.Iterations;
+                    candidate_point = result.FunctionInfoAtMinimum;
+                }
+                else
+                {
+                    candidate_point = objective.Evaluate(candidate_point.Point + search_direction);
+                }                
+                
+                tmp_line_search = false;
+
+                iterations += 1;
+            }
+
+            if (iterations == this.MaximumIterations)
+                throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", this.MaximumIterations));
+
+            return new MinimizationWithLineSearchOutput(candidate_point, iterations, total_line_search_steps, iterations_with_nontrivial_line_search);
+        }
+
+        private bool ExitCriteriaSatisfied(Vector<double> candidate_point, Vector<double> gradient)
+        {
+            return gradient.Norm(2.0) < this.GradientTolerance;
+        }
+
+        private void ValidateGradient(Vector<double> gradient, Vector<double> input)
+        {
+            foreach (var x in gradient)
+            {
+                if (Double.IsNaN(x) || Double.IsInfinity(x))
+                    throw new EvaluationException("Non-finite gradient returned.");
+            }
+        }
+
+        private void ValidateObjective(double objective, Vector<double> input)
+        {
+            if (Double.IsNaN(objective) || Double.IsInfinity(objective))
+                throw new EvaluationException("Non-finite objective function returned.");
+        }
+
+        private void ValidateHessian(Matrix<double> hessian, Vector<double> input)
+        {
+            for (int ii = 0; ii < hessian.RowCount; ++ii)
+            {
+                for (int jj = 0; jj < hessian.ColumnCount; ++jj)
+                {
+                    if (Double.IsNaN(hessian[ii,jj]) || Double.IsInfinity(hessian[ii,jj]))
+                        throw new EvaluationException("Non-finite Hessian returned.");
+                }
+            }
+        }
+    }
+}

src/Numerics/Optimization/WeakWolfeLineSearch.cs

@@ -23,6 +23,10 @@ namespace MathNet.Numerics.Optimization
         // Implemented following http://www.math.washington.edu/~burke/crs/408/lectures/L9-weak-Wolfe.pdf
         public LineSearchOutput FindConformingStep(IObjectiveFunction objective, IEvaluation starting_point, Vector<double> search_direction, double initial_step)
         {
+
+            if (!(objective is ObjectiveChecker))
+                objective = new ObjectiveChecker(objective, this.ValidateValue, this.ValidateGradient, null);
+
             double lower_bound = 0.0;
             double upper_bound = Double.PositiveInfinity;
             double step = initial_step;
@@ -71,5 +75,24 @@ namespace MathNet.Numerics.Optimization
             return step > 0 && sufficient_decrease && not_too_steep;
         }
 
+        private void ValidateValue(double value, Vector<double> input)
+        {
+            if (!this.IsFinite(value))
+                throw new EvaluationException(String.Format("Non-finite value returned by objective function: {0}", value));
+        }
+
+        private void ValidateGradient(Vector<double> gradient, Vector<double> input)
+        {
+            foreach (double x in gradient)
+                if (!this.IsFinite(x))
+                {
+                    throw new EvaluationException(String.Format("Non-finite value returned by gradient: {0}", x));
+                }
+        }
+
+        private bool IsFinite(double x)
+        {
+            return !(Double.IsNaN(x) || Double.IsInfinity(x));
+        }
     }
 }

src/UnitTests/OptimizationTests/RosenbrockFunction.cs

@@ -17,9 +17,20 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
         public static Vector<double> Gradient(Vector<double> input)
         {
             Vector<double> output = new MathNet.Numerics.LinearAlgebra.Double.DenseVector(2);
-            output[0] = -2 * (1 - input[0]) + 2 * 100 * (input[1] - input[0] * input[0]) * (-2 * input[0]);
+            output[0] = -2 * (1 - input[0]) + 200 * (input[1] - input[0] * input[0]) * (-2 * input[0]);
             output[1] = 2 * 100 * (input[1] - input[0] * input[0]);
             return output;
         }
+
+        public static Matrix<double> Hessian(Vector<double> input)
+        {
+
+            Matrix<double> output = new MathNet.Numerics.LinearAlgebra.Double.DenseMatrix(2,2);
+            output[0, 0] = 2 - 400 * input[1] + 1200 * input[0] * input[0];
+            output[1, 1] = 200;
+            output[0, 1] = -400 * input[0];
+            output[1, 0] = output[0, 1];
+            return output;
+        }
     }
 }

src/UnitTests/OptimizationTests/TestBfgsMinimizer.cs

@@ -0,0 +1,50 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+
+using NUnit.Framework;
+
+using MathNet.Numerics.Optimization;
+
+namespace MathNet.Numerics.UnitTests.OptimizationTests
+{
+    [TestFixture]
+    public class TestBfgsMinimizer
+    {
+
+        [Test]
+        public void FindMinimum_Rosenbrock_Easy()
+        {
+            var obj = new SimpleObjectiveFunction(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
+            var solver = new BfgsMinimizer(1e-5, 1000);
+            var result = solver.FindMinimum(obj, new MathNet.Numerics.LinearAlgebra.Double.DenseVector(new double[] { 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 = new SimpleObjectiveFunction(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
+            var solver = new BfgsMinimizer(1e-5, 1000);
+            var result = solver.FindMinimum(obj, new MathNet.Numerics.LinearAlgebra.Double.DenseVector(new double[] { -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 = new SimpleObjectiveFunction(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
+            var solver = new BfgsMinimizer(1e-5, 1000);
+            var result = solver.FindMinimum(obj, new MathNet.Numerics.LinearAlgebra.Double.DenseVector(new double[] { -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/TestNewtonMinimizer.cs

@@ -0,0 +1,82 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+
+using NUnit.Framework;
+
+using MathNet.Numerics.Optimization;
+
+namespace MathNet.Numerics.UnitTests.OptimizationTests
+{
+    [TestFixture]
+    public class TestNewtonMinimizer
+    {
+
+        [Test]
+        public void FindMinimum_Rosenbrock_Easy()
+        {
+            var obj = new SimpleObjectiveFunction(RosenbrockFunction.Value, RosenbrockFunction.Gradient, RosenbrockFunction.Hessian);
+            var solver = new NewtonMinimizer(1e-5, 1000);
+            var result = solver.FindMinimum(obj, new MathNet.Numerics.LinearAlgebra.Double.DenseVector(new double[] { 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 = new SimpleObjectiveFunction(RosenbrockFunction.Value, RosenbrockFunction.Gradient, RosenbrockFunction.Hessian);
+            var solver = new NewtonMinimizer(1e-5, 1000);
+            var result = solver.FindMinimum(obj, new MathNet.Numerics.LinearAlgebra.Double.DenseVector(new double[] { -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 = new SimpleObjectiveFunction(RosenbrockFunction.Value, RosenbrockFunction.Gradient, RosenbrockFunction.Hessian);
+            var solver = new NewtonMinimizer(1e-5, 1000);
+            var result = solver.FindMinimum(obj, new MathNet.Numerics.LinearAlgebra.Double.DenseVector(new double[] { -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));
+        }
+
+        [Test]
+        public void FindMinimum_Linesearch_Rosenbrock_Easy()
+        {
+            var obj = new SimpleObjectiveFunction(RosenbrockFunction.Value, RosenbrockFunction.Gradient, RosenbrockFunction.Hessian);
+            var solver = new NewtonMinimizer(1e-5, 1000, true);
+            var result = solver.FindMinimum(obj, new MathNet.Numerics.LinearAlgebra.Double.DenseVector(new double[] { 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_Linesearch_Rosenbrock_Hard()
+        {
+            var obj = new SimpleObjectiveFunction(RosenbrockFunction.Value, RosenbrockFunction.Gradient, RosenbrockFunction.Hessian);
+            var solver = new NewtonMinimizer(1e-5, 1000, true);
+            var result = solver.FindMinimum(obj, new MathNet.Numerics.LinearAlgebra.Double.DenseVector(new double[] { -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_Linesearch_Rosenbrock_Overton()
+        {
+            var obj = new SimpleObjectiveFunction(RosenbrockFunction.Value, RosenbrockFunction.Gradient, RosenbrockFunction.Hessian);
+            var solver = new NewtonMinimizer(1e-5, 1000, true);
+            var result = solver.FindMinimum(obj, new MathNet.Numerics.LinearAlgebra.Double.DenseVector(new double[] { -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,59 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using NUnit.Framework;
+namespace MathNet.Numerics.UnitTests.OptimizationTests
+{
+    [TestFixture]
+    class TestRosenbrockFunction
+    {
+        [Test]
+        public void TestGradient()
+        {
+            var input = new LinearAlgebra.Double.DenseVector(new double[]{ -0.9, -0.5 } );
+
+            var v1 = RosenbrockFunction.Value(input);
+            var g = RosenbrockFunction.Gradient(input);
+
+            var eps = 1e-5;
+            var eps0 = (new LinearAlgebra.Double.DenseVector(new double[] { 1.0, 0.0 })) * eps;
+            var eps1 = (new LinearAlgebra.Double.DenseVector(new double[] { 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 LinearAlgebra.Double.DenseVector(new double[] { -0.9, -0.5 });
+
+            var v1 = RosenbrockFunction.Value(input);
+            var h = RosenbrockFunction.Hessian(input);
+
+            var eps = 1e-5;
+
+            var eps0 = (new LinearAlgebra.Double.DenseVector(new double[] { 1.0, 0.0 })) * eps;
+            var eps1 = (new LinearAlgebra.Double.DenseVector(new double[] { 0.0, 1.0 })) * eps;
+
+            var epsuu = (new LinearAlgebra.Double.DenseVector(new double[] { 1.0, 1.0 })) * eps;
+            var epsud = (new LinearAlgebra.Double.DenseVector(new double[] { 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

@@ -361,8 +361,11 @@
     <Compile Include="Random\SystemRandomSourceTests.cs" />
     <Compile Include="RootFindingTests\BisectionTest.cs" />
     <Compile Include="OptimizationTests\RosenbrockFunction.cs" />
+    <Compile Include="OptimizationTests\TestBfgsMinimizer.cs" />
     <Compile Include="OptimizationTests\TestBisectionRootFinder.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" />