diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj
index 27b0246a..2ee74519 100644
--- a/src/Numerics/Numerics.csproj
+++ b/src/Numerics/Numerics.csproj
@@ -262,7 +262,7 @@
-
+
diff --git a/src/Numerics/Optimization/BFGS.cs b/src/Numerics/Optimization/BFGS.cs
deleted file mode 100644
index 7b9942d2..00000000
--- a/src/Numerics/Optimization/BFGS.cs
+++ /dev/null
@@ -1,11 +0,0 @@
-using System;
-using System.Collections.Generic;
-using System.Linq;
-using System.Text;
-
-namespace MathNet.Numerics.Optimization
-{
- class BFGS
- {
- }
-}
diff --git a/src/Numerics/Optimization/BfgsMinimizer.cs b/src/Numerics/Optimization/BfgsMinimizer.cs
new file mode 100644
index 00000000..09b2264f
--- /dev/null
+++ b/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 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(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(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 candidatePoint, Vector 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);
+ }
+ }
+}
diff --git a/src/Numerics/Optimization/ConjugateGradientMinimizer.cs b/src/Numerics/Optimization/ConjugateGradientMinimizer.cs
index 707459d6..f33a9043 100644
--- a/src/Numerics/Optimization/ConjugateGradientMinimizer.cs
+++ b/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 candidatePoint, Vector gradient)
diff --git a/src/UnitTests/OptimizationTests/TestBfgsMinimizer.cs b/src/UnitTests/OptimizationTests/TestBfgsMinimizer.cs
new file mode 100644
index 00000000..851d30ee
--- /dev/null
+++ b/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));
+ }
+ }
+}
diff --git a/src/UnitTests/OptimizationTests/TestRosenbrockFunction.cs b/src/UnitTests/OptimizationTests/TestRosenbrockFunction.cs
new file mode 100644
index 00000000..14a0f9f7
--- /dev/null
+++ b/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);
+ }
+ }
+}
diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj
index 3f1241dd..16a7a82f 100644
--- a/src/UnitTests/UnitTests.csproj
+++ b/src/UnitTests/UnitTests.csproj
@@ -372,7 +372,10 @@
+
+
+