diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 6b37bd2e..a478aab5 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -115,6 +115,19 @@ + + + + + + + + + + + + + @@ -487,5 +500,6 @@ Resources.Designer.cs + \ No newline at end of file diff --git a/src/Numerics/Optimization/BaseEvaluation.cs b/src/Numerics/Optimization/BaseEvaluation.cs new file mode 100644 index 00000000..fa19f5e8 --- /dev/null +++ b/src/Numerics/Optimization/BaseEvaluation.cs @@ -0,0 +1,70 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization +{ + public abstract class BaseEvaluation : IEvaluation + { + public EvaluationStatus Status { get; set; } + public Vector Point { get; set; } + public double ValueRaw { get; set; } + public Vector GradientRaw { get; set; } + public Matrix HessianRaw { get; set; } + + protected BaseEvaluation() + { + Status = EvaluationStatus.None; + } + + public double Value + { + get + { + if (!Status.HasFlag(EvaluationStatus.Value)) + { + setValue(); + Status |= EvaluationStatus.Value; + } + return ValueRaw; + } + } + public Vector Gradient + { + get + { + if (!Status.HasFlag(EvaluationStatus.Gradient)) + { + setGradient(); + Status |= EvaluationStatus.Gradient; + } + return GradientRaw; + } + } + public Matrix Hessian + { + get + { + if (!Status.HasFlag(EvaluationStatus.Hessian)) + { + setHessian(); + Status |= EvaluationStatus.Hessian; + } + return HessianRaw; + } + } + + public void Reset(Vector new_point) + { + this.Point = new_point; + this.Status = EvaluationStatus.None; + } + + protected abstract void setValue(); + protected abstract void setGradient(); + protected abstract void setHessian(); + } +} diff --git a/src/Numerics/Optimization/BaseObjectiveFunction.cs b/src/Numerics/Optimization/BaseObjectiveFunction.cs new file mode 100644 index 00000000..5cd5f68f --- /dev/null +++ b/src/Numerics/Optimization/BaseObjectiveFunction.cs @@ -0,0 +1,39 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +namespace MathNet.Numerics.Optimization +{ + public class BaseObjectiveFunction : IObjectiveFunction where T : IEvaluation + { + public BaseObjectiveFunction(bool gradient_supported, bool hessian_supported) + { + _gradient_supported = gradient_supported; + _hessian_supported = hessian_supported; + } + + private bool _gradient_supported; + private bool _hessian_supported; + + public bool GradientSupported + { + get { return _gradient_supported; } + } + + public bool HessianSupported + { + get { return _hessian_supported; } + } + + public void Evaluate(LinearAlgebra.Vector point, IEvaluation output) + { + output.Reset(point); + } + + public virtual IEvaluation CreateEvaluationObject() + { + return default(T); + } + } +} diff --git a/src/Numerics/Optimization/Exceptions.cs b/src/Numerics/Optimization/Exceptions.cs new file mode 100644 index 00000000..7999b093 --- /dev/null +++ b/src/Numerics/Optimization/Exceptions.cs @@ -0,0 +1,67 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +namespace MathNet.Numerics.Optimization +{ + public class OptimizationException : Exception + { + public OptimizationException(string message) + : base(message) { } + + public OptimizationException(string message, Exception inner_exception) + : base(message, inner_exception) { } + } + + public class MaximumIterationsException : OptimizationException + { + public MaximumIterationsException(string message) + : base(message) { } + } + + public class EvaluationException : OptimizationException + { + public IEvaluation Evaluation { get; private set; } + + public EvaluationException(string message, IEvaluation eval) + : base(message) + { + this.Evaluation = eval; + } + + public EvaluationException(string message, IEvaluation eval, Exception inner_exception) + : base(message, inner_exception) + { + this.Evaluation = eval; + } + + //public EvaluationException(string message, IEvaluation1D eval) + // : base(message) + //{ + // this.Evaluation = new OneDEvaluationExpander(eval); + //} + + //public EvaluationException(string message, IEvaluation1D eval, Exception inner_exception) + // : base(message, inner_exception) + //{ + // this.Evaluation = new OneDEvaluationExpander(eval); + //} + + } + + public class InnerOptimizationException : OptimizationException + { + public InnerOptimizationException(string message) + : base(message) { } + + public InnerOptimizationException(string message, Exception inner_exception) + : base(message, inner_exception) { } + } + + public class IncompatibleObjectiveException : OptimizationException + { + public IncompatibleObjectiveException(string message) + : base(message) { } + } +} diff --git a/src/Numerics/Optimization/IEvaluation.cs b/src/Numerics/Optimization/IEvaluation.cs new file mode 100644 index 00000000..9150d5fd --- /dev/null +++ b/src/Numerics/Optimization/IEvaluation.cs @@ -0,0 +1,28 @@ +using MathNet.Numerics.LinearAlgebra; +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +namespace MathNet.Numerics.Optimization +{ + [Flags] + public enum EvaluationStatus { None = 0, Value = 1, Gradient = 2, Hessian = 4 } + + public interface IEvaluation + { + Vector Point { get; set; } + EvaluationStatus Status { get; set; } + + // Used by algorithm + double Value { get; } + Vector Gradient { get; } + Matrix Hessian { get; } + + // Used by ObjectiveFunction + void Reset(Vector new_point); + double ValueRaw { get; set; } + Vector GradientRaw { get; set; } + Matrix HessianRaw { get; set; } + } +} diff --git a/src/Numerics/Optimization/IObjectiveFunction.cs b/src/Numerics/Optimization/IObjectiveFunction.cs new file mode 100644 index 00000000..cd553ac2 --- /dev/null +++ b/src/Numerics/Optimization/IObjectiveFunction.cs @@ -0,0 +1,18 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization +{ + public interface IObjectiveFunction + { + bool GradientSupported { get; } + bool HessianSupported { get; } + + IEvaluation CreateEvaluationObject(); + void Evaluate(Vector point, IEvaluation output); + } +} diff --git a/src/Numerics/Optimization/IUnconstrainedMinimizer.cs b/src/Numerics/Optimization/IUnconstrainedMinimizer.cs new file mode 100644 index 00000000..d7200b01 --- /dev/null +++ b/src/Numerics/Optimization/IUnconstrainedMinimizer.cs @@ -0,0 +1,14 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization +{ + public interface IUnconstrainedMinimizer + { + MinimizationOutput FindMinimum(IObjectiveFunction objective, Vector initial_guess); + } + +} diff --git a/src/Numerics/Optimization/Implementation/LineSearchOutput.cs b/src/Numerics/Optimization/Implementation/LineSearchOutput.cs new file mode 100644 index 00000000..01a28276 --- /dev/null +++ b/src/Numerics/Optimization/Implementation/LineSearchOutput.cs @@ -0,0 +1,18 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +namespace MathNet.Numerics.Optimization.Implementation +{ + public class LineSearchOutput : MinimizationOutput + { + 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) + { + this.FinalStep = final_step; + } + } +} diff --git a/src/Numerics/Optimization/Implementation/NullEvaluation.cs b/src/Numerics/Optimization/Implementation/NullEvaluation.cs new file mode 100644 index 00000000..fcc31c4b --- /dev/null +++ b/src/Numerics/Optimization/Implementation/NullEvaluation.cs @@ -0,0 +1,32 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization.Implementation +{ + public class NullEvaluation : BaseEvaluation + { + public NullEvaluation(Vector point) + : base() + { + this.Point = point; + } + protected override void setValue() + { + throw new NotImplementedException(); + } + + protected override void setGradient() + { + throw new NotImplementedException(); + } + + protected override void setHessian() + { + throw new NotImplementedException(); + } + } +} diff --git a/src/Numerics/Optimization/Implementation/ObjectiveChecker.cs b/src/Numerics/Optimization/Implementation/ObjectiveChecker.cs new file mode 100644 index 00000000..3aa0d052 --- /dev/null +++ b/src/Numerics/Optimization/Implementation/ObjectiveChecker.cs @@ -0,0 +1,183 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization.Implementation +{ + public class CheckedEvaluation : IEvaluation + { + private ObjectiveChecker Checker; + public IEvaluation InnerEvaluation { get; private set; } + private bool ValueChecked; + private bool GradientChecked; + private bool HessianChecked; + + public CheckedEvaluation(ObjectiveChecker checker, IEvaluation evaluation) + { + this.Checker = checker; + this.InnerEvaluation = evaluation; + } + + public Vector Point + { + get { return this.InnerEvaluation.Point; } + set { this.InnerEvaluation.Point = value; } + } + + public EvaluationStatus Status + { + get + { + return this.InnerEvaluation.Status; + } + set + { + this.InnerEvaluation.Status = value; + } + } + + public double ValueRaw + { + get + { + return this.InnerEvaluation.Value; + } + set + { + this.InnerEvaluation.ValueRaw = value; + } + } + + public Vector GradientRaw + { + get { return this.InnerEvaluation.GradientRaw; } + set { this.InnerEvaluation.GradientRaw = value; } + } + + public Matrix HessianRaw + { + get { return this.InnerEvaluation.HessianRaw; } + set { this.InnerEvaluation.HessianRaw = value; } + } + + public double Value + { + get + { + + if (!this.ValueChecked) + { + double tmp; + try + { + tmp = this.InnerEvaluation.Value; + } + catch (Exception e) + { + throw new EvaluationException("Objective function evaluation failed.", this.InnerEvaluation, e); + } + this.Checker.ValueChecker(this.InnerEvaluation); + this.ValueChecked = true; + } + return this.InnerEvaluation.Value; + } + } + + public Vector Gradient + { + get + { + + if (!this.GradientChecked) + { + Vector tmp; + try + { + tmp = this.InnerEvaluation.Gradient; + } + catch (Exception e) + { + throw new EvaluationException("Objective gradient evaluation failed.", this.InnerEvaluation, e); + } + this.Checker.GradientChecker(this.InnerEvaluation); + this.GradientChecked = true; + } + return this.InnerEvaluation.Gradient; + } + } + + public Matrix Hessian + { + get + { + + if (!this.HessianChecked) + { + Matrix tmp; + try + { + tmp = this.InnerEvaluation.Hessian; + } + catch (Exception e) + { + throw new EvaluationException("Objective hessian evaluation failed.", this.InnerEvaluation, e); + } + this.Checker.HessianChecker(InnerEvaluation); + this.HessianChecked = true; + } + return this.InnerEvaluation.Hessian; + } + } + + public void Reset(Vector new_point) + { + this.InnerEvaluation.Reset(new_point); + } + } + + public class ObjectiveChecker : IObjectiveFunction + { + public IObjectiveFunction InnerObjective { get; private set; } + public Action ValueChecker { get; private set; } + 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) + { + this.InnerObjective = objective; + this.ValueChecker = value_checker; + this.GradientChecker = gradient_checker; + this.HessianChecker = hessian_checker; + } + + public bool GradientSupported + { + get { return this.InnerObjective.GradientSupported; } + } + + public bool HessianSupported + { + get { return this.InnerObjective.HessianSupported; } + } + + public void Evaluate(Vector point, IEvaluation output) + { + try + { + this.InnerObjective.Evaluate(point, output); + } + catch (Exception e) + { + throw new EvaluationException("Objective evaluation failed.", new NullEvaluation(point), e); + } + } + + + public IEvaluation CreateEvaluationObject() + { + return this.InnerObjective.CreateEvaluationObject(); + } + } +} diff --git a/src/Numerics/Optimization/Implementation/WeakWolfeLineSearch.cs b/src/Numerics/Optimization/Implementation/WeakWolfeLineSearch.cs new file mode 100644 index 00000000..605c424b --- /dev/null +++ b/src/Numerics/Optimization/Implementation/WeakWolfeLineSearch.cs @@ -0,0 +1,119 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization.Implementation +{ + public class WeakWolfeLineSearch + { + public double C1 { get; set; } + public double C2 { 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) + { + this.C1 = c1; + this.C2 = c2; + this.ParameterTolerance = parameter_tolerance; + this.MaximumIterations = max_iterations; + } + + // 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) + { + + if (!(objective is ObjectiveChecker)) + objective = new ObjectiveChecker(objective, this.ValidateValue, this.ValidateGradient, null); + + double lower_bound = 0.0; + double upper_bound = Double.PositiveInfinity; + double step = initial_step; + + double initial_value = starting_point.Value; + Vector initial_gradient = starting_point.Gradient; + + double initial_dd = search_direction * initial_gradient; + + int ii; + IEvaluation candidate_eval = objective.CreateEvaluationObject(); + MinimizationOutput.ExitCondition reason_for_exit = MinimizationOutput.ExitCondition.None; + for (ii = 0; ii < this.MaximumIterations; ++ii) + { + objective.Evaluate(starting_point.Point + search_direction * step, candidate_eval); + + double step_dd = search_direction * candidate_eval.Gradient; + + if (candidate_eval.Value > initial_value + this.C1 * step * initial_dd) + { + upper_bound = step; + step = 0.5 * (lower_bound + upper_bound); + } + else if (step_dd < this.C2 * initial_dd) + { + lower_bound = step; + step = Double.IsPositiveInfinity(upper_bound) ? 2 * lower_bound : 0.5 * (lower_bound + upper_bound); + } + else + { + reason_for_exit = MinimizationOutput.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 = MinimizationOutput.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); + } + + private bool Conforms(IEvaluation starting_point, Vector search_direction, double step, IEvaluation ending_point) + { + + 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; + } + + private 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); + } + + private void ValidateGradient(IEvaluation eval) + { + foreach (double x in eval.Gradient) + if (!this.IsFinite(x)) + { + throw new EvaluationException(String.Format("Non-finite value returned by gradient: {0}", x), eval); + } + } + + private bool IsFinite(double x) + { + return !(Double.IsNaN(x) || Double.IsInfinity(x)); + } + } +} diff --git a/src/Numerics/Optimization/MinimizationOutput.cs b/src/Numerics/Optimization/MinimizationOutput.cs new file mode 100644 index 00000000..d80e4e80 --- /dev/null +++ b/src/Numerics/Optimization/MinimizationOutput.cs @@ -0,0 +1,25 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization +{ + public class MinimizationOutput + { + public enum ExitCondition { None, RelativeGradient, LackOfProgress, AbsoluteGradient, WeakWolfeCriteria, BoundTolerance, StrongWolfeCriteria, LackOfFunctionImprovement } + + 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 MinimizationOutput(IEvaluation function_info, int iterations, ExitCondition reason_for_exit) + { + this.FunctionInfoAtMinimum = function_info; + this.Iterations = iterations; + this.ReasonForExit = reason_for_exit; + } + } +} diff --git a/src/Numerics/Optimization/MinimizationWithLineSearchOutput.cs b/src/Numerics/Optimization/MinimizationWithLineSearchOutput.cs new file mode 100644 index 00000000..2cd00ad9 --- /dev/null +++ b/src/Numerics/Optimization/MinimizationWithLineSearchOutput.cs @@ -0,0 +1,20 @@ +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/NewtonMinimizer.cs b/src/Numerics/Optimization/NewtonMinimizer.cs new file mode 100644 index 00000000..810c859f --- /dev/null +++ b/src/Numerics/Optimization/NewtonMinimizer.cs @@ -0,0 +1,130 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using MathNet.Numerics.LinearAlgebra; +using LU = MathNet.Numerics.LinearAlgebra.Factorization.LU; +using MathNet.Numerics.Optimization.Implementation; + +namespace MathNet.Numerics.Optimization +{ + public class NewtonMinimizer + { + public double GradientTolerance { get; set; } + public int MaximumIterations { get; set; } + public bool UseLineSearch { get; set; } + + public NewtonMinimizer(double gradient_tolerance, int maximum_iterations, bool use_line_search = false) + { + this.GradientTolerance = gradient_tolerance; + this.MaximumIterations = maximum_iterations; + this.UseLineSearch = use_line_search; + } + + public MinimizationOutput FindMinimum(IObjectiveFunction objective, Vector initial_guess) + { + if (!objective.GradientSupported) + throw new IncompatibleObjectiveException("Gradient not supported in objective function, but required for Newton minimization."); + + if (!objective.HessianSupported) + 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); + + IEvaluation initial_eval = objective.CreateEvaluationObject(); + objective.Evaluate(initial_guess, initial_eval); + + // Check that we're not already done + if (this.ExitCriteriaSatisfied(initial_guess, initial_eval.Gradient)) + return new MinimizationOutput(initial_eval, 0, MinimizationOutput.ExitCondition.AbsoluteGradient); + + // Set up line search algorithm + var line_searcher = new WeakWolfeLineSearch(1e-4, 0.9, 1e-4, max_iterations: 1000); + + // Declare state variables + IEvaluation candidate_point = initial_eval; + Vector search_direction; + LineSearchOutput result; + + // 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) + { + + search_direction = candidate_point.Hessian.LU().Solve(-candidate_point.Gradient); + + if (search_direction * candidate_point.Gradient >= 0) + { + search_direction = -candidate_point.Gradient; + steepest_descent_resets += 1; + tmp_line_search = true; + } + + if (this.UseLineSearch || tmp_line_search) + { + try + { + result = line_searcher.FindConformingStep(objective, candidate_point, search_direction, 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; + } + else + { + objective.Evaluate(candidate_point.Point + search_direction, candidate_point); + } + + tmp_line_search = false; + + iterations += 1; + } + + if (iterations == this.MaximumIterations) + throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", this.MaximumIterations)); + + return new MinimizationWithLineSearchOutput(candidate_point, iterations, MinimizationOutput.ExitCondition.AbsoluteGradient, total_line_search_steps, iterations_with_nontrivial_line_search); + } + + private bool ExitCriteriaSatisfied(Vector candidate_point, Vector gradient) + { + return gradient.Norm(2.0) < this.GradientTolerance; + } + + private void ValidateGradient(IEvaluation eval) + { + foreach (var x in eval.Gradient) + { + if (Double.IsNaN(x) || Double.IsInfinity(x)) + throw new EvaluationException("Non-finite gradient returned.", eval); + } + } + + private 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) + { + 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])) + throw new EvaluationException("Non-finite Hessian returned.", eval); + } + } + } + } +} diff --git a/src/UnitTests/OptimizationTests/RosenbrockFunction.cs b/src/UnitTests/OptimizationTests/RosenbrockFunction.cs new file mode 100644 index 00000000..e27b6147 --- /dev/null +++ b/src/UnitTests/OptimizationTests/RosenbrockFunction.cs @@ -0,0 +1,55 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.UnitTests.OptimizationTests +{ + public static class RosenbrockFunction + { + public static double Value(Vector input) + { + 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]); + 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]; + output[1, 1] = 200; + output[0, 1] = -400 * input[0]; + output[1, 0] = output[0, 1]; + return output; + } + } + + public static class BigRosenbrockFunction + { + public static double Value(Vector input) + { + return 1000.0 + 100.0 * RosenbrockFunction.Value(input / 100.0); + } + + public static Vector Gradient(Vector input) + { + return 100.0 * RosenbrockFunction.Gradient(input / 100.0); + } + + public static Matrix Hessian(Vector input) + { + return 100.0 * RosenbrockFunction.Hessian(input / 100.0); + } + + } +} diff --git a/src/UnitTests/OptimizationTests/TestNewtonMinimizer.cs b/src/UnitTests/OptimizationTests/TestNewtonMinimizer.cs new file mode 100644 index 00000000..fd6b6779 --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestNewtonMinimizer.cs @@ -0,0 +1,104 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; + +using NUnit.Framework; +using MathNet.Numerics.Optimization; + +namespace MathNet.Numerics.UnitTests.OptimizationTests +{ + public class RosenbrockEvaluation : BaseEvaluation + { + public const bool SupportsGradient = true; + public const bool SupportsHessian = true; + + protected override void setValue() + { + this.ValueRaw = RosenbrockFunction.Value(this.Point); + } + + protected override void setGradient() + { + this.GradientRaw = RosenbrockFunction.Gradient(this.Point); + } + + protected override void setHessian() + { + this.HessianRaw = RosenbrockFunction.Hessian(this.Point); + } + } + + [TestFixture] + public class TestNewtonMinimizer + { + + [Test] + public void FindMinimum_Rosenbrock_Easy() + { + var obj = new BaseObjectiveFunction(RosenbrockEvaluation.SupportsGradient, RosenbrockEvaluation.SupportsHessian); + + var solver = new NewtonMinimizer(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)); + Assert.That(Math.Abs(result.MinimizingPoint[1] - 1.0), Is.LessThan(1e-3)); + } + + [Test] + public void FindMinimum_Rosenbrock_Hard() + { + var obj = new BaseObjectiveFunction(RosenbrockEvaluation.SupportsGradient, RosenbrockEvaluation.SupportsHessian); + var solver = new NewtonMinimizer(1e-5, 1000); + var result = solver.FindMinimum(obj, new MathNet.Numerics.LinearAlgebra.Double.DenseVector(new double[] { -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 BaseObjectiveFunction(RosenbrockEvaluation.SupportsGradient, RosenbrockEvaluation.SupportsHessian); + var solver = new NewtonMinimizer(1e-5, 1000); + var result = solver.FindMinimum(obj, new MathNet.Numerics.LinearAlgebra.Double.DenseVector(new double[] { -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)); + } + + [Test] + public void FindMinimum_Linesearch_Rosenbrock_Easy() + { + var obj = new BaseObjectiveFunction(RosenbrockEvaluation.SupportsGradient, RosenbrockEvaluation.SupportsHessian); + 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 })); + + 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_Linesearch_Rosenbrock_Hard() + { + var obj = new BaseObjectiveFunction(RosenbrockEvaluation.SupportsGradient, RosenbrockEvaluation.SupportsHessian); + 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 })); + + 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_Linesearch_Rosenbrock_Overton() + { + var obj = new BaseObjectiveFunction(RosenbrockEvaluation.SupportsGradient, RosenbrockEvaluation.SupportsHessian); + 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 })); + + 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/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 814e64c2..61d1ae4f 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -366,6 +366,8 @@ + +