From 94a9e26f8c3690d914ae4542f245d88208979cc0 Mon Sep 17 00:00:00 2001 From: strMikhailPotapenko Date: Wed, 28 Feb 2024 15:59:33 -0500 Subject: [PATCH 1/8] Created unit test for NewtonMinimizer on Beale's function which fails to converge in 1 under 1000 steps --- .../OptimizationTests/NewtonMinimizerTests.cs | 55 +++++++++++ .../TestFunctions/BealeFunction2D.cs | 98 +++++++++++++++++++ 2 files changed, 153 insertions(+) create mode 100644 src/Numerics.Tests/OptimizationTests/TestFunctions/BealeFunction2D.cs diff --git a/src/Numerics.Tests/OptimizationTests/NewtonMinimizerTests.cs b/src/Numerics.Tests/OptimizationTests/NewtonMinimizerTests.cs index 602f2d66c..8480b6428 100644 --- a/src/Numerics.Tests/OptimizationTests/NewtonMinimizerTests.cs +++ b/src/Numerics.Tests/OptimizationTests/NewtonMinimizerTests.cs @@ -84,6 +84,50 @@ protected override void Evaluate() } } + public class LazyBealebjectiveFunction : LazyObjectiveFunctionBase + { + public LazyBealebjectiveFunction() : base(true, true) { } + + public override IObjectiveFunction CreateNew() + { + return new LazyBealebjectiveFunction(); + } + + protected override void EvaluateValue() + { + Value = BealeFunction2D.Value(Point); + } + + protected override void EvaluateGradient() + { + Gradient = BealeFunction2D.Gradient(Point); + } + + protected override void EvaluateHessian() + { + Hessian = BealeFunction2D.Hessian(Point); + } + } + + public class BealeObjectiveFunction : ObjectiveFunctionBase + { + public BealeObjectiveFunction() : base(true, true) { } + + public override IObjectiveFunction CreateNew() + { + return new BealeObjectiveFunction(); + } + + protected override void Evaluate() + { + // here we could directly overwrite the existing matrix cells instead. + // note: values must then be initialized manually first, if null. + Value = BealeFunction2D.Value(Point); ; + Gradient = BealeFunction2D.Gradient(Point); + Hessian = BealeFunction2D.Hessian(Point); + } + } + [TestFixture] public class NewtonMinimizerTests { @@ -153,6 +197,17 @@ public void FindMinimum_Linesearch_Rosenbrock_Overton() Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3)); } + [Test] + public void FindMinimum_Beale_IndefiniteHessian() + { + var obj = new LazyBealebjectiveFunction(); + var solver = new NewtonMinimizer(1e-5, 1000, true); + var result = solver.FindMinimum(obj, new DenseVector(new double[] { 4, 0.9 })); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - 3.0), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - 0.5), Is.LessThan(1e-3)); + } + private class MghTestCaseEnumerator : IEnumerable { private static readonly string[] _ignore_list = diff --git a/src/Numerics.Tests/OptimizationTests/TestFunctions/BealeFunction2D.cs b/src/Numerics.Tests/OptimizationTests/TestFunctions/BealeFunction2D.cs new file mode 100644 index 000000000..28b9f1897 --- /dev/null +++ b/src/Numerics.Tests/OptimizationTests/TestFunctions/BealeFunction2D.cs @@ -0,0 +1,98 @@ +// +// 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 MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.LinearAlgebra.Double; + +namespace MathNet.Numerics.Tests.OptimizationTests.TestFunctions +{ + public static class BealeFunction2D + { + public static double Value(Vector input) + { + double x = input[0]; + double y = input[1]; + + double a = 1.5 - x + x * y; + double b = 2.25 - x + x * y * y; + double c = 2.625 - x + x * Math.Pow(y, 3); + + return Math.Pow(a, 2) + Math.Pow(b, 2) + Math.Pow(c, 2); + } + + public static Vector Gradient(Vector input) + { + double x = input[0]; + double y = input[1]; + + double a = 1.5 - x + x * y; + double b = 2.25 - x + x * y * y; + double c = 2.625 - x + x * Math.Pow(y, 3); + + Vector output = new DenseVector(2); + + output[0] = 2 * (y - 1) * a + 2 * (y * y - 1) * b + 2 * (Math.Pow(y, 3) - 1) * c; + output[1] = 2 * x * a + 4 * x * y * b + 6 * x * y * y * c; + return output; + } + + public static Matrix Hessian(Vector input) + { + double x = input[0]; + double y = input[1]; + + double a = 1.5 - x + x * y; + double axprime = y - 1; + + double b = 2.25 - x + x * y * y; + double bxprime = y * y - 1; + + double c = 2.625 - x + x * Math.Pow(y, 3); + double cxprime = Math.Pow(y, 3) - 1; + + + Matrix output = new DenseMatrix(2, 2); + output[0, 0] = 2 * (Math.Pow(axprime, 2) + Math.Pow(bxprime, 2) + Math.Pow(cxprime, 2)); + output[1, 1] = 2 * Math.Pow(x, 2) + 8 * Math.Pow(x, 2) * Math.Pow(y, 2) + 18 * Math.Pow(x, 2) * Math.Pow(y, 4) + 4 * x * b + 12 * x * y * c; + output[0, 1] = 2 * x * axprime + 2 * a + 4 * x * y * bxprime + 4 * y * b + 6 * x * Math.Pow(y, 2) * cxprime + 6 * Math.Pow(y, 2) * c; + output[1, 0] = output[0, 1]; + return output; + } + + public static Vector Minimum + { + get + { + return new DenseVector(new double[] { 3, 0.5 }); + } + } + } +} + From c091e71fe0102ccbbac69bd60ad307f144f92220 Mon Sep 17 00:00:00 2001 From: strMikhailPotapenko Date: Thu, 29 Feb 2024 15:03:36 -0500 Subject: [PATCH 2/8] Added Six Hump Camel Function to set of test functions, removed 2D Beale Function --- ...eFunction2D.cs => SixHumpCamelFunction.cs} | 44 ++++++++----------- 1 file changed, 19 insertions(+), 25 deletions(-) rename src/Numerics.Tests/OptimizationTests/TestFunctions/{BealeFunction2D.cs => SixHumpCamelFunction.cs} (62%) diff --git a/src/Numerics.Tests/OptimizationTests/TestFunctions/BealeFunction2D.cs b/src/Numerics.Tests/OptimizationTests/TestFunctions/SixHumpCamelFunction.cs similarity index 62% rename from src/Numerics.Tests/OptimizationTests/TestFunctions/BealeFunction2D.cs rename to src/Numerics.Tests/OptimizationTests/TestFunctions/SixHumpCamelFunction.cs index 28b9f1897..9bc223261 100644 --- a/src/Numerics.Tests/OptimizationTests/TestFunctions/BealeFunction2D.cs +++ b/src/Numerics.Tests/OptimizationTests/TestFunctions/SixHumpCamelFunction.cs @@ -1,4 +1,4 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // https://numerics.mathdotnet.com // https://github.com/mathnet/mathnet-numerics @@ -27,24 +27,23 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using System; using MathNet.Numerics.LinearAlgebra; using MathNet.Numerics.LinearAlgebra.Double; namespace MathNet.Numerics.Tests.OptimizationTests.TestFunctions { - public static class BealeFunction2D + public static class SixHumpCamelFunction { public static double Value(Vector input) { double x = input[0]; double y = input[1]; - double a = 1.5 - x + x * y; - double b = 2.25 - x + x * y * y; - double c = 2.625 - x + x * Math.Pow(y, 3); + double x2 = x * x; + double x4 = x * x * x * x; + double y2 = y * y; - return Math.Pow(a, 2) + Math.Pow(b, 2) + Math.Pow(c, 2); + return x2 * (4 - 2.1 * x2 + x4 / 3) + x * y + 4 * y2 * (y2 - 1); } public static Vector Gradient(Vector input) @@ -52,14 +51,15 @@ public static Vector Gradient(Vector input) double x = input[0]; double y = input[1]; - double a = 1.5 - x + x * y; - double b = 2.25 - x + x * y * y; - double c = 2.625 - x + x * Math.Pow(y, 3); + 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 * (y - 1) * a + 2 * (y * y - 1) * b + 2 * (Math.Pow(y, 3) - 1) * c; - output[1] = 2 * x * a + 4 * x * y * b + 6 * x * y * y * c; + output[0] = 2 * (4 * x - 4.2 * x3 + x5 + 0.5 * y); + output[1] = x - 8 * y + 16 * y3; return output; } @@ -68,20 +68,15 @@ public static Matrix Hessian(Vector input) double x = input[0]; double y = input[1]; - double a = 1.5 - x + x * y; - double axprime = y - 1; - - double b = 2.25 - x + x * y * y; - double bxprime = y * y - 1; - - double c = 2.625 - x + x * Math.Pow(y, 3); - double cxprime = Math.Pow(y, 3) - 1; + double x2 = x * x; + double x4 = x * x * x * x; + double y2 = y * y; Matrix output = new DenseMatrix(2, 2); - output[0, 0] = 2 * (Math.Pow(axprime, 2) + Math.Pow(bxprime, 2) + Math.Pow(cxprime, 2)); - output[1, 1] = 2 * Math.Pow(x, 2) + 8 * Math.Pow(x, 2) * Math.Pow(y, 2) + 18 * Math.Pow(x, 2) * Math.Pow(y, 4) + 4 * x * b + 12 * x * y * c; - output[0, 1] = 2 * x * axprime + 2 * a + 4 * x * y * bxprime + 4 * y * b + 6 * x * Math.Pow(y, 2) * cxprime + 6 * Math.Pow(y, 2) * c; + 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; } @@ -90,9 +85,8 @@ public static Vector Minimum { get { - return new DenseVector(new double[] { 3, 0.5 }); + return new DenseVector(new double[] { 0.0898, -0.7126 }); } } } } - From 8ff2f03735edc5493c9c3e623c52edd6e0d46f7f Mon Sep 17 00:00:00 2001 From: strMikhailPotapenko Date: Thu, 29 Feb 2024 15:07:23 -0500 Subject: [PATCH 3/8] Changed test case to use six hump camel test function and set initial guess so that the following conditions are all met: 1. Hessian is indefinite 2. Directional Derivative is negative 3. The initial guess is between a local maximum and a global minimum This causes Newton's method with and without line search to converge to the local maximum --- .../OptimizationTests/NewtonMinimizerTests.cs | 40 +++++-------------- 1 file changed, 11 insertions(+), 29 deletions(-) diff --git a/src/Numerics.Tests/OptimizationTests/NewtonMinimizerTests.cs b/src/Numerics.Tests/OptimizationTests/NewtonMinimizerTests.cs index 8480b6428..77431a6c4 100644 --- a/src/Numerics.Tests/OptimizationTests/NewtonMinimizerTests.cs +++ b/src/Numerics.Tests/OptimizationTests/NewtonMinimizerTests.cs @@ -84,48 +84,30 @@ protected override void Evaluate() } } - public class LazyBealebjectiveFunction : LazyObjectiveFunctionBase + public class LazySixHumpCamelObjectiveFunction : LazyObjectiveFunctionBase { - public LazyBealebjectiveFunction() : base(true, true) { } + public LazySixHumpCamelObjectiveFunction() : base(true, true) { } public override IObjectiveFunction CreateNew() { - return new LazyBealebjectiveFunction(); + return new LazySixHumpCamelObjectiveFunction(); } protected override void EvaluateValue() { - Value = BealeFunction2D.Value(Point); + Value = SixHumpCamelFunction.Value(Point); } protected override void EvaluateGradient() { - Gradient = BealeFunction2D.Gradient(Point); + Gradient = SixHumpCamelFunction.Gradient(Point); } protected override void EvaluateHessian() { - Hessian = BealeFunction2D.Hessian(Point); - } - } - - public class BealeObjectiveFunction : ObjectiveFunctionBase - { - public BealeObjectiveFunction() : base(true, true) { } - - public override IObjectiveFunction CreateNew() - { - return new BealeObjectiveFunction(); + Hessian = SixHumpCamelFunction.Hessian(Point); } - protected override void Evaluate() - { - // here we could directly overwrite the existing matrix cells instead. - // note: values must then be initialized manually first, if null. - Value = BealeFunction2D.Value(Point); ; - Gradient = BealeFunction2D.Gradient(Point); - Hessian = BealeFunction2D.Hessian(Point); - } } [TestFixture] @@ -198,14 +180,14 @@ public void FindMinimum_Linesearch_Rosenbrock_Overton() } [Test] - public void FindMinimum_Beale_IndefiniteHessian() + public void FindMinimum_SixHumpCamel_IndefiniteHessian() { - var obj = new LazyBealebjectiveFunction(); + var obj = new LazySixHumpCamelObjectiveFunction(); var solver = new NewtonMinimizer(1e-5, 1000, true); - var result = solver.FindMinimum(obj, new DenseVector(new double[] { 4, 0.9 })); + var result = solver.FindMinimum(obj, new DenseVector(new double[] { 1.0, -0.6 })); - Assert.That(Math.Abs(result.MinimizingPoint[0] - 3.0), Is.LessThan(1e-3)); - Assert.That(Math.Abs(result.MinimizingPoint[1] - 0.5), Is.LessThan(1e-3)); + 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 From 86920fb8ab4471f139fadf8b1f547cf76a8a8fca Mon Sep 17 00:00:00 2001 From: strMikhailPotapenko Date: Thu, 29 Feb 2024 15:31:25 -0500 Subject: [PATCH 4/8] Implemented modification of Hessian to force it to be positive definite by reversing negative eigenvalues. Added optional enum parameter to NewtonMinimizer class to allow for future implementations of other methods of Hessian modification. --- src/Numerics/Optimization/HessianModifiers.cs | 43 ++++++++++++++++++ src/Numerics/Optimization/NewtonMinimizer.cs | 45 +++++++++++++++++-- 2 files changed, 84 insertions(+), 4 deletions(-) create mode 100644 src/Numerics/Optimization/HessianModifiers.cs diff --git a/src/Numerics/Optimization/HessianModifiers.cs b/src/Numerics/Optimization/HessianModifiers.cs new file mode 100644 index 000000000..4ee950d2a --- /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..7b9336034 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,40 @@ 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); + } + + static Vector ReverseNegativeEigenValuesAndSolve(IObjectiveFunction objective) + { + Evd oEVD = objective.Hessian.Evd(Symmetricity.Symmetric); + for (int i = 0; i < oEVD.EigenValues.Count; i++) + { + if (oEVD.EigenValues[i].Real < double.Epsilon) + { + oEVD.EigenValues[i] = Math.Max(-oEVD.EigenValues[i].Real, double.Epsilon); + } + } + return oEVD.Solve(-objective.Gradient); + } + static void ValidateGradient(IObjectiveFunctionEvaluation eval) { foreach (var x in eval.Gradient) From ccd10c7bd4ade3d6400557ed0bc817f2e29e9ef1 Mon Sep 17 00:00:00 2001 From: strMikhailPotapenko Date: Thu, 29 Feb 2024 15:32:03 -0500 Subject: [PATCH 5/8] Changed test case to use eigenvalue reversal of Hessian. This test now passes. --- src/Numerics.Tests/OptimizationTests/NewtonMinimizerTests.cs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Numerics.Tests/OptimizationTests/NewtonMinimizerTests.cs b/src/Numerics.Tests/OptimizationTests/NewtonMinimizerTests.cs index 77431a6c4..813c233b0 100644 --- a/src/Numerics.Tests/OptimizationTests/NewtonMinimizerTests.cs +++ b/src/Numerics.Tests/OptimizationTests/NewtonMinimizerTests.cs @@ -183,7 +183,7 @@ public void FindMinimum_Linesearch_Rosenbrock_Overton() public void FindMinimum_SixHumpCamel_IndefiniteHessian() { var obj = new LazySixHumpCamelObjectiveFunction(); - var solver = new NewtonMinimizer(1e-5, 1000, true); + 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)); From 040a9debcd2993e674b73a4ca156c7df2ff3d2ba Mon Sep 17 00:00:00 2001 From: strMikhailPotapenko Date: Fri, 1 Mar 2024 08:46:23 -0500 Subject: [PATCH 6/8] Changed to Camel Case --- src/Numerics/Optimization/NewtonMinimizer.cs | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/Numerics/Optimization/NewtonMinimizer.cs b/src/Numerics/Optimization/NewtonMinimizer.cs index 7b9336034..7984b8edb 100644 --- a/src/Numerics/Optimization/NewtonMinimizer.cs +++ b/src/Numerics/Optimization/NewtonMinimizer.cs @@ -151,15 +151,15 @@ static Vector SolveLU(IObjectiveFunction objective) static Vector ReverseNegativeEigenValuesAndSolve(IObjectiveFunction objective) { - Evd oEVD = objective.Hessian.Evd(Symmetricity.Symmetric); - for (int i = 0; i < oEVD.EigenValues.Count; i++) + Evd evd = objective.Hessian.Evd(Symmetricity.Symmetric); + for (int i = 0; i < evd.EigenValues.Count; i++) { - if (oEVD.EigenValues[i].Real < double.Epsilon) + if (evd.EigenValues[i].Real < double.Epsilon) { - oEVD.EigenValues[i] = Math.Max(-oEVD.EigenValues[i].Real, double.Epsilon); + evd.EigenValues[i] = Math.Max(-evd.EigenValues[i].Real, double.Epsilon); } } - return oEVD.Solve(-objective.Gradient); + return evd.Solve(-objective.Gradient); } static void ValidateGradient(IObjectiveFunctionEvaluation eval) From 59c96d6527cae7526c25f65044c5057f3b7c7e6c Mon Sep 17 00:00:00 2001 From: strMikhailPotapenko Date: Fri, 1 Mar 2024 09:08:28 -0500 Subject: [PATCH 7/8] Added References and links --- .../OptimizationTests/TestFunctions/SixHumpCamelFunction.cs | 3 +++ src/Numerics/Optimization/NewtonMinimizer.cs | 6 ++++++ 2 files changed, 9 insertions(+) diff --git a/src/Numerics.Tests/OptimizationTests/TestFunctions/SixHumpCamelFunction.cs b/src/Numerics.Tests/OptimizationTests/TestFunctions/SixHumpCamelFunction.cs index 9bc223261..1268d02f7 100644 --- a/src/Numerics.Tests/OptimizationTests/TestFunctions/SixHumpCamelFunction.cs +++ b/src/Numerics.Tests/OptimizationTests/TestFunctions/SixHumpCamelFunction.cs @@ -32,6 +32,9 @@ 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) diff --git a/src/Numerics/Optimization/NewtonMinimizer.cs b/src/Numerics/Optimization/NewtonMinimizer.cs index 7984b8edb..e61ebe020 100644 --- a/src/Numerics/Optimization/NewtonMinimizer.cs +++ b/src/Numerics/Optimization/NewtonMinimizer.cs @@ -149,6 +149,12 @@ 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); From 8ff2a6a2b0e2a7af361ba2086a86431cbb58d505 Mon Sep 17 00:00:00 2001 From: strMikhailPotapenko Date: Tue, 2 Jul 2024 11:04:07 -0400 Subject: [PATCH 8/8] Fixed spelling error of Eigenvalue, simplified logic that forces postive definiteness. --- .../OptimizationTests/NewtonMinimizerTests.cs | 2 +- src/Numerics/Optimization/HessianModifiers.cs | 2 +- src/Numerics/Optimization/NewtonMinimizer.cs | 11 ++++------- 3 files changed, 6 insertions(+), 9 deletions(-) diff --git a/src/Numerics.Tests/OptimizationTests/NewtonMinimizerTests.cs b/src/Numerics.Tests/OptimizationTests/NewtonMinimizerTests.cs index 813c233b0..bb20a2927 100644 --- a/src/Numerics.Tests/OptimizationTests/NewtonMinimizerTests.cs +++ b/src/Numerics.Tests/OptimizationTests/NewtonMinimizerTests.cs @@ -183,7 +183,7 @@ public void FindMinimum_Linesearch_Rosenbrock_Overton() public void FindMinimum_SixHumpCamel_IndefiniteHessian() { var obj = new LazySixHumpCamelObjectiveFunction(); - var solver = new NewtonMinimizer(1e-5, 1000, true, HessianModifiers.ReverseNegativeEigenValues); + 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)); diff --git a/src/Numerics/Optimization/HessianModifiers.cs b/src/Numerics/Optimization/HessianModifiers.cs index 4ee950d2a..1458a8d76 100644 --- a/src/Numerics/Optimization/HessianModifiers.cs +++ b/src/Numerics/Optimization/HessianModifiers.cs @@ -38,6 +38,6 @@ namespace MathNet.Numerics.Optimization public enum HessianModifiers { None, - ReverseNegativeEigenValues + ReverseNegativeEigenvalues } } diff --git a/src/Numerics/Optimization/NewtonMinimizer.cs b/src/Numerics/Optimization/NewtonMinimizer.cs index e61ebe020..ac836ebf2 100644 --- a/src/Numerics/Optimization/NewtonMinimizer.cs +++ b/src/Numerics/Optimization/NewtonMinimizer.cs @@ -136,8 +136,8 @@ static Vector CalculateSearchDirection(IObjectiveFunction objective, Hes case HessianModifiers.None: searchDirection = SolveLU(objective); break; - case HessianModifiers.ReverseNegativeEigenValues: - searchDirection = ReverseNegativeEigenValuesAndSolve(objective); + case HessianModifiers.ReverseNegativeEigenvalues: + searchDirection = ReverseNegativeEigenvaluesAndSolve(objective); break; } @@ -155,15 +155,12 @@ static Vector SolveLU(IObjectiveFunction objective) /// /// /// - static Vector ReverseNegativeEigenValuesAndSolve(IObjectiveFunction objective) + static Vector ReverseNegativeEigenvaluesAndSolve(IObjectiveFunction objective) { Evd evd = objective.Hessian.Evd(Symmetricity.Symmetric); for (int i = 0; i < evd.EigenValues.Count; i++) { - if (evd.EigenValues[i].Real < double.Epsilon) - { - evd.EigenValues[i] = Math.Max(-evd.EigenValues[i].Real, double.Epsilon); - } + evd.EigenValues[i] = Math.Max(Math.Abs(evd.EigenValues[i].Real), double.Epsilon); } return evd.Solve(-objective.Gradient); }