diff --git a/src/Numerics.Tests/OptimizationTests/NewtonMinimizerTests.cs b/src/Numerics.Tests/OptimizationTests/NewtonMinimizerTests.cs index 602f2d66c..bb20a2927 100644 --- a/src/Numerics.Tests/OptimizationTests/NewtonMinimizerTests.cs +++ b/src/Numerics.Tests/OptimizationTests/NewtonMinimizerTests.cs @@ -84,6 +84,32 @@ protected override void Evaluate() } } + public class LazySixHumpCamelObjectiveFunction : LazyObjectiveFunctionBase + { + public LazySixHumpCamelObjectiveFunction() : base(true, true) { } + + public override IObjectiveFunction CreateNew() + { + return new LazySixHumpCamelObjectiveFunction(); + } + + protected override void EvaluateValue() + { + Value = SixHumpCamelFunction.Value(Point); + } + + protected override void EvaluateGradient() + { + Gradient = SixHumpCamelFunction.Gradient(Point); + } + + protected override void EvaluateHessian() + { + Hessian = SixHumpCamelFunction.Hessian(Point); + } + + } + [TestFixture] public class NewtonMinimizerTests { @@ -153,6 +179,17 @@ public void FindMinimum_Linesearch_Rosenbrock_Overton() Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3)); } + [Test] + public void FindMinimum_SixHumpCamel_IndefiniteHessian() + { + var obj = new LazySixHumpCamelObjectiveFunction(); + var solver = new NewtonMinimizer(1e-5, 1000, true, HessianModifiers.ReverseNegativeEigenvalues); + var result = solver.FindMinimum(obj, new DenseVector(new double[] { 1.0, -0.6 })); + + Assert.That(result.MinimizingPoint[0], Is.EqualTo(0.0898).Within(1e-3)); + Assert.That(result.MinimizingPoint[1], Is.EqualTo(-0.7126).Within(1e-3)); + } + private class MghTestCaseEnumerator : IEnumerable { private static readonly string[] _ignore_list = diff --git a/src/Numerics.Tests/OptimizationTests/TestFunctions/SixHumpCamelFunction.cs b/src/Numerics.Tests/OptimizationTests/TestFunctions/SixHumpCamelFunction.cs new file mode 100644 index 000000000..1268d02f7 --- /dev/null +++ b/src/Numerics.Tests/OptimizationTests/TestFunctions/SixHumpCamelFunction.cs @@ -0,0 +1,95 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// https://numerics.mathdotnet.com +// https://github.com/mathnet/mathnet-numerics +// +// Copyright (c) 2009-$CURRENT_YEAR$ Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.LinearAlgebra.Double; + +namespace MathNet.Numerics.Tests.OptimizationTests.TestFunctions +{ + /// + /// Six-Hump Camel Function, see http://www.sfu.ca/~ssurjano/camel6.html for formula and global minimum locations + /// + public static class SixHumpCamelFunction + { + public static double Value(Vector input) + { + double x = input[0]; + double y = input[1]; + + double x2 = x * x; + double x4 = x * x * x * x; + double y2 = y * y; + + return x2 * (4 - 2.1 * x2 + x4 / 3) + x * y + 4 * y2 * (y2 - 1); + } + + public static Vector Gradient(Vector input) + { + double x = input[0]; + double y = input[1]; + + double x3 = x * x * x; + double x5 = x * x * x * x * x; + double y2 = y * y; + double y3 = y * y * y; + + Vector output = new DenseVector(2); + + output[0] = 2 * (4 * x - 4.2 * x3 + x5 + 0.5 * y); + output[1] = x - 8 * y + 16 * y3; + return output; + } + + public static Matrix Hessian(Vector input) + { + double x = input[0]; + double y = input[1]; + + double x2 = x * x; + double x4 = x * x * x * x; + double y2 = y * y; + + + Matrix output = new DenseMatrix(2, 2); + output[0, 0] = 10 * (0.8 - 2.52 * x2 + x4); + output[1, 1] = 48 * y2 - 8; + output[0, 1] = 1; + output[1, 0] = output[0, 1]; + return output; + } + + public static Vector Minimum + { + get + { + return new DenseVector(new double[] { 0.0898, -0.7126 }); + } + } + } +} diff --git a/src/Numerics/Optimization/HessianModifiers.cs b/src/Numerics/Optimization/HessianModifiers.cs new file mode 100644 index 000000000..1458a8d76 --- /dev/null +++ b/src/Numerics/Optimization/HessianModifiers.cs @@ -0,0 +1,43 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// https://numerics.mathdotnet.com +// https://github.com/mathnet/mathnet-numerics +// +// Copyright (c) 2009-$CURRENT_YEAR$ Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; + +namespace MathNet.Numerics.Optimization +{ + public enum HessianModifiers + { + None, + ReverseNegativeEigenvalues + } +} diff --git a/src/Numerics/Optimization/NewtonMinimizer.cs b/src/Numerics/Optimization/NewtonMinimizer.cs index e8fd26fa7..ac836ebf2 100644 --- a/src/Numerics/Optimization/NewtonMinimizer.cs +++ b/src/Numerics/Optimization/NewtonMinimizer.cs @@ -30,6 +30,7 @@ using System; using MathNet.Numerics.LinearAlgebra; using MathNet.Numerics.Optimization.LineSearch; +using MathNet.Numerics.LinearAlgebra.Factorization; namespace MathNet.Numerics.Optimization { @@ -38,20 +39,22 @@ public sealed class NewtonMinimizer : IUnconstrainedMinimizer public double GradientTolerance { get; set; } public int MaximumIterations { get; set; } public bool UseLineSearch { get; set; } + public HessianModifiers Modifier { get; set; } - public NewtonMinimizer(double gradientTolerance, int maximumIterations, bool useLineSearch = false) + public NewtonMinimizer(double gradientTolerance, int maximumIterations, bool useLineSearch = false, HessianModifiers modifier = HessianModifiers.None) { GradientTolerance = gradientTolerance; MaximumIterations = maximumIterations; UseLineSearch = useLineSearch; + Modifier = modifier; } public MinimizationResult FindMinimum(IObjectiveFunction objective, Vector initialGuess) { - return Minimum(objective, initialGuess, GradientTolerance, MaximumIterations, UseLineSearch); + return Minimum(objective, initialGuess, GradientTolerance, MaximumIterations, UseLineSearch, Modifier); } - public static MinimizationResult Minimum(IObjectiveFunction objective, Vector initialGuess, double gradientTolerance=1e-8, int maxIterations=1000, bool useLineSearch=false) + public static MinimizationResult Minimum(IObjectiveFunction objective, Vector initialGuess, double gradientTolerance=1e-8, int maxIterations=1000, bool useLineSearch=false, HessianModifiers modifier=HessianModifiers.None) { if (!objective.IsGradientSupported) { @@ -83,7 +86,7 @@ public static MinimizationResult Minimum(IObjectiveFunction objective, Vector= 0) { searchDirection = -objective.Gradient; @@ -125,6 +128,43 @@ public static MinimizationResult Minimum(IObjectiveFunction objective, Vector CalculateSearchDirection(IObjectiveFunction objective, HessianModifiers modifier) + { + Vector searchDirection = null; + switch (modifier) + { + case HessianModifiers.None: + searchDirection = SolveLU(objective); + break; + case HessianModifiers.ReverseNegativeEigenvalues: + searchDirection = ReverseNegativeEigenvaluesAndSolve(objective); + break; + } + + return searchDirection; + } + + static Vector SolveLU(IObjectiveFunction objective) + { + return objective.Hessian.LU().Solve(-objective.Gradient); + } + + /// + /// Force the Hessian to be positive definite by reversing the negative eigenvalues. Use the EVD decomposition to then solve for the search direction. + /// Implementation based on Philip E. Gill, Walter Murray, and Margaret H. Wright, Practical Optimization, 1981, 107–8. + /// + /// + /// + static Vector ReverseNegativeEigenvaluesAndSolve(IObjectiveFunction objective) + { + Evd evd = objective.Hessian.Evd(Symmetricity.Symmetric); + for (int i = 0; i < evd.EigenValues.Count; i++) + { + evd.EigenValues[i] = Math.Max(Math.Abs(evd.EigenValues[i].Real), double.Epsilon); + } + return evd.Solve(-objective.Gradient); + } + static void ValidateGradient(IObjectiveFunctionEvaluation eval) { foreach (var x in eval.Gradient)