forked from tsai/mathnet-numerics
Scott Stephens
13 years ago
committed by
Erik Ovegard
Erik Ovegard
8 changed files with 628 additions and 11 deletions
@ -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<double> lowerBound, Vector<double> upperBound, Vector<double> 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<double> reducedSolution1, reducedGradient, reducedInitialPoint, reducedCauchyPoint, solution1; |
|||
Matrix<double> reducedHessian; |
|||
List<int> reducedMap; |
|||
|
|||
// First step
|
|||
var pseudoHessian = CreateMatrix.DiagonalIdentity<double>(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<int>(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<int>(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<double> ReducedToFull(List<int> reducedMap, Vector<double> reducedVector, Vector<double> fullVector) |
|||
{ |
|||
var output = fullVector.Clone(); |
|||
for (int ii = 0; ii < reducedMap.Count; ++ii) |
|||
output[reducedMap[ii]] = reducedVector[ii]; |
|||
return output; |
|||
} |
|||
|
|||
static double FindMaxStep(Vector<double> startingPoint, Vector<double> searchDirection, Vector<double> lowerBound, Vector<double> 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<double> initialPoint, Vector<double> cauchyPoint, List<bool> isFixed, Vector<double> lowerBound, Vector<double> upperBound, Vector<double> gradient, Matrix<double> pseudoHessian, Vector<double> reducedInitialPoint, Vector<double> reducedCauchyPoint, Vector<double> reducedGradient, Matrix<double> reducedHessian, List<int> 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<double> lowerBound, Vector<double> upperBound, int iterations) |
|||
{ |
|||
Vector<double> 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); |
|||
} |
|||
} |
|||
} |
|||
@ -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<double> searchDirection, double initialStep, double upperBound = Double.PositiveInfinity) |
|||
{ |
|||
double lowerBound = 0.0; |
|||
double step = initialStep; |
|||
|
|||
double initialValue = objective.Value; |
|||
Vector<double> 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<double> 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)); |
|||
} |
|||
} |
|||
} |
|||
|
|||
@ -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<Vector<double>,int,List<bool>> Search(Vector<double> x0, Vector<double> gradient, Matrix<double> hessian, Vector<double> lowerBound, Vector<double> upperBound) |
|||
{ |
|||
List<bool> isFixed = new List<bool>(x0.Count); |
|||
List<double> breakpoint = new List<double>(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<double>(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); |
|||
} |
|||
} |
|||
} |
|||
} |
|||
} |
|||
@ -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)); |
|||
} |
|||
} |
|||
} |
|||
|
|||
Loading…
Reference in new issue