From e7f2894c7ab4a6af575fdad47802016d43dec554 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 4 Jan 2014 13:33:54 +0100 Subject: [PATCH] Optimization: cosmetics --- src/Numerics/Numerics.csproj | 4 +- src/Numerics/Optimization/BfgsMinimizer.cs | 205 ++++++++++-------- .../Optimization/BisectionRootFinder.cs | 122 +++++++---- .../ConjugateGradientMinimizer.cs | 130 ++++++----- src/Numerics/Optimization/Exceptions.cs | 85 +++++--- src/Numerics/Optimization/ExitCondition.cs | 41 +++- .../Optimization/GoldenSectionMinimizer.cs | 75 +++++-- src/Numerics/Optimization/LineSearchOutput.cs | 39 +++- .../LineSearchingMinimizerOutput.cs | 20 -- .../Optimization/MinimizationOutput.cs | 51 ++++- .../Optimization/MinimizationOutput1D.cs | 47 +++- .../MinimizationWithLineSearchOutput.cs | 45 ++++ src/Numerics/Optimization/NewtonMinimizer.cs | 121 ++++++----- src/Numerics/Optimization/ObjectiveChecker.cs | 104 +++++---- .../Optimization/ObjectiveChecker1D.cs | 106 +++++---- .../Optimization/ObjectiveFunction.cs | 167 ++++++++------ .../Optimization/ObjectiveFunction1D.cs | 108 ++++++--- .../Optimization/OptimizationResult.cs | 8 - .../Optimization/StrongWolfeLineSearch.cs | 11 - .../Optimization/WeakWolfeLineSearch.cs | 160 ++++++++------ .../OptimizationTests/RosenbrockFunction.cs | 48 +++- .../OptimizationTests/TestBfgsMinimizer.cs | 45 +++- .../TestBisectionRootFinder.cs | 42 +++- .../TestConjugateGradientMinimizer.cs | 44 +++- .../OptimizationTests/TestNewtonMinimizer.cs | 48 +++- .../TestRosenbrockFunction.cs | 75 +++++-- 26 files changed, 1267 insertions(+), 684 deletions(-) delete mode 100644 src/Numerics/Optimization/LineSearchingMinimizerOutput.cs create mode 100644 src/Numerics/Optimization/MinimizationWithLineSearchOutput.cs delete mode 100644 src/Numerics/Optimization/OptimizationResult.cs delete mode 100644 src/Numerics/Optimization/StrongWolfeLineSearch.cs diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index dfb48a3a..beee4f21 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -200,7 +200,7 @@ - + @@ -209,8 +209,6 @@ - - diff --git a/src/Numerics/Optimization/BfgsMinimizer.cs b/src/Numerics/Optimization/BfgsMinimizer.cs index 9e82f8e7..89ac26c6 100644 --- a/src/Numerics/Optimization/BfgsMinimizer.cs +++ b/src/Numerics/Optimization/BfgsMinimizer.cs @@ -1,152 +1,173 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; +// +// 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-2013 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 { - public class BfgsMinimizer { public double GradientTolerance { get; set; } - public double ParameterTolerance { get; set; } + public double ParameterTolerance { get; set; } public int MaximumIterations { get; set; } - public BfgsMinimizer(double gradient_tolerance, double parameter_tolerance, int maximum_iterations=1000) + public BfgsMinimizer(double gradientTolerance, double parameterTolerance, int maximumIterations = 1000) { - this.GradientTolerance = gradient_tolerance; - this.ParameterTolerance = parameter_tolerance; - this.MaximumIterations = maximum_iterations; + GradientTolerance = gradientTolerance; + ParameterTolerance = parameterTolerance; + MaximumIterations = maximumIterations; } - public MinimizationOutput FindMinimum(IObjectiveFunction objective, Vector initial_guess) + public MinimizationOutput FindMinimum(IObjectiveFunction objective, Vector initialGuess) { if (!objective.GradientSupported) throw new IncompatibleObjectiveException("Gradient not supported in objective function, but required for BFGS minimization."); if (!(objective is ObjectiveChecker)) - objective = new ObjectiveChecker(objective, this.ValidateObjective, this.ValidateGradient, null); + objective = new ObjectiveChecker(objective, ValidateObjective, ValidateGradient, null); + + IEvaluation initialEval = objective.Evaluate(initialGuess); - IEvaluation initial_eval = objective.Evaluate(initial_guess); - // Check that we're not already done - ExitCondition current_exit_condition = this.ExitCriteriaSatisfied(initial_eval, null); - if (current_exit_condition != ExitCondition.None) - return new MinimizationOutput(initial_eval, 0, current_exit_condition); - + ExitCondition currentExitCondition = ExitCriteriaSatisfied(initialEval, null); + if (currentExitCondition != ExitCondition.None) + return new MinimizationOutput(initialEval, 0, currentExitCondition); + // Set up line search algorithm - var line_searcher = new WeakWolfeLineSearch(1e-4, 0.9,this.ParameterTolerance, max_iterations:1000); + var lineSearcher = new WeakWolfeLineSearch(1e-4, 0.9, ParameterTolerance, maxIterations: 1000); // Declare state variables - IEvaluation candidate_point, previous_point; - double step_size; - Vector gradient, step, search_direction; - Matrix inverse_pseudo_hessian; + Vector gradient; // First step - inverse_pseudo_hessian = Matrix.Build.DiagonalIdentity(initial_guess.Count); - search_direction = -initial_eval.Gradient; - step_size = 100 * this.GradientTolerance / (search_direction * search_direction); - + Matrix inversePseudoHessian = Matrix.Build.DiagonalIdentity(initialGuess.Count); + Vector searchDirection = -initialEval.Gradient; + double stepSize = 100*GradientTolerance/(searchDirection*searchDirection); + LineSearchOutput result; - try + try { - result = line_searcher.FindConformingStep(objective, initial_eval, search_direction, step_size); - } - catch (Exception e) + result = lineSearcher.FindConformingStep(objective, initialEval, searchDirection, stepSize); + } + catch (Exception e) { throw new InnerOptimizationException("Line search failed.", e); } - previous_point = initial_eval; - candidate_point = result.FunctionInfoAtMinimum; - gradient = candidate_point.Gradient; - step = candidate_point.Point - initial_guess; - step_size = result.FinalStep; + IEvaluation previousPoint = initialEval; + IEvaluation candidatePoint = result.FunctionInfoAtMinimum; + gradient = candidatePoint.Gradient; + Vector step = candidatePoint.Point - initialGuess; + stepSize = result.FinalStep; // Subsequent steps int iterations; - int total_line_search_steps = result.Iterations; - int iterations_with_nontrivial_line_search = result.Iterations > 0 ? 0 : 1; - int steepest_descent_resets = 0; - for (iterations = 1; iterations < this.MaximumIterations; ++iterations) + int totalLineSearchSteps = result.Iterations; + int iterationsWithNontrivialLineSearch = result.Iterations > 0 ? 0 : 1; + for (iterations = 1; iterations < MaximumIterations; ++iterations) { - var y = candidate_point.Gradient - previous_point.Gradient; + var y = candidatePoint.Gradient - previousPoint.Gradient; - double sy = step * y; - 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); + 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); - search_direction = -inverse_pseudo_hessian * candidate_point.Gradient; + searchDirection = -inversePseudoHessian*candidatePoint.Gradient; - if (search_direction * candidate_point.Gradient >= -this.GradientTolerance*this.GradientTolerance) + if (searchDirection*candidatePoint.Gradient >= -GradientTolerance*GradientTolerance) { - search_direction = -candidate_point.Gradient; - inverse_pseudo_hessian = Matrix.Build.DiagonalIdentity(initial_guess.Count); - steepest_descent_resets += 1; + searchDirection = -candidatePoint.Gradient; + inversePseudoHessian = Matrix.Build.DiagonalIdentity(initialGuess.Count); } try { - result = line_searcher.FindConformingStep(objective, candidate_point, search_direction, 1.0); + result = lineSearcher.FindConformingStep(objective, candidatePoint, searchDirection, 1.0); } catch (Exception e) { throw new InnerOptimizationException("Line search failed.", e); } - iterations_with_nontrivial_line_search += result.Iterations > 0 ? 1 : 0; - total_line_search_steps += result.Iterations; + iterationsWithNontrivialLineSearch += result.Iterations > 0 ? 1 : 0; + totalLineSearchSteps += result.Iterations; - step_size = result.FinalStep; - step = result.FunctionInfoAtMinimum.Point - candidate_point.Point; - previous_point = candidate_point; - candidate_point = result.FunctionInfoAtMinimum; + stepSize = result.FinalStep; + step = result.FunctionInfoAtMinimum.Point - candidatePoint.Point; + previousPoint = candidatePoint; + candidatePoint = result.FunctionInfoAtMinimum; - current_exit_condition = this.ExitCriteriaSatisfied(candidate_point, previous_point); - if (current_exit_condition != ExitCondition.None) - break; + currentExitCondition = ExitCriteriaSatisfied(candidatePoint, previousPoint); + if (currentExitCondition != ExitCondition.None) + break; } - if (iterations == this.MaximumIterations) - throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", this.MaximumIterations)); + if (iterations == MaximumIterations) + throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", MaximumIterations)); - return new MinimizationWithLineSearchOutput(candidate_point, iterations, current_exit_condition, total_line_search_steps, iterations_with_nontrivial_line_search); + return new MinimizationWithLineSearchOutput(candidatePoint, iterations, currentExitCondition, totalLineSearchSteps, iterationsWithNontrivialLineSearch); } - private ExitCondition ExitCriteriaSatisfied(IEvaluation candidate_point, IEvaluation last_point) + ExitCondition ExitCriteriaSatisfied(IEvaluation candidatePoint, IEvaluation lastPoint) { - Vector rel_grad = new MathNet.Numerics.LinearAlgebra.Double.DenseVector(candidate_point.Point.Count); - double relative_gradient = 0.0; - double normalizer = Math.Max(Math.Abs(candidate_point.Value),1.0); - for (int ii = 0; ii < rel_grad.Count; ++ii) - { - double tmp = candidate_point.Gradient[ii]*Math.Max(Math.Abs(candidate_point.Point[ii]), 1.0) / normalizer; - relative_gradient = Math.Max(relative_gradient, Math.Abs(tmp)); - } - if (relative_gradient < this.GradientTolerance) - { - return ExitCondition.RelativeGradient; - } - - if (last_point != null) - { - double most_progress = 0.0; - for (int ii = 0; ii < candidate_point.Point.Count; ++ii) - { - var tmp = Math.Abs(candidate_point.Point[ii] - last_point.Point[ii])/Math.Max(Math.Abs(last_point.Point[ii]),1.0); - most_progress = Math.Max(most_progress, tmp); - } - if ( most_progress < this.ParameterTolerance ) - { - return ExitCondition.LackOfProgress; - } - } - - return ExitCondition.None; + 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 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 ExitCondition.LackOfProgress; + } + } + + return ExitCondition.None; } - private void ValidateGradient(IEvaluation eval) + void ValidateGradient(IEvaluation eval) { foreach (var x in eval.Gradient) { @@ -155,7 +176,7 @@ namespace MathNet.Numerics.Optimization } } - private void ValidateObjective(IEvaluation eval) + void ValidateObjective(IEvaluation 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/BisectionRootFinder.cs b/src/Numerics/Optimization/BisectionRootFinder.cs index d3c196cb..0e0d2911 100644 --- a/src/Numerics/Optimization/BisectionRootFinder.cs +++ b/src/Numerics/Optimization/BisectionRootFinder.cs @@ -1,82 +1,108 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; +// +// 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-2013 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 { public class BisectionRootFinder { - public double ObjectiveTolerance { get; set; } public double XTolerance { get; set; } public double LowerExpansionFactor { get; set; } public double UpperExpansionFactor { get; set; } public int MaxExpansionSteps { get; set; } - public BisectionRootFinder(double objective_tolerance=1e-5, double x_tolerance=1e-5, double lower_expansion_factor=-1.0, double upper_expansion_factor=-1.0, int max_expansion_steps=10) + public BisectionRootFinder(double objectiveTolerance = 1e-5, double xTolerance = 1e-5, double lowerExpansionFactor = -1.0, double upperExpansionFactor = -1.0, int maxExpansionSteps = 10) { - this.ObjectiveTolerance = objective_tolerance; - this.XTolerance = x_tolerance; - this.LowerExpansionFactor = lower_expansion_factor; - this.UpperExpansionFactor = upper_expansion_factor; - this.MaxExpansionSteps = max_expansion_steps; - } + ObjectiveTolerance = objectiveTolerance; + XTolerance = xTolerance; + LowerExpansionFactor = lowerExpansionFactor; + UpperExpansionFactor = upperExpansionFactor; + MaxExpansionSteps = maxExpansionSteps; + } - public double FindRoot(Func objective_function, double lower_bound, double upper_bound) + public double FindRoot(Func objectiveFunction, double lowerBound, double upperBound) { - double lower_val = objective_function(lower_bound); - double upper_val = objective_function(upper_bound); + double lowerVal = objectiveFunction(lowerBound); + double upperVal = objectiveFunction(upperBound); - if (lower_val == 0.0) - return lower_bound; - if (upper_val == 0.0) - return upper_bound; + if (lowerVal == 0.0) + return lowerBound; + if (upperVal == 0.0) + return upperBound; - this.ValidateEvaluation(lower_val, lower_bound); - this.ValidateEvaluation(upper_val, upper_bound); + ValidateEvaluation(lowerVal, lowerBound); + ValidateEvaluation(upperVal, upperBound); - if (Math.Sign(lower_val) == Math.Sign(upper_val) && this.LowerExpansionFactor <= 1.0 && this.UpperExpansionFactor <= 1.0) + if (Math.Sign(lowerVal) == Math.Sign(upperVal) && LowerExpansionFactor <= 1.0 && UpperExpansionFactor <= 1.0) throw new Exception("Bounds do not necessarily span a root, and StepExpansionFactor is not set to expand the interval in this case."); - int expansion_steps = 0; - while (Math.Sign(lower_val) == Math.Sign(upper_val) && expansion_steps < this.MaxExpansionSteps) + int expansionSteps = 0; + while (Math.Sign(lowerVal) == Math.Sign(upperVal) && expansionSteps < MaxExpansionSteps) { - double midpoint = 0.5 * (upper_bound + lower_bound); - double range = upper_bound - lower_bound; - if (this.UpperExpansionFactor <= 0.0 || (this.LowerExpansionFactor > 0.0 && Math.Abs(lower_val) < Math.Abs(upper_val)) ) + double midpoint = 0.5*(upperBound + lowerBound); + double range = upperBound - lowerBound; + if (UpperExpansionFactor <= 0.0 || (LowerExpansionFactor > 0.0 && Math.Abs(lowerVal) < Math.Abs(upperVal))) { - lower_bound = upper_bound - this.LowerExpansionFactor * range; - lower_val = objective_function(lower_bound); - this.ValidateEvaluation(lower_val, lower_bound); + lowerBound = upperBound - LowerExpansionFactor*range; + lowerVal = objectiveFunction(lowerBound); + ValidateEvaluation(lowerVal, lowerBound); } else { - upper_bound = lower_bound + this.UpperExpansionFactor * range; - upper_val = objective_function(upper_bound); - this.ValidateEvaluation(upper_val, upper_bound); + upperBound = lowerBound + UpperExpansionFactor*range; + upperVal = objectiveFunction(upperBound); + ValidateEvaluation(upperVal, upperBound); } - expansion_steps += 1; + expansionSteps += 1; } - if (Math.Sign(lower_val) == Math.Sign(upper_val) && expansion_steps == this.MaxExpansionSteps) + if (Math.Sign(lowerVal) == Math.Sign(upperVal) && expansionSteps == MaxExpansionSteps) throw new MaximumIterationsException("Could not bound root in maximum expansion iterations."); - while (Math.Abs(upper_val - lower_val) > 0.5 * this.ObjectiveTolerance || Math.Abs(upper_bound - lower_bound) > 0.5 * this.XTolerance) + while (Math.Abs(upperVal - lowerVal) > 0.5*ObjectiveTolerance || Math.Abs(upperBound - lowerBound) > 0.5*XTolerance) { - double midpoint = 0.5 * (upper_bound + lower_bound); - double midval = objective_function(midpoint); - this.ValidateEvaluation(midval, midpoint); + double midpoint = 0.5*(upperBound + lowerBound); + double midval = objectiveFunction(midpoint); + ValidateEvaluation(midval, midpoint); - if (Math.Sign(midval) == Math.Sign(lower_val)) + if (Math.Sign(midval) == Math.Sign(lowerVal)) { - lower_bound = midpoint; - lower_val = midval; + lowerBound = midpoint; + lowerVal = midval; } - else if (Math.Sign(midval) == Math.Sign(upper_val)) + else if (Math.Sign(midval) == Math.Sign(upperVal)) { - upper_bound = midpoint; - upper_val = midval; + upperBound = midpoint; + upperVal = midval; } else { @@ -84,16 +110,16 @@ namespace MathNet.Numerics.Optimization } } - return 0.5 * (lower_bound + upper_bound); + return 0.5*(lowerBound + upperBound); } - private void ValidateEvaluation(double output, double input) + void ValidateEvaluation(double output, double input) { if (!IsFinite(output)) throw new Exception(String.Format("Objective function returned non-finite result: f({0}) = {1}", input, output)); } - private static bool IsFinite(double x) + static bool IsFinite(double x) { return !(Double.IsInfinity(x) || Double.IsNaN(x)); } diff --git a/src/Numerics/Optimization/ConjugateGradientMinimizer.cs b/src/Numerics/Optimization/ConjugateGradientMinimizer.cs index 3bcc509f..2bf778a8 100644 --- a/src/Numerics/Optimization/ConjugateGradientMinimizer.cs +++ b/src/Numerics/Optimization/ConjugateGradientMinimizer.cs @@ -1,8 +1,34 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; - +// +// 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-2013 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 @@ -12,109 +38,105 @@ namespace MathNet.Numerics.Optimization public double GradientTolerance { get; set; } public int MaximumIterations { get; set; } - public ConjugateGradientMinimizer(double gradient_tolerance, int maximum_iterations) + public ConjugateGradientMinimizer(double gradientTolerance, int maximumIterations) { - this.GradientTolerance = gradient_tolerance; - this.MaximumIterations = maximum_iterations; + GradientTolerance = gradientTolerance; + MaximumIterations = maximumIterations; } - public MinimizationOutput FindMinimum(IObjectiveFunction objective, Vector initial_guess) + public MinimizationOutput FindMinimum(IObjectiveFunction objective, Vector initialGuess) { if (!objective.GradientSupported) throw new IncompatibleObjectiveException("Gradient not supported in objective function, but required for ConjugateGradient minimization."); if (!(objective is ObjectiveChecker)) - objective = new ObjectiveChecker(objective, this.ValidateObjective, this.ValidateGradient, null); + objective = new ObjectiveChecker(objective, ValidateObjective, ValidateGradient, null); + + IEvaluation initialEval = objective.Evaluate(initialGuess); + var gradient = initialEval.Gradient; - IEvaluation initial_eval = objective.Evaluate(initial_guess); - var gradient = initial_eval.Gradient; - // Check that we're not already done - if (this.ExitCriteriaSatisfied(initial_guess, gradient)) - return new MinimizationOutput(initial_eval, 0, ExitCondition.AbsoluteGradient); - + if (ExitCriteriaSatisfied(initialGuess, gradient)) + return new MinimizationOutput(initialEval, 0, ExitCondition.AbsoluteGradient); + // Set up line search algorithm - var line_searcher = new WeakWolfeLineSearch(1e-4, 0.1, 1e-4, max_iterations:1000); + var lineSearcher = new WeakWolfeLineSearch(1e-4, 0.1, 1e-4, maxIterations: 1000); // Declare state variables - IEvaluation candidate_point; - Vector steepest_direction, previous_steepest_direction, search_direction; // First step - steepest_direction = -gradient; - search_direction = steepest_direction; - double initial_step_size = 100 * this.GradientTolerance / (gradient * gradient); + Vector steepestDirection = -gradient; + Vector searchDirection = steepestDirection; + double initialStepSize = 100*GradientTolerance/(gradient*gradient); LineSearchOutput result; - try + try { - result = line_searcher.FindConformingStep(objective, initial_eval, search_direction, initial_step_size); - } - catch (Exception e) + result = lineSearcher.FindConformingStep(objective, initialEval, searchDirection, initialStepSize); + } + catch (Exception e) { throw new InnerOptimizationException("Line search failed.", e); } - candidate_point = result.FunctionInfoAtMinimum; + IEvaluation candidatePoint = result.FunctionInfoAtMinimum; + + double stepSize = result.FinalStep; - double step_size = result.FinalStep; - // Subsequent steps int iterations = 1; - int total_line_search_steps = result.Iterations; - int iterations_with_nontrivial_line_search = result.Iterations > 0 ? 0 : 1; - int steepest_descent_resets = 0; - while (!this.ExitCriteriaSatisfied(candidate_point.Point, candidate_point.Gradient) && iterations < this.MaximumIterations) + int totalLineSearchSteps = result.Iterations; + int iterationsWithNontrivialLineSearch = result.Iterations > 0 ? 0 : 1; + while (!ExitCriteriaSatisfied(candidatePoint.Point, candidatePoint.Gradient) && iterations < MaximumIterations) { - previous_steepest_direction = steepest_direction; - steepest_direction = -candidate_point.Gradient; - var search_direction_adjuster = Math.Max(0,steepest_direction * (steepest_direction - previous_steepest_direction) / (previous_steepest_direction * previous_steepest_direction)); - + Vector previousSteepestDirection = steepestDirection; + steepestDirection = -candidatePoint.Gradient; + var searchDirectionAdjuster = Math.Max(0, steepestDirection*(steepestDirection - previousSteepestDirection)/(previousSteepestDirection*previousSteepestDirection)); + //double prev_grad_mag = previous_steepest_direction*previous_steepest_direction; //double grad_overlap = steepest_direction*previous_steepest_direction; //double search_grad_overlap = candidate_point.Gradient*search_direction; //if (iterations % initial_guess.Count == 0 || (Math.Abs(grad_overlap) >= 0.2 * prev_grad_mag) || (-2 * prev_grad_mag >= search_grad_overlap) || (search_grad_overlap >= -0.2 * prev_grad_mag)) // search_direction = steepest_direction; - //else + //else // search_direction = steepest_direction + search_direction_adjuster * search_direction; - search_direction = steepest_direction + search_direction_adjuster * search_direction; - if (search_direction * candidate_point.Gradient >= 0) + searchDirection = steepestDirection + searchDirectionAdjuster*searchDirection; + if (searchDirection*candidatePoint.Gradient >= 0) { - search_direction = steepest_direction; - steepest_descent_resets += 1; + searchDirection = steepestDirection; } try { - result = line_searcher.FindConformingStep(objective, candidate_point, search_direction, step_size); + result = lineSearcher.FindConformingStep(objective, candidatePoint, searchDirection, stepSize); } catch (Exception e) { throw new InnerOptimizationException("Line search failed.", e); } - iterations_with_nontrivial_line_search += result.Iterations > 0 ? 1 : 0; - total_line_search_steps += result.Iterations; + iterationsWithNontrivialLineSearch += result.Iterations > 0 ? 1 : 0; + totalLineSearchSteps += result.Iterations; - step_size = result.FinalStep; - candidate_point = result.FunctionInfoAtMinimum; + stepSize = result.FinalStep; + candidatePoint = result.FunctionInfoAtMinimum; iterations += 1; } - if (iterations == this.MaximumIterations) - throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", this.MaximumIterations)); + if (iterations == MaximumIterations) + throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", MaximumIterations)); - return new MinimizationWithLineSearchOutput(candidate_point, iterations, ExitCondition.AbsoluteGradient, total_line_search_steps, iterations_with_nontrivial_line_search); + return new MinimizationWithLineSearchOutput(candidatePoint, iterations, ExitCondition.AbsoluteGradient, totalLineSearchSteps, iterationsWithNontrivialLineSearch); } - private bool ExitCriteriaSatisfied(Vector candidate_point, Vector gradient) + bool ExitCriteriaSatisfied(Vector candidatePoint, Vector gradient) { - return gradient.Norm(2.0) < this.GradientTolerance; + return gradient.Norm(2.0) < GradientTolerance; } - private void ValidateGradient(IEvaluation eval) + void ValidateGradient(IEvaluation eval) { foreach (var x in eval.Gradient) { @@ -123,7 +145,7 @@ namespace MathNet.Numerics.Optimization } } - private void ValidateObjective(IEvaluation eval) + void ValidateObjective(IEvaluation 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/Exceptions.cs b/src/Numerics/Optimization/Exceptions.cs index aab4f690..216004d8 100644 --- a/src/Numerics/Optimization/Exceptions.cs +++ b/src/Numerics/Optimization/Exceptions.cs @@ -1,25 +1,56 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; -using MathNet.Numerics.LinearAlgebra; -using MathNet.Numerics.LinearAlgebra.Double; +// +// 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-2013 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 { public class OptimizationException : Exception { public OptimizationException(string message) - : base(message) {} + : base(message) + { + } - public OptimizationException(string message, Exception inner_exception) - : base(message, inner_exception) { } + public OptimizationException(string message, Exception innerException) + : base(message, innerException) + { + } } public class MaximumIterationsException : OptimizationException { - public MaximumIterationsException(string message) - : base(message) {} + public MaximumIterationsException(string message) + : base(message) + { + } } public class EvaluationException : OptimizationException @@ -29,25 +60,25 @@ namespace MathNet.Numerics.Optimization public EvaluationException(string message, IEvaluation eval) : base(message) { - this.Evaluation = eval; + Evaluation = eval; } - public EvaluationException(string message, IEvaluation eval, Exception inner_exception) - : base(message, inner_exception) + public EvaluationException(string message, IEvaluation eval, Exception innerException) + : base(message, innerException) { - this.Evaluation = eval; + Evaluation = eval; } public EvaluationException(string message, IEvaluation1D eval) : base(message) { - this.Evaluation = new OneDEvaluationExpander(eval); + Evaluation = new OneDEvaluationExpander(eval); } - public EvaluationException(string message, IEvaluation1D eval, Exception inner_exception) - : base(message, inner_exception) + public EvaluationException(string message, IEvaluation1D eval, Exception innerException) + : base(message, innerException) { - this.Evaluation = new OneDEvaluationExpander(eval); + Evaluation = new OneDEvaluationExpander(eval); } } @@ -55,17 +86,21 @@ namespace MathNet.Numerics.Optimization public class InnerOptimizationException : OptimizationException { public InnerOptimizationException(string message) - : base(message) {} + : base(message) + { + } - public InnerOptimizationException(string message, Exception inner_exception) - : base(message, inner_exception) { } + public InnerOptimizationException(string message, Exception innerException) + : base(message, innerException) + { + } } public class IncompatibleObjectiveException : OptimizationException { public IncompatibleObjectiveException(string message) - : base(message) {} + : base(message) + { + } } - - } diff --git a/src/Numerics/Optimization/ExitCondition.cs b/src/Numerics/Optimization/ExitCondition.cs index e17e14bd..8a16dee1 100644 --- a/src/Numerics/Optimization/ExitCondition.cs +++ b/src/Numerics/Optimization/ExitCondition.cs @@ -1,7 +1,42 @@ -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-2013 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 { - public enum ExitCondition { None, RelativeGradient, LackOfProgress, AbsoluteGradient, WeakWolfeCriteria, BoundTolerance } + public enum ExitCondition + { + None, + RelativeGradient, + LackOfProgress, + AbsoluteGradient, + WeakWolfeCriteria, + BoundTolerance + } } - diff --git a/src/Numerics/Optimization/GoldenSectionMinimizer.cs b/src/Numerics/Optimization/GoldenSectionMinimizer.cs index 538610ff..20cf7cfa 100644 --- a/src/Numerics/Optimization/GoldenSectionMinimizer.cs +++ b/src/Numerics/Optimization/GoldenSectionMinimizer.cs @@ -1,7 +1,34 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; +// +// 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-2013 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 { @@ -10,33 +37,33 @@ namespace MathNet.Numerics.Optimization public double XTolerance { get; set; } public int MaximumIterations { get; set; } - public GoldenSectionMinimizer(double x_tolerance=1e-5, int max_iterations=1000) + public GoldenSectionMinimizer(double xTolerance = 1e-5, int maxIterations = 1000) { - this.XTolerance = x_tolerance; - this.MaximumIterations = max_iterations; + XTolerance = xTolerance; + MaximumIterations = maxIterations; } - public MinimizationOutput1D FindMinimum(IObjectiveFunction1D objective, double lower_bound, double upper_bound) + public MinimizationOutput1D FindMinimum(IObjectiveFunction1D objective, double lowerBound, double upperBound) { if (!(objective is ObjectiveChecker1D)) - objective = new ObjectiveChecker1D(objective, this.ValueChecker, null, null); + objective = new ObjectiveChecker1D(objective, ValueChecker, null, null); + + double middlePointX = lowerBound + (upperBound - lowerBound)/(1 + GoldenRatio); + IEvaluation1D lower = objective.Evaluate(lowerBound); + IEvaluation1D middle = objective.Evaluate(middlePointX); + IEvaluation1D upper = objective.Evaluate(upperBound); - double middle_point_x = lower_bound + (upper_bound - lower_bound) / (1 + _golden_ratio); - IEvaluation1D lower = objective.Evaluate(lower_bound); - IEvaluation1D middle = objective.Evaluate(middle_point_x); - IEvaluation1D upper = objective.Evaluate(upper_bound); - - if (upper_bound <= lower_bound) + if (upperBound <= lowerBound) throw new OptimizationException("Lower bound must be lower than upper bound."); if (upper.Value < middle.Value || lower.Value < middle.Value) throw new OptimizationException("Lower and upper bounds do not necessarily bound a minimum."); int iterations = 0; - while (Math.Abs(upper.Point - lower.Point) > this.XTolerance && iterations < this.MaximumIterations) + while (Math.Abs(upper.Point - lower.Point) > XTolerance && iterations < MaximumIterations) { - double test_x = lower.Point + (upper.Point - middle.Point); - var test = objective.Evaluate(test_x); + double testX = lower.Point + (upper.Point - middle.Point); + var test = objective.Evaluate(testX); if (test.Point < middle.Point) { @@ -66,18 +93,18 @@ namespace MathNet.Numerics.Optimization iterations += 1; } - if (iterations == this.MaximumIterations) + if (iterations == MaximumIterations) throw new MaximumIterationsException("Max iterations reached."); - else - return new MinimizationOutput1D(middle, iterations, ExitCondition.BoundTolerance); - + + return new MinimizationOutput1D(middle, iterations, ExitCondition.BoundTolerance); } - private void ValueChecker(IEvaluation1D eval) + void ValueChecker(IEvaluation1D eval) { if (Double.IsNaN(eval.Value) || Double.IsInfinity(eval.Value)) throw new EvaluationException("Objective function returned non-finite value.", eval); } - private static double _golden_ratio = (1.0 + Math.Sqrt(5)) / 2.0; + + static readonly double GoldenRatio = (1.0 + Math.Sqrt(5))/2.0; } } diff --git a/src/Numerics/Optimization/LineSearchOutput.cs b/src/Numerics/Optimization/LineSearchOutput.cs index 0659f60e..274e897a 100644 --- a/src/Numerics/Optimization/LineSearchOutput.cs +++ b/src/Numerics/Optimization/LineSearchOutput.cs @@ -1,7 +1,32 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; +// +// 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-2013 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 { @@ -9,10 +34,10 @@ namespace MathNet.Numerics.Optimization { public double FinalStep { get; private set; } - public LineSearchOutput(IEvaluation function_info, int iterations, double final_step, ExitCondition reason_for_exit) - : base(function_info, iterations, reason_for_exit) + public LineSearchOutput(IEvaluation functionInfo, int iterations, double finalStep, ExitCondition reasonForExit) + : base(functionInfo, iterations, reasonForExit) { - this.FinalStep = final_step; + FinalStep = finalStep; } } } diff --git a/src/Numerics/Optimization/LineSearchingMinimizerOutput.cs b/src/Numerics/Optimization/LineSearchingMinimizerOutput.cs deleted file mode 100644 index 2cd00ad9..00000000 --- a/src/Numerics/Optimization/LineSearchingMinimizerOutput.cs +++ /dev/null @@ -1,20 +0,0 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; - -namespace MathNet.Numerics.Optimization -{ - public class MinimizationWithLineSearchOutput : MinimizationOutput - { - public int TotalLineSearchIterations { get; private set; } - public int IterationsWithNonTrivialLineSearch { get; private set; } - - public MinimizationWithLineSearchOutput(IEvaluation function_info, int iterations, ExitCondition reason_for_exit, int total_line_search_iterations, int iterations_with_non_trivial_line_search) - : base(function_info, iterations, reason_for_exit) - { - this.TotalLineSearchIterations = total_line_search_iterations; - this.IterationsWithNonTrivialLineSearch = iterations_with_non_trivial_line_search; - } - } -} diff --git a/src/Numerics/Optimization/MinimizationOutput.cs b/src/Numerics/Optimization/MinimizationOutput.cs index d4ef5104..a02ce30a 100644 --- a/src/Numerics/Optimization/MinimizationOutput.cs +++ b/src/Numerics/Optimization/MinimizationOutput.cs @@ -1,7 +1,32 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; +// +// 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-2013 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; @@ -9,16 +34,20 @@ namespace MathNet.Numerics.Optimization { public class MinimizationOutput { - public Vector MinimizingPoint { get { return FunctionInfoAtMinimum.Point; } } + public Vector MinimizingPoint + { + get { return FunctionInfoAtMinimum.Point; } + } + public IEvaluation FunctionInfoAtMinimum { get; private set; } - public int Iterations { get; private set; } - public ExitCondition ReasonForExit { get; private set; } + public int Iterations { get; private set; } + public ExitCondition ReasonForExit { get; private set; } - public MinimizationOutput(IEvaluation function_info, int iterations, ExitCondition reason_for_exit) + public MinimizationOutput(IEvaluation functionInfo, int iterations, ExitCondition reasonForExit) { - this.FunctionInfoAtMinimum = function_info; - this.Iterations = iterations; - this.ReasonForExit = reason_for_exit; + FunctionInfoAtMinimum = functionInfo; + Iterations = iterations; + ReasonForExit = reasonForExit; } } } diff --git a/src/Numerics/Optimization/MinimizationOutput1D.cs b/src/Numerics/Optimization/MinimizationOutput1D.cs index dd68bbfc..4a4068e0 100644 --- a/src/Numerics/Optimization/MinimizationOutput1D.cs +++ b/src/Numerics/Optimization/MinimizationOutput1D.cs @@ -1,22 +1,51 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; +// +// 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-2013 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 { public class MinimizationOutput1D { - public double MinimizingPoint { get { return FunctionInfoAtMinimum.Point; } } + public double MinimizingPoint + { + get { return FunctionInfoAtMinimum.Point; } + } + public IEvaluation1D FunctionInfoAtMinimum { get; private set; } public int Iterations { get; private set; } public ExitCondition ReasonForExit { get; private set; } - public MinimizationOutput1D(IEvaluation1D function_info, int iterations, ExitCondition reason_for_exit) + public MinimizationOutput1D(IEvaluation1D functionInfo, int iterations, ExitCondition reasonForExit) { - this.FunctionInfoAtMinimum = function_info; - this.Iterations = iterations; - this.ReasonForExit = reason_for_exit; + FunctionInfoAtMinimum = functionInfo; + Iterations = iterations; + ReasonForExit = reasonForExit; } } } diff --git a/src/Numerics/Optimization/MinimizationWithLineSearchOutput.cs b/src/Numerics/Optimization/MinimizationWithLineSearchOutput.cs new file mode 100644 index 00000000..000bd171 --- /dev/null +++ b/src/Numerics/Optimization/MinimizationWithLineSearchOutput.cs @@ -0,0 +1,45 @@ +// +// 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-2013 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 +{ + public class MinimizationWithLineSearchOutput : MinimizationOutput + { + public int TotalLineSearchIterations { get; private set; } + public int IterationsWithNonTrivialLineSearch { get; private set; } + + public MinimizationWithLineSearchOutput(IEvaluation functionInfo, int iterations, ExitCondition reasonForExit, int totalLineSearchIterations, int iterationsWithNonTrivialLineSearch) + : base(functionInfo, iterations, reasonForExit) + { + TotalLineSearchIterations = totalLineSearchIterations; + IterationsWithNonTrivialLineSearch = iterationsWithNonTrivialLineSearch; + } + } +} diff --git a/src/Numerics/Optimization/NewtonMinimizer.cs b/src/Numerics/Optimization/NewtonMinimizer.cs index c6b329ad..843289bb 100644 --- a/src/Numerics/Optimization/NewtonMinimizer.cs +++ b/src/Numerics/Optimization/NewtonMinimizer.cs @@ -1,9 +1,35 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; +// +// 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-2013 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 LU = MathNet.Numerics.LinearAlgebra.Factorization.LU; namespace MathNet.Numerics.Optimization { @@ -13,14 +39,14 @@ namespace MathNet.Numerics.Optimization public int MaximumIterations { get; set; } public bool UseLineSearch { get; set; } - public NewtonMinimizer(double gradient_tolerance, int maximum_iterations, bool use_line_search=false) + public NewtonMinimizer(double gradientTolerance, int maximumIterations, bool useLineSearch = false) { - this.GradientTolerance = gradient_tolerance; - this.MaximumIterations = maximum_iterations; - this.UseLineSearch = use_line_search; + GradientTolerance = gradientTolerance; + MaximumIterations = maximumIterations; + UseLineSearch = useLineSearch; } - public MinimizationOutput FindMinimum(IObjectiveFunction objective, Vector initial_guess) + public MinimizationOutput FindMinimum(IObjectiveFunction objective, Vector initialGuess) { if (!objective.GradientSupported) throw new IncompatibleObjectiveException("Gradient not supported in objective function, but required for Newton minimization."); @@ -29,76 +55,73 @@ namespace MathNet.Numerics.Optimization throw new IncompatibleObjectiveException("Hessian not supported in objective function, but required for Newton minimization."); if (!(objective is ObjectiveChecker)) - objective = new ObjectiveChecker(objective, this.ValidateObjective, this.ValidateGradient, this.ValidateHessian); + objective = new ObjectiveChecker(objective, ValidateObjective, ValidateGradient, ValidateHessian); + + IEvaluation initialEval = objective.Evaluate(initialGuess); - IEvaluation initial_eval = objective.Evaluate(initial_guess); - // Check that we're not already done - if (this.ExitCriteriaSatisfied(initial_guess, initial_eval.Gradient)) - return new MinimizationOutput(initial_eval, 0, ExitCondition.AbsoluteGradient); - - // Set up line search algorithm - var line_searcher = new WeakWolfeLineSearch(1e-4, 0.9, 1e-4, max_iterations:1000); + if (ExitCriteriaSatisfied(initialGuess, initialEval.Gradient)) + return new MinimizationOutput(initialEval, 0, ExitCondition.AbsoluteGradient); + + // Set up line search algorithm + var lineSearcher = new WeakWolfeLineSearch(1e-4, 0.9, 1e-4, maxIterations: 1000); // Declare state variables - IEvaluation candidate_point = initial_eval; - Vector search_direction; - LineSearchOutput result; + IEvaluation candidatePoint = initialEval; // Subsequent steps int iterations = 0; - int total_line_search_steps = 0; - int iterations_with_nontrivial_line_search = 0; - int steepest_descent_resets = 0; - bool tmp_line_search = false; - while (!this.ExitCriteriaSatisfied(candidate_point.Point, candidate_point.Gradient) && iterations < this.MaximumIterations) + int totalLineSearchSteps = 0; + int iterationsWithNontrivialLineSearch = 0; + bool tmpLineSearch = false; + while (!ExitCriteriaSatisfied(candidatePoint.Point, candidatePoint.Gradient) && iterations < MaximumIterations) { - search_direction = candidate_point.Hessian.LU().Solve(-candidate_point.Gradient); + Vector searchDirection = candidatePoint.Hessian.LU().Solve(-candidatePoint.Gradient); - if (search_direction * candidate_point.Gradient >= 0) + if (searchDirection*candidatePoint.Gradient >= 0) { - search_direction = -candidate_point.Gradient; - steepest_descent_resets += 1; - tmp_line_search = true; + searchDirection = -candidatePoint.Gradient; + tmpLineSearch = true; } - if (this.UseLineSearch || tmp_line_search) + if (UseLineSearch || tmpLineSearch) { + LineSearchOutput result; try { - result = line_searcher.FindConformingStep(objective, candidate_point, search_direction, 1.0); + result = lineSearcher.FindConformingStep(objective, candidatePoint, searchDirection, 1.0); } catch (Exception e) { throw new InnerOptimizationException("Line search failed.", e); } - iterations_with_nontrivial_line_search += result.Iterations > 0 ? 1 : 0; - total_line_search_steps += result.Iterations; - candidate_point = result.FunctionInfoAtMinimum; + iterationsWithNontrivialLineSearch += result.Iterations > 0 ? 1 : 0; + totalLineSearchSteps += result.Iterations; + candidatePoint = result.FunctionInfoAtMinimum; } else { - candidate_point = objective.Evaluate(candidate_point.Point + search_direction); - } - - tmp_line_search = false; + candidatePoint = objective.Evaluate(candidatePoint.Point + searchDirection); + } + + tmpLineSearch = false; iterations += 1; } - if (iterations == this.MaximumIterations) - throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", this.MaximumIterations)); + if (iterations == MaximumIterations) + throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", MaximumIterations)); - return new MinimizationWithLineSearchOutput(candidate_point, iterations, ExitCondition.AbsoluteGradient, total_line_search_steps, iterations_with_nontrivial_line_search); + return new MinimizationWithLineSearchOutput(candidatePoint, iterations, ExitCondition.AbsoluteGradient, totalLineSearchSteps, iterationsWithNontrivialLineSearch); } - private bool ExitCriteriaSatisfied(Vector candidate_point, Vector gradient) + bool ExitCriteriaSatisfied(Vector candidatePoint, Vector gradient) { - return gradient.Norm(2.0) < this.GradientTolerance; + return gradient.Norm(2.0) < GradientTolerance; } - private void ValidateGradient(IEvaluation eval) + void ValidateGradient(IEvaluation eval) { foreach (var x in eval.Gradient) { @@ -107,19 +130,19 @@ namespace MathNet.Numerics.Optimization } } - private void ValidateObjective(IEvaluation eval) + void ValidateObjective(IEvaluation eval) { if (Double.IsNaN(eval.Value) || Double.IsInfinity(eval.Value)) throw new EvaluationException("Non-finite objective function returned.", eval); } - private void ValidateHessian(IEvaluation eval) + void ValidateHessian(IEvaluation eval) { for (int ii = 0; ii < eval.Hessian.RowCount; ++ii) { for (int jj = 0; jj < eval.Hessian.ColumnCount; ++jj) { - if (Double.IsNaN(eval.Hessian[ii,jj]) || Double.IsInfinity(eval.Hessian[ii,jj])) + if (Double.IsNaN(eval.Hessian[ii, jj]) || Double.IsInfinity(eval.Hessian[ii, jj])) throw new EvaluationException("Non-finite Hessian returned.", eval); } } diff --git a/src/Numerics/Optimization/ObjectiveChecker.cs b/src/Numerics/Optimization/ObjectiveChecker.cs index 53e9a011..6bcc2fc4 100644 --- a/src/Numerics/Optimization/ObjectiveChecker.cs +++ b/src/Numerics/Optimization/ObjectiveChecker.cs @@ -1,50 +1,80 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; +// +// 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-2013 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 { public class CheckedEvaluation : IEvaluation { - private ObjectiveChecker Checker; + readonly ObjectiveChecker Checker; public IEvaluation InnerEvaluation { get; private set; } - private bool ValueChecked; - private bool GradientChecked; - private bool HessianChecked; + bool ValueChecked; + bool GradientChecked; + bool HessianChecked; public CheckedEvaluation(ObjectiveChecker checker, IEvaluation evaluation) { - this.Checker = checker; - this.InnerEvaluation = evaluation; + Checker = checker; + InnerEvaluation = evaluation; } public Vector Point { - get { return this.InnerEvaluation.Point; } + get { return InnerEvaluation.Point; } + } + + public EvaluationStatus Status + { + get { return InnerEvaluation.Status; } } - public EvaluationStatus Status { get { return this.InnerEvaluation.Status; } } public double Value { get { - - if (!this.ValueChecked) + if (!ValueChecked) { double tmp; try { - tmp = this.InnerEvaluation.Value; + tmp = InnerEvaluation.Value; } catch (Exception e) { - throw new EvaluationException("Objective function evaluation failed.", this.InnerEvaluation, e); + throw new EvaluationException("Objective function evaluation failed.", InnerEvaluation, e); } - this.Checker.ValueChecker(this.InnerEvaluation); + Checker.ValueChecker(InnerEvaluation); } - return this.InnerEvaluation.Value; + return InnerEvaluation.Value; } } @@ -52,21 +82,20 @@ namespace MathNet.Numerics.Optimization { get { - - if (!this.GradientChecked) + if (!GradientChecked) { Vector tmp; try { - tmp = this.InnerEvaluation.Gradient; + tmp = InnerEvaluation.Gradient; } catch (Exception e) { - throw new EvaluationException("Objective gradient evaluation failed.", this.InnerEvaluation, e); + throw new EvaluationException("Objective gradient evaluation failed.", InnerEvaluation, e); } - this.Checker.GradientChecker(this.InnerEvaluation); + Checker.GradientChecker(InnerEvaluation); } - return this.InnerEvaluation.Gradient; + return InnerEvaluation.Gradient; } } @@ -74,21 +103,20 @@ namespace MathNet.Numerics.Optimization { get { - - if (!this.HessianChecked) + if (!HessianChecked) { Matrix tmp; try { - tmp = this.InnerEvaluation.Hessian; + tmp = InnerEvaluation.Hessian; } catch (Exception e) { - throw new EvaluationException("Objective hessian evaluation failed.", this.InnerEvaluation, e); + throw new EvaluationException("Objective hessian evaluation failed.", InnerEvaluation, e); } - this.Checker.HessianChecker(InnerEvaluation); + Checker.HessianChecker(InnerEvaluation); } - return this.InnerEvaluation.Hessian; + return InnerEvaluation.Hessian; } } } @@ -100,29 +128,29 @@ namespace MathNet.Numerics.Optimization public Action GradientChecker { get; private set; } public Action HessianChecker { get; private set; } - public ObjectiveChecker(IObjectiveFunction objective, Action value_checker, Action gradient_checker, Action hessian_checker) + public ObjectiveChecker(IObjectiveFunction objective, Action valueChecker, Action gradientChecker, Action hessianChecker) { - this.InnerObjective = objective; - this.ValueChecker = value_checker; - this.GradientChecker = gradient_checker; - this.HessianChecker = hessian_checker; + InnerObjective = objective; + ValueChecker = valueChecker; + GradientChecker = gradientChecker; + HessianChecker = hessianChecker; } public bool GradientSupported { - get { return this.InnerObjective.GradientSupported; } + get { return InnerObjective.GradientSupported; } } public bool HessianSupported { - get { return this.InnerObjective.HessianSupported; } + get { return InnerObjective.HessianSupported; } } public IEvaluation Evaluate(Vector point) { try { - return new CheckedEvaluation(this, this.InnerObjective.Evaluate(point)); + return new CheckedEvaluation(this, InnerObjective.Evaluate(point)); } catch (Exception e) { diff --git a/src/Numerics/Optimization/ObjectiveChecker1D.cs b/src/Numerics/Optimization/ObjectiveChecker1D.cs index 49c489cd..4b6a98d1 100644 --- a/src/Numerics/Optimization/ObjectiveChecker1D.cs +++ b/src/Numerics/Optimization/ObjectiveChecker1D.cs @@ -1,51 +1,79 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; -using MathNet.Numerics.LinearAlgebra.Double; +// +// 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-2013 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 { public class CheckedEvaluation1D : IEvaluation1D { - private ObjectiveChecker1D Checker; - private IEvaluation1D InnerEvaluation; - private bool ValueChecked; - private bool DerivativeChecked; - private bool SecondDerivativeChecked; + readonly ObjectiveChecker1D Checker; + readonly IEvaluation1D InnerEvaluation; + bool ValueChecked; + bool DerivativeChecked; + bool SecondDerivativeChecked; public CheckedEvaluation1D(ObjectiveChecker1D checker, IEvaluation1D evaluation) { - this.Checker = checker; - this.InnerEvaluation = evaluation; + Checker = checker; + InnerEvaluation = evaluation; } public double Point { - get { return this.InnerEvaluation.Point; } + get { return InnerEvaluation.Point; } } - public EvaluationStatus Status { get { return this.InnerEvaluation.Status; } } + public EvaluationStatus Status + { + get { return InnerEvaluation.Status; } + } public double Value { get { - - if (!this.ValueChecked) + if (!ValueChecked) { double tmp; try { - tmp = this.InnerEvaluation.Value; + tmp = InnerEvaluation.Value; } catch (Exception e) { - throw new EvaluationException("Objective function evaluation failed.", this.InnerEvaluation, e); + throw new EvaluationException("Objective function evaluation failed.", InnerEvaluation, e); } - this.Checker.ValueChecker(this.InnerEvaluation); + Checker.ValueChecker(InnerEvaluation); } - return this.InnerEvaluation.Value; + return InnerEvaluation.Value; } } @@ -53,21 +81,20 @@ namespace MathNet.Numerics.Optimization { get { - - if (!this.DerivativeChecked) + if (!DerivativeChecked) { double tmp; try { - tmp = this.InnerEvaluation.Derivative; + tmp = InnerEvaluation.Derivative; } catch (Exception e) { - throw new EvaluationException("Objective derivative evaluation failed.", this.InnerEvaluation, e); + throw new EvaluationException("Objective derivative evaluation failed.", InnerEvaluation, e); } - this.Checker.DerivativeChecker(this.InnerEvaluation); + Checker.DerivativeChecker(InnerEvaluation); } - return this.InnerEvaluation.Derivative; + return InnerEvaluation.Derivative; } } @@ -75,21 +102,20 @@ namespace MathNet.Numerics.Optimization { get { - - if (!this.SecondDerivativeChecked) + if (!SecondDerivativeChecked) { double tmp; try { - tmp = this.InnerEvaluation.SecondDerivative; + tmp = InnerEvaluation.SecondDerivative; } catch (Exception e) { - throw new EvaluationException("Objective second derivative evaluation failed.", this.InnerEvaluation, e); + throw new EvaluationException("Objective second derivative evaluation failed.", InnerEvaluation, e); } - this.Checker.SecondDerivativeChecker(this.InnerEvaluation); + Checker.SecondDerivativeChecker(InnerEvaluation); } - return this.InnerEvaluation.SecondDerivative; + return InnerEvaluation.SecondDerivative; } } } @@ -101,29 +127,29 @@ namespace MathNet.Numerics.Optimization public Action DerivativeChecker { get; private set; } public Action SecondDerivativeChecker { get; private set; } - public ObjectiveChecker1D(IObjectiveFunction1D objective, Action value_checker, Action gradient_checker, Action hessian_checker) + public ObjectiveChecker1D(IObjectiveFunction1D objective, Action valueChecker, Action gradientChecker, Action hessianChecker) { - this.InnerObjective = objective; - this.ValueChecker = value_checker; - this.DerivativeChecker = gradient_checker; - this.SecondDerivativeChecker = hessian_checker; + InnerObjective = objective; + ValueChecker = valueChecker; + DerivativeChecker = gradientChecker; + SecondDerivativeChecker = hessianChecker; } public bool DerivativeSupported { - get { return this.InnerObjective.DerivativeSupported; } + get { return InnerObjective.DerivativeSupported; } } public bool SecondDerivativeSupported { - get { return this.InnerObjective.SecondDerivativeSupported; } + get { return InnerObjective.SecondDerivativeSupported; } } public IEvaluation1D Evaluate(double point) { try { - return new CheckedEvaluation1D(this, this.InnerObjective.Evaluate(point)); + return new CheckedEvaluation1D(this, InnerObjective.Evaluate(point)); } catch (Exception e) { diff --git a/src/Numerics/Optimization/ObjectiveFunction.cs b/src/Numerics/Optimization/ObjectiveFunction.cs index 6c91c694..7de8c782 100644 --- a/src/Numerics/Optimization/ObjectiveFunction.cs +++ b/src/Numerics/Optimization/ObjectiveFunction.cs @@ -1,15 +1,47 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; - +// +// 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-2013 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.LinearAlgebra.Double; namespace MathNet.Numerics.Optimization { [Flags] - public enum EvaluationStatus { None = 0, Value = 1, Gradient = 2, Hessian = 4 } + public enum EvaluationStatus + { + None = 0, + Value = 1, + Gradient = 2, + Hessian = 4 + } public interface IEvaluation { @@ -34,56 +66,64 @@ namespace MathNet.Numerics.Optimization protected Vector _gradient; protected Matrix _hessian; protected Vector _point; - - public Vector Point { get { return _point; } } + + public Vector Point + { + get { return _point; } + } protected BaseEvaluation() { _status = EvaluationStatus.None; } - public EvaluationStatus Status { get { return _status; } } + public EvaluationStatus Status + { + get { return _status; } + } - public double Value - { - get + public double Value + { + get { if (!_value.HasValue) { - setValue(); + SetValue(); _status |= EvaluationStatus.Value; } - return _value.Value; - } + return _value.Value; + } } - public Vector Gradient - { - get + + public Vector Gradient + { + get { if (_gradient == null) { - setGradient(); + SetGradient(); _status |= EvaluationStatus.Gradient; } - return _gradient; - } + return _gradient; + } } - public Matrix Hessian - { + + public Matrix Hessian + { get { if (_hessian == null) { - setHessian(); + SetHessian(); _status |= EvaluationStatus.Hessian; } - return _hessian; - } + return _hessian; + } } - protected abstract void setValue(); - protected abstract void setGradient(); - protected abstract void setHessian(); + protected abstract void SetValue(); + protected abstract void SetGradient(); + protected abstract void SetHessian(); } @@ -102,17 +142,17 @@ namespace MathNet.Numerics.Optimization _status = EvaluationStatus.None; } - protected override void setValue() + protected override void SetValue() { throw new NotImplementedException(); } - protected override void setGradient() + protected override void SetGradient() { throw new NotImplementedException(); } - protected override void setHessian() + protected override void SetHessian() { throw new NotImplementedException(); } @@ -121,99 +161,100 @@ namespace MathNet.Numerics.Optimization public class OneDEvaluationExpander : IEvaluation { public IEvaluation1D InnerEval { get; private set; } + public OneDEvaluationExpander(IEvaluation1D eval) { - this.InnerEval = eval; + InnerEval = eval; } public Vector Point { - get { return DenseVector.Create(1,this.InnerEval.Point); } + get { return DenseVector.Create(1, InnerEval.Point); } } public EvaluationStatus Status { - get { return this.InnerEval.Status; } + get { return InnerEval.Status; } } public double Value { - get { return this.InnerEval.Value; } + get { return InnerEval.Value; } } public Vector Gradient { - get { return DenseVector.Create(1,this.InnerEval.Derivative) ; } + get { return DenseVector.Create(1, InnerEval.Derivative); } } public Matrix Hessian { - get { return DenseMatrix.Create(1,1,this.InnerEval.SecondDerivative); } + get { return DenseMatrix.Create(1, 1, InnerEval.SecondDerivative); } } } public class CachedEvaluation : BaseEvaluation { - private SimpleObjectiveFunction _objective_object; + readonly SimpleObjectiveFunction _objectiveObject; public CachedEvaluation(SimpleObjectiveFunction f, Vector point) { - _objective_object = f; - _point = point; + _objectiveObject = f; + _point = point; } - protected override void setValue() + protected override void SetValue() { - _value = _objective_object.Objective(_point); + _value = _objectiveObject.Objective(_point); } - protected override void setGradient() + protected override void SetGradient() { - _gradient = _objective_object.Gradient(_point); + _gradient = _objectiveObject.Gradient(_point); } - protected override void setHessian() + protected override void SetHessian() { - _hessian = _objective_object.Hessian(_point); + _hessian = _objectiveObject.Hessian(_point); } } public class SimpleObjectiveFunction : IObjectiveFunction { - public Func,double> Objective { get; private set; } - public Func,Vector> Gradient { get; private set; } + public Func, double> Objective { get; private set; } + public Func, Vector> Gradient { get; private set; } public Func, Matrix> Hessian { get; private set; } - public SimpleObjectiveFunction(Func,double> objective) + public SimpleObjectiveFunction(Func, double> objective) { - this.Objective = objective; - this.Gradient = null; - this.Hessian = null; + Objective = objective; + Gradient = null; + Hessian = null; } - public SimpleObjectiveFunction(Func,double> objective, Func,Vector> gradient) + public SimpleObjectiveFunction(Func, double> objective, Func, Vector> gradient) { - this.Objective = objective; - this.Gradient = gradient; - this.Hessian = null; + Objective = objective; + Gradient = gradient; + Hessian = null; } - public SimpleObjectiveFunction(Func,double> objective, Func,Vector> gradient, Func, Matrix> hessian) + public SimpleObjectiveFunction(Func, double> objective, Func, Vector> gradient, Func, Matrix> hessian) { - this.Objective = objective; - this.Gradient = gradient; - this.Hessian = hessian; + Objective = objective; + Gradient = gradient; + Hessian = hessian; } public bool GradientSupported { - get { return this.Gradient != null; } + get { return Gradient != null; } } public bool HessianSupported { - get { return this.Hessian != null; } + get { return Hessian != null; } } public IEvaluation Evaluate(Vector point) diff --git a/src/Numerics/Optimization/ObjectiveFunction1D.cs b/src/Numerics/Optimization/ObjectiveFunction1D.cs index 339bf1b9..99a82b37 100644 --- a/src/Numerics/Optimization/ObjectiveFunction1D.cs +++ b/src/Numerics/Optimization/ObjectiveFunction1D.cs @@ -1,7 +1,34 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; +// +// 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-2013 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 { @@ -27,16 +54,22 @@ namespace MathNet.Numerics.Optimization protected EvaluationStatus _status; protected double? _value; protected double? _derivative; - protected double? _second_derivative; + protected double? _secondDerivative; - public double Point { get { return _point; } } + public double Point + { + get { return _point; } + } protected BaseEvaluation1D() { _status = EvaluationStatus.None; } - public EvaluationStatus Status { get { return _status; } } + public EvaluationStatus Status + { + get { return _status; } + } public double Value { @@ -44,62 +77,67 @@ namespace MathNet.Numerics.Optimization { if (!_value.HasValue) { - setValue(); + SetValue(); _status |= EvaluationStatus.Value; } return _value.Value; } } + public double Derivative { get { if (_derivative == null) { - setDerivative(); + SetDerivative(); _status |= EvaluationStatus.Gradient; } return _derivative.Value; } } + public double SecondDerivative { get { - if (_second_derivative == null) + if (_secondDerivative == null) { - setSecondDerivative(); + SetSecondDerivative(); _status |= EvaluationStatus.Hessian; } - return _second_derivative.Value; + return _secondDerivative.Value; } } - protected abstract void setValue(); - protected abstract void setDerivative(); - protected abstract void setSecondDerivative(); + protected abstract void SetValue(); + protected abstract void SetDerivative(); + protected abstract void SetSecondDerivative(); } public class CachedEvaluation1D : BaseEvaluation1D { - private SimpleObjectiveFunction1D _objective_object; + readonly SimpleObjectiveFunction1D _objectiveObject; public CachedEvaluation1D(SimpleObjectiveFunction1D f, double point) { - _objective_object = f; + _objectiveObject = f; _point = point; } - protected override void setValue() + + protected override void SetValue() { - _value = _objective_object.Objective(_point); + _value = _objectiveObject.Objective(_point); } - protected override void setDerivative() + + protected override void SetDerivative() { - _derivative = _objective_object.Derivative(_point); + _derivative = _objectiveObject.Derivative(_point); } - protected override void setSecondDerivative() + + protected override void SetSecondDerivative() { - _second_derivative = _objective_object.SecondDerivative(_point); + _secondDerivative = _objectiveObject.SecondDerivative(_point); } } @@ -111,33 +149,33 @@ namespace MathNet.Numerics.Optimization public SimpleObjectiveFunction1D(Func objective) { - this.Objective = objective; - this.Derivative = null; - this.SecondDerivative = null; + Objective = objective; + Derivative = null; + SecondDerivative = null; } public SimpleObjectiveFunction1D(Func objective, Func derivative) { - this.Objective = objective; - this.Derivative = derivative; - this.SecondDerivative = null; + Objective = objective; + Derivative = derivative; + SecondDerivative = null; } - public SimpleObjectiveFunction1D(Func objective, Func derivative, Func second_derivative) + public SimpleObjectiveFunction1D(Func objective, Func derivative, Func secondDerivative) { - this.Objective = objective; - this.Derivative = derivative; - this.SecondDerivative = second_derivative; + Objective = objective; + Derivative = derivative; + SecondDerivative = secondDerivative; } public bool DerivativeSupported { - get { return this.Derivative != null; } + get { return Derivative != null; } } public bool SecondDerivativeSupported { - get { return this.SecondDerivative != null; } + get { return SecondDerivative != null; } } public IEvaluation1D Evaluate(double point) diff --git a/src/Numerics/Optimization/OptimizationResult.cs b/src/Numerics/Optimization/OptimizationResult.cs deleted file mode 100644 index 44017dec..00000000 --- a/src/Numerics/Optimization/OptimizationResult.cs +++ /dev/null @@ -1,8 +0,0 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; - -namespace MathNet.Numerics.Optimization -{ -} diff --git a/src/Numerics/Optimization/StrongWolfeLineSearch.cs b/src/Numerics/Optimization/StrongWolfeLineSearch.cs deleted file mode 100644 index 236c63c9..00000000 --- a/src/Numerics/Optimization/StrongWolfeLineSearch.cs +++ /dev/null @@ -1,11 +0,0 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; - -namespace MathNet.Numerics.Optimization -{ - class StrongWolfeLineSearch - { - } -} diff --git a/src/Numerics/Optimization/WeakWolfeLineSearch.cs b/src/Numerics/Optimization/WeakWolfeLineSearch.cs index f7481154..1bbcdf9c 100644 --- a/src/Numerics/Optimization/WeakWolfeLineSearch.cs +++ b/src/Numerics/Optimization/WeakWolfeLineSearch.cs @@ -1,8 +1,34 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; - +// +// 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-2013 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 @@ -11,107 +37,105 @@ namespace MathNet.Numerics.Optimization { public double C1 { get; set; } public double C2 { get; set; } - public double ParameterTolerance { get; set; } - public int MaximumIterations { get; set; } + public double ParameterTolerance { get; set; } + public int MaximumIterations { get; set; } - public WeakWolfeLineSearch(double c1, double c2, double parameter_tolerance, int max_iterations=10) + public WeakWolfeLineSearch(double c1, double c2, double parameterTolerance, int maxIterations = 10) { - this.C1 = c1; - this.C2 = c2; - this.ParameterTolerance = parameter_tolerance; - this.MaximumIterations = max_iterations; + C1 = c1; + C2 = c2; + ParameterTolerance = parameterTolerance; + MaximumIterations = maxIterations; } // Implemented following http://www.math.washington.edu/~burke/crs/408/lectures/L9-weak-Wolfe.pdf - public LineSearchOutput FindConformingStep(IObjectiveFunction objective, IEvaluation starting_point, Vector search_direction, double initial_step) + public LineSearchOutput FindConformingStep(IObjectiveFunction objective, IEvaluation startingPoint, Vector searchDirection, double initialStep) { if (!(objective is ObjectiveChecker)) - objective = new ObjectiveChecker(objective, this.ValidateValue, this.ValidateGradient, null); + objective = new ObjectiveChecker(objective, ValidateValue, ValidateGradient, null); + + double lowerBound = 0.0; + double upperBound = Double.PositiveInfinity; + double step = initialStep; - double lower_bound = 0.0; - double upper_bound = Double.PositiveInfinity; - double step = initial_step; + double initialValue = startingPoint.Value; + Vector initialGradient = startingPoint.Gradient; - double initial_value = starting_point.Value; - Vector initial_gradient = starting_point.Gradient; - - double initial_dd = search_direction*initial_gradient; + double initialDd = searchDirection*initialGradient; int ii; - IEvaluation candidate_eval = null; - ExitCondition reason_for_exit = ExitCondition.None; - for (ii = 0; ii < this.MaximumIterations; ++ii) + IEvaluation candidateEval = null; + var reasonForExit = ExitCondition.None; + for (ii = 0; ii < MaximumIterations; ++ii) { - candidate_eval = objective.Evaluate(starting_point.Point + search_direction * step); + candidateEval = objective.Evaluate(startingPoint.Point + searchDirection*step); - double step_dd = search_direction * candidate_eval.Gradient; + double stepDd = searchDirection*candidateEval.Gradient; - if (candidate_eval.Value > initial_value + this.C1 * step * initial_dd) + if (candidateEval.Value > initialValue + C1*step*initialDd) { - upper_bound = step; - step = 0.5 * (lower_bound + upper_bound); - } - else if (step_dd < this.C2*initial_dd) + upperBound = step; + step = 0.5*(lowerBound + upperBound); + } + else if (stepDd < C2*initialDd) { - lower_bound = step; - step = Double.IsPositiveInfinity(upper_bound) ? 2 * lower_bound : 0.5 * (lower_bound + upper_bound); + lowerBound = step; + step = Double.IsPositiveInfinity(upperBound) ? 2*lowerBound : 0.5*(lowerBound + upperBound); } - else + else { - reason_for_exit = ExitCondition.WeakWolfeCriteria; + reasonForExit = ExitCondition.WeakWolfeCriteria; break; } - if (!Double.IsInfinity(upper_bound)) - { - double max_rel_change = 0.0; - for (int jj = 0; jj < candidate_eval.Point.Count; ++jj) - { - double tmp = Math.Abs (search_direction[jj]*(upper_bound - lower_bound)) / Math.Max(Math.Abs(candidate_eval.Point[jj]),1.0); - max_rel_change = Math.Max(max_rel_change, tmp); - } - if (max_rel_change < this.ParameterTolerance) - { - reason_for_exit = ExitCondition.LackOfProgress; - 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 = ExitCondition.LackOfProgress; + break; + } + } } - - if (ii == this.MaximumIterations && Double.IsPositiveInfinity(upper_bound)) - throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached. Function appears to be unbounded in search direction.",this.MaximumIterations)); - else if (ii == this.MaximumIterations) - throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.",this.MaximumIterations)); - else - return new LineSearchOutput(candidate_eval, ii, step, reason_for_exit); + + 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 LineSearchOutput(candidateEval, ii, step, reasonForExit); } - private bool Conforms(IEvaluation starting_point, Vector search_direction, double step, IEvaluation ending_point) + bool Conforms(IEvaluation startingPoint, Vector searchDirection, double step, IEvaluation endingPoint) { + bool sufficientDecrease = endingPoint.Value <= startingPoint.Value + C1*step*(startingPoint.Gradient*searchDirection); + bool notTooSteep = endingPoint.Gradient*searchDirection >= C2*startingPoint.Gradient*searchDirection; - bool sufficient_decrease = ending_point.Value <= starting_point.Value + this.C1 * step * (starting_point.Gradient * search_direction); - bool not_too_steep = ending_point.Gradient * search_direction >= this.C2 * starting_point.Gradient * search_direction; - - return step > 0 && sufficient_decrease && not_too_steep; + return step > 0 && sufficientDecrease && notTooSteep; } - private void ValidateValue(IEvaluation eval) + void ValidateValue(IEvaluation eval) { - if (!this.IsFinite(eval.Value)) - throw new EvaluationException(String.Format("Non-finite value returned by objective function: {0}", eval.Value),eval); + if (!IsFinite(eval.Value)) + throw new EvaluationException(String.Format("Non-finite value returned by objective function: {0}", eval.Value), eval); } - private void ValidateGradient(IEvaluation eval) + void ValidateGradient(IEvaluation eval) { foreach (double x in eval.Gradient) - if (!this.IsFinite(x)) + if (!IsFinite(x)) { - throw new EvaluationException(String.Format("Non-finite value returned by gradient: {0}", x),eval); + throw new EvaluationException(String.Format("Non-finite value returned by gradient: {0}", x), eval); } } - private bool IsFinite(double x) + bool IsFinite(double x) { return !(Double.IsNaN(x) || Double.IsInfinity(x)); } diff --git a/src/UnitTests/OptimizationTests/RosenbrockFunction.cs b/src/UnitTests/OptimizationTests/RosenbrockFunction.cs index 7dea0f51..4ca8a02d 100644 --- a/src/UnitTests/OptimizationTests/RosenbrockFunction.cs +++ b/src/UnitTests/OptimizationTests/RosenbrockFunction.cs @@ -1,8 +1,34 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; +// +// 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-2013 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.UnitTests.OptimizationTests @@ -11,24 +37,24 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests { public static double Value(Vector input) { - return Math.Pow((1 - input[0]), 2) + 100 * Math.Pow((input[1] - input[0] * input[0]), 2); + return Math.Pow((1 - input[0]), 2) + 100*Math.Pow((input[1] - input[0]*input[0]), 2); } public static Vector Gradient(Vector input) { - Vector output = new MathNet.Numerics.LinearAlgebra.Double.DenseVector(2); - output[0] = -2 * (1 - input[0]) + 200 * (input[1] - input[0] * input[0]) * (-2 * input[0]); - output[1] = 2 * 100 * (input[1] - input[0] * input[0]); + Vector output = new LinearAlgebra.Double.DenseVector(2); + output[0] = -2*(1 - input[0]) + 200*(input[1] - input[0]*input[0])*(-2*input[0]); + output[1] = 2*100*(input[1] - input[0]*input[0]); return output; } public static Matrix Hessian(Vector input) { - Matrix output = new MathNet.Numerics.LinearAlgebra.Double.DenseMatrix(2,2); - output[0, 0] = 2 - 400 * input[1] + 1200 * input[0] * input[0]; + Matrix output = new LinearAlgebra.Double.DenseMatrix(2, 2); + output[0, 0] = 2 - 400*input[1] + 1200*input[0]*input[0]; output[1, 1] = 200; - output[0, 1] = -400 * input[0]; + output[0, 1] = -400*input[0]; output[1, 0] = output[0, 1]; return output; } diff --git a/src/UnitTests/OptimizationTests/TestBfgsMinimizer.cs b/src/UnitTests/OptimizationTests/TestBfgsMinimizer.cs index 742f7b59..21a02cb5 100644 --- a/src/UnitTests/OptimizationTests/TestBfgsMinimizer.cs +++ b/src/UnitTests/OptimizationTests/TestBfgsMinimizer.cs @@ -1,10 +1,35 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; +// +// 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-2013 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 NUnit.Framework; - using MathNet.Numerics.Optimization; namespace MathNet.Numerics.UnitTests.OptimizationTests @@ -12,13 +37,12 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests [TestFixture] public class TestBfgsMinimizer { - [Test] public void FindMinimum_Rosenbrock_Easy() { var obj = new SimpleObjectiveFunction(RosenbrockFunction.Value, RosenbrockFunction.Gradient); var solver = new BfgsMinimizer(1e-5, 1000); - var result = solver.FindMinimum(obj, new MathNet.Numerics.LinearAlgebra.Double.DenseVector(new double[] { 1.2, 1.2 })); + var result = solver.FindMinimum(obj, new LinearAlgebra.Double.DenseVector(new[] { 1.2, 1.2 })); 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)); @@ -29,22 +53,21 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests { var obj = new SimpleObjectiveFunction(RosenbrockFunction.Value, RosenbrockFunction.Gradient); var solver = new BfgsMinimizer(1e-5, 1000); - var result = solver.FindMinimum(obj, new MathNet.Numerics.LinearAlgebra.Double.DenseVector(new double[] { -1.2, 1.0 })); + var result = solver.FindMinimum(obj, new LinearAlgebra.Double.DenseVector(new[] { -1.2, 1.0 })); 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 = new SimpleObjectiveFunction(RosenbrockFunction.Value, RosenbrockFunction.Gradient); var solver = new BfgsMinimizer(1e-5, 1000); - var result = solver.FindMinimum(obj, new MathNet.Numerics.LinearAlgebra.Double.DenseVector(new double[] { -0.9, -0.5 })); + var result = solver.FindMinimum(obj, new LinearAlgebra.Double.DenseVector(new[] { -0.9, -0.5 })); 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/TestBisectionRootFinder.cs b/src/UnitTests/OptimizationTests/TestBisectionRootFinder.cs index 561e9965..4d1149db 100644 --- a/src/UnitTests/OptimizationTests/TestBisectionRootFinder.cs +++ b/src/UnitTests/OptimizationTests/TestBisectionRootFinder.cs @@ -1,26 +1,52 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; -using NUnit.Framework; +// +// 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-2013 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 NUnit.Framework; using MathNet.Numerics.Optimization; namespace MathNet.Numerics.UnitTests.OptimizationTests { [TestFixture] - class TestBisectionRootFinder + internal class TestBisectionRootFinder { [Test] public void FindRoot_Works() { var algorithm = new BisectionRootFinder(0.001, 0.001); - var f1 = new Func((x) => (x - 3) * (x - 4)); + var f1 = new Func(x => (x - 3)*(x - 4)); var r1 = algorithm.FindRoot(f1, 2.1, 3.9); Assert.That(Math.Abs(f1(r1)), Is.LessThan(0.001)); Assert.That(Math.Abs(r1 - 3.0), Is.LessThan(0.001)); - var f2 = new Func((x) => (x - 3) * (x - 4)); + var f2 = new Func(x => (x - 3)*(x - 4)); var r2 = algorithm.FindRoot(f1, 2.1, 3.4); Assert.That(Math.Abs(f2(r2)), Is.LessThan(0.001)); Assert.That(Math.Abs(r2 - 3.0), Is.LessThan(0.001)); diff --git a/src/UnitTests/OptimizationTests/TestConjugateGradientMinimizer.cs b/src/UnitTests/OptimizationTests/TestConjugateGradientMinimizer.cs index 79b12b45..12aa6283 100644 --- a/src/UnitTests/OptimizationTests/TestConjugateGradientMinimizer.cs +++ b/src/UnitTests/OptimizationTests/TestConjugateGradientMinimizer.cs @@ -1,10 +1,35 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; +// +// 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-2013 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 NUnit.Framework; - using MathNet.Numerics.Optimization; namespace MathNet.Numerics.UnitTests.OptimizationTests @@ -12,15 +37,14 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests [TestFixture] public class TestConjugateGradientMinimizer { - [Test] public void FindMinimum_Rosenbrock_Easy() { var obj = new SimpleObjectiveFunction(RosenbrockFunction.Value, RosenbrockFunction.Gradient); var solver = new ConjugateGradientMinimizer(1e-5, 1000); - var result = solver.FindMinimum(obj, new MathNet.Numerics.LinearAlgebra.Double.DenseVector(new double[]{1.2,1.2})); - - Assert.That(Math.Abs(result.MinimizingPoint[0]-1.0), Is.LessThan(1e-3)); + var result = solver.FindMinimum(obj, new LinearAlgebra.Double.DenseVector(new[] { 1.2, 1.2 })); + + 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)); } @@ -29,7 +53,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests { var obj = new SimpleObjectiveFunction(RosenbrockFunction.Value, RosenbrockFunction.Gradient); var solver = new ConjugateGradientMinimizer(1e-5, 1000); - var result = solver.FindMinimum(obj, new MathNet.Numerics.LinearAlgebra.Double.DenseVector(new double[] { -1.2, 1.0 })); + var result = solver.FindMinimum(obj, new LinearAlgebra.Double.DenseVector(new[] { -1.2, 1.0 })); 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/TestNewtonMinimizer.cs b/src/UnitTests/OptimizationTests/TestNewtonMinimizer.cs index e97ffd86..a5645b4d 100644 --- a/src/UnitTests/OptimizationTests/TestNewtonMinimizer.cs +++ b/src/UnitTests/OptimizationTests/TestNewtonMinimizer.cs @@ -1,10 +1,35 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; +// +// 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-2013 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 NUnit.Framework; - using MathNet.Numerics.Optimization; namespace MathNet.Numerics.UnitTests.OptimizationTests @@ -12,13 +37,12 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests [TestFixture] public class TestNewtonMinimizer { - [Test] public void FindMinimum_Rosenbrock_Easy() { var obj = new SimpleObjectiveFunction(RosenbrockFunction.Value, RosenbrockFunction.Gradient, RosenbrockFunction.Hessian); var solver = new NewtonMinimizer(1e-5, 1000); - var result = solver.FindMinimum(obj, new MathNet.Numerics.LinearAlgebra.Double.DenseVector(new double[] { 1.2, 1.2 })); + var result = solver.FindMinimum(obj, new LinearAlgebra.Double.DenseVector(new[] { 1.2, 1.2 })); 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)); @@ -29,7 +53,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests { var obj = new SimpleObjectiveFunction(RosenbrockFunction.Value, RosenbrockFunction.Gradient, RosenbrockFunction.Hessian); var solver = new NewtonMinimizer(1e-5, 1000); - var result = solver.FindMinimum(obj, new MathNet.Numerics.LinearAlgebra.Double.DenseVector(new double[] { -1.2, 1.0 })); + var result = solver.FindMinimum(obj, new LinearAlgebra.Double.DenseVector(new[] { -1.2, 1.0 })); 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)); @@ -40,7 +64,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests { var obj = new SimpleObjectiveFunction(RosenbrockFunction.Value, RosenbrockFunction.Gradient, RosenbrockFunction.Hessian); var solver = new NewtonMinimizer(1e-5, 1000); - var result = solver.FindMinimum(obj, new MathNet.Numerics.LinearAlgebra.Double.DenseVector(new double[] { -0.9, -0.5 })); + var result = solver.FindMinimum(obj, new LinearAlgebra.Double.DenseVector(new[] { -0.9, -0.5 })); 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)); @@ -51,7 +75,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests { var obj = new SimpleObjectiveFunction(RosenbrockFunction.Value, RosenbrockFunction.Gradient, RosenbrockFunction.Hessian); var solver = new NewtonMinimizer(1e-5, 1000, true); - var result = solver.FindMinimum(obj, new MathNet.Numerics.LinearAlgebra.Double.DenseVector(new double[] { 1.2, 1.2 })); + var result = solver.FindMinimum(obj, new LinearAlgebra.Double.DenseVector(new[] { 1.2, 1.2 })); 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)); @@ -62,7 +86,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests { var obj = new SimpleObjectiveFunction(RosenbrockFunction.Value, RosenbrockFunction.Gradient, RosenbrockFunction.Hessian); var solver = new NewtonMinimizer(1e-5, 1000, true); - var result = solver.FindMinimum(obj, new MathNet.Numerics.LinearAlgebra.Double.DenseVector(new double[] { -1.2, 1.0 })); + var result = solver.FindMinimum(obj, new LinearAlgebra.Double.DenseVector(new[] { -1.2, 1.0 })); 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)); @@ -73,7 +97,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests { var obj = new SimpleObjectiveFunction(RosenbrockFunction.Value, RosenbrockFunction.Gradient, RosenbrockFunction.Hessian); var solver = new NewtonMinimizer(1e-5, 1000, true); - var result = solver.FindMinimum(obj, new MathNet.Numerics.LinearAlgebra.Double.DenseVector(new double[] { -0.9, -0.5 })); + var result = solver.FindMinimum(obj, new LinearAlgebra.Double.DenseVector(new[] { -0.9, -0.5 })); 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/TestRosenbrockFunction.cs b/src/UnitTests/OptimizationTests/TestRosenbrockFunction.cs index 28098b39..fe964e1e 100644 --- a/src/UnitTests/OptimizationTests/TestRosenbrockFunction.cs +++ b/src/UnitTests/OptimizationTests/TestRosenbrockFunction.cs @@ -1,27 +1,55 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; +// +// 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-2013 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 NUnit.Framework; + namespace MathNet.Numerics.UnitTests.OptimizationTests { [TestFixture] - class TestRosenbrockFunction + internal class TestRosenbrockFunction { [Test] public void TestGradient() { - var input = new LinearAlgebra.Double.DenseVector(new double[]{ -0.9, -0.5 } ); + var input = new LinearAlgebra.Double.DenseVector(new[] { -0.9, -0.5 }); var v1 = RosenbrockFunction.Value(input); var g = RosenbrockFunction.Gradient(input); - var eps = 1e-5; - var eps0 = (new LinearAlgebra.Double.DenseVector(new double[] { 1.0, 0.0 })) * eps; - var eps1 = (new LinearAlgebra.Double.DenseVector(new double[] { 0.0, 1.0 })) * eps; + const double eps = 1e-5; + var eps0 = (new LinearAlgebra.Double.DenseVector(new[] { 1.0, 0.0 }))*eps; + var eps1 = (new LinearAlgebra.Double.DenseVector(new[] { 0.0, 1.0 }))*eps; - var g0 = (RosenbrockFunction.Value(input + eps0) - RosenbrockFunction.Value(input - eps0)) / (2 * eps); - var g1 = (RosenbrockFunction.Value(input + eps1) - RosenbrockFunction.Value(input - eps1)) / (2 * eps); + var g0 = (RosenbrockFunction.Value(input + eps0) - RosenbrockFunction.Value(input - eps0))/(2*eps); + var g1 = (RosenbrockFunction.Value(input + eps1) - RosenbrockFunction.Value(input - eps1))/(2*eps); Assert.That(Math.Abs(g0 - g[0]) < 1e-3); Assert.That(Math.Abs(g1 - g[1]) < 1e-3); @@ -30,28 +58,27 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests [Test] public void TestHessian() { - var input = new LinearAlgebra.Double.DenseVector(new double[] { -0.9, -0.5 }); + var input = new LinearAlgebra.Double.DenseVector(new[] { -0.9, -0.5 }); var v1 = RosenbrockFunction.Value(input); var h = RosenbrockFunction.Hessian(input); - var eps = 1e-5; + const double eps = 1e-5; + + var eps0 = (new LinearAlgebra.Double.DenseVector(new[] { 1.0, 0.0 }))*eps; + var eps1 = (new LinearAlgebra.Double.DenseVector(new[] { 0.0, 1.0 }))*eps; - var eps0 = (new LinearAlgebra.Double.DenseVector(new double[] { 1.0, 0.0 })) * eps; - var eps1 = (new LinearAlgebra.Double.DenseVector(new double[] { 0.0, 1.0 })) * eps; + var epsuu = (new LinearAlgebra.Double.DenseVector(new[] { 1.0, 1.0 }))*eps; + var epsud = (new LinearAlgebra.Double.DenseVector(new[] { 1.0, -1.0 }))*eps; - var epsuu = (new LinearAlgebra.Double.DenseVector(new double[] { 1.0, 1.0 })) * eps; - var epsud = (new LinearAlgebra.Double.DenseVector(new double[] { 1.0, -1.0 })) * eps; - - - var h00 = (RosenbrockFunction.Value(input + eps0) - 2*RosenbrockFunction.Value(input) + RosenbrockFunction.Value(input - eps0)) / (eps*eps); - var h11 = (RosenbrockFunction.Value(input + eps1) - 2 * RosenbrockFunction.Value(input) + RosenbrockFunction.Value(input - eps1)) / (eps * eps); - var h01 = (RosenbrockFunction.Value(input + epsuu) - RosenbrockFunction.Value(input + epsud) - RosenbrockFunction.Value(input - epsud) + RosenbrockFunction.Value(input - epsuu)) / (4*eps * eps); + var h00 = (RosenbrockFunction.Value(input + eps0) - 2*RosenbrockFunction.Value(input) + RosenbrockFunction.Value(input - eps0))/(eps*eps); + var h11 = (RosenbrockFunction.Value(input + eps1) - 2*RosenbrockFunction.Value(input) + RosenbrockFunction.Value(input - eps1))/(eps*eps); + var h01 = (RosenbrockFunction.Value(input + epsuu) - RosenbrockFunction.Value(input + epsud) - RosenbrockFunction.Value(input - epsud) + RosenbrockFunction.Value(input - epsuu))/(4*eps*eps); - Assert.That(Math.Abs(h00 - h[0,0]) < 1e-3); - Assert.That(Math.Abs(h11 - h[1,1]) < 1e-3); + Assert.That(Math.Abs(h00 - h[0, 0]) < 1e-3); + Assert.That(Math.Abs(h11 - h[1, 1]) < 1e-3); Assert.That(Math.Abs(h01 - h[0, 1]) < 1e-3); Assert.That(Math.Abs(h01 - h[1, 0]) < 1e-3); }