From b04465e6ec9cc2ba029c76189b342d548d7e2a7a Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Thu, 14 May 2015 22:47:17 +0200 Subject: [PATCH] Optimization: clarify objective function forking and mutation, fixes unit tests --- .../Optimization/IObjectiveFunction.cs | 8 +++- .../Implementation/WeakWolfeLineSearch.cs | 17 ++++---- src/Numerics/Optimization/NewtonMinimizer.cs | 39 +++++++++++-------- 3 files changed, 37 insertions(+), 27 deletions(-) diff --git a/src/Numerics/Optimization/IObjectiveFunction.cs b/src/Numerics/Optimization/IObjectiveFunction.cs index d78bee72..f257db48 100644 --- a/src/Numerics/Optimization/IObjectiveFunction.cs +++ b/src/Numerics/Optimization/IObjectiveFunction.cs @@ -2,9 +2,8 @@ namespace MathNet.Numerics.Optimization { - public interface IObjectiveFunction + public interface IObjectiveFunctionEvaluation { - void EvaluateAt(Vector point); IObjectiveFunction Fork(); Vector Point { get; } @@ -16,4 +15,9 @@ namespace MathNet.Numerics.Optimization bool IsHessianSupported { get; } Matrix Hessian { get; } } + + public interface IObjectiveFunction : IObjectiveFunctionEvaluation + { + void EvaluateAt(Vector point); + } } diff --git a/src/Numerics/Optimization/Implementation/WeakWolfeLineSearch.cs b/src/Numerics/Optimization/Implementation/WeakWolfeLineSearch.cs index 8d291ab1..cc6e3c44 100644 --- a/src/Numerics/Optimization/Implementation/WeakWolfeLineSearch.cs +++ b/src/Numerics/Optimization/Implementation/WeakWolfeLineSearch.cs @@ -19,8 +19,9 @@ namespace MathNet.Numerics.Optimization.Implementation } // Implemented following http://www.math.washington.edu/~burke/crs/408/lectures/L9-weak-Wolfe.pdf - public LineSearchOutput FindConformingStep(IObjectiveFunction objective, IObjectiveFunction startingPoint, Vector searchDirection, double initialStep) + public LineSearchOutput FindConformingStep(IObjectiveFunctionEvaluation startingPoint, Vector searchDirection, double initialStep) { + var objective = startingPoint.Fork(); if (!(objective is CheckedObjectiveFunction)) { objective = new CheckedObjectiveFunction(objective, ValidateValue, ValidateGradient, null); @@ -30,21 +31,21 @@ namespace MathNet.Numerics.Optimization.Implementation double upperBound = Double.PositiveInfinity; double step = initialStep; + Vector initialPoint = startingPoint.Point; double initialValue = startingPoint.Value; Vector initialGradient = startingPoint.Gradient; double initialDd = searchDirection * initialGradient; int ii; - IObjectiveFunction candidateEval = objective.Fork(); MinimizationOutput.ExitCondition reasonForExit = MinimizationOutput.ExitCondition.None; for (ii = 0; ii < _maximumIterations; ++ii) { - candidateEval.EvaluateAt(startingPoint.Point + searchDirection * step); + objective.EvaluateAt(initialPoint + searchDirection * step); - double stepDd = searchDirection * candidateEval.Gradient; + double stepDd = searchDirection * objective.Gradient; - if (candidateEval.Value > initialValue + _c1 * step * initialDd) + if (objective.Value > initialValue + _c1 * step * initialDd) { upperBound = step; step = 0.5 * (lowerBound + upperBound); @@ -63,9 +64,9 @@ namespace MathNet.Numerics.Optimization.Implementation if (!Double.IsInfinity(upperBound)) { double maxRelChange = 0.0; - for (int jj = 0; jj < candidateEval.Point.Count; ++jj) + for (int jj = 0; jj < objective.Point.Count; ++jj) { - double tmp = Math.Abs(searchDirection[jj] * (upperBound - lowerBound)) / Math.Max(Math.Abs(candidateEval.Point[jj]), 1.0); + double tmp = Math.Abs(searchDirection[jj] * (upperBound - lowerBound)) / Math.Max(Math.Abs(objective.Point[jj]), 1.0); maxRelChange = Math.Max(maxRelChange, tmp); } if (maxRelChange < _parameterTolerance) @@ -86,7 +87,7 @@ namespace MathNet.Numerics.Optimization.Implementation throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", _maximumIterations)); } - return new LineSearchOutput(candidateEval, ii, step, reasonForExit); + return new LineSearchOutput(objective, ii, step, reasonForExit); } bool Conforms(IObjectiveFunction startingPoint, Vector searchDirection, double step, IObjectiveFunction endingPoint) diff --git a/src/Numerics/Optimization/NewtonMinimizer.cs b/src/Numerics/Optimization/NewtonMinimizer.cs index 69cb7193..b419b7c7 100644 --- a/src/Numerics/Optimization/NewtonMinimizer.cs +++ b/src/Numerics/Optimization/NewtonMinimizer.cs @@ -21,38 +21,41 @@ namespace MathNet.Numerics.Optimization public MinimizationOutput FindMinimum(IObjectiveFunction objective, Vector initialGuess) { if (!objective.IsGradientSupported) + { throw new IncompatibleObjectiveException("Gradient not supported in objective function, but required for Newton minimization."); + } if (!objective.IsHessianSupported) + { throw new IncompatibleObjectiveException("Hessian not supported in objective function, but required for Newton minimization."); + } if (!(objective is CheckedObjectiveFunction)) + { objective = new CheckedObjectiveFunction(objective, ValidateObjective, ValidateGradient, ValidateHessian); - - IObjectiveFunction initialEval = objective.Fork(); - initialEval.EvaluateAt(initialGuess); + } // Check that we're not already done - if (ExitCriteriaSatisfied(initialGuess, initialEval.Gradient)) - return new MinimizationOutput(initialEval, 0, MinimizationOutput.ExitCondition.AbsoluteGradient); + objective.EvaluateAt(initialGuess); + if (ExitCriteriaSatisfied(objective.Gradient)) + { + return new MinimizationOutput(objective, 0, MinimizationOutput.ExitCondition.AbsoluteGradient); + } // Set up line search algorithm var lineSearcher = new WeakWolfeLineSearch(1e-4, 0.9, 1e-4, maxIterations: 1000); - // Declare state variables - IObjectiveFunction candidatePoint = initialEval; - // Subsequent steps int iterations = 0; int totalLineSearchSteps = 0; int iterationsWithNontrivialLineSearch = 0; bool tmpLineSearch = false; - while (!ExitCriteriaSatisfied(candidatePoint.Point, candidatePoint.Gradient) && iterations < MaximumIterations) + while (!ExitCriteriaSatisfied(objective.Gradient) && iterations < MaximumIterations) { - var searchDirection = candidatePoint.Hessian.LU().Solve(-candidatePoint.Gradient); - if (searchDirection * candidatePoint.Gradient >= 0) + var searchDirection = objective.Hessian.LU().Solve(-objective.Gradient); + if (searchDirection * objective.Gradient >= 0) { - searchDirection = -candidatePoint.Gradient; + searchDirection = -objective.Gradient; tmpLineSearch = true; } @@ -61,7 +64,7 @@ namespace MathNet.Numerics.Optimization LineSearchOutput result; try { - result = lineSearcher.FindConformingStep(objective, candidatePoint, searchDirection, 1.0); + result = lineSearcher.FindConformingStep(objective.Fork(), searchDirection, 1.0); } catch (Exception e) { @@ -69,11 +72,11 @@ namespace MathNet.Numerics.Optimization } iterationsWithNontrivialLineSearch += result.Iterations > 0 ? 1 : 0; totalLineSearchSteps += result.Iterations; - candidatePoint = result.FunctionInfoAtMinimum; + objective = result.FunctionInfoAtMinimum; } else { - candidatePoint.EvaluateAt(candidatePoint.Point + searchDirection); + objective.EvaluateAt(objective.Point + searchDirection); } tmpLineSearch = false; @@ -82,12 +85,14 @@ namespace MathNet.Numerics.Optimization } if (iterations == MaximumIterations) + { throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", MaximumIterations)); + } - return new MinimizationWithLineSearchOutput(candidatePoint, iterations, MinimizationOutput.ExitCondition.AbsoluteGradient, totalLineSearchSteps, iterationsWithNontrivialLineSearch); + return new MinimizationWithLineSearchOutput(objective, iterations, MinimizationOutput.ExitCondition.AbsoluteGradient, totalLineSearchSteps, iterationsWithNontrivialLineSearch); } - private bool ExitCriteriaSatisfied(Vector candidatePoint, Vector gradient) + private bool ExitCriteriaSatisfied(Vector gradient) { return gradient.Norm(2.0) < GradientTolerance; }