diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj
index 78035dd0..ecfa8180 100644
--- a/src/Numerics/Numerics.csproj
+++ b/src/Numerics/Numerics.csproj
@@ -261,6 +261,7 @@
+
@@ -270,6 +271,7 @@
+
diff --git a/src/Numerics/Optimization/BfgsBMinimizer.cs b/src/Numerics/Optimization/BfgsBMinimizer.cs
new file mode 100644
index 00000000..457b3ce9
--- /dev/null
+++ b/src/Numerics/Optimization/BfgsBMinimizer.cs
@@ -0,0 +1,333 @@
+using System;
+using System.Collections.Generic;
+using MathNet.Numerics.LinearAlgebra;
+using MathNet.Numerics.LinearAlgebra.Double;
+using MathNet.Numerics.Optimization.LineSearch;
+
+namespace MathNet.Numerics.Optimization
+{
+ public class BfgsBMinimizer
+ {
+ 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)
+ {
+ GradientTolerance = gradientTolerance;
+ ParameterTolerance = parameterTolerance;
+ MaximumIterations = maximumIterations;
+ FunctionProgressTolerance = functionProgressTolerance;
+ }
+
+ public MinimizationResult FindMinimum(IObjectiveFunction objective, Vector lowerBound, Vector upperBound, Vector initialGuess)
+ {
+ if (!objective.IsGradientSupported)
+ throw new IncompatibleObjectiveException("Gradient not supported in objective function, but required for BFGS minimization.");
+
+ // Check that dimensions match
+ if (lowerBound.Count != upperBound.Count || lowerBound.Count != initialGuess.Count)
+ throw new ArgumentException("Dimensions of bounds and/or initial guess do not match.");
+
+ // Check that initial guess is feasible
+ for (int ii = 0; ii < initialGuess.Count; ++ii)
+ if (initialGuess[ii] < lowerBound[ii] || initialGuess[ii] > upperBound[ii])
+ throw new ArgumentException("Initial guess is not in the feasible region");
+
+ objective.EvaluateAt(initialGuess);
+
+ // Check that we're not already done
+ MinimizationResult.ExitCondition currentExitCondition = ExitCriteriaSatisfied(objective, null, lowerBound, upperBound, 0);
+ if (currentExitCondition != MinimizationResult.ExitCondition.None)
+ return new MinimizationResult(objective, 0, currentExitCondition);
+
+ // Set up line search algorithm
+ var lineSearcher = new StrongWolfeLineSearch(1e-4, 0.9, Math.Max(ParameterTolerance, 1e-5), maxIterations: 1000);
+
+ // Declare state variables
+ Vector reducedSolution1, reducedGradient, reducedInitialPoint, reducedCauchyPoint, solution1;
+ Matrix reducedHessian;
+ List reducedMap;
+
+ // First step
+ var pseudoHessian = CreateMatrix.DiagonalIdentity(initialGuess.Count);
+
+ // 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 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);
+
+ CreateReducedData(objective.Point, cauchyPoint, isFixed, lowerBound, upperBound, objective.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;
+ }
+
+ var directionFromCauchy = solution1 - cauchyPoint;
+ var maxStepFromCauchyPoint = FindMaxStep(cauchyPoint, directionFromCauchy, lowerBound, upperBound);
+
+ var solution2 = cauchyPoint + Math.Min(maxStepFromCauchyPoint, 1.0)*directionFromCauchy;
+
+ var lineSearchDirection = solution2 - objective.Point;
+ var maxLineSearchStep = FindMaxStep(objective.Point, lineSearchDirection, lowerBound, upperBound);
+ var estStepSize = -objective.Gradient*lineSearchDirection/(lineSearchDirection*pseudoHessian*lineSearchDirection);
+
+ var startingStepSize = Math.Min(Math.Max(estStepSize, 1.0), maxLineSearchStep);
+
+ // Line search
+ LineSearchResult result;
+ try
+ {
+ result = lineSearcher.FindConformingStep(objective, lineSearchDirection, startingStepSize, upperBound: maxLineSearchStep);
+ }
+ catch (Exception e)
+ {
+ throw new InnerOptimizationException("Line search failed.", e);
+ }
+
+ var previousPoint = objective.Fork();
+ var candidatePoint = result.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);
+ }
+
+ // 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);
+
+ 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;
+ }
+
+ directionFromCauchy = solution1 - cauchyPoint;
+ maxStepFromCauchyPoint = FindMaxStep(cauchyPoint, directionFromCauchy, lowerBound, upperBound);
+ //var cauchy_eval = objective.Evaluate(cauchy_point);
+
+ solution2 = cauchyPoint + Math.Min(maxStepFromCauchyPoint, 1.0)*directionFromCauchy;
+
+ lineSearchDirection = solution2 - candidatePoint.Point;
+ maxLineSearchStep = FindMaxStep(candidatePoint.Point, lineSearchDirection, lowerBound, upperBound);
+
+ //line_search_direction = solution1 - candidate_point.Point;
+ //max_line_search_step = FindMaxStep(candidate_point.Point, line_search_direction, lower_bound, upper_bound);
+
+ if (maxLineSearchStep == 0.0)
+ {
+ lineSearchDirection = cauchyPoint - candidatePoint.Point;
+ maxLineSearchStep = FindMaxStep(candidatePoint.Point, lineSearchDirection, lowerBound, upperBound);
+ }
+
+ estStepSize = -candidatePoint.Gradient*lineSearchDirection/(lineSearchDirection*pseudoHessian*lineSearchDirection);
+
+ startingStepSize = Math.Min(Math.Max(estStepSize, 1.0), maxLineSearchStep);
+
+ // 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);
+ }
+
+ iterationsWithNontrivialLineSearch += result.Iterations > 0 ? 1 : 0;
+ totalLineSearchSteps += result.Iterations;
+
+ step = result.FunctionInfoAtMinimum.Point - candidatePoint.Point;
+ previousPoint = candidatePoint;
+ candidatePoint = result.FunctionInfoAtMinimum;
+
+ currentExitCondition = ExitCriteriaSatisfied(candidatePoint, previousPoint, lowerBound, upperBound, iterations);
+ if (currentExitCondition != MinimizationResult.ExitCondition.None)
+ break;
+ }
+
+ if (iterations == MaximumIterations && currentExitCondition == MinimizationResult.ExitCondition.None)
+ throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", MaximumIterations));
+
+ return new MinimizationWithLineSearchResult(candidatePoint, iterations, currentExitCondition, totalLineSearchSteps, iterationsWithNontrivialLineSearch);
+ }
+
+ static Vector ReducedToFull(List reducedMap, Vector reducedVector, Vector fullVector)
+ {
+ var output = fullVector.Clone();
+ for (int ii = 0; ii < reducedMap.Count; ++ii)
+ output[reducedMap[ii]] = reducedVector[ii];
+ return output;
+ }
+
+ static double FindMaxStep(Vector startingPoint, Vector searchDirection, Vector lowerBound, Vector upperBound)
+ {
+ double maxStep = Double.PositiveInfinity;
+ for (int ii = 0; ii < startingPoint.Count; ++ii)
+ {
+ double paramMaxStep;
+ if (searchDirection[ii] > 0)
+ paramMaxStep = (upperBound[ii] - startingPoint[ii])/searchDirection[ii];
+ else if (searchDirection[ii] < 0)
+ paramMaxStep = (startingPoint[ii] - lowerBound[ii])/-searchDirection[ii];
+ else
+ paramMaxStep = Double.PositiveInfinity;
+
+ if (paramMaxStep < maxStep)
+ maxStep = paramMaxStep;
+ }
+ 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)
+ {
+ int ll = 0;
+ for (int ii = 0; ii < lowerBound.Count; ++ii)
+ {
+ if (!isFixed[ii])
+ {
+ // hessian
+ int mm = 0;
+ for (int jj = 0; jj < lowerBound.Count; ++jj)
+ {
+ if (!isFixed[jj])
+ {
+ reducedHessian[ll, mm++] = pseudoHessian[ii, jj];
+ }
+ }
+
+ // gradient
+ reducedInitialPoint[ll] = initialPoint[ii];
+ reducedCauchyPoint[ll] = cauchyPoint[ii];
+ reducedGradient[ll] = gradient[ii];
+ ll += 1;
+ reducedMap.Add(ii);
+
+ }
+ }
+ }
+
+ 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)
+ {
+ 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);
+ }
+ }
+}
diff --git a/src/Numerics/Optimization/BfgsMinimizer.cs b/src/Numerics/Optimization/BfgsMinimizer.cs
index 0258c50f..e240fbcc 100644
--- a/src/Numerics/Optimization/BfgsMinimizer.cs
+++ b/src/Numerics/Optimization/BfgsMinimizer.cs
@@ -33,7 +33,7 @@ namespace MathNet.Numerics.Optimization
return new MinimizationResult(objective, 0, currentExitCondition);
// Set up line search algorithm
- var lineSearcher = new WeakWolfeLineSearch(1e-4, 0.9, ParameterTolerance, 1000);
+ var lineSearcher = new WeakWolfeLineSearch(1e-4, 0.9, Math.Max(ParameterTolerance, 1e-10), 1000);
// First step
var inversePseudoHessian = CreateMatrix.DenseIdentity(initialGuess.Count);
@@ -61,6 +61,7 @@ namespace MathNet.Numerics.Optimization
stepSize = result.FinalStep;
// Subsequent steps
+ Matrix I = CreateMatrix.DiagonalIdentity(initialGuess.Count);
int iterations;
int totalLineSearchSteps = result.Iterations;
int iterationsWithNontrivialLineSearch = result.Iterations > 0 ? 0 : 1;
@@ -70,14 +71,18 @@ namespace MathNet.Numerics.Optimization
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 >= -GradientTolerance*GradientTolerance)
+ if (searchDirection * objective.Gradient >= 0.0)
{
searchDirection = -objective.Gradient;
inversePseudoHessian = CreateMatrix.DenseIdentity(initialGuess.Count);
}
+ //else if (searchDirection * objective.Gradient >= -GradientTolerance*GradientTolerance)
+ //{
+ // searchDirection = -objective.Gradient;
+ // inversePseudoHessian = CreateMatrix.DenseIdentity(initialGuess.Count);
+ //}
previousGradient = objective.Gradient;
previousPoint = objective.Point;
diff --git a/src/Numerics/Optimization/LineSearch/StrongWolfeLineSearch.cs b/src/Numerics/Optimization/LineSearch/StrongWolfeLineSearch.cs
index 236c63c9..b543ba48 100644
--- a/src/Numerics/Optimization/LineSearch/StrongWolfeLineSearch.cs
+++ b/src/Numerics/Optimization/LineSearch/StrongWolfeLineSearch.cs
@@ -1,11 +1,111 @@
using System;
-using System.Collections.Generic;
-using System.Linq;
-using System.Text;
+using MathNet.Numerics.LinearAlgebra;
-namespace MathNet.Numerics.Optimization
+namespace MathNet.Numerics.Optimization.LineSearch
{
- class StrongWolfeLineSearch
+ public class StrongWolfeLineSearch
{
+ 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)
+ {
+ C1 = c1;
+ C2 = c2;
+ ParameterTolerance = parameterTolerance;
+ MaximumIterations = maxIterations;
+ }
+
+ // 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);
+ }
+ }
+ }
+
+ bool IsFinite(double x)
+ {
+ return !(Double.IsNaN(x) || Double.IsInfinity(x));
+ }
}
}
diff --git a/src/Numerics/Optimization/QuadraticGradientProjectionSearch.cs b/src/Numerics/Optimization/QuadraticGradientProjectionSearch.cs
new file mode 100644
index 00000000..69d30c48
--- /dev/null
+++ b/src/Numerics/Optimization/QuadraticGradientProjectionSearch.cs
@@ -0,0 +1,86 @@
+using System;
+using System.Collections.Generic;
+using MathNet.Numerics.LinearAlgebra;
+
+namespace MathNet.Numerics.Optimization
+{
+ public static class QuadraticGradientProjectionSearch
+ {
+ public static Tuple,int,List> Search(Vector x0, Vector gradient, Matrix hessian, Vector lowerBound, Vector upperBound)
+ {
+ List isFixed = new List(x0.Count);
+ List breakpoint = new List(x0.Count);
+ for (int ii = 0; ii < x0.Count; ++ii)
+ {
+ breakpoint.Add(0.0);
+ isFixed.Add(false);
+ if (gradient[ii] < 0)
+ breakpoint[ii] = (x0[ii] - upperBound[ii]) / gradient[ii];
+ else if (gradient[ii] > 0)
+ breakpoint[ii] = (x0[ii] - lowerBound[ii]) / gradient[ii];
+ else
+ {
+ if (Math.Abs(x0[ii] - upperBound[ii]) < 100 * Double.Epsilon || Math.Abs(x0[ii] - lowerBound[ii]) < 100 * Double.Epsilon)
+ breakpoint[ii] = 0.0;
+ else
+ breakpoint[ii] = Double.PositiveInfinity;
+ }
+ }
+
+ var orderedBreakpoint = new List(x0.Count);
+ orderedBreakpoint.AddRange(breakpoint);
+ orderedBreakpoint.Sort();
+
+ // Compute initial state variables
+ var d = -gradient;
+ for (int ii = 0; ii < d.Count; ++ii)
+ if (breakpoint[ii] <= 0.0)
+ d[ii] *= 0.0;
+
+
+ int jj = -1;
+ var x = x0;
+ var f1 = gradient * d;
+ var f2 = 0.5 * d * hessian * d;
+ var sMin = -f1 / f2;
+ var maxS = orderedBreakpoint[0];
+
+ if (sMin < maxS)
+ return Tuple.Create(x + sMin * d, 0,isFixed);
+
+ // while minimum of the last quadratic piece observed is beyond the interval searched
+ while (true)
+ {
+ // update data to the beginning of the interval we're searching
+ jj += 1;
+ x = x + d * maxS;
+ maxS = orderedBreakpoint[jj+1] - orderedBreakpoint[jj];
+
+ int fixedCount = 0;
+ for (int ii = 0; ii < d.Count; ++ii)
+ if (orderedBreakpoint[jj] >= breakpoint[ii])
+ {
+ d[ii] *= 0.0;
+ isFixed[ii] = true;
+ fixedCount += 1;
+ }
+
+ if (Double.IsPositiveInfinity(orderedBreakpoint[jj + 1]))
+ return Tuple.Create(x, fixedCount, isFixed);
+
+ f1 = gradient * d + (x - x0) * hessian * d;
+ f2 = d * hessian * d;
+
+ sMin = -f1 / f2;
+
+ if (sMin < maxS)
+ return Tuple.Create(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);
+ }
+ }
+ }
+ }
+}
diff --git a/src/UnitTests/OptimizationTests/TestBfgsBMinimizer.cs b/src/UnitTests/OptimizationTests/TestBfgsBMinimizer.cs
new file mode 100644
index 00000000..eac38e8f
--- /dev/null
+++ b/src/UnitTests/OptimizationTests/TestBfgsBMinimizer.cs
@@ -0,0 +1,90 @@
+using System;
+using MathNet.Numerics.LinearAlgebra.Double;
+using MathNet.Numerics.Optimization;
+using NUnit.Framework;
+
+namespace MathNet.Numerics.UnitTests.OptimizationTests
+{
+ [TestFixture]
+ public class TestBfgsBMinimizer
+ {
+ [Test]
+ public void FindMinimum_Rosenbrock_Easy()
+ {
+ 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[]{ -5.0, -5.0 });
+ var upperBound = new DenseVector(new[]{ 5.0, 5.0 });
+ var initialGuess = new DenseVector(new[] { 1.2, 1.2 });
+
+ 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));
+ }
+
+ [Test]
+ public void FindMinimum_Rosenbrock_Hard()
+ {
+ 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[]{ -5.0, -5.0 });
+ var upperBound = new DenseVector(new[]{ 5.0, 5.0 });
+
+ var initialGuess = new DenseVector (new[]{ -1.2, 1.0 });
+
+ 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));
+ }
+
+ [Test]
+ public void FindMinimum_Rosenbrock_Overton()
+ {
+ 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[]{ -5.0, -5.0 });
+ var upperBound = new DenseVector(new[]{ 5.0, 5.0 });
+ var initialGuess = new DenseVector (new[]{ -0.9, -0.5 });
+
+ 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));
+ }
+
+ [Test]
+ public void FindMinimum_Rosenbrock_Easy_OneBoundary()
+ {
+ 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[]{ 1.0, -5.0 });
+ var upperBound = new DenseVector(new[]{ 5.0, 5.0 });
+ var initialGuess = new DenseVector(new[] { 1.2, 1.2 });
+
+ 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));
+ }
+
+ [Test]
+ public void FindMinimum_Rosenbrock_Easy_TwoBoundaries()
+ {
+ 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[]{ 1.0, 1.0 });
+ var upperBound = new DenseVector(new[]{ 5.0, 5.0 });
+ var initialGuess = new DenseVector(new[] { 1.2, 1.2 });
+
+ 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));
+ }
+ }
+}
+
diff --git a/src/UnitTests/OptimizationTests/TestBfgsMinimizer.cs b/src/UnitTests/OptimizationTests/TestBfgsMinimizer.cs
index 281a5b7a..5486f4b1 100644
--- a/src/UnitTests/OptimizationTests/TestBfgsMinimizer.cs
+++ b/src/UnitTests/OptimizationTests/TestBfgsMinimizer.cs
@@ -44,7 +44,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
[Test]
public void FindMinimum_BigRosenbrock_Easy()
{
- var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
+ var obj = ObjectiveFunction.Gradient(BigRosenbrockFunction.Value, BigRosenbrockFunction.Gradient);
var solver = new BfgsMinimizer(1e-10, 1e-5, 1000);
var result = solver.FindMinimum(obj, new DenseVector(new[] { 1.2*100.0, 1.2*100.0 }));
@@ -55,7 +55,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
[Test]
public void FindMinimum_BigRosenbrock_Hard()
{
- var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
+ var obj = ObjectiveFunction.Gradient(BigRosenbrockFunction.Value, BigRosenbrockFunction.Gradient);
var solver = new BfgsMinimizer(1e-5, 1e-5, 1000);
var result = solver.FindMinimum(obj, new DenseVector(new[] { -1.2*100.0, 1.0*100.0 }));
@@ -66,7 +66,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
[Test]
public void FindMinimum_BigRosenbrock_Overton()
{
- var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
+ var obj = ObjectiveFunction.Gradient(BigRosenbrockFunction.Value, BigRosenbrockFunction.Gradient);
var solver = new BfgsMinimizer(1e-5, 1e-5, 1000);
var result = solver.FindMinimum(obj, new DenseVector(new[] { -0.9*100.0, -0.5*100.0 }));
diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj
index 608647a0..1ca11826 100644
--- a/src/UnitTests/UnitTests.csproj
+++ b/src/UnitTests/UnitTests.csproj
@@ -425,6 +425,7 @@
+