From b9b0686a75e4c5ce88d72bf624ce1b16eaefb758 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 15 Jul 2017 20:29:18 +0200 Subject: [PATCH] Optimization: Simple to use FindMinimum facade; static impl rename --- src/Numerics/FindMinimum.cs | 149 ++++++++++++++++++ src/Numerics/Numerics.csproj | 1 + .../Optimization/GoldenSectionMinimizer.cs | 25 ++- .../Optimization/NelderMeadSimplex.cs | 10 +- src/Numerics/Optimization/NewtonMinimizer.cs | 4 +- .../GoldenSectionMinimizerTests.cs | 2 +- .../NelderMeadSimplexTests.cs | 2 +- .../OptimizationTests/NewtonMinimizerTests.cs | 2 +- 8 files changed, 179 insertions(+), 16 deletions(-) create mode 100644 src/Numerics/FindMinimum.cs diff --git a/src/Numerics/FindMinimum.cs b/src/Numerics/FindMinimum.cs new file mode 100644 index 00000000..b32a59d2 --- /dev/null +++ b/src/Numerics/FindMinimum.cs @@ -0,0 +1,149 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// +// Copyright (c) 2009-2017 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.Optimization; +using MathNet.Numerics.Optimization.ObjectiveFunctions; + +namespace MathNet.Numerics +{ + public static class FindMinimum + { + /// + /// Find value x that minimizes the scalar function f(x), constrained within bounds, using the Golden Section algorithm. + /// For more options and diagnostics consider to use directly. + /// + public static double OfScalarFunctionConstrained(Func function, double lowerBound, double upperBound, double tolerance=1e-5, int maxIterations=1000) + { + var objective = new SimpleObjectiveFunction1D(function); + var algorithm = new GoldenSectionMinimizer(tolerance, maxIterations); + var result = algorithm.FindMinimum(objective, lowerBound, upperBound); + return result.MinimizingPoint; + } + + /// + /// Find vector x that minimizes the function f(x) using the Nelder-Mead Simplex algorithm. + /// For more options and diagnostics consider to use directly. + /// + public static Vector OfFunction(Func, double> function, Vector initialGuess, double tolerance=1e-8, int maxIterations=1000) + { + var objective = ObjectiveFunction.Value(function); + var algorithm = new NelderMeadSimplex(tolerance, maxIterations); + var result = algorithm.FindMinimum(objective, initialGuess); + return result.MinimizingPoint; + } + + /// + /// Find vector x that minimizes the function f(x), constrained within bounds, using the Broyden–Fletcher–Goldfarb–Shanno Bounded (BFGS-B) algorithm. + /// The missing gradient is evaluated numerically (forward difference). + /// For more options and diagnostics consider to use directly. + /// + public static Vector OfFunctionConstrained(Func, double> function, Vector lowerBound, Vector upperBound, Vector initialGuess, double gradientTolerance=1e-5, double parameterTolerance=1e-5, double functionProgressTolerance=1e-5, int maxIterations=1000) + { + var objective = ObjectiveFunction.Value(function); + var objectiveWithGradient = new ForwardDifferenceGradientObjectiveFunction(objective, lowerBound, upperBound); + var algorithm = new BfgsBMinimizer(gradientTolerance, parameterTolerance, functionProgressTolerance, maxIterations); + var result = algorithm.FindMinimum(objectiveWithGradient, lowerBound, upperBound, initialGuess); + return result.MinimizingPoint; + } + + /// + /// Find vector x that minimizes the function f(x) using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm. + /// For more options and diagnostics consider to use directly. + /// + public static Vector OfFunctionGradient(Func, double> function, Func, Vector> gradient, Vector initialGuess, double gradientTolerance=1e-5, double parameterTolerance=1e-5, double functionProgressTolerance=1e-5, int maxIterations=1000) + { + var objective = ObjectiveFunction.Gradient(function, gradient); + var algorithm = new BfgsMinimizer(gradientTolerance, parameterTolerance, functionProgressTolerance, maxIterations); + var result = algorithm.FindMinimum(objective, initialGuess); + return result.MinimizingPoint; + } + + /// + /// Find vector x that minimizes the function f(x) using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm. + /// For more options and diagnostics consider to use directly. + /// + public static Vector OfFunctionGradient(Func, Tuple>> functionGradient, Vector initialGuess, double gradientTolerance=1e-5, double parameterTolerance=1e-5, double functionProgressTolerance=1e-5, int maxIterations=1000) + { + var objective = ObjectiveFunction.Gradient(functionGradient); + var algorithm = new BfgsMinimizer(gradientTolerance, parameterTolerance, functionProgressTolerance, maxIterations); + var result = algorithm.FindMinimum(objective, initialGuess); + return result.MinimizingPoint; + } + + /// + /// Find vector x that minimizes the function f(x), constrained within bounds, using the Broyden–Fletcher–Goldfarb–Shanno Bounded (BFGS-B) algorithm. + /// For more options and diagnostics consider to use directly. + /// + public static Vector OfFunctionGradientConstrained(Func, double> function, Func, Vector> gradient, Vector lowerBound, Vector upperBound, Vector initialGuess, double gradientTolerance=1e-5, double parameterTolerance=1e-5, double functionProgressTolerance=1e-5, int maxIterations=1000) + { + var objective = ObjectiveFunction.Gradient(function, gradient); + var algorithm = new BfgsBMinimizer(gradientTolerance, parameterTolerance, functionProgressTolerance, maxIterations); + var result = algorithm.FindMinimum(objective, lowerBound, upperBound, initialGuess); + return result.MinimizingPoint; + } + + /// + /// Find vector x that minimizes the function f(x), constrained within bounds, using the Broyden–Fletcher–Goldfarb–Shanno Bounded (BFGS-B) algorithm. + /// For more options and diagnostics consider to use directly. + /// + public static Vector OfFunctionGradientConstrained(Func, Tuple>> functionGradient, Vector lowerBound, Vector upperBound, Vector initialGuess, double gradientTolerance=1e-5, double parameterTolerance=1e-5, double functionProgressTolerance=1e-5, int maxIterations=1000) + { + var objective = ObjectiveFunction.Gradient(functionGradient); + var algorithm = new BfgsBMinimizer(gradientTolerance, parameterTolerance, functionProgressTolerance, maxIterations); + var result = algorithm.FindMinimum(objective, lowerBound, upperBound, initialGuess); + return result.MinimizingPoint; + } + + /// + /// Find vector x that minimizes the function f(x) using the Newton algorithm. + /// For more options and diagnostics consider to use directly. + /// + public static Vector OfFunctionGradientHessian(Func, double> function, Func, Vector> gradient, Func, Matrix> hessian, Vector initialGuess, double tolerance=1e-8, int maxIterations=1000) + { + var objective = ObjectiveFunction.GradientHessian(function, gradient, hessian); + var algorithm = new NewtonMinimizer(tolerance, maxIterations); + var result = algorithm.FindMinimum(objective, initialGuess); + return result.MinimizingPoint; + } + + /// + /// Find vector x that minimizes the function f(x) using the Newton algorithm. + /// For more options and diagnostics consider to use directly. + /// + public static Vector OfFunctionGradientHessian(Func, Tuple, Matrix>> functionGradientHessian, Vector initialGuess, double tolerance=1e-8, int maxIterations=1000) + { + var objective = ObjectiveFunction.GradientHessian(functionGradientHessian); + var algorithm = new NewtonMinimizer(tolerance, maxIterations); + var result = algorithm.FindMinimum(objective, initialGuess); + return result.MinimizingPoint; + } + } +} diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 588fd070..8d8a30ad 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -92,6 +92,7 @@ + diff --git a/src/Numerics/Optimization/GoldenSectionMinimizer.cs b/src/Numerics/Optimization/GoldenSectionMinimizer.cs index f156d7a7..43d43016 100644 --- a/src/Numerics/Optimization/GoldenSectionMinimizer.cs +++ b/src/Numerics/Optimization/GoldenSectionMinimizer.cs @@ -49,9 +49,16 @@ namespace MathNet.Numerics.Optimization } public MinimizationResult1D FindMinimum(IObjectiveFunction1D objective, double lowerBound, double upperBound) + { + return Minimum(objective, lowerBound, upperBound, XTolerance, MaximumIterations, MaximumExpansionSteps, LowerExpansionFactor, UpperExpansionFactor); + } + + public static MinimizationResult1D Minimum(IObjectiveFunction1D objective, double lowerBound, double upperBound, double xTolerance = 1e-5, int maxIterations = 1000, int maxExpansionSteps = 10, double lowerExpansionFactor = 2.0, double upperExpansionFactor = 2.0) { if (upperBound <= lowerBound) + { throw new OptimizationException("Lower bound must be lower than upper bound."); + } double middlePointX = lowerBound + (upperBound - lowerBound)/(1 + Constants.GoldenRatio); IEvaluation1D lower = objective.Evaluate(lowerBound); @@ -63,17 +70,17 @@ namespace MathNet.Numerics.Optimization ValueChecker(upper.Value, upperBound); int expansion_steps = 0; - while ((expansion_steps < this.MaximumExpansionSteps) && (upper.Value < middle.Value || lower.Value < middle.Value)) + while ((expansion_steps < maxExpansionSteps) && (upper.Value < middle.Value || lower.Value < middle.Value)) { if (lower.Value < middle.Value) { - lowerBound = 0.5*(upperBound + lowerBound) - this.LowerExpansionFactor*0.5*(upperBound - lowerBound); + lowerBound = 0.5*(upperBound + lowerBound) - lowerExpansionFactor*0.5*(upperBound - lowerBound); lower = objective.Evaluate(lowerBound); } if (upper.Value < middle.Value) { - upperBound = 0.5*(upperBound + lowerBound) + this.UpperExpansionFactor*0.5*(upperBound - lowerBound); + upperBound = 0.5*(upperBound + lowerBound) + upperExpansionFactor*0.5*(upperBound - lowerBound); upper = objective.Evaluate(upperBound); } @@ -84,10 +91,12 @@ namespace MathNet.Numerics.Optimization } if (upper.Value < middle.Value || lower.Value < middle.Value) + { throw new OptimizationException("Lower and upper bounds do not necessarily bound a minimum."); + } int iterations = 0; - while (Math.Abs(upper.Point - lower.Point) > XTolerance && iterations < MaximumIterations) + while (Math.Abs(upper.Point - lower.Point) > xTolerance && iterations < maxIterations) { double testX = lower.Point + (upper.Point - middle.Point); var test = objective.Evaluate(testX); @@ -121,16 +130,20 @@ namespace MathNet.Numerics.Optimization iterations += 1; } - if (iterations == MaximumIterations) + if (iterations == maxIterations) + { throw new MaximumIterationsException("Max iterations reached."); + } return new MinimizationResult1D(middle, iterations, ExitCondition.BoundTolerance); } - void ValueChecker(double value, double point) + static void ValueChecker(double value, double point) { if (Double.IsNaN(value) || Double.IsInfinity(value)) + { throw new Exception("Objective function returned non-finite value."); + } } } } diff --git a/src/Numerics/Optimization/NelderMeadSimplex.cs b/src/Numerics/Optimization/NelderMeadSimplex.cs index 50fed8e2..8f1c2d58 100644 --- a/src/Numerics/Optimization/NelderMeadSimplex.cs +++ b/src/Numerics/Optimization/NelderMeadSimplex.cs @@ -64,7 +64,7 @@ namespace MathNet.Numerics.Optimization /// The minimum point public MinimizationResult FindMinimum(IObjectiveFunction objectiveFunction, Vector initialGuess) { - return FindMinimum(objectiveFunction, initialGuess, ConvergenceTolerance, MaximumIterations); + return Minimum(objectiveFunction, initialGuess, ConvergenceTolerance, MaximumIterations); } /// @@ -76,7 +76,7 @@ namespace MathNet.Numerics.Optimization /// The minimum point public MinimizationResult FindMinimum(IObjectiveFunction objectiveFunction, Vector initialGuess, Vector initalPertubation) { - return FindMinimum(objectiveFunction, initialGuess, initalPertubation, ConvergenceTolerance, MaximumIterations); + return Minimum(objectiveFunction, initialGuess, initalPertubation, ConvergenceTolerance, MaximumIterations); } /// @@ -87,14 +87,14 @@ namespace MathNet.Numerics.Optimization /// The objective function, no gradient or hessian needed /// The intial guess /// The minimum point - public static MinimizationResult FindMinimum(IObjectiveFunction objectiveFunction, Vector initialGuess, double convergenceTolerance, int maximumIterations=1000) + public static MinimizationResult Minimum(IObjectiveFunction objectiveFunction, Vector initialGuess, double convergenceTolerance, int maximumIterations=1000) { var initalPertubation = new LinearAlgebra.Double.DenseVector(initialGuess.Count); for (int i = 0; i < initialGuess.Count; i++) { initalPertubation[i] = initialGuess[i] == 0.0 ? 0.00025 : initialGuess[i] * 0.05; } - return FindMinimum(objectiveFunction, initialGuess, initalPertubation, convergenceTolerance, maximumIterations); + return Minimum(objectiveFunction, initialGuess, initalPertubation, convergenceTolerance, maximumIterations); } /// @@ -104,7 +104,7 @@ namespace MathNet.Numerics.Optimization /// The intial guess /// The inital pertubation /// The minimum point - public static MinimizationResult FindMinimum(IObjectiveFunction objectiveFunction, Vector initialGuess, Vector initalPertubation, double convergenceTolerance, int maximumIterations=1000) + public static MinimizationResult Minimum(IObjectiveFunction objectiveFunction, Vector initialGuess, Vector initalPertubation, double convergenceTolerance, int maximumIterations=1000) { // confirm that we are in a position to commence if (objectiveFunction == null) diff --git a/src/Numerics/Optimization/NewtonMinimizer.cs b/src/Numerics/Optimization/NewtonMinimizer.cs index 6d57086d..eac9a003 100644 --- a/src/Numerics/Optimization/NewtonMinimizer.cs +++ b/src/Numerics/Optimization/NewtonMinimizer.cs @@ -48,10 +48,10 @@ namespace MathNet.Numerics.Optimization public MinimizationResult FindMinimum(IObjectiveFunction objective, Vector initialGuess) { - return FindMinimum(objective, initialGuess, GradientTolerance, MaximumIterations, UseLineSearch); + return Minimum(objective, initialGuess, GradientTolerance, MaximumIterations, UseLineSearch); } - public static MinimizationResult FindMinimum(IObjectiveFunction objective, Vector initialGuess, double gradientTolerance, int maxIterations=1000, bool useLineSearch = false) + public static MinimizationResult Minimum(IObjectiveFunction objective, Vector initialGuess, double gradientTolerance, int maxIterations=1000, bool useLineSearch = false) { if (!objective.IsGradientSupported) { diff --git a/src/UnitTests/OptimizationTests/GoldenSectionMinimizerTests.cs b/src/UnitTests/OptimizationTests/GoldenSectionMinimizerTests.cs index 293940a4..6ff308ae 100644 --- a/src/UnitTests/OptimizationTests/GoldenSectionMinimizerTests.cs +++ b/src/UnitTests/OptimizationTests/GoldenSectionMinimizerTests.cs @@ -42,7 +42,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests var algorithm = new GoldenSectionMinimizer(1e-5, 1000); var f1 = new Func(x => (x - 3)*(x - 3)); var obj = new SimpleObjectiveFunction1D(f1); - var r1 = algorithm.FindMinimum(obj, -100, 100); + var r1 = GoldenSectionMinimizer.Minimum(obj, -100, 100); Assert.That(Math.Abs(r1.MinimizingPoint - 3.0), Is.LessThan(1e-4)); } diff --git a/src/UnitTests/OptimizationTests/NelderMeadSimplexTests.cs b/src/UnitTests/OptimizationTests/NelderMeadSimplexTests.cs index 54fc124a..30e9eac2 100644 --- a/src/UnitTests/OptimizationTests/NelderMeadSimplexTests.cs +++ b/src/UnitTests/OptimizationTests/NelderMeadSimplexTests.cs @@ -112,7 +112,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests { var obj = new MghObjectiveFunction(test_case.Function, true, true); - var result = NelderMeadSimplex.FindMinimum(obj, test_case.InitialGuess, 1e-8, 1000); + var result = NelderMeadSimplex.Minimum(obj, test_case.InitialGuess, 1e-8, 1000); if (test_case.MinimizingPoint != null) { diff --git a/src/UnitTests/OptimizationTests/NewtonMinimizerTests.cs b/src/UnitTests/OptimizationTests/NewtonMinimizerTests.cs index 7cdf65f5..3366ccd1 100644 --- a/src/UnitTests/OptimizationTests/NewtonMinimizerTests.cs +++ b/src/UnitTests/OptimizationTests/NewtonMinimizerTests.cs @@ -199,7 +199,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests { var obj = new MghObjectiveFunction(test_case.Function, true, true); - var result = NewtonMinimizer.FindMinimum(obj, test_case.InitialGuess, 1e-8, 1000, useLineSearch: false); + var result = NewtonMinimizer.Minimum(obj, test_case.InitialGuess, 1e-8, 1000, useLineSearch: false); if (test_case.MinimizingPoint != null) {