From dafc80100cafa3de81c72961e83614d9cacfe935 Mon Sep 17 00:00:00 2001 From: Scott Stephens Date: Thu, 31 Jan 2013 09:23:50 -0600 Subject: [PATCH] First commit working on optimization tools (Updated by cdrnet) --- src/Numerics/Numerics.csproj | 9 ++ src/Numerics/Optimization/BFGS.cs | 11 ++ .../Optimization/BisectionRootFinder.cs | 78 +++++++++++++ .../ConjugateGradientMinimizer.cs | 86 +++++++++++++++ .../Optimization/GoldenSectionMinimizer.cs | 11 ++ .../Optimization/MinimizationOutput.cs | 22 ++++ .../Optimization/ObjectiveFunction.cs | 103 ++++++++++++++++++ .../Optimization/OptimizationResult.cs | 8 ++ .../Optimization/StrongWolfeLineSearch.cs | 11 ++ .../Optimization/WeakWolfeLineSearch.cs | 65 +++++++++++ .../OptimizationTests/RosenbrockFunction.cs | 25 +++++ .../TestBisectionRootFinder.cs | 30 +++++ .../TestConjugateGradientMinimizer.cs | 32 ++++++ src/UnitTests/UnitTests.csproj | 4 + 14 files changed, 495 insertions(+) create mode 100644 src/Numerics/Optimization/BFGS.cs create mode 100644 src/Numerics/Optimization/BisectionRootFinder.cs create mode 100644 src/Numerics/Optimization/ConjugateGradientMinimizer.cs create mode 100644 src/Numerics/Optimization/GoldenSectionMinimizer.cs create mode 100644 src/Numerics/Optimization/MinimizationOutput.cs create mode 100644 src/Numerics/Optimization/ObjectiveFunction.cs create mode 100644 src/Numerics/Optimization/OptimizationResult.cs create mode 100644 src/Numerics/Optimization/StrongWolfeLineSearch.cs create mode 100644 src/Numerics/Optimization/WeakWolfeLineSearch.cs create mode 100644 src/UnitTests/OptimizationTests/RosenbrockFunction.cs create mode 100644 src/UnitTests/OptimizationTests/TestBisectionRootFinder.cs create mode 100644 src/UnitTests/OptimizationTests/TestConjugateGradientMinimizer.cs diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 31785ea8..499dfdae 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -195,6 +195,15 @@ + + + + + + + + + diff --git a/src/Numerics/Optimization/BFGS.cs b/src/Numerics/Optimization/BFGS.cs new file mode 100644 index 00000000..7b9942d2 --- /dev/null +++ b/src/Numerics/Optimization/BFGS.cs @@ -0,0 +1,11 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +namespace MathNet.Numerics.Optimization +{ + class BFGS + { + } +} diff --git a/src/Numerics/Optimization/BisectionRootFinder.cs b/src/Numerics/Optimization/BisectionRootFinder.cs new file mode 100644 index 00000000..f19444f1 --- /dev/null +++ b/src/Numerics/Optimization/BisectionRootFinder.cs @@ -0,0 +1,78 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +namespace MathNet.Numerics.Optimization +{ + public class BisectionRootFinder + { + + public double ObjectiveTolerance { get; set; } + public double XTolerance { get; set; } + public double LowerExpansionFactor { get; set; } + public double UpperExapansionFactor { 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) + { + this.ObjectiveTolerance = objective_tolerance; + this.XTolerance = x_tolerance; + this.LowerExpansionFactor = lower_expansion_factor; + this.UpperExapansionFactor = upper_expansion_factor; + this.MaxExpansionSteps = max_expansion_steps; + } + + public double FindRoot(Func objective_function, double lower_bound, double upper_bound) + { + double lower_val = objective_function(lower_bound); + double upper_val = objective_function(upper_bound); + + if (lower_val == 0.0) + return lower_bound; + if (upper_val == 0.0) + return upper_bound; + + this.ValidateEvaluation(lower_val, lower_bound); + this.ValidateEvaluation(upper_val, upper_bound); + + if (Math.Sign(lower_val) == Math.Sign(upper_val) && this.LowerExpansionFactor <= 1.0 && this.UpperExapansionFactor <= 1.0) + throw new Exception("Bounds do not necessarily span a root, and StepExpansionFactor is not set to expand the interval in this case."); + + while (Math.Abs(upper_val - lower_val) > 0.5 * this.ObjectiveTolerance || Math.Abs(upper_bound - lower_bound) > 0.5 * this.XTolerance) + { + double midpoint = 0.5 * (upper_bound + lower_bound); + double midval = objective_function(midpoint); + this.ValidateEvaluation(midval, midpoint); + + if (Math.Sign(midval) == Math.Sign(lower_val)) + { + lower_bound = midpoint; + lower_val = midval; + } + else if (Math.Sign(midval) == Math.Sign(upper_val)) + { + upper_bound = midpoint; + upper_val = midval; + } + else + { + return midpoint; + } + } + + return 0.5 * (lower_bound + upper_bound); + } + + private 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) + { + return !(Double.IsInfinity(x) || Double.IsNaN(x)); + } + } +} diff --git a/src/Numerics/Optimization/ConjugateGradientMinimizer.cs b/src/Numerics/Optimization/ConjugateGradientMinimizer.cs new file mode 100644 index 00000000..9b539433 --- /dev/null +++ b/src/Numerics/Optimization/ConjugateGradientMinimizer.cs @@ -0,0 +1,86 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization +{ + public class ConjugateGradientMinimizer + { + public double GradientTolerance { get; set; } + public int MaximumIterations { get; set; } + + public ConjugateGradientMinimizer(double gradient_tolerance, int maximum_iterations) + { + this.GradientTolerance = gradient_tolerance; + this.MaximumIterations = maximum_iterations; + } + + public MinimizationOutput FindMinimum(IObjectiveFunction objective, Vector initial_guess) + { + if (!objective.GradientSupported) + throw new Exception("Gradient not supported in objective function, but required for ConjugateGradient minimization."); + + IEvaluation initial_eval = objective.Evaluate(initial_guess); + var gradient = initial_eval.Gradient; + this.ValidateGradient(gradient, initial_guess); + + // Check that we're not already done + if (this.ExitCriteriaSatisfied(initial_guess, gradient)) + return new MinimizationOutput(initial_eval, 0); + + // Set up line search algorithm + var line_searcher = new WeakWolfeLineSearch(1e-4, 0.1); + + // Declare state variables + IEvaluation candidate_point; + Vector steepest_direction, previous_steepest_direction, search_direction; + + // First step + steepest_direction = -gradient; + search_direction = steepest_direction; + var result = line_searcher.FindConformingStep(objective, initial_eval, search_direction, 1.0); + candidate_point = result.FunctionInfoAtMinimum; + this.ValidateGradient(candidate_point.Gradient, candidate_point.Point); + + // Subsequent steps + int iterations = 1; + while (!this.ExitCriteriaSatisfied(candidate_point.Point, candidate_point.Gradient) && iterations < this.MaximumIterations) + { + previous_steepest_direction = steepest_direction; + steepest_direction = -candidate_point.Gradient; + var search_direction_adjuster = steepest_direction * (steepest_direction - previous_steepest_direction) / (previous_steepest_direction * previous_steepest_direction); + search_direction = steepest_direction + search_direction_adjuster * previous_steepest_direction; + result = line_searcher.FindConformingStep(objective, candidate_point, search_direction, 1.0); + + candidate_point = result.FunctionInfoAtMinimum; + + iterations += 1; + } + + return new MinimizationOutput(candidate_point, iterations); + } + + private bool ExitCriteriaSatisfied(Vector candidate_point, Vector gradient) + { + return gradient.Norm(2.0) < this.GradientTolerance; + } + + private void ValidateGradient(Vector gradient, Vector input) + { + foreach (var x in gradient) + { + if (Double.IsNaN(x) || Double.IsInfinity(x)) + throw new Exception("Non-finite gradient returned."); + } + } + + private void ValidateObjective(double objective, Vector input) + { + if (Double.IsNaN(objective) || Double.IsInfinity(objective)) + throw new Exception("Non-finite objective function returned."); + } + } +} diff --git a/src/Numerics/Optimization/GoldenSectionMinimizer.cs b/src/Numerics/Optimization/GoldenSectionMinimizer.cs new file mode 100644 index 00000000..afa35d26 --- /dev/null +++ b/src/Numerics/Optimization/GoldenSectionMinimizer.cs @@ -0,0 +1,11 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +namespace MathNet.Numerics.Optimization +{ + class GoldenSectionMinimizer + { + } +} diff --git a/src/Numerics/Optimization/MinimizationOutput.cs b/src/Numerics/Optimization/MinimizationOutput.cs new file mode 100644 index 00000000..25a09cc4 --- /dev/null +++ b/src/Numerics/Optimization/MinimizationOutput.cs @@ -0,0 +1,22 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization +{ + public class MinimizationOutput + { + public Vector MinimizingPoint { get { return FunctionInfoAtMinimum.Point; } } + public IEvaluation FunctionInfoAtMinimum { get; private set; } + public int Iterations { get; private set; } + + public MinimizationOutput(IEvaluation function_info, int iterations) + { + this.FunctionInfoAtMinimum = function_info; + this.Iterations = iterations; + } + } +} diff --git a/src/Numerics/Optimization/ObjectiveFunction.cs b/src/Numerics/Optimization/ObjectiveFunction.cs new file mode 100644 index 00000000..969d6858 --- /dev/null +++ b/src/Numerics/Optimization/ObjectiveFunction.cs @@ -0,0 +1,103 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization +{ + public interface IEvaluation + { + Vector Point { get; } + double Value { get; } + Vector Gradient { get; } + Matrix Hessian { get; } + } + + public interface IObjectiveFunction + { + bool GradientSupported { get; } + bool HessianSupported { get; } + IEvaluation Evaluate(Vector point); + } + + public class CachedEvaluation : IEvaluation + { + private double? _value; + private Vector _gradient; + private Matrix _hessian; + private SimpleObjectiveFunction _objective_object; + private Vector _point; + + public CachedEvaluation(SimpleObjectiveFunction f, Vector point) + { + _objective_object = f; + _point = point; + } + private double setValue() + { + _value = _objective_object.Objective(_point); + return _value.Value; + } + private Vector setGradient() + { + _gradient = _objective_object.Gradient(_point); + return _gradient; + } + private Matrix setHessian() + { + _hessian = _objective_object.Hessian(_point); + return _hessian; + } + + public Vector Point { get { return _point; } } + public double Value { get { return _value ?? setValue(); } } + public Vector Gradient { get { return _gradient ?? setGradient(); } } + public Matrix Hessian { get { return _hessian ?? setHessian(); } } + + } + + public class SimpleObjectiveFunction : IObjectiveFunction + { + 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) + { + this.Objective = objective; + this.Gradient = null; + this.Hessian = null; + } + + public SimpleObjectiveFunction(Func,double> objective, Func,Vector> gradient) + { + this.Objective = objective; + this.Gradient = gradient; + this.Hessian = null; + } + + public SimpleObjectiveFunction(Func,double> objective, Func,Vector> gradient, Func, Matrix> hessian) + { + this.Objective = objective; + this.Gradient = gradient; + this.Hessian = hessian; + } + + public bool GradientSupported + { + get { return this.Gradient != null; } + } + + public bool HessianSupported + { + get { return this.Hessian != null; } + } + + public IEvaluation Evaluate(Vector point) + { + return new CachedEvaluation(this, point); + } + } +} diff --git a/src/Numerics/Optimization/OptimizationResult.cs b/src/Numerics/Optimization/OptimizationResult.cs new file mode 100644 index 00000000..44017dec --- /dev/null +++ b/src/Numerics/Optimization/OptimizationResult.cs @@ -0,0 +1,8 @@ +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 new file mode 100644 index 00000000..236c63c9 --- /dev/null +++ b/src/Numerics/Optimization/StrongWolfeLineSearch.cs @@ -0,0 +1,11 @@ +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 new file mode 100644 index 00000000..12e53aef --- /dev/null +++ b/src/Numerics/Optimization/WeakWolfeLineSearch.cs @@ -0,0 +1,65 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization +{ + public class WeakWolfeLineSearch + { + public double C1 { get; set; } + public double C2 { get; set; } + + public WeakWolfeLineSearch(double c1, double c2) + { + this.C1 = c1; + this.C2 = c2; + } + + // Implemented following http://www.math.washington.edu/~burke/crs/408/lectures/L9-weak-Wolfe.pdf + public MinimizationOutput FindConformingStep(IObjectiveFunction objective, IEvaluation starting_point, Vector search_direction, double initial_step) + { + 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 = null; + for (ii = 0; ii < 10; ++ii) + { + candidate_eval = objective.Evaluate(starting_point.Point + search_direction * step); + + 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 + { + break; + } + } + + if (ii == 10) + throw new Exception("Line search failed with max iterations. Function is likely unbounded in search direction."); + else + return new MinimizationOutput(candidate_eval, ii); + } + + + } +} diff --git a/src/UnitTests/OptimizationTests/RosenbrockFunction.cs b/src/UnitTests/OptimizationTests/RosenbrockFunction.cs new file mode 100644 index 00000000..1e82c22a --- /dev/null +++ b/src/UnitTests/OptimizationTests/RosenbrockFunction.cs @@ -0,0 +1,25 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +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]) + 2 * 100 * (input[1] - input[0] * input[0]) * (-2 * input[0]); + output[1] = 2 * 100 * (input[1] - input[0] * input[0]); + return output; + } + } +} diff --git a/src/UnitTests/OptimizationTests/TestBisectionRootFinder.cs b/src/UnitTests/OptimizationTests/TestBisectionRootFinder.cs new file mode 100644 index 00000000..561e9965 --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestBisectionRootFinder.cs @@ -0,0 +1,30 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using NUnit.Framework; + +using MathNet.Numerics.Optimization; + +namespace MathNet.Numerics.UnitTests.OptimizationTests +{ + [TestFixture] + 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 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 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 new file mode 100644 index 00000000..b06f410c --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestConjugateGradientMinimizer.cs @@ -0,0 +1,32 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +using NUnit.Framework; + +using MathNet.Numerics.Optimization; + +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, 100); + var result = solver.FindMinimum(obj, new MathNet.Numerics.LinearAlgebra.Double.DenseVector(new double[]{1.2,1.2})); + Assert.That(result.MinimizingPoint[0], Is.EqualTo(1.0)); + Assert.That(result.MinimizingPoint[1], Is.EqualTo(1.0)); + } + + [Test] + public void FindMinimum_Rosenbrock_Hard() + { + + } + } +} diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 54104d8d..2da443ab 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -360,6 +360,9 @@ + + + @@ -682,5 +685,6 @@ + \ No newline at end of file