From fc76933ceae566aad01de649ba397ee15ad641a0 Mon Sep 17 00:00:00 2001 From: Scott Stephens Date: Wed, 12 Oct 2016 18:04:28 -0500 Subject: [PATCH] Optimization: add finite difference objective function adapter --- src/Numerics/Numerics.csproj | 1 + ...wardDifferenceGradientObjectiveFunction.cs | 151 ++++++++++++++++++ .../OptimizationTests/BfgsBMinimizerTests.cs | 49 +++++- 3 files changed, 199 insertions(+), 2 deletions(-) create mode 100644 src/Numerics/Optimization/ObjectiveFunctions/ForwardDifferenceGradientObjectiveFunction.cs diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 9c89a636..00b170d5 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -120,6 +120,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") { } + } } }