diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj
index b310fd3c..c572888b 100644
--- a/src/Numerics/Numerics.csproj
+++ b/src/Numerics/Numerics.csproj
@@ -105,6 +105,7 @@
+
diff --git a/src/Numerics/Optimization/ObjectiveFunctions/ForwardDifferenceGradientObjectiveFunction.cs b/src/Numerics/Optimization/ObjectiveFunctions/ForwardDifferenceGradientObjectiveFunction.cs
new file mode 100644
index 00000000..1a3733bf
--- /dev/null
+++ b/src/Numerics/Optimization/ObjectiveFunctions/ForwardDifferenceGradientObjectiveFunction.cs
@@ -0,0 +1,151 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using MathNet.Numerics.LinearAlgebra;
+
+namespace MathNet.Numerics.Optimization.ObjectiveFunctions
+{
+ ///
+ /// Adapts an objective function with only value implemented
+ /// to provide a gradient as well. Gradient calculation is
+ /// done using the finite difference method, specifically
+ /// forward differences.
+ ///
+ /// For each gradient computed, the algorithm requires an
+ /// additional number of function evaluations equal to the
+ /// functions's number of input parameters.
+ ///
+ public class ForwardDifferenceGradientObjectiveFunction : IObjectiveFunction
+ {
+ public IObjectiveFunction InnerObjectiveFunction { get; protected set; }
+ protected Vector _lower_bound;
+ protected Vector _upper_bound;
+
+ protected Vector _point;
+ protected bool _value_evaluated = false;
+ protected bool _gradient_evaluated = false;
+ protected Vector _gradient;
+
+ public double MinimumIncrement;
+ public double RelativeIncrement;
+
+ public ForwardDifferenceGradientObjectiveFunction(IObjectiveFunction valueOnlyObj, Vector lowerBound, Vector upperBound, double relativeIncrement=1e-5, double minimumIncrement=1e-8)
+ {
+ this.InnerObjectiveFunction = valueOnlyObj;
+ _lower_bound = lowerBound;
+ _upper_bound = upperBound;
+ _gradient = new LinearAlgebra.Double.DenseVector(_lower_bound.Count);
+ this.RelativeIncrement = relativeIncrement;
+ this.MinimumIncrement = minimumIncrement;
+ }
+
+ protected void EvaluateValue()
+ {
+ _value_evaluated = true;
+ }
+
+ protected void EvaluateGradient()
+ {
+ if (!_value_evaluated)
+ this.EvaluateValue();
+
+ var tmp_point = _point.Clone();
+ var tmp_obj = this.InnerObjectiveFunction.CreateNew();
+ for (int ii = 0; ii < _gradient.Count; ++ii)
+ {
+ var orig_point = tmp_point[ii];
+ var rel_incr = orig_point * this.RelativeIncrement;
+ var h = Math.Max(rel_incr, this.MinimumIncrement);
+ var mult = 1;
+ if (orig_point + h > _upper_bound[ii])
+ mult = -1;
+
+ tmp_point[ii] = orig_point + mult*h;
+ tmp_obj.EvaluateAt(tmp_point);
+ double bumped_value = tmp_obj.Value;
+ _gradient[ii] = (mult * bumped_value - mult * this.InnerObjectiveFunction.Value) / h;
+
+ tmp_point[ii] = orig_point;
+ }
+ _gradient_evaluated = true;
+ }
+
+ public Vector Gradient
+ {
+ get
+ {
+ if (!_gradient_evaluated)
+ this.EvaluateGradient();
+ return _gradient;
+ }
+ }
+
+ public Matrix Hessian
+ {
+ get
+ {
+ throw new NotImplementedException();
+ }
+ }
+
+ public bool IsGradientSupported
+ {
+ get
+ {
+ return true;
+ }
+ }
+
+ public bool IsHessianSupported
+ {
+ get
+ {
+ return false;
+ }
+ }
+
+ public Vector Point
+ {
+ get
+ {
+ return _point;
+ }
+ }
+
+ public double Value
+ {
+ get
+ {
+ if (!_value_evaluated)
+ this.EvaluateValue();
+ return this.InnerObjectiveFunction.Value;
+ }
+ }
+
+ public IObjectiveFunction CreateNew()
+ {
+ var tmp = new ForwardDifferenceGradientObjectiveFunction(this.InnerObjectiveFunction.CreateNew(), _lower_bound, _upper_bound, this.RelativeIncrement, this.MinimumIncrement);
+ return tmp;
+ }
+
+ public void EvaluateAt(Vector point)
+ {
+ _point = point;
+ _value_evaluated = false;
+ _gradient_evaluated = false;
+ this.InnerObjectiveFunction.EvaluateAt(point);
+ }
+
+ public IObjectiveFunction Fork()
+ {
+ var tmp = new ForwardDifferenceGradientObjectiveFunction(this.InnerObjectiveFunction.Fork(), _lower_bound, _upper_bound, this.RelativeIncrement, this.MinimumIncrement);
+ tmp._point = _point?.Clone();
+ tmp._gradient_evaluated = _gradient_evaluated;
+ tmp._value_evaluated = _value_evaluated;
+ tmp._gradient = _gradient?.Clone();
+
+ return tmp;
+ }
+ }
+}
diff --git a/src/UnitTests/OptimizationTests/BfgsBMinimizerTests.cs b/src/UnitTests/OptimizationTests/BfgsBMinimizerTests.cs
index 785ddb81..33125613 100644
--- a/src/UnitTests/OptimizationTests/BfgsBMinimizerTests.cs
+++ b/src/UnitTests/OptimizationTests/BfgsBMinimizerTests.cs
@@ -37,6 +37,7 @@ using System.Text;
using System.Collections.Generic;
using MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions;
using System.Collections;
+using MathNet.Numerics.Optimization.ObjectiveFunctions;
namespace MathNet.Numerics.UnitTests.OptimizationTests
{
@@ -178,8 +179,42 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
Assert.That(success, "Minimal function value is not as expected.");
}
- private class MghTestCaseEnumerator : IEnumerable
+ [Test]
+ [TestCaseSource(typeof(FdMghTestCaseEnumerator))]
+ public void Mgh_FiniteDifference_Tests(TestFunctions.TestCase test_case)
+ {
+ var obj1 = new MghObjectiveFunction(test_case.Function, true, true);
+ var obj = new ForwardDifferenceGradientObjectiveFunction(obj1, test_case.LowerBound, test_case.UpperBound, 1e-10, 1e-10);
+ var solver = new BfgsBMinimizer(1e-8, 1e-8, 1e-8, 1000);
+
+ var result = solver.FindMinimum(obj, test_case.LowerBound, test_case.UpperBound, test_case.InitialGuess);
+
+ if (test_case.MinimizingPoint != null)
+ {
+ Assert.That((result.MinimizingPoint - test_case.MinimizingPoint).L2Norm(), Is.LessThan(1e-3));
+ }
+
+ var val1 = result.FunctionInfoAtMinimum.Value;
+ var val2 = test_case.MinimalValue;
+ var abs_min = Math.Min(Math.Abs(val1), Math.Abs(val2));
+ var abs_err = Math.Abs(val1 - val2);
+ var rel_err = abs_err / abs_min;
+ var success = (abs_min <= 1 && abs_err < 1e-3) || (abs_min > 1 && rel_err < 1e-3);
+ Assert.That(success, "Minimal function value is not as expected.");
+ }
+
+
+ private class BaseMghTestCaseEnumerator : IEnumerable
{
+ private string _prefix = "";
+
+ public BaseMghTestCaseEnumerator(string prefix)
+ {
+ if (prefix.EndsWith(" "))
+ _prefix = prefix;
+ else
+ _prefix = prefix + " ";
+ }
public IEnumerator GetEnumerator()
{
return
@@ -192,7 +227,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
.Concat(BrownAndDennisFunction.TestCases)
.Where(x => x.IsBounded)
.Select(x => new TestCaseData(x)
- .SetName(x.FullName)
+ .SetName(_prefix + x.FullName)
)
.GetEnumerator();
}
@@ -202,6 +237,16 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
return this.GetEnumerator();
}
}
+
+ private class MghTestCaseEnumerator : BaseMghTestCaseEnumerator
+ {
+ public MghTestCaseEnumerator() : base("") { }
+ }
+
+ private class FdMghTestCaseEnumerator : BaseMghTestCaseEnumerator
+ {
+ public FdMghTestCaseEnumerator() : base("FD") { }
+ }
}
}