From decf9fe5be8b180f5e49c135c658de74e486e9dd Mon Sep 17 00:00:00 2001 From: Scott Stephens Date: Thu, 31 Jan 2013 14:52:54 -0600 Subject: [PATCH] Very minimally tested ConjugateGradient minimizer is complete. --- src/Numerics/Numerics.csproj | 2 + .../ConjugateGradientMinimizer.cs | 39 +++++++++++++++---- src/Numerics/Optimization/LineSearchOutput.cs | 18 +++++++++ .../LineSearchingMinimizerOutput.cs | 20 ++++++++++ .../Optimization/WeakWolfeLineSearch.cs | 22 ++++++++--- .../TestConjugateGradientMinimizer.cs | 12 ++++-- 6 files changed, 97 insertions(+), 16 deletions(-) create mode 100644 src/Numerics/Optimization/LineSearchOutput.cs create mode 100644 src/Numerics/Optimization/LineSearchingMinimizerOutput.cs diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 499dfdae..cadd0bf1 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -199,6 +199,8 @@ + + diff --git a/src/Numerics/Optimization/ConjugateGradientMinimizer.cs b/src/Numerics/Optimization/ConjugateGradientMinimizer.cs index 9b539433..c16e2173 100644 --- a/src/Numerics/Optimization/ConjugateGradientMinimizer.cs +++ b/src/Numerics/Optimization/ConjugateGradientMinimizer.cs @@ -32,7 +32,7 @@ namespace MathNet.Numerics.Optimization return new MinimizationOutput(initial_eval, 0); // Set up line search algorithm - var line_searcher = new WeakWolfeLineSearch(1e-4, 0.1); + var line_searcher = new WeakWolfeLineSearch(1e-4, 0.1,1000); // Declare state variables IEvaluation candidate_point; @@ -41,26 +41,51 @@ namespace MathNet.Numerics.Optimization // First step steepest_direction = -gradient; search_direction = steepest_direction; - var result = line_searcher.FindConformingStep(objective, initial_eval, search_direction, 1.0); + double initial_step_size = 100 * this.GradientTolerance / (gradient * gradient); + var result = line_searcher.FindConformingStep(objective, initial_eval, search_direction, initial_step_size); candidate_point = result.FunctionInfoAtMinimum; this.ValidateGradient(candidate_point.Gradient, candidate_point.Point); + double step_size = (candidate_point.Point - initial_guess).Norm(2.0); // Subsequent steps int iterations = 1; + int total_line_search_steps = result.Iterations; + int no_line_search_iterations = result.Iterations > 0 ? 0 : 1; + int steepest_descent_resets = 0; 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); + var search_direction_adjuster = Math.Max(0,steepest_direction * (steepest_direction - previous_steepest_direction) / (previous_steepest_direction * previous_steepest_direction)); + + //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; - candidate_point = result.FunctionInfoAtMinimum; + //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 + // 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) + { + search_direction = steepest_direction; + steepest_descent_resets += 1; + } + + result = line_searcher.FindConformingStep(objective, candidate_point, search_direction, step_size); + + no_line_search_iterations += result.Iterations == 0 ? 1 : 0; + total_line_search_steps += result.Iterations; + + step_size = (result.FunctionInfoAtMinimum.Point - candidate_point.Point).Norm (2.0); + candidate_point = result.FunctionInfoAtMinimum; iterations += 1; } - return new MinimizationOutput(candidate_point, iterations); + return new MinimizationWithLineSearchOutput(candidate_point, iterations, total_line_search_steps, no_line_search_iterations); } private bool ExitCriteriaSatisfied(Vector candidate_point, Vector gradient) diff --git a/src/Numerics/Optimization/LineSearchOutput.cs b/src/Numerics/Optimization/LineSearchOutput.cs new file mode 100644 index 00000000..bf2db97d --- /dev/null +++ b/src/Numerics/Optimization/LineSearchOutput.cs @@ -0,0 +1,18 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +namespace MathNet.Numerics.Optimization +{ + public class LineSearchOutput : MinimizationOutput + { + public double FinalStep { get; private set; } + + public LineSearchOutput(IEvaluation function_info, int iterations, double final_step) + : base(function_info, iterations) + { + this.FinalStep = final_step; + } + } +} diff --git a/src/Numerics/Optimization/LineSearchingMinimizerOutput.cs b/src/Numerics/Optimization/LineSearchingMinimizerOutput.cs new file mode 100644 index 00000000..e0c50add --- /dev/null +++ b/src/Numerics/Optimization/LineSearchingMinimizerOutput.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 IterationsWithNoLineSearch { get; private set; } + + public MinimizationWithLineSearchOutput(IEvaluation function_info, int iterations, int total_line_search_iterations, int iterations_with_no_line_search) + : base(function_info, iterations) + { + this.TotalLineSearchIterations = total_line_search_iterations; + this.IterationsWithNoLineSearch = iterations_with_no_line_search; + } + } +} diff --git a/src/Numerics/Optimization/WeakWolfeLineSearch.cs b/src/Numerics/Optimization/WeakWolfeLineSearch.cs index 12e53aef..7286a79c 100644 --- a/src/Numerics/Optimization/WeakWolfeLineSearch.cs +++ b/src/Numerics/Optimization/WeakWolfeLineSearch.cs @@ -11,15 +11,17 @@ namespace MathNet.Numerics.Optimization { public double C1 { get; set; } public double C2 { get; set; } + public int MaximumIterations { get; set; } - public WeakWolfeLineSearch(double c1, double c2) + public WeakWolfeLineSearch(double c1, double c2, int max_iterations=10) { this.C1 = c1; this.C2 = c2; + this.MaximumIterations = max_iterations; } // 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) + public LineSearchOutput FindConformingStep(IObjectiveFunction objective, IEvaluation starting_point, Vector search_direction, double initial_step) { double lower_bound = 0.0; double upper_bound = Double.PositiveInfinity; @@ -32,7 +34,7 @@ namespace MathNet.Numerics.Optimization int ii; IEvaluation candidate_eval = null; - for (ii = 0; ii < 10; ++ii) + for (ii = 0; ii < this.MaximumIterations; ++ii) { candidate_eval = objective.Evaluate(starting_point.Point + search_direction * step); @@ -53,13 +55,21 @@ namespace MathNet.Numerics.Optimization break; } } - - if (ii == 10) + + if (ii == this.MaximumIterations) throw new Exception("Line search failed with max iterations. Function is likely unbounded in search direction."); else - return new MinimizationOutput(candidate_eval, ii); + return new LineSearchOutput(candidate_eval, ii, step); } + 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; + } } } diff --git a/src/UnitTests/OptimizationTests/TestConjugateGradientMinimizer.cs b/src/UnitTests/OptimizationTests/TestConjugateGradientMinimizer.cs index b06f410c..79b12b45 100644 --- a/src/UnitTests/OptimizationTests/TestConjugateGradientMinimizer.cs +++ b/src/UnitTests/OptimizationTests/TestConjugateGradientMinimizer.cs @@ -17,16 +17,22 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests public void FindMinimum_Rosenbrock_Easy() { var obj = new SimpleObjectiveFunction(RosenbrockFunction.Value, RosenbrockFunction.Gradient); - var solver = new ConjugateGradientMinimizer(1e-5, 100); + 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(result.MinimizingPoint[0], Is.EqualTo(1.0)); - Assert.That(result.MinimizingPoint[1], Is.EqualTo(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_Hard() { + 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 })); + 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)); } } }