diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj
index 4ef1934c..9c89a636 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 eda05765..9281a589 100644
--- a/src/Numerics/Optimization/BfgsBMinimizer.cs
+++ b/src/Numerics/Optimization/BfgsBMinimizer.cs
@@ -1,4 +1,34 @@
-using System;
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2016 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 MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearAlgebra.Double;
@@ -10,23 +40,25 @@ 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 class BfgsBMinimizer : BfgsMinimizerBase
{
- public double GradientTolerance { get; set; }
- public double ParameterTolerance { get; set; }
- public int MaximumIterations { get; set; }
- public double FunctionProgressTolerance { get; set; }
-
public BfgsBMinimizer(double gradientTolerance, double parameterTolerance, double functionProgressTolerance, int maximumIterations = 1000)
+ : base(gradientTolerance,parameterTolerance,functionProgressTolerance,maximumIterations)
{
- GradientTolerance = gradientTolerance;
- ParameterTolerance = parameterTolerance;
- MaximumIterations = maximumIterations;
- FunctionProgressTolerance = functionProgressTolerance;
}
+ ///
+ /// Find the minimum of the objective function given lower and upper bounds
+ ///
+ /// The objective function, must support a gradient
+ /// The lower bound
+ /// The upper bound
+ /// The initial guess
+ /// The MinimizationResult which contains the minimum and the ExitCondition
public MinimizationResult FindMinimum(IObjectiveFunction objective, Vector lowerBound, Vector upperBound, Vector initialGuess)
{
+ _lowerBound = lowerBound;
+ _upperBound = upperBound;
if (!objective.IsGradientSupported)
throw new IncompatibleObjectiveException("Gradient not supported in objective function, but required for BFGS minimization.");
@@ -40,9 +72,10 @@ namespace MathNet.Numerics.Optimization
throw new ArgumentException("Initial guess is not in the feasible region");
objective.EvaluateAt(initialGuess);
+ ValidateGradientAndObjective(objective);
// Check that we're not already done
- MinimizationResult.ExitCondition currentExitCondition = ExitCriteriaSatisfied(objective, null, lowerBound, upperBound, 0);
+ MinimizationResult.ExitCondition currentExitCondition = ExitCriteriaSatisfied(objective, null, 0);
if (currentExitCondition != MinimizationResult.ExitCondition.None)
return new MinimizationResult(objective, 0, currentExitCondition);
@@ -59,9 +92,9 @@ namespace MathNet.Numerics.Optimization
// Determine active set
var gradientProjectionResult = QuadraticGradientProjectionSearch.Search(objective.Point, objective.Gradient, pseudoHessian, lowerBound, upperBound);
- var cauchyPoint = gradientProjectionResult.Item1;
- var fixedCount = gradientProjectionResult.Item2;
- var isFixed = gradientProjectionResult.Item3;
+ var cauchyPoint = gradientProjectionResult.CauchyPoint;
+ var fixedCount = gradientProjectionResult.FixedCount;
+ var isFixed = gradientProjectionResult.IsFixed;
var freeCount = lowerBound.Count - fixedCount;
if (freeCount > 0)
@@ -96,10 +129,10 @@ namespace MathNet.Numerics.Optimization
var startingStepSize = Math.Min(Math.Max(estStepSize, 1.0), maxLineSearchStep);
// Line search
- LineSearchResult result;
+ LineSearchResult lineSearchResult;
try
{
- result = lineSearcher.FindConformingStep(objective, lineSearchDirection, startingStepSize, upperBound: maxLineSearchStep);
+ lineSearchResult = lineSearcher.FindConformingStep(objective, lineSearchDirection, startingStepSize, upperBound: maxLineSearchStep);
}
catch (Exception e)
{
@@ -107,112 +140,94 @@ namespace MathNet.Numerics.Optimization
}
var previousPoint = objective.Fork();
- var candidatePoint = result.FunctionInfoAtMinimum;
+ var candidatePoint = lineSearchResult.FunctionInfoAtMinimum;
var gradient = candidatePoint.Gradient;
var step = candidatePoint.Point - initialGuess;
// Subsequent steps
- int iterations;
- int totalLineSearchSteps = result.Iterations;
- int iterationsWithNontrivialLineSearch = result.Iterations > 0 ? 0 : 1;
- for (iterations = 1; iterations < MaximumIterations; ++iterations)
- {
- // Do BFGS update
- var y = candidatePoint.Gradient - previousPoint.Gradient;
-
- double sy = step*y;
- if (sy > 0.0) // only do update if it will create a positive definite matrix
- {
- double sts = step*step;
- //inverse_pseudo_hessian = inverse_pseudo_hessian + ((sy + y * inverse_pseudo_hessian * y) / Math.Pow(sy, 2.0)) * step.OuterProduct(step) - ((inverse_pseudo_hessian * y.ToColumnMatrix()) * step.ToRowMatrix() + step.ToColumnMatrix() * (y.ToRowMatrix() * inverse_pseudo_hessian)) * (1.0 / sy);
- var Hs = pseudoHessian*step;
- var sHs = step*pseudoHessian*step;
- pseudoHessian = pseudoHessian + y.OuterProduct(y)*(1.0/sy) - Hs.OuterProduct(Hs)*(1.0/sHs);
- }
- else
- {
- //pseudo_hessian = LinearAlgebra.Double.DiagonalMatrix.Identity(initial_guess.Count);
- }
+ int totalLineSearchSteps = lineSearchResult.Iterations;
+ int iterationsWithNontrivialLineSearch = lineSearchResult.Iterations > 0 ? 0 : 1;
- // Determine active set
- gradientProjectionResult = QuadraticGradientProjectionSearch.Search(candidatePoint.Point, candidatePoint.Gradient, pseudoHessian, lowerBound, upperBound);
- cauchyPoint = gradientProjectionResult.Item1;
- fixedCount = gradientProjectionResult.Item2;
- isFixed = gradientProjectionResult.Item3;
- freeCount = lowerBound.Count - fixedCount;
-
- if (freeCount > 0)
- {
- reducedGradient = new DenseVector(freeCount);
- reducedHessian = new DenseMatrix(freeCount, freeCount);
- reducedMap = new List(freeCount);
- reducedInitialPoint = new DenseVector(freeCount);
- reducedCauchyPoint = new DenseVector(freeCount);
+ int iterations = DoBfgsUpdate(ref currentExitCondition, lineSearcher, ref pseudoHessian, ref lineSearchDirection, ref previousPoint, ref lineSearchResult, ref candidatePoint, ref step, ref totalLineSearchSteps, ref iterationsWithNontrivialLineSearch);
- CreateReducedData(candidatePoint.Point, cauchyPoint, isFixed, lowerBound, upperBound, candidatePoint.Gradient, pseudoHessian, reducedInitialPoint, reducedCauchyPoint, reducedGradient, reducedHessian, reducedMap);
-
- // Determine search direction and maximum step size
- reducedSolution1 = reducedInitialPoint + reducedHessian.Cholesky().Solve(-reducedGradient);
-
- solution1 = ReducedToFull(reducedMap, reducedSolution1, cauchyPoint);
- }
- else
- {
- solution1 = cauchyPoint;
- }
+ if (iterations == MaximumIterations && currentExitCondition == MinimizationResult.ExitCondition.None)
+ throw new MaximumIterationsException(string.Format("Maximum iterations ({0}) reached.", MaximumIterations));
- directionFromCauchy = solution1 - cauchyPoint;
- maxStepFromCauchyPoint = FindMaxStep(cauchyPoint, directionFromCauchy, lowerBound, upperBound);
- //var cauchy_eval = objective.Evaluate(cauchy_point);
+ return new MinimizationWithLineSearchResult(candidatePoint, iterations, currentExitCondition, totalLineSearchSteps, iterationsWithNontrivialLineSearch);
+ }
- solution2 = cauchyPoint + Math.Min(maxStepFromCauchyPoint, 1.0)*directionFromCauchy;
+ protected override Vector CalculateSearchDirection(ref Matrix pseudoHessian,
+ out double maxLineSearchStep,
+ out double startingStepSize,
+ IObjectiveFunction previousPoint,
+ IObjectiveFunction candidatePoint,
+ Vector step)
+ {
+ Vector lineSearchDirection;
+ var y = candidatePoint.Gradient - previousPoint.Gradient;
- lineSearchDirection = solution2 - candidatePoint.Point;
- maxLineSearchStep = FindMaxStep(candidatePoint.Point, lineSearchDirection, lowerBound, upperBound);
+ double sy = step * y;
+ if (sy > 0.0) // only do update if it will create a positive definite matrix
+ {
+ double sts = step * step;
+
+ var Hs = pseudoHessian * step;
+ var sHs = step * pseudoHessian * step;
+ pseudoHessian = pseudoHessian + y.OuterProduct(y) * (1.0 / sy) - Hs.OuterProduct(Hs) * (1.0 / sHs);
+ }
+ else
+ {
+ //pseudo_hessian = LinearAlgebra.Double.DiagonalMatrix.Identity(initial_guess.Count);
+ }
- //line_search_direction = solution1 - candidate_point.Point;
- //max_line_search_step = FindMaxStep(candidate_point.Point, line_search_direction, lower_bound, upper_bound);
+ // Determine active set
+ var gradientProjectionResult = QuadraticGradientProjectionSearch.Search(candidatePoint.Point, candidatePoint.Gradient, pseudoHessian, _lowerBound, _upperBound);
+ var cauchyPoint = gradientProjectionResult.CauchyPoint;
+ var fixedCount = gradientProjectionResult.FixedCount;
+ var isFixed = gradientProjectionResult.IsFixed;
+ var freeCount = _lowerBound.Count - fixedCount;
+ Vector solution1;
+ if (freeCount > 0)
+ {
+ var reducedGradient = new DenseVector(freeCount);
+ var reducedHessian = new DenseMatrix(freeCount, freeCount);
+ var reducedMap = new List(freeCount);
+ var reducedInitialPoint = new DenseVector(freeCount);
+ var reducedCauchyPoint = new DenseVector(freeCount);
- if (maxLineSearchStep == 0.0)
- {
- lineSearchDirection = cauchyPoint - candidatePoint.Point;
- maxLineSearchStep = FindMaxStep(candidatePoint.Point, lineSearchDirection, lowerBound, upperBound);
- }
+ CreateReducedData(candidatePoint.Point, cauchyPoint, isFixed, _lowerBound, _upperBound, candidatePoint.Gradient, pseudoHessian, reducedInitialPoint, reducedCauchyPoint, reducedGradient, reducedHessian, reducedMap);
- estStepSize = -candidatePoint.Gradient*lineSearchDirection/(lineSearchDirection*pseudoHessian*lineSearchDirection);
+ // Determine search direction and maximum step size
+ Vector reducedSolution1 = reducedInitialPoint + reducedHessian.Cholesky().Solve(-reducedGradient);
- startingStepSize = Math.Min(Math.Max(estStepSize, 1.0), maxLineSearchStep);
+ solution1 = ReducedToFull(reducedMap, reducedSolution1, cauchyPoint);
+ }
+ else
+ {
+ solution1 = cauchyPoint;
+ }
- // Line search
- try
- {
- result = lineSearcher.FindConformingStep(candidatePoint, lineSearchDirection, startingStepSize, upperBound: maxLineSearchStep);
- //result = line_searcher.FindConformingStep(objective, cauchy_eval, direction_from_cauchy, Math.Min(1.0, max_step_from_cauchy_point), upper_bound: max_step_from_cauchy_point);
- }
- catch (Exception e)
- {
- throw new InnerOptimizationException("Line search failed.", e);
- }
+ var directionFromCauchy = solution1 - cauchyPoint;
+ var maxStepFromCauchyPoint = FindMaxStep(cauchyPoint, directionFromCauchy, _lowerBound, _upperBound);
- iterationsWithNontrivialLineSearch += result.Iterations > 0 ? 1 : 0;
- totalLineSearchSteps += result.Iterations;
+ var solution2 = cauchyPoint + Math.Min(maxStepFromCauchyPoint, 1.0) * directionFromCauchy;
- step = result.FunctionInfoAtMinimum.Point - candidatePoint.Point;
- previousPoint = candidatePoint;
- candidatePoint = result.FunctionInfoAtMinimum;
+ lineSearchDirection = solution2 - candidatePoint.Point;
+ maxLineSearchStep = FindMaxStep(candidatePoint.Point, lineSearchDirection, _lowerBound, _upperBound);
- currentExitCondition = ExitCriteriaSatisfied(candidatePoint, previousPoint, lowerBound, upperBound, iterations);
- if (currentExitCondition != MinimizationResult.ExitCondition.None)
- break;
+ if (maxLineSearchStep == 0.0)
+ {
+ lineSearchDirection = cauchyPoint - candidatePoint.Point;
+ maxLineSearchStep = FindMaxStep(candidatePoint.Point, lineSearchDirection, _lowerBound, _upperBound);
}
- if (iterations == MaximumIterations && currentExitCondition == MinimizationResult.ExitCondition.None)
- throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", MaximumIterations));
+ double estStepSize = -candidatePoint.Gradient * lineSearchDirection / (lineSearchDirection * pseudoHessian * lineSearchDirection);
- return new MinimizationWithLineSearchResult(candidatePoint, iterations, currentExitCondition, totalLineSearchSteps, iterationsWithNontrivialLineSearch);
+ startingStepSize = Math.Min(Math.Max(estStepSize, 1.0), maxLineSearchStep);
+ return lineSearchDirection;
}
- static Vector ReducedToFull(List reducedMap, Vector reducedVector, Vector fullVector)
+ private static Vector ReducedToFull(List reducedMap, Vector reducedVector, Vector fullVector)
{
var output = fullVector.Clone();
for (int ii = 0; ii < reducedMap.Count; ++ii)
@@ -220,7 +235,10 @@ namespace MathNet.Numerics.Optimization
return output;
}
- static double FindMaxStep(Vector startingPoint, Vector searchDirection, Vector lowerBound, Vector upperBound)
+ private Vector _lowerBound;
+ private Vector _upperBound;
+
+ private static double FindMaxStep(Vector startingPoint, Vector searchDirection, Vector lowerBound, Vector upperBound)
{
double maxStep = Double.PositiveInfinity;
for (int ii = 0; ii < startingPoint.Count; ++ii)
@@ -239,7 +257,7 @@ namespace MathNet.Numerics.Optimization
return maxStep;
}
- static void CreateReducedData(Vector initialPoint, Vector cauchyPoint, List isFixed, Vector lowerBound, Vector upperBound, Vector gradient, Matrix pseudoHessian, Vector reducedInitialPoint, Vector reducedCauchyPoint, Vector reducedGradient, Matrix reducedHessian, List reducedMap)
+ private static void CreateReducedData(Vector initialPoint, Vector cauchyPoint, List isFixed, Vector lowerBound, Vector upperBound, Vector gradient, Matrix pseudoHessian, Vector reducedInitialPoint, Vector reducedCauchyPoint, Vector reducedGradient, Matrix reducedHessian, List reducedMap)
{
int ll = 0;
for (int ii = 0; ii < lowerBound.Count; ++ii)
@@ -267,71 +285,21 @@ namespace MathNet.Numerics.Optimization
}
}
- const double VerySmall = 1e-15;
-
- MinimizationResult.ExitCondition ExitCriteriaSatisfied(IObjectiveFunction candidatePoint, IObjectiveFunction lastPoint, Vector lowerBound, Vector upperBound, 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;
-
- bool atLowerBound = candidatePoint.Point[ii] - lowerBound[ii] < VerySmall;
- bool atUpperBound = upperBound[ii] - candidatePoint.Point[ii] < VerySmall;
-
- if (atLowerBound && atUpperBound)
- projectedGradient = 0.0;
- else if (atLowerBound)
- projectedGradient = Math.Min(candidatePoint.Gradient[ii], 0.0);
- else if (atUpperBound)
- projectedGradient = Math.Max(candidatePoint.Gradient[ii], 0.0);
- else
- projectedGradient = candidatePoint.Gradient[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 MinimizationResult.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 MinimizationResult.ExitCondition.LackOfProgress;
- }
-
- double functionChange = candidatePoint.Value - lastPoint.Value;
- if (iterations > 500 && functionChange < 0 && Math.Abs(functionChange) < FunctionProgressTolerance)
- return MinimizationResult.ExitCondition.LackOfProgress;
- }
-
- return MinimizationResult.ExitCondition.None;
- }
-
- void ValidateGradient(IObjectiveFunction eval)
+ protected override double GetProjectedGradient(IObjectiveFunction candidatePoint, int ii)
{
- foreach (var x in eval.Gradient)
- {
- if (Double.IsNaN(x) || Double.IsInfinity(x))
- throw new EvaluationException("Non-finite gradient returned.", eval);
- }
- }
-
- void ValidateObjective(IObjectiveFunction eval)
- {
- if (Double.IsNaN(eval.Value) || Double.IsInfinity(eval.Value))
- throw new EvaluationException("Non-finite objective function returned.", eval);
+ double projectedGradient;
+ bool atLowerBound = candidatePoint.Point[ii] - _lowerBound[ii] < VerySmall;
+ bool atUpperBound = _upperBound[ii] - candidatePoint.Point[ii] < VerySmall;
+
+ if (atLowerBound && atUpperBound)
+ projectedGradient = 0.0;
+ else if (atLowerBound)
+ projectedGradient = Math.Min(candidatePoint.Gradient[ii], 0.0);
+ else if (atUpperBound)
+ projectedGradient = Math.Max(candidatePoint.Gradient[ii], 0.0);
+ else
+ projectedGradient = base.GetProjectedGradient(candidatePoint, ii);
+ return projectedGradient;
}
}
}
diff --git a/src/Numerics/Optimization/BfgsMinimizer.cs b/src/Numerics/Optimization/BfgsMinimizer.cs
index 96d76518..afac00f7 100644
--- a/src/Numerics/Optimization/BfgsMinimizer.cs
+++ b/src/Numerics/Optimization/BfgsMinimizer.cs
@@ -1,4 +1,34 @@
-using System;
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2016 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.LineSearch;
@@ -7,31 +37,36 @@ namespace MathNet.Numerics.Optimization
///
/// Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm is an iterative method for solving unconstrained nonlinear optimization problems
///
- public class BfgsMinimizer
+ public class BfgsMinimizer : BfgsMinimizerBase
{
- public double GradientTolerance { get; set; }
- public double ParameterTolerance { get; set; }
- public int MaximumIterations { get; set; }
-
- public BfgsMinimizer(double gradientTolerance, double parameterTolerance, int maximumIterations=1000)
+ ///
+ /// Creates BFGS minimizer
+ ///
+ /// The gradient tolerance
+ /// The parameter tolerance
+ /// The funciton progress tolerance
+ /// The maximum number of iterations
+ public BfgsMinimizer(double gradientTolerance, double parameterTolerance, double functionProgressTolerance, int maximumIterations=1000)
+ :base(gradientTolerance,parameterTolerance,functionProgressTolerance,maximumIterations)
{
- GradientTolerance = gradientTolerance;
- ParameterTolerance = parameterTolerance;
- MaximumIterations = maximumIterations;
}
+ ///
+ /// 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 BFGS minimization.");
objective.EvaluateAt(initialGuess);
- ValidateGradient(objective);
-
- var initial = objective.Fork();
+ ValidateGradientAndObjective(objective);
// Check that we're not already done
- MinimizationResult.ExitCondition currentExitCondition = ExitCriteriaSatisfied(objective, null);
+ MinimizationResult.ExitCondition currentExitCondition = ExitCriteriaSatisfied(objective, null, 0);
if (currentExitCondition != MinimizationResult.ExitCondition.None)
return new MinimizationResult(objective, 0, currentExitCondition);
@@ -40,122 +75,68 @@ namespace MathNet.Numerics.Optimization
// First step
var inversePseudoHessian = CreateMatrix.DenseIdentity(initialGuess.Count);
- var searchDirection = -objective.Gradient;
- var stepSize = 100 * GradientTolerance / (searchDirection * searchDirection);
+ var lineSearchDirection = -objective.Gradient;
+ var stepSize = 100 * GradientTolerance / (lineSearchDirection * lineSearchDirection);
- var previousPoint = objective.Point;
- var previousGradient = objective.Gradient;
+ var previousPoint = objective;
- LineSearchResult result;
+ LineSearchResult lineSearchResult;
try
{
- result = lineSearcher.FindConformingStep(objective, searchDirection, stepSize);
+ lineSearchResult = lineSearcher.FindConformingStep(objective, lineSearchDirection, stepSize);
}
- catch (Exception e)
+ catch (OptimizationException e)
+ {
+ throw new InnerOptimizationException("Line search failed.", e);
+ }
+ catch (ArgumentException e)
{
throw new InnerOptimizationException("Line search failed.", e);
}
- objective = result.FunctionInfoAtMinimum;
- ValidateGradient(objective);
+ var candidate = lineSearchResult.FunctionInfoAtMinimum;
+ ValidateGradientAndObjective(candidate);
- var gradient = objective.Gradient;
- var step = objective.Point - initialGuess;
- stepSize = result.FinalStep;
+ var gradient = candidate.Gradient;
+ var step = candidate.Point - initialGuess;
// Subsequent steps
Matrix I = CreateMatrix.DiagonalIdentity(initialGuess.Count);
int iterations;
- int totalLineSearchSteps = result.Iterations;
- int iterationsWithNontrivialLineSearch = result.Iterations > 0 ? 0 : 1;
- for (iterations = 1; iterations < MaximumIterations; ++iterations)
- {
- var y = objective.Gradient - previousGradient;
-
- double sy = step * y;
- inversePseudoHessian = inversePseudoHessian + ((sy + y * inversePseudoHessian * y) / Math.Pow(sy, 2.0)) * step.OuterProduct(step) - ( (inversePseudoHessian * y.ToColumnMatrix())*step.ToRowMatrix() + step.ToColumnMatrix()*(y.ToRowMatrix() * inversePseudoHessian)) * (1.0 / sy);
- searchDirection = -inversePseudoHessian * objective.Gradient;
-
- if (searchDirection * objective.Gradient >= 0.0)
- {
- searchDirection = -objective.Gradient;
- inversePseudoHessian = CreateMatrix.DenseIdentity(initialGuess.Count);
- }
-
- previousGradient = objective.Gradient;
- previousPoint = objective.Point;
-
- try
- {
- result = lineSearcher.FindConformingStep(objective, searchDirection, 1.0);
- }
- catch (Exception e)
- {
- throw new InnerOptimizationException("Line search failed.", e);
- }
-
- iterationsWithNontrivialLineSearch += result.Iterations > 0 ? 1 : 0;
- totalLineSearchSteps += result.Iterations;
- stepSize = result.FinalStep;
- step = result.FunctionInfoAtMinimum.Point - previousPoint;
- objective = result.FunctionInfoAtMinimum;
-
- currentExitCondition = ExitCriteriaSatisfied(objective, previousPoint);
- if (currentExitCondition != MinimizationResult.ExitCondition.None)
- break;
- }
+ int totalLineSearchSteps = lineSearchResult.Iterations;
+ int iterationsWithNontrivialLineSearch = lineSearchResult.Iterations > 0 ? 0 : 1;
+ iterations = DoBfgsUpdate(ref currentExitCondition, lineSearcher, ref inversePseudoHessian, ref lineSearchDirection, ref previousPoint, ref lineSearchResult, ref candidate, ref step, ref totalLineSearchSteps, ref iterationsWithNontrivialLineSearch);
if (iterations == MaximumIterations && currentExitCondition == MinimizationResult.ExitCondition.None)
throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", MaximumIterations));
- return new MinimizationWithLineSearchResult(objective, iterations, MinimizationResult.ExitCondition.AbsoluteGradient, totalLineSearchSteps, iterationsWithNontrivialLineSearch);
+ return new MinimizationWithLineSearchResult(candidate, iterations, MinimizationResult.ExitCondition.AbsoluteGradient, totalLineSearchSteps, iterationsWithNontrivialLineSearch);
}
- private MinimizationResult.ExitCondition ExitCriteriaSatisfied(IObjectiveFunction candidatePoint, Vector lastPoint)
+ protected override Vector CalculateSearchDirection(ref Matrix inversePseudoHessian,
+ out double maxLineSearchStep,
+ out double startingStepSize,
+ IObjectiveFunction previousPoint,
+ IObjectiveFunction candidate,
+ Vector step)
{
- Vector relGrad = new LinearAlgebra.Double.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 tmp = candidatePoint.Gradient[ii]*Math.Max(Math.Abs(candidatePoint.Point[ii]), 1.0) / normalizer;
- relativeGradient = Math.Max(relativeGradient, Math.Abs(tmp));
- }
- if (relativeGradient < GradientTolerance)
- {
- return MinimizationResult.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[ii])/Math.Max(Math.Abs(lastPoint[ii]), 1.0);
- mostProgress = Math.Max(mostProgress, tmp);
- }
- if ( mostProgress < ParameterTolerance )
- {
- return MinimizationResult.ExitCondition.LackOfProgress;
- }
- }
-
- return MinimizationResult.ExitCondition.None;
- }
+ startingStepSize = 1.0;
+ maxLineSearchStep = double.PositiveInfinity;
- private void ValidateGradient(IObjectiveFunction objective)
- {
- foreach (var x in objective.Gradient)
+ Vector lineSearchDirection;
+ var y = candidate.Gradient - previousPoint.Gradient;
+
+ double sy = step * y;
+ inversePseudoHessian = inversePseudoHessian + ((sy + y * inversePseudoHessian * y) / Math.Pow(sy, 2.0)) * step.OuterProduct(step) - ((inversePseudoHessian * y.ToColumnMatrix()) * step.ToRowMatrix() + step.ToColumnMatrix() * (y.ToRowMatrix() * inversePseudoHessian)) * (1.0 / sy);
+ lineSearchDirection = -inversePseudoHessian * candidate.Gradient;
+
+ if (lineSearchDirection * candidate.Gradient >= 0.0)
{
- if (Double.IsNaN(x) || Double.IsInfinity(x))
- throw new EvaluationException("Non-finite gradient returned.", objective);
+ lineSearchDirection = -candidate.Gradient;
+ inversePseudoHessian = CreateMatrix.DenseIdentity(candidate.Point.Count);
}
- }
- private void ValidateObjective(IObjectiveFunction objective)
- {
- if (Double.IsNaN(objective.Value) || Double.IsInfinity(objective.Value))
- throw new EvaluationException("Non-finite objective function returned.", objective);
+ return lineSearchDirection;
}
}
}
diff --git a/src/Numerics/Optimization/BfgsMinimizerBase.cs b/src/Numerics/Optimization/BfgsMinimizerBase.cs
new file mode 100644
index 00000000..16b06490
--- /dev/null
+++ b/src/Numerics/Optimization/BfgsMinimizerBase.cs
@@ -0,0 +1,158 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2016 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 BfgsMinimizerBase
+ {
+ 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
+ public BfgsMinimizerBase(double gradientTolerance, double parameterTolerance, double functionProgressTolerance, int maximumIterations)
+ {
+ GradientTolerance = gradientTolerance;
+ ParameterTolerance = parameterTolerance;
+ FunctionProgressTolerance = functionProgressTolerance;
+ MaximumIterations = maximumIterations;
+ }
+
+ protected MinimizationResult.ExitCondition ExitCriteriaSatisfied(IObjectiveFunction candidatePoint, IObjectiveFunction 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 MinimizationResult.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 MinimizationResult.ExitCondition.LackOfProgress;
+ }
+
+ double functionChange = candidatePoint.Value - lastPoint.Value;
+ if (iterations > 500 && functionChange < 0 && Math.Abs(functionChange) < FunctionProgressTolerance)
+ return MinimizationResult.ExitCondition.LackOfProgress;
+ }
+
+ return MinimizationResult.ExitCondition.None;
+ }
+
+ protected virtual double GetProjectedGradient(IObjectiveFunction candidatePoint, int ii)
+ {
+ return candidatePoint.Gradient[ii];
+ }
+
+ protected void ValidateGradientAndObjective(IObjectiveFunction 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 MinimizationResult.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)
+ {
+ int iterations;
+ for (iterations = 1; iterations < MaximumIterations; ++iterations)
+ {
+ double startingStepSize;
+ double maxLineSearchStep;
+ lineSearchDirection = CalculateSearchDirection(ref inversePseudoHessian, out maxLineSearchStep, out startingStepSize, previousPoint, candidate, step);
+
+ try
+ {
+ lineSearchResult = lineSearcher.FindConformingStep(candidate, lineSearchDirection, startingStepSize, maxLineSearchStep);
+ }
+ catch (Exception e)
+ {
+ throw new InnerOptimizationException("Line search failed.", e);
+ }
+
+ iterationsWithNontrivialLineSearch += lineSearchResult.Iterations > 0 ? 1 : 0;
+ totalLineSearchSteps += lineSearchResult.Iterations;
+
+ step = lineSearchResult.FunctionInfoAtMinimum.Point - candidate.Point;
+ previousPoint = candidate;
+ candidate = lineSearchResult.FunctionInfoAtMinimum;
+
+ currentExitCondition = ExitCriteriaSatisfied(candidate, previousPoint, iterations);
+ if (currentExitCondition != MinimizationResult.ExitCondition.None)
+ break;
+ }
+
+ return iterations;
+ }
+
+ protected abstract Vector CalculateSearchDirection(ref Matrix inversePseudoHessian,
+ out double maxLineSearchStep,
+ out double startingStepSize,
+ IObjectiveFunction previousPoint,
+ IObjectiveFunction candidate,
+ Vector step);
+ }
+}
diff --git a/src/Numerics/Optimization/BfgsSolver.cs b/src/Numerics/Optimization/BfgsSolver.cs
index e457ac67..d02c7583 100644
--- a/src/Numerics/Optimization/BfgsSolver.cs
+++ b/src/Numerics/Optimization/BfgsSolver.cs
@@ -28,7 +28,6 @@
// OTHER DEALINGS IN THE SOFTWARE.
//
-
using System;
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearAlgebra.Double;
@@ -56,6 +55,9 @@ namespace MathNet.Numerics.Optimization
/// The minimum found
public static Vector Solve(Vector initialGuess, Func, double> functionValue, Func, Vector> functionGradient)
{
+ var objectiveFunction = ObjectiveFunction.Gradient(functionValue, functionGradient);
+ objectiveFunction.EvaluateAt(initialGuess);
+
int dim = initialGuess.Count;
int iter = 0;
// H represents the approximation of the inverse hessian matrix
@@ -65,18 +67,20 @@ namespace MathNet.Numerics.Optimization
Vector x = initialGuess;
Vector x_old = x;
Vector grad;
-
+ WolfeLineSearch wolfeLineSearch = new WeakWolfeLineSearch(1e-4, 0.9, 1e-5, 200);
do
{
// search along the direction of the gradient
- grad = functionGradient(x);
+ grad = objectiveFunction.Gradient;
Vector p = -1 * H * grad;
- double rate = WolfeRule.LineSearch(x, p, functionValue, functionGradient);
+ var lineSearchResult = wolfeLineSearch.FindConformingStep(objectiveFunction, p, 1.0);
+ double rate = lineSearchResult.FinalStep;
x = x + rate * p;
Vector grad_old = grad;
// update the gradient
- grad = functionGradient(x);
+ objectiveFunction.EvaluateAt(x);
+ grad = objectiveFunction.Gradient;// functionGradient(x);
Vector s = x - x_old;
Vector y = grad - grad_old;
diff --git a/src/Numerics/Optimization/LineSearch/LineSearchResult.cs b/src/Numerics/Optimization/LineSearch/LineSearchResult.cs
index df263be0..2473e891 100644
--- a/src/Numerics/Optimization/LineSearch/LineSearchResult.cs
+++ b/src/Numerics/Optimization/LineSearch/LineSearchResult.cs
@@ -1,4 +1,34 @@
-namespace MathNet.Numerics.Optimization.LineSearch
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2016 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.
+//
+
+namespace MathNet.Numerics.Optimization.LineSearch
{
public class LineSearchResult : MinimizationResult
{
diff --git a/src/Numerics/Optimization/LineSearch/StrongWolfeLineSearch.cs b/src/Numerics/Optimization/LineSearch/StrongWolfeLineSearch.cs
index a2bf0da7..69a8675a 100644
--- a/src/Numerics/Optimization/LineSearch/StrongWolfeLineSearch.cs
+++ b/src/Numerics/Optimization/LineSearch/StrongWolfeLineSearch.cs
@@ -1,5 +1,34 @@
-using System;
-using MathNet.Numerics.LinearAlgebra;
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2016 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;
namespace MathNet.Numerics.Optimization.LineSearch
{
@@ -8,7 +37,7 @@ namespace MathNet.Numerics.Optimization.LineSearch
public StrongWolfeLineSearch(double c1, double c2, double parameterTolerance, int maxIterations = 10)
: base(c1, c2, parameterTolerance, maxIterations)
{
- // Validation in base class
+ // Argument validation in base class
}
protected override MinimizationResult.ExitCondition WolfeExitCondition { get { return MinimizationResult.ExitCondition.StrongWolfeCriteria; } }
diff --git a/src/Numerics/Optimization/LineSearch/WeakWolfeLineSearch.cs b/src/Numerics/Optimization/LineSearch/WeakWolfeLineSearch.cs
index ff5ed3ab..870aaf4e 100644
--- a/src/Numerics/Optimization/LineSearch/WeakWolfeLineSearch.cs
+++ b/src/Numerics/Optimization/LineSearch/WeakWolfeLineSearch.cs
@@ -1,4 +1,34 @@
-using System;
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2016 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;
namespace MathNet.Numerics.Optimization.LineSearch
diff --git a/src/Numerics/Optimization/LineSearch/WolfeLineSearch.cs b/src/Numerics/Optimization/LineSearch/WolfeLineSearch.cs
index e8ff4313..60ab6562 100644
--- a/src/Numerics/Optimization/LineSearch/WolfeLineSearch.cs
+++ b/src/Numerics/Optimization/LineSearch/WolfeLineSearch.cs
@@ -1,8 +1,35 @@
-using MathNet.Numerics.LinearAlgebra;
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2016 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 System;
-using System.Collections.Generic;
-using System.Linq;
-using System.Text;
namespace MathNet.Numerics.Optimization.LineSearch
{
diff --git a/src/Numerics/Optimization/LineSearch/WolfeRule.cs b/src/Numerics/Optimization/LineSearch/WolfeRule.cs
deleted file mode 100644
index 1f30565d..00000000
--- a/src/Numerics/Optimization/LineSearch/WolfeRule.cs
+++ /dev/null
@@ -1,111 +0,0 @@
-//
-// Math.NET Numerics, part of the Math.NET Project
-// http://numerics.mathdotnet.com
-// http://github.com/mathnet/mathnet-numerics
-// http://mathnetnumerics.codeplex.com
-//
-// Copyright (c) 2009-2015 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;
-
-namespace MathNet.Numerics.Optimization.LineSearch
-{
- ///
- /// Performs an inexact line search. This is used as a part of quasi-Newton optimization methods to figure
- /// out how far to move along a certain gradient.
- /// See http://en.wikipedia.org/wiki/Wolfe_conditions
- /// Inspired by implementation: https://github.com/PatWie/CppNumericalSolvers/blob/master/src/linesearch/WolfeRule.h
- ///
- internal static class WolfeRule
- {
- ///
- /// Searches along a line to satisfy the Wolfe conditions (inexact search for minimum)
- ///
- /// Starting point of search
- /// Search direction
- /// Evaluates the function being minimized
- /// Evaluates the gradient of the function
- /// Initial value for the coefficient of z (distance to travel in z direction)
- ///
- public static double LineSearch(
- Vector x0,
- Vector z,
- Func, double> functionValue,
- Func, Vector> functionGradient,
- float alphaInit = 1)
- {
- Vector x = x0;
-
- // evaluate phi(0)
- double phi0 = functionValue(x0);
-
- // evaluate phi'(0)
- Vector grad = functionGradient(x);
- double phi0_dash = z * grad;
-
- double alpha = alphaInit;
-
- bool decrease_direction = true;
-
- // 200 guesses max
- for (int iter = 0; iter < 200; ++iter) {
-
- // new guess for phi(alpha)
- Vector x_candidate = x + alpha * z;
- double phi = functionValue(x_candidate);
-
- // decrease condition invalid --> shrink interval
- if (phi > phi0 + 0.0001 * alpha * phi0_dash)
- {
- alpha *= 0.5;
- decrease_direction = false;
- }
- else
- {
- // valid decrease --> test strong wolfe condition
- Vector grad2 = functionGradient(x_candidate);
- double phi_dash = z * grad2;
-
- // curvature condition invalid ?
- if ((phi_dash < 0.9 * phi0_dash) || !decrease_direction) {
- // increase interval
- alpha *= 4.0;
- }
- else {
- // both condition are valid --> we are happy
- x = x_candidate;
- grad = grad2;
- phi0 = phi;
- return alpha;
- }
- }
- }
-
-
- return alpha;
- }
- }
-}
diff --git a/src/Numerics/Optimization/QuadraticGradientProjectionSearch.cs b/src/Numerics/Optimization/QuadraticGradientProjectionSearch.cs
index 69d30c48..abf97f1d 100644
--- a/src/Numerics/Optimization/QuadraticGradientProjectionSearch.cs
+++ b/src/Numerics/Optimization/QuadraticGradientProjectionSearch.cs
@@ -6,7 +6,7 @@ namespace MathNet.Numerics.Optimization
{
public static class QuadraticGradientProjectionSearch
{
- public static Tuple,int,List> Search(Vector x0, Vector gradient, Matrix hessian, Vector lowerBound, Vector upperBound)
+ public static GradientProjectionResult Search(Vector x0, Vector gradient, Matrix hessian, Vector lowerBound, Vector upperBound)
{
List isFixed = new List(x0.Count);
List breakpoint = new List(x0.Count);
@@ -46,7 +46,7 @@ namespace MathNet.Numerics.Optimization
var maxS = orderedBreakpoint[0];
if (sMin < maxS)
- return Tuple.Create(x + sMin * d, 0,isFixed);
+ return new GradientProjectionResult(x + sMin * d, 0,isFixed);
// while minimum of the last quadratic piece observed is beyond the interval searched
while (true)
@@ -66,7 +66,7 @@ namespace MathNet.Numerics.Optimization
}
if (Double.IsPositiveInfinity(orderedBreakpoint[jj + 1]))
- return Tuple.Create(x, fixedCount, isFixed);
+ return new GradientProjectionResult(x, fixedCount, isFixed);
f1 = gradient * d + (x - x0) * hessian * d;
f2 = d * hessian * d;
@@ -74,13 +74,26 @@ namespace MathNet.Numerics.Optimization
sMin = -f1 / f2;
if (sMin < maxS)
- return Tuple.Create(x + sMin * d, fixedCount, isFixed);
+ return new GradientProjectionResult(x + sMin * d, fixedCount, isFixed);
else if (jj + 1 >= orderedBreakpoint.Count - 1)
{
isFixed[isFixed.Count - 1] = true;
- return Tuple.Create(x + maxS * d, lowerBound.Count, isFixed);
+ return new GradientProjectionResult(x + maxS * d, lowerBound.Count, isFixed);
}
}
}
+
+ public struct GradientProjectionResult
+ {
+ public GradientProjectionResult(Vector cauchyPoint, int fixedCount, List isFixed)
+ {
+ CauchyPoint = cauchyPoint;
+ FixedCount = fixedCount;
+ IsFixed = isFixed;
+ }
+ public Vector CauchyPoint { get; }
+ public int FixedCount { get; }
+ public List IsFixed { get; }
+ }
}
}
diff --git a/src/UnitTests/OptimizationTests/BfgsTest.cs b/src/UnitTests/OptimizationTests/BfgsTest.cs
index 7f612529..b2c17916 100644
--- a/src/UnitTests/OptimizationTests/BfgsTest.cs
+++ b/src/UnitTests/OptimizationTests/BfgsTest.cs
@@ -4,7 +4,7 @@
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
-// Copyright (c) 2009-2015 Math.NET
+// Copyright (c) 2009-2016 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
@@ -28,7 +28,6 @@
// OTHER DEALINGS IN THE SOFTWARE.
//
-using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearAlgebra.Double;
using MathNet.Numerics.Optimization;
using NUnit.Framework;
@@ -50,23 +49,8 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
private static void CheckRosenbrock(double a, double b, double expectedMin)
{
- var x = BfgsSolver.Solve(new DenseVector(new[] { a, b }), Rosenbrock, RosenbrockGradient);
- Numerics.Precision.AlmostEqual(expectedMin, Rosenbrock(x), Precision);
- }
-
- private static double Rosenbrock(Vector x)
- {
- double t1 = (1 - x[0]);
- double t2 = (x[1] - x[0] * x[0]);
- return t1 * t1 + 100 * t2 * t2;
- }
-
- private static Vector RosenbrockGradient(Vector x)
- {
- var grad = new DenseVector(2);
- grad[0] = -2 * (1 - x[0]) + 200 * (x[1] - x[0] * x[0]) * (-2 * x[0]);
- grad[1] = 200 * (x[1] - x[0] * x[0]);
- return grad;
+ var x = BfgsSolver.Solve(new DenseVector(new[] { a, b }), RosenbrockFunction.Value, RosenbrockFunction.Gradient);
+ Numerics.Precision.AlmostEqual(expectedMin, RosenbrockFunction.Value(x), Precision);
}
}
}
diff --git a/src/UnitTests/OptimizationTests/TestBfgsBMinimizer.cs b/src/UnitTests/OptimizationTests/TestBfgsBMinimizer.cs
index eac38e8f..da70b9d1 100644
--- a/src/UnitTests/OptimizationTests/TestBfgsBMinimizer.cs
+++ b/src/UnitTests/OptimizationTests/TestBfgsBMinimizer.cs
@@ -1,3 +1,33 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2016 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.Double;
using MathNet.Numerics.Optimization;
@@ -19,8 +49,8 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
var result = solver.FindMinimum(obj, lowerBound, upperBound, initialGuess);
- Assert.That(Math.Abs(result.MinimizingPoint[0] - 1.0), Is.LessThan(1e-3));
- Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3));
+ 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]
@@ -36,8 +66,8 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
var result = solver.FindMinimum(obj, lowerBound, upperBound, initialGuess);
- Assert.That(Math.Abs(result.MinimizingPoint[0] - 1.0), Is.LessThan(1e-3));
- Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3));
+ 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]
@@ -52,8 +82,8 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
var result = solver.FindMinimum (obj, lowerBound, upperBound, initialGuess);
- Assert.That(Math.Abs(result.MinimizingPoint[0] - 1.0), Is.LessThan(1e-3));
- Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3));
+ 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]
@@ -67,8 +97,8 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
var result = solver.FindMinimum(obj, lowerBound, upperBound, initialGuess);
- Assert.That(Math.Abs(result.MinimizingPoint[0] - 1.0), Is.LessThan(1e-3));
- Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3));
+ 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]
@@ -82,9 +112,43 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
var result = solver.FindMinimum(obj, lowerBound, upperBound, initialGuess);
- Assert.That(Math.Abs(result.MinimizingPoint[0] - 1.0), Is.LessThan(1e-3));
- Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3));
+ 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_MinimumGreateerOrEqualToLowerBoundary()
+ {
+ var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
+ var solver = new BfgsBMinimizer(1e-5, 1e-5, 1e-5, maximumIterations: 1000);
+
+ var lowerBound = new DenseVector(new[] { 2, 2.0 });
+ var upperBound = new DenseVector(new[] { 5.0, 5.0 });
+
+ var initialGuess = new DenseVector(new[] { 2.5, 2.5 });
+
+ var result = solver.FindMinimum(obj, lowerBound, upperBound, initialGuess);
+
+ Assert.GreaterOrEqual(result.MinimizingPoint[0],lowerBound[0]);
+ Assert.GreaterOrEqual(result.MinimizingPoint[1], lowerBound[1]);
+ }
+
+ [Test]
+ public void FindMinimum_Rosenbrock_MinimumLesserOrEqualToUpperBoundary()
+ {
+ var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
+ var solver = new BfgsBMinimizer(1e-5, 1e-5, 1e-5, maximumIterations: 1000);
+
+ var lowerBound = new DenseVector(new[] { -2.0, -2.0 });
+ var upperBound = new DenseVector(new[] { 0.5, 0.5 });
+
+ var initialGuess = new DenseVector(new[] { -0.9, -0.5 });
+
+ var result = solver.FindMinimum(obj, lowerBound, upperBound, initialGuess);
+
+ Assert.LessOrEqual(result.MinimizingPoint[0],upperBound[0]);
+ Assert.LessOrEqual(result.MinimizingPoint[1],upperBound[1]);
+ }
+ }
}
diff --git a/src/UnitTests/OptimizationTests/TestBfgsMinimizer.cs b/src/UnitTests/OptimizationTests/TestBfgsMinimizer.cs
index f321449c..b2300d35 100644
--- a/src/UnitTests/OptimizationTests/TestBfgsMinimizer.cs
+++ b/src/UnitTests/OptimizationTests/TestBfgsMinimizer.cs
@@ -12,7 +12,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
public void FindMinimum_Rosenbrock_Easy()
{
var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
- var solver = new BfgsMinimizer(1e-5, 1e-5, 1000);
+ var solver = new BfgsMinimizer(1e-5, 1e-5, 1e-5, 1000);
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));
@@ -23,7 +23,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
public void FindMinimum_Rosenbrock_Hard()
{
var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
- var solver = new BfgsMinimizer(1e-5, 1e-5, 1000);
+ var solver = new BfgsMinimizer(1e-5, 1e-5, 1e-5, 1000);
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));
@@ -34,7 +34,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
public void FindMinimum_Rosenbrock_Overton()
{
var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
- var solver = new BfgsMinimizer(1e-5, 1e-5, 1000);
+ var solver = new BfgsMinimizer(1e-5, 1e-5, 1e-5, 1000);
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));
@@ -45,7 +45,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
public void FindMinimum_BigRosenbrock_Easy()
{
var obj = ObjectiveFunction.Gradient(BigRosenbrockFunction.Value, BigRosenbrockFunction.Gradient);
- var solver = new BfgsMinimizer(1e-10, 1e-5, 1000);
+ var solver = new BfgsMinimizer(1e-10, 1e-5, 1e-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));
@@ -56,7 +56,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
public void FindMinimum_BigRosenbrock_Hard()
{
var obj = ObjectiveFunction.Gradient(BigRosenbrockFunction.Value, BigRosenbrockFunction.Gradient);
- var solver = new BfgsMinimizer(1e-5, 1e-5, 1000);
+ var solver = new BfgsMinimizer(1e-5, 1e-5, 1e-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));
@@ -67,7 +67,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
public void FindMinimum_BigRosenbrock_Overton()
{
var obj = ObjectiveFunction.Gradient(BigRosenbrockFunction.Value, BigRosenbrockFunction.Gradient);
- var solver = new BfgsMinimizer(1e-5, 1e-5, 1000);
+ var solver = new BfgsMinimizer(1e-5, 1e-5, 1e-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));