From 415b2e9ac68424c3a43345dd228c80a28a9911db Mon Sep 17 00:00:00 2001 From: Erik Ovegard Date: Mon, 23 May 2016 15:48:49 +0200 Subject: [PATCH] Switching computer --- src/Numerics/Numerics.csproj | 1 + src/Numerics/Optimization/BfgsBMinimizer.cs | 4 + src/Numerics/Optimization/BfgsMinimizer.cs | 3 + .../LineSearch/StrongWolfeLineSearch.cs | 102 +------------- .../LineSearch/WeakWolfeLineSearch.cs | 115 +++------------ .../LineSearch/WolfeLineSearch.cs | 131 ++++++++++++++++++ .../NelderMeadSimplexTests.cs | 2 +- 7 files changed, 162 insertions(+), 196 deletions(-) create mode 100644 src/Numerics/Optimization/LineSearch/WolfeLineSearch.cs diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 763b6278..4ef1934c 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -117,6 +117,7 @@ + diff --git a/src/Numerics/Optimization/BfgsBMinimizer.cs b/src/Numerics/Optimization/BfgsBMinimizer.cs index 457b3ce9..eda05765 100644 --- a/src/Numerics/Optimization/BfgsBMinimizer.cs +++ b/src/Numerics/Optimization/BfgsBMinimizer.cs @@ -6,6 +6,10 @@ using MathNet.Numerics.Optimization.LineSearch; namespace MathNet.Numerics.Optimization { + /// + /// Broyden–Fletcher–Goldfarb–Shanno Bounded (BFGS-B) algorithm is an iterative method for solving box-constrained nonlinear optimization problems + /// http://www.ece.northwestern.edu/~nocedal/PSfiles/limited.ps.gz + /// public class BfgsBMinimizer { public double GradientTolerance { get; set; } diff --git a/src/Numerics/Optimization/BfgsMinimizer.cs b/src/Numerics/Optimization/BfgsMinimizer.cs index a0791a4a..96d76518 100644 --- a/src/Numerics/Optimization/BfgsMinimizer.cs +++ b/src/Numerics/Optimization/BfgsMinimizer.cs @@ -4,6 +4,9 @@ using MathNet.Numerics.Optimization.LineSearch; namespace MathNet.Numerics.Optimization { + /// + /// Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm is an iterative method for solving unconstrained nonlinear optimization problems + /// public class BfgsMinimizer { public double GradientTolerance { get; set; } diff --git a/src/Numerics/Optimization/LineSearch/StrongWolfeLineSearch.cs b/src/Numerics/Optimization/LineSearch/StrongWolfeLineSearch.cs index b543ba48..a2bf0da7 100644 --- a/src/Numerics/Optimization/LineSearch/StrongWolfeLineSearch.cs +++ b/src/Numerics/Optimization/LineSearch/StrongWolfeLineSearch.cs @@ -3,109 +3,19 @@ using MathNet.Numerics.LinearAlgebra; namespace MathNet.Numerics.Optimization.LineSearch { - public class StrongWolfeLineSearch + public class StrongWolfeLineSearch : WolfeLineSearch { - public double C1 { get; set; } - public double C2 { get; set; } - public double ParameterTolerance { get; set; } - public int MaximumIterations { get; set; } - public StrongWolfeLineSearch(double c1, double c2, double parameterTolerance, int maxIterations = 10) + : base(c1, c2, parameterTolerance, maxIterations) { - C1 = c1; - C2 = c2; - ParameterTolerance = parameterTolerance; - MaximumIterations = maxIterations; + // Validation in base class } - // Implemented following http://www.math.washington.edu/~burke/crs/408/lectures/L9-weak-Wolfe.pdf - public LineSearchResult FindConformingStep(IObjectiveFunctionEvaluation objective, Vector searchDirection, double initialStep, double upperBound = Double.PositiveInfinity) - { - double lowerBound = 0.0; - double step = initialStep; - - double initialValue = objective.Value; - Vector initialGradient = objective.Gradient; - - double initialDd = searchDirection*initialGradient; - - int ii; - IObjectiveFunction candidateEval = objective.CreateNew(); - MinimizationResult.ExitCondition reasonForExit = MinimizationResult.ExitCondition.None; - for (ii = 0; ii < this.MaximumIterations; ++ii) - { - candidateEval.EvaluateAt(objective.Point + searchDirection*step); - - double stepDd = searchDirection*candidateEval.Gradient; - - if (candidateEval.Value > initialValue + C1*step*initialDd) - { - upperBound = step; - step = 0.5*(lowerBound + upperBound); - } - else if (Math.Abs(stepDd) > C2*Math.Abs(initialDd)) - { - lowerBound = step; - step = Double.IsPositiveInfinity(upperBound) ? 2*lowerBound : 0.5*(lowerBound + upperBound); - } - else - { - reasonForExit = MinimizationResult.ExitCondition.StrongWolfeCriteria; - break; - } - - if (!Double.IsInfinity(upperBound)) - { - double maxRelChange = 0.0; - for (int jj = 0; jj < candidateEval.Point.Count; ++jj) - { - double tmp = Math.Abs(searchDirection[jj]*(upperBound - lowerBound))/Math.Max(Math.Abs(candidateEval.Point[jj]), 1.0); - maxRelChange = Math.Max(maxRelChange, tmp); - } - if (maxRelChange < ParameterTolerance) - { - reasonForExit = MinimizationResult.ExitCondition.LackOfProgress; - break; - } - } - } - - if (ii == MaximumIterations && Double.IsPositiveInfinity(upperBound)) - throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached. Function appears to be unbounded in search direction.", MaximumIterations)); - if (ii == MaximumIterations) - throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", MaximumIterations)); - - return new LineSearchResult(candidateEval, ii, step, reasonForExit); - } - - bool Conforms(IObjectiveFunction startingPoint, Vector searchDirection, double step, IObjectiveFunction endingPoint) - { - bool sufficientDecrease = endingPoint.Value <= startingPoint.Value + C1*step*(startingPoint.Gradient*searchDirection); - bool notTooSteep = endingPoint.Gradient*searchDirection >= C2*startingPoint.Gradient*searchDirection; - - return step > 0 && sufficientDecrease && notTooSteep; - } - - void ValidateValue(IObjectiveFunction eval) - { - if (!IsFinite(eval.Value)) - throw new EvaluationException(String.Format("Non-finite value returned by objective function: {0}", eval.Value), eval); - } - - void ValidateGradient(IObjectiveFunction eval) - { - foreach (double x in eval.Gradient) - { - if (!IsFinite(x)) - { - throw new EvaluationException(String.Format("Non-finite value returned by gradient: {0}", x), eval); - } - } - } + protected override MinimizationResult.ExitCondition WolfeExitCondition { get { return MinimizationResult.ExitCondition.StrongWolfeCriteria; } } - bool IsFinite(double x) + protected override bool WolfeCondition(double stepDd, double initialDd) { - return !(Double.IsNaN(x) || Double.IsInfinity(x)); + return Math.Abs(stepDd) > C2 * Math.Abs(initialDd); } } } diff --git a/src/Numerics/Optimization/LineSearch/WeakWolfeLineSearch.cs b/src/Numerics/Optimization/LineSearch/WeakWolfeLineSearch.cs index bf12b146..ff5ed3ab 100644 --- a/src/Numerics/Optimization/LineSearch/WeakWolfeLineSearch.cs +++ b/src/Numerics/Optimization/LineSearch/WeakWolfeLineSearch.cs @@ -16,111 +16,22 @@ namespace MathNet.Numerics.Optimization.LineSearch /// http://en.wikipedia.org/wiki/Wolfe_conditions /// http://www.math.washington.edu/~burke/crs/408/lectures/L9-weak-Wolfe.pdf /// - public class WeakWolfeLineSearch + public class WeakWolfeLineSearch : WolfeLineSearch { - readonly double _c1; - readonly double _c2; - readonly double _parameterTolerance; - readonly int _maximumIterations; - public WeakWolfeLineSearch(double c1, double c2, double parameterTolerance, int maxIterations = 10) + : base(c1,c2,parameterTolerance,maxIterations) { - if (c1 <= 0) - throw new ArgumentException(string.Format("c1 {0} should be greater than 0", c1)); - if (c2 <= c1) - throw new ArgumentException(string.Format("c1 {0} should be less than c2 {1}", c1, c2)); - if (c2 >= 1) - throw new ArgumentException(string.Format("c2 {0} should be less than 1", c2)); - - _c1 = c1; - _c2 = c2; - _parameterTolerance = parameterTolerance; - _maximumIterations = maxIterations; + // Validation in base class } - /// The objective function being optimized, evaluated at the starting point of the search - /// Search direction - /// Initial size of the step in the search direction - public LineSearchResult FindConformingStep(IObjectiveFunctionEvaluation startingPoint, Vector searchDirection, double initialStep) - { - if (!startingPoint.IsGradientSupported) - throw new ArgumentException("objective function does not support gradient"); - - double lowerBound = 0.0; - double upperBound = Double.PositiveInfinity; - double step = initialStep; - - Vector initialPoint = startingPoint.Point; - double initialValue = startingPoint.Value; - Vector initialGradient = startingPoint.Gradient; - - double initialDd = searchDirection * initialGradient; - - var objective = startingPoint.CreateNew(); - int ii; - MinimizationResult.ExitCondition reasonForExit = MinimizationResult.ExitCondition.None; - for (ii = 0; ii < _maximumIterations; ++ii) - { - objective.EvaluateAt(initialPoint + searchDirection * step); - ValidateGradient(objective); - ValidateValue(objective); - - double stepDd = searchDirection * objective.Gradient; - - if (objective.Value > initialValue + _c1 * step * initialDd) - { - upperBound = step; - step = 0.5 * (lowerBound + upperBound); - } - else if (stepDd < _c2 * initialDd) - { - lowerBound = step; - step = Double.IsPositiveInfinity(upperBound) ? 2 * lowerBound : 0.5 * (lowerBound + upperBound); - } - else - { - reasonForExit = MinimizationResult.ExitCondition.WeakWolfeCriteria; - break; - } - - if (!Double.IsInfinity(upperBound)) - { - double maxRelChange = 0.0; - for (int jj = 0; jj < objective.Point.Count; ++jj) - { - double tmp = Math.Abs(searchDirection[jj] * (upperBound - lowerBound)) / Math.Max(Math.Abs(objective.Point[jj]), 1.0); - maxRelChange = Math.Max(maxRelChange, tmp); - } - if (maxRelChange < _parameterTolerance) - { - reasonForExit = MinimizationResult.ExitCondition.LackOfProgress; - break; - } - } - } - - if (ii == _maximumIterations && Double.IsPositiveInfinity(upperBound)) - { - throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached. Function appears to be unbounded in search direction.", _maximumIterations)); - } - - if (ii == _maximumIterations) - { - throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", _maximumIterations)); - } - - return new LineSearchResult(objective, ii, step, reasonForExit); - } + protected override MinimizationResult.ExitCondition WolfeExitCondition { get { return MinimizationResult.ExitCondition.WeakWolfeCriteria; } } - bool Conforms(IObjectiveFunction startingPoint, Vector searchDirection, double step, IObjectiveFunction endingPoint) + protected override bool WolfeCondition(double stepDd, double initialDd) { - bool sufficientDecrease = endingPoint.Value <= startingPoint.Value + _c1 * step * (startingPoint.Gradient * searchDirection); - bool notTooSteep = endingPoint.Gradient * searchDirection >= _c2 * startingPoint.Gradient * searchDirection; - - return step > 0 && sufficientDecrease && notTooSteep; + return stepDd < C2 * initialDd; } - static void ValidateValue(IObjectiveFunction eval) + protected override void ValidateValue(IObjectiveFunction eval) { if (!IsFinite(eval.Value)) { @@ -128,20 +39,26 @@ namespace MathNet.Numerics.Optimization.LineSearch } } - static void ValidateGradient(IObjectiveFunction eval) + protected override void ValidateInputArguments(IObjectiveFunctionEvaluation startingPoint, Vector searchDirection, double initialStep, double upperBound) + { + if (!startingPoint.IsGradientSupported) + throw new ArgumentException("objective function does not support gradient"); + } + + protected override void ValidateGradient(IObjectiveFunction eval) { foreach (double x in eval.Gradient) { if (!IsFinite(x)) { - throw new EvaluationException(String.Format("Non-finite value returned by gradient: {0}", x), eval); + throw new EvaluationException(string.Format("Non-finite value returned by gradient: {0}", x), eval); } } } static bool IsFinite(double x) { - return !(Double.IsNaN(x) || Double.IsInfinity(x)); + return !(double.IsNaN(x) || double.IsInfinity(x)); } } } diff --git a/src/Numerics/Optimization/LineSearch/WolfeLineSearch.cs b/src/Numerics/Optimization/LineSearch/WolfeLineSearch.cs new file mode 100644 index 00000000..e8ff4313 --- /dev/null +++ b/src/Numerics/Optimization/LineSearch/WolfeLineSearch.cs @@ -0,0 +1,131 @@ +using MathNet.Numerics.LinearAlgebra; +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +namespace MathNet.Numerics.Optimization.LineSearch +{ + public abstract class WolfeLineSearch + { + protected double C1 { get; } + protected double C2 { get; } + protected double ParameterTolerance { get; } + protected int MaximumIterations { get; } + + public WolfeLineSearch(double c1, double c2, double parameterTolerance, int maxIterations = 10) + { + if (c1 <= 0) + throw new ArgumentException(string.Format("c1 {0} should be greater than 0", c1)); + if (c2 <= c1) + throw new ArgumentException(string.Format("c1 {0} should be less than c2 {1}", c1, c2)); + if (c2 >= 1) + throw new ArgumentException(string.Format("c2 {0} should be less than 1", c2)); + + C1 = c1; + C2 = c2; + ParameterTolerance = parameterTolerance; + MaximumIterations = maxIterations; + } + + /// Implemented following http://www.math.washington.edu/~burke/crs/408/lectures/L9-weak-Wolfe.pdf + /// The objective function being optimized, evaluated at the starting point of the search + /// Search direction + /// Initial size of the step in the search direction + public LineSearchResult FindConformingStep(IObjectiveFunctionEvaluation startingPoint, Vector searchDirection, double initialStep) + { + return FindConformingStep(startingPoint, searchDirection, initialStep, double.PositiveInfinity); + } + + /// + /// The objective function being optimized, evaluated at the starting point of the search + /// Search direction + /// Initial size of the step in the search direction + /// The upper bound + public LineSearchResult FindConformingStep(IObjectiveFunctionEvaluation startingPoint, Vector searchDirection, double initialStep, double upperBound) + { + ValidateInputArguments(startingPoint, searchDirection, initialStep, upperBound); + + double lowerBound = 0.0; + double step = initialStep; + + double initialValue = startingPoint.Value; + Vector initialGradient = startingPoint.Gradient; + + double initialDd = searchDirection * initialGradient; + + IObjectiveFunction objective = startingPoint.CreateNew(); + int ii; + MinimizationResult.ExitCondition reasonForExit = MinimizationResult.ExitCondition.None; + for (ii = 0; ii < MaximumIterations; ++ii) + { + objective.EvaluateAt(startingPoint.Point + searchDirection * step); + ValidateGradient(objective); // Differ! (added) + ValidateValue(objective); // Differ! (added) + + double stepDd = searchDirection * objective.Gradient; + + if (objective.Value > initialValue + C1 * step * initialDd) + { + upperBound = step; + step = 0.5 * (lowerBound + upperBound); + } + else if (WolfeCondition(stepDd,initialDd)) // Differ, weak Wolfe + { + lowerBound = step; + step = double.IsPositiveInfinity(upperBound) ? 2 * lowerBound : 0.5 * (lowerBound + upperBound); + } + else + { + reasonForExit = WolfeExitCondition; + break; + } + + if (!double.IsInfinity(upperBound)) + { + double maxRelChange = 0.0; + for (int jj = 0; jj < objective.Point.Count; ++jj) + { + double tmp = Math.Abs(searchDirection[jj] * (upperBound - lowerBound)) / Math.Max(Math.Abs(objective.Point[jj]), 1.0); + maxRelChange = Math.Max(maxRelChange, tmp); + } + if (maxRelChange < ParameterTolerance) + { + reasonForExit = MinimizationResult.ExitCondition.LackOfProgress; + break; + } + } + } + + if (ii == MaximumIterations && Double.IsPositiveInfinity(upperBound)) + { + throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached. Function appears to be unbounded in search direction.", MaximumIterations)); + } + + if (ii == MaximumIterations) + { + throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", MaximumIterations)); + } + + return new LineSearchResult(objective, ii, step, reasonForExit); + } + protected abstract MinimizationResult.ExitCondition WolfeExitCondition { get; } + + protected abstract bool WolfeCondition(double stepDd, double initialDd); + + protected virtual void ValidateGradient(IObjectiveFunction objective) + { + } + protected virtual void ValidateValue(IObjectiveFunction objective) + { + } + + protected virtual void ValidateInputArguments(IObjectiveFunctionEvaluation startingPoint, Vector searchDirection, double initialStep, double upperBound) + { + + } + } + + + +} diff --git a/src/UnitTests/OptimizationTests/NelderMeadSimplexTests.cs b/src/UnitTests/OptimizationTests/NelderMeadSimplexTests.cs index f1a2e493..aed9409c 100644 --- a/src/UnitTests/OptimizationTests/NelderMeadSimplexTests.cs +++ b/src/UnitTests/OptimizationTests/NelderMeadSimplexTests.cs @@ -52,10 +52,10 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests var nms = new NelderMeadSimplex(1e-6, 1000); double a = 5; double b = 10; + System.Random r = new System.Random(2); IObjectiveFunction objFun = ObjectiveFunction.Value((constants)=> { double ssq = 0; - System.Random r = new System.Random(); for (double x = -10; x < 10; x += .1) { double yTrue = a * x * x + b * x + r.NextDouble();