Optimization: add finite difference objective function adapter

src/Numerics/Numerics.csproj

@@ -120,6 +120,7 @@
     <Compile Include="Optimization\BfgsMinimizerBase.cs" />
     <Compile Include="Optimization\LineSearch\WolfeLineSearch.cs" />
     <Compile Include="Optimization\NelderMeadSimplex.cs" />
+    <Compile Include="Optimization\ObjectiveFunctions\ForwardDifferenceGradientObjectiveFunction.cs" />
     <Compile Include="Optimization\ObjectiveFunctions\LazyObjectiveFunctionBase.cs" />
     <Compile Include="Optimization\Exceptions.cs" />
     <Compile Include="Optimization\ObjectiveFunctions\ObjectiveFunctionBase.cs" />

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
+{
+    /// <summary>
+    /// 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.
+    /// </summary>
+    public class ForwardDifferenceGradientObjectiveFunction : IObjectiveFunction
+    {
+        public IObjectiveFunction InnerObjectiveFunction { get; protected set; }
+        protected Vector<double> _lower_bound;
+        protected Vector<double> _upper_bound;
+
+        protected Vector<double> _point;
+        protected bool _value_evaluated = false;
+        protected bool _gradient_evaluated = false;
+        protected Vector<double> _gradient;
+
+        public double MinimumIncrement;
+        public double RelativeIncrement;
+        
+        public ForwardDifferenceGradientObjectiveFunction(IObjectiveFunction valueOnlyObj, Vector<double> lowerBound, Vector<double> 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<double> Gradient
+        {
+            get
+            {
+                if (!_gradient_evaluated)
+                    this.EvaluateGradient();
+                return _gradient;
+            }
+        }
+
+        public Matrix<double> Hessian
+        {
+            get
+            {
+                throw new NotImplementedException();
+            }
+        }
+
+        public bool IsGradientSupported
+        {
+            get
+            {
+                return true;
+            }
+        }
+
+        public bool IsHessianSupported
+        {
+            get
+            {
+                return false;
+            }
+        }
+
+        public Vector<double> 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<double> 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;
+        }
+    }
+}

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<ITestCaseData>
+        [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<ITestCaseData>
         {
+            private string _prefix = "";
+
+            public BaseMghTestCaseEnumerator(string prefix)
+            {
+                if (prefix.EndsWith(" "))
+                    _prefix = prefix;
+                else
+                    _prefix = prefix + " ";
+            }
             public IEnumerator<ITestCaseData> 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") { }
+        }
     }
 }