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();