diff --git a/src/Numerics/Optimization/BfgsMinimizerBase.cs b/src/Numerics/Optimization/BfgsMinimizerBase.cs index 1917960e..9e9abe94 100644 --- a/src/Numerics/Optimization/BfgsMinimizerBase.cs +++ b/src/Numerics/Optimization/BfgsMinimizerBase.cs @@ -34,85 +34,16 @@ using System; namespace MathNet.Numerics.Optimization { - public abstract class BfgsMinimizerBase + public abstract class BfgsMinimizerBase : MinimizerBase { - public double GradientTolerance { get; set; } - public double ParameterTolerance { get; set; } - public double FunctionProgressTolerance { get; set; } - public int MaximumIterations { get; set; } - - protected const double VerySmall = 1e-15; - + /// /// /// Creates a base class for BFGS minimization /// - /// The gradient tolerance - /// The parameter tolerance - /// The funciton progress tolerance - /// The maximum number of iterations - protected BfgsMinimizerBase(double gradientTolerance, double parameterTolerance, double functionProgressTolerance, int maximumIterations) - { - GradientTolerance = gradientTolerance; - ParameterTolerance = parameterTolerance; - FunctionProgressTolerance = functionProgressTolerance; - MaximumIterations = maximumIterations; - } - - protected ExitCondition ExitCriteriaSatisfied(IObjectiveFunctionEvaluation candidatePoint, IObjectiveFunctionEvaluation lastPoint, int iterations) - { - Vector relGrad = new DenseVector(candidatePoint.Point.Count); - double relativeGradient = 0.0; - double normalizer = Math.Max(Math.Abs(candidatePoint.Value), 1.0); - for (int ii = 0; ii < relGrad.Count; ++ii) - { - double projectedGradient = GetProjectedGradient(candidatePoint, ii); - - double tmp = projectedGradient * - Math.Max(Math.Abs(candidatePoint.Point[ii]), 1.0) / normalizer; - relativeGradient = Math.Max(relativeGradient, Math.Abs(tmp)); - } - if (relativeGradient < GradientTolerance) - { - return ExitCondition.RelativeGradient; - } - - if (lastPoint != null) - { - double mostProgress = 0.0; - for (int ii = 0; ii < candidatePoint.Point.Count; ++ii) - { - var tmp = Math.Abs(candidatePoint.Point[ii] - lastPoint.Point[ii]) / - Math.Max(Math.Abs(lastPoint.Point[ii]), 1.0); - mostProgress = Math.Max(mostProgress, tmp); - } - if (mostProgress < ParameterTolerance) - { - return ExitCondition.LackOfProgress; - } - - double functionChange = candidatePoint.Value - lastPoint.Value; - if (iterations > 500 && functionChange < 0 && Math.Abs(functionChange) < FunctionProgressTolerance) - return ExitCondition.LackOfProgress; - } - - return ExitCondition.None; - } - - protected virtual double GetProjectedGradient(IObjectiveFunctionEvaluation candidatePoint, int ii) + protected BfgsMinimizerBase(double gradientTolerance, double parameterTolerance, double functionProgressTolerance, int maximumIterations) : base(gradientTolerance, parameterTolerance, functionProgressTolerance, maximumIterations) { - return candidatePoint.Gradient[ii]; } - protected void ValidateGradientAndObjective(IObjectiveFunctionEvaluation eval) - { - foreach (var x in eval.Gradient) - { - if (Double.IsNaN(x) || Double.IsInfinity(x)) - throw new EvaluationException("Non-finite gradient returned.", eval); - } - if (Double.IsNaN(eval.Value) || Double.IsInfinity(eval.Value)) - throw new EvaluationException("Non-finite objective function returned.", eval); - } protected int DoBfgsUpdate(ref ExitCondition currentExitCondition, WolfeLineSearch lineSearcher, ref Matrix inversePseudoHessian, ref Vector lineSearchDirection, ref IObjectiveFunction previousPoint, ref LineSearchResult lineSearchResult, ref IObjectiveFunction candidate, ref Vector step, ref int totalLineSearchSteps, ref int iterationsWithNontrivialLineSearch) { diff --git a/src/Numerics/Optimization/LimitedMemoryBfgsMinimizer.cs b/src/Numerics/Optimization/LimitedMemoryBfgsMinimizer.cs new file mode 100644 index 00000000..0c1ac53e --- /dev/null +++ b/src/Numerics/Optimization/LimitedMemoryBfgsMinimizer.cs @@ -0,0 +1,179 @@ +// +// 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 System.Collections.Generic; +using System.Linq; +using MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.Optimization.LineSearch; + +namespace MathNet.Numerics.Optimization +{ + /// + /// Limited Memory version of Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm + /// + public class LimitedMemoryBfgsMinimizer : MinimizerBase, IUnconstrainedMinimizer + { + public int Memory { get; set; } + + /// + /// + /// Creates L-BFGS minimizer + /// + /// Numbers of gradients and steps to store. + public LimitedMemoryBfgsMinimizer(double gradientTolerance, double parameterTolerance, double functionProgressTolerance, int memory, int maximumIterations=1000) : base(gradientTolerance, parameterTolerance, functionProgressTolerance, maximumIterations) + { + Memory = memory; + } + + /// + /// Find the minimum of the objective function given lower and upper bounds + /// + /// The objective function, must support a gradient + /// The initial guess + /// The MinimizationResult which contains the minimum and the ExitCondition + public MinimizationResult FindMinimum(IObjectiveFunction objective, Vector initialGuess) + { + if (!objective.IsGradientSupported) + throw new IncompatibleObjectiveException("Gradient not supported in objective function, but required for L-BFGS minimization."); + + objective.EvaluateAt(initialGuess); + ValidateGradientAndObjective(objective); + + // Check that we're not already done + ExitCondition currentExitCondition = ExitCriteriaSatisfied(objective, null, 0); + if (currentExitCondition != ExitCondition.None) + return new MinimizationResult(objective, 0, currentExitCondition); + + // Set up line search algorithm + var lineSearcher = new WeakWolfeLineSearch(1e-4, 0.9, Math.Max(ParameterTolerance, 1e-10), 1000); + + // First step + + var lineSearchDirection = -objective.Gradient; + var stepSize = 100 * GradientTolerance / (lineSearchDirection * lineSearchDirection); + + var previousPoint = objective; + + LineSearchResult lineSearchResult; + try + { + lineSearchResult = lineSearcher.FindConformingStep(objective, lineSearchDirection, stepSize); + } + catch (OptimizationException e) + { + throw new InnerOptimizationException("Line search failed.", e); + } + catch (ArgumentException e) + { + throw new InnerOptimizationException("Line search failed.", e); + } + + var candidate = lineSearchResult.FunctionInfoAtMinimum; + ValidateGradientAndObjective(candidate); + + var gradient = candidate.Gradient; + var step = candidate.Point - initialGuess; + var yk = candidate.Gradient - previousPoint.Gradient; + var ykhistory = new List>() {yk}; + var skhistory = new List>() {step}; + var rhokhistory = new List() {1.0/yk.DotProduct(step)}; + + // Subsequent steps + int iterations = 1; + int totalLineSearchSteps = lineSearchResult.Iterations; + int iterationsWithNontrivialLineSearch = lineSearchResult.Iterations > 0 ? 0 : 1; + previousPoint = candidate; + while (iterations++ < MaximumIterations && previousPoint.Gradient.Norm(2) >= GradientTolerance) + { + lineSearchDirection = -ApplyLbfgsUpdate(previousPoint, ykhistory, skhistory, rhokhistory); + var directionalDerivative = previousPoint.Gradient.DotProduct(lineSearchDirection); + if (directionalDerivative > 0) + throw new InnerOptimizationException("Direction is not a descent direction."); + try + { + lineSearchResult = lineSearcher.FindConformingStep(previousPoint, lineSearchDirection, 1.0); + } + catch (OptimizationException e) + { + throw new InnerOptimizationException("Line search failed.", e); + } + catch (ArgumentException e) + { + throw new InnerOptimizationException("Line search failed.", e); + } + iterationsWithNontrivialLineSearch += lineSearchResult.Iterations > 0 ? 1 : 0; + totalLineSearchSteps += lineSearchResult.Iterations; + + candidate = lineSearchResult.FunctionInfoAtMinimum; + currentExitCondition = ExitCriteriaSatisfied(candidate, previousPoint, iterations); + if (currentExitCondition != ExitCondition.None) + break; + step = candidate.Point - previousPoint.Point; + yk = candidate.Gradient - previousPoint.Gradient; + ykhistory.Add(yk); + skhistory.Add(step); + rhokhistory.Add(1.0/yk.DotProduct(step)); + previousPoint = candidate; + if (ykhistory.Count > Memory) + { + ykhistory.RemoveAt(0); + skhistory.RemoveAt(0); + rhokhistory.RemoveAt(0); + } + } + + if (iterations == MaximumIterations && currentExitCondition == ExitCondition.None) + throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", MaximumIterations)); + + return new MinimizationWithLineSearchResult(candidate, iterations, ExitCondition.AbsoluteGradient, totalLineSearchSteps, iterationsWithNontrivialLineSearch); + } + + private Vector ApplyLbfgsUpdate(IObjectiveFunction previousPoint, List> ykhistory, List> skhistory, List rhokhistory) + { + var q = previousPoint.Gradient.Clone(); + var alphas = new Stack(); + for (int k = ykhistory.Count - 1; k >= 0; k--) + { + var alpha = rhokhistory[k]*q.DotProduct(skhistory[k]); + alphas.Push(alpha); + q -= alpha*ykhistory[k]; + } + var yk = ykhistory.Last(); + var sk = skhistory.Last(); + q *= yk.DotProduct(sk)/yk.DotProduct(yk); + for (int k = 0; k < ykhistory.Count; k++) + { + var beta = rhokhistory[k]*ykhistory[k].DotProduct(q); + q += skhistory[k]*(alphas.Pop() - beta); + } + return q; + } + } +} diff --git a/src/Numerics/Optimization/MinimizerBase.cs b/src/Numerics/Optimization/MinimizerBase.cs new file mode 100644 index 00000000..e7f27c7f --- /dev/null +++ b/src/Numerics/Optimization/MinimizerBase.cs @@ -0,0 +1,117 @@ +// +// 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 MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.LinearAlgebra.Double; +using MathNet.Numerics.Optimization.LineSearch; +using System; + +namespace MathNet.Numerics.Optimization +{ + public abstract class MinimizerBase + { + public double GradientTolerance { get; set; } + public double ParameterTolerance { get; set; } + public double FunctionProgressTolerance { get; set; } + public int MaximumIterations { get; set; } + + protected const double VerySmall = 1e-15; + + /// + /// Creates a base class for minimization + /// + /// The gradient tolerance + /// The parameter tolerance + /// The funciton progress tolerance + /// The maximum number of iterations + protected MinimizerBase(double gradientTolerance, double parameterTolerance, double functionProgressTolerance, int maximumIterations) + { + GradientTolerance = gradientTolerance; + ParameterTolerance = parameterTolerance; + FunctionProgressTolerance = functionProgressTolerance; + MaximumIterations = maximumIterations; + } + + protected ExitCondition ExitCriteriaSatisfied(IObjectiveFunctionEvaluation candidatePoint, IObjectiveFunctionEvaluation lastPoint, int iterations) + { + Vector relGrad = new DenseVector(candidatePoint.Point.Count); + double relativeGradient = 0.0; + double normalizer = Math.Max(Math.Abs(candidatePoint.Value), 1.0); + for (int ii = 0; ii < relGrad.Count; ++ii) + { + double projectedGradient = GetProjectedGradient(candidatePoint, ii); + + double tmp = projectedGradient * + Math.Max(Math.Abs(candidatePoint.Point[ii]), 1.0) / normalizer; + relativeGradient = Math.Max(relativeGradient, Math.Abs(tmp)); + } + if (relativeGradient < GradientTolerance) + { + return ExitCondition.RelativeGradient; + } + + if (lastPoint != null) + { + double mostProgress = 0.0; + for (int ii = 0; ii < candidatePoint.Point.Count; ++ii) + { + var tmp = Math.Abs(candidatePoint.Point[ii] - lastPoint.Point[ii]) / + Math.Max(Math.Abs(lastPoint.Point[ii]), 1.0); + mostProgress = Math.Max(mostProgress, tmp); + } + if (mostProgress < ParameterTolerance) + { + return ExitCondition.LackOfProgress; + } + + double functionChange = candidatePoint.Value - lastPoint.Value; + if (iterations > 500 && functionChange < 0 && Math.Abs(functionChange) < FunctionProgressTolerance) + return ExitCondition.LackOfProgress; + } + + return ExitCondition.None; + } + + protected virtual double GetProjectedGradient(IObjectiveFunctionEvaluation candidatePoint, int ii) + { + return candidatePoint.Gradient[ii]; + } + + protected void ValidateGradientAndObjective(IObjectiveFunctionEvaluation eval) + { + foreach (var x in eval.Gradient) + { + if (Double.IsNaN(x) || Double.IsInfinity(x)) + throw new EvaluationException("Non-finite gradient returned.", eval); + } + if (Double.IsNaN(eval.Value) || Double.IsInfinity(eval.Value)) + throw new EvaluationException("Non-finite objective function returned.", eval); + } + } +} diff --git a/src/UnitTests/OptimizationTests/LBfgsMinimizerTests.cs b/src/UnitTests/OptimizationTests/LBfgsMinimizerTests.cs new file mode 100644 index 00000000..efbe2385 --- /dev/null +++ b/src/UnitTests/OptimizationTests/LBfgsMinimizerTests.cs @@ -0,0 +1,159 @@ +// +// 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 System.Linq; +using MathNet.Numerics.LinearAlgebra.Double; +using MathNet.Numerics.Optimization; +using NUnit.Framework; +using MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions; +using System.Collections.Generic; +using System.Collections; +using NUnit.Framework.Interfaces; + +namespace MathNet.Numerics.UnitTests.OptimizationTests +{ + [TestFixture] + public class LBfgsMinimizerTests + { + [Test] + public void FindMinimum_Rosenbrock_Easy() + { + var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient); + var solver = new LimitedMemoryBfgsMinimizer(1e-5, 1e-5, 1e-5, 5, 100); + var result = solver.FindMinimum(obj, new DenseVector(new[] { 1.2, 1.2 })); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - RosenbrockFunction.Minimum[0]), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - RosenbrockFunction.Minimum[1]), Is.LessThan(1e-3)); + } + + [Test] + public void FindMinimum_Rosenbrock_Hard() + { + var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient); + var solver = new LimitedMemoryBfgsMinimizer(1e-5, 1e-5, 1e-5, 5, 100); + var result = solver.FindMinimum(obj, new DenseVector(new[] { -1.2, 1.0 })); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - RosenbrockFunction.Minimum[0]), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - RosenbrockFunction.Minimum[1]), Is.LessThan(1e-3)); + } + + [Test] + public void FindMinimum_Rosenbrock_Overton() + { + var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient); + var solver = new LimitedMemoryBfgsMinimizer(1e-5, 1e-5, 1e-5, 5, 100); + var result = solver.FindMinimum(obj, new DenseVector(new[] { -0.9, -0.5 })); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - RosenbrockFunction.Minimum[0]), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - RosenbrockFunction.Minimum[1]), Is.LessThan(1e-3)); + } + + [Test] + public void FindMinimum_BigRosenbrock_Easy() + { + var obj = ObjectiveFunction.Gradient(BigRosenbrockFunction.Value, BigRosenbrockFunction.Gradient); + var solver = new LimitedMemoryBfgsMinimizer(1e-10, 1e-5, 1e-5, 5, 1000); + var result = solver.FindMinimum(obj, new DenseVector(new[] { 1.2 * 100.0, 1.2 * 100.0 })); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - BigRosenbrockFunction.Minimum[0]), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - BigRosenbrockFunction.Minimum[1]), Is.LessThan(1e-3)); + } + + [Test] + public void FindMinimum_BigRosenbrock_Hard() + { + var obj = ObjectiveFunction.Gradient(BigRosenbrockFunction.Value, BigRosenbrockFunction.Gradient); + var solver = new LimitedMemoryBfgsMinimizer(1e-5, 1e-5, 1e-5, 5, 1000); + var result = solver.FindMinimum(obj, new DenseVector(new[] { -1.2 * 100.0, 1.0 * 100.0 })); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - BigRosenbrockFunction.Minimum[0]), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - BigRosenbrockFunction.Minimum[1]), Is.LessThan(1e-3)); + } + + [Test] + public void FindMinimum_BigRosenbrock_Overton() + { + var obj = ObjectiveFunction.Gradient(BigRosenbrockFunction.Value, BigRosenbrockFunction.Gradient); + var solver = new LimitedMemoryBfgsMinimizer(1e-5, 1e-5, 1e-5, 5, 1000); + var result = solver.FindMinimum(obj, new DenseVector(new[] { -0.9 * 100.0, -0.5 * 100.0 })); + + Assert.That(Math.Abs(result.MinimizingPoint[0] - BigRosenbrockFunction.Minimum[0]), Is.LessThan(1e-3)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - BigRosenbrockFunction.Minimum[1]), Is.LessThan(1e-3)); + } + + private class MghTestCaseEnumerator : IEnumerable + { + public IEnumerator GetEnumerator() + { + return + RosenbrockFunction2.TestCases + .Concat(BealeFunction.TestCases) + .Concat(HelicalValleyFunction.TestCases) + .Concat(MeyerFunction.TestCases) + .Concat(PowellSingularFunction.TestCases) + .Concat(WoodFunction.TestCases) + .Concat(BrownAndDennisFunction.TestCases) + .Where(x => x.IsUnbounded) + .Select(x => new TestCaseData(x) + .SetName(x.FullName) + ) + .GetEnumerator(); + } + + IEnumerator IEnumerable.GetEnumerator() + { + return this.GetEnumerator(); + } + } + + [Test] + [TestCaseSource(typeof(MghTestCaseEnumerator))] + public void Mgh_Tests(TestFunctions.TestCase test_case) + { + var obj = new MghObjectiveFunction(test_case.Function, true, true); + var solver = new LimitedMemoryBfgsMinimizer(1e-8, 1e-8, 1e-8, 5, 1000); + + var result = solver.FindMinimum(obj, 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."); + } + } +}