diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 163c728e..aab028a8 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -195,7 +195,7 @@ - + @@ -203,6 +203,7 @@ + diff --git a/src/Numerics/Optimization/BFGS.cs b/src/Numerics/Optimization/BFGS.cs deleted file mode 100644 index 7b9942d2..00000000 --- a/src/Numerics/Optimization/BFGS.cs +++ /dev/null @@ -1,11 +0,0 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; - -namespace MathNet.Numerics.Optimization -{ - class BFGS - { - } -} diff --git a/src/Numerics/Optimization/BfgsMinimizer.cs b/src/Numerics/Optimization/BfgsMinimizer.cs new file mode 100644 index 00000000..8068a3cf --- /dev/null +++ b/src/Numerics/Optimization/BfgsMinimizer.cs @@ -0,0 +1,131 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization +{ + public class BfgsMinimizer + { + public double GradientTolerance { get; set; } + public int MaximumIterations { get; set; } + + public BfgsMinimizer(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 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); + + 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); + + // Set up line search algorithm + var line_searcher = new WeakWolfeLineSearch(1e-4, 0.9, 1000); + + // Declare state variables + IEvaluation candidate_point; + double step_size; + Vector gradient, previous_gradient, step, search_direction; + Matrix inverse_pseudo_hessian; + + // 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); + + LineSearchOutput result; + try + { + result = line_searcher.FindConformingStep(objective, initial_eval, search_direction, step_size); + } + catch (Exception e) + { + throw new InnerOptimizationException("Line search failed.", e); + } + + candidate_point = result.FunctionInfoAtMinimum; + gradient = candidate_point.Gradient; + previous_gradient = initial_eval.Gradient; + step = candidate_point.Point - initial_guess; + 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) + { + var y = candidate_point.Gradient - previous_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); + + search_direction = -inverse_pseudo_hessian * candidate_point.Gradient; + + if (search_direction * candidate_point.Gradient >= 0) + { + search_direction = -candidate_point.Gradient; + inverse_pseudo_hessian = Matrix.Build.DiagonalIdentity(initial_guess.Count); + steepest_descent_resets += 1; + } + + 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; + + step_size = result.FinalStep; + step = result.FunctionInfoAtMinimum.Point - candidate_point.Point; + previous_gradient = candidate_point.Gradient; + candidate_point = result.FunctionInfoAtMinimum; + + iterations += 1; + } + + if (iterations == this.MaximumIterations) + throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", this.MaximumIterations)); + + return new MinimizationWithLineSearchOutput(candidate_point, iterations, 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(Vector gradient, Vector input) + { + foreach (var x in gradient) + { + if (Double.IsNaN(x) || Double.IsInfinity(x)) + throw new EvaluationException("Non-finite gradient returned."); + } + } + + private void ValidateObjective(double objective, Vector input) + { + if (Double.IsNaN(objective) || Double.IsInfinity(objective)) + throw new EvaluationException("Non-finite objective function returned."); + } + } +} diff --git a/src/Numerics/Optimization/ConjugateGradientMinimizer.cs b/src/Numerics/Optimization/ConjugateGradientMinimizer.cs index 4fe75e83..34758ad6 100644 --- a/src/Numerics/Optimization/ConjugateGradientMinimizer.cs +++ b/src/Numerics/Optimization/ConjugateGradientMinimizer.cs @@ -23,7 +23,8 @@ namespace MathNet.Numerics.Optimization if (!objective.GradientSupported) throw new IncompatibleObjectiveException("Gradient not supported in objective function, but required for ConjugateGradient minimization."); - objective = new ObjectiveChecker(objective, this.ValidateObjective, this.ValidateGradient, null); + if (!(objective is ObjectiveChecker)) + objective = new ObjectiveChecker(objective, this.ValidateObjective, this.ValidateGradient, null); IEvaluation initial_eval = objective.Evaluate(initial_guess); var gradient = initial_eval.Gradient; @@ -53,14 +54,14 @@ namespace MathNet.Numerics.Optimization throw new InnerOptimizationException("Line search failed.", e); } - candidate_point = result.FunctionInfoAtMinimum; - - double step_size = (candidate_point.Point - initial_guess).Norm(2.0); + candidate_point = result.FunctionInfoAtMinimum; + + double step_size = result.FinalStep; // Subsequent steps int iterations = 1; int total_line_search_steps = result.Iterations; - int no_line_search_iterations = result.Iterations > 0 ? 0 : 1; + 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) { @@ -93,10 +94,10 @@ namespace MathNet.Numerics.Optimization throw new InnerOptimizationException("Line search failed.", e); } - no_line_search_iterations += result.Iterations == 0 ? 1 : 0; + iterations_with_nontrivial_line_search += result.Iterations > 0 ? 1 : 0; total_line_search_steps += result.Iterations; - step_size = (result.FunctionInfoAtMinimum.Point - candidate_point.Point).Norm (2.0); + step_size = result.FinalStep; candidate_point = result.FunctionInfoAtMinimum; iterations += 1; @@ -105,7 +106,7 @@ namespace MathNet.Numerics.Optimization if (iterations == this.MaximumIterations) throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", this.MaximumIterations)); - return new MinimizationWithLineSearchOutput(candidate_point, iterations, total_line_search_steps, no_line_search_iterations); + return new MinimizationWithLineSearchOutput(candidate_point, iterations, total_line_search_steps, iterations_with_nontrivial_line_search); } private bool ExitCriteriaSatisfied(Vector candidate_point, Vector gradient) diff --git a/src/Numerics/Optimization/LineSearchingMinimizerOutput.cs b/src/Numerics/Optimization/LineSearchingMinimizerOutput.cs index e0c50add..a6890a1d 100644 --- a/src/Numerics/Optimization/LineSearchingMinimizerOutput.cs +++ b/src/Numerics/Optimization/LineSearchingMinimizerOutput.cs @@ -8,13 +8,13 @@ namespace MathNet.Numerics.Optimization public class MinimizationWithLineSearchOutput : MinimizationOutput { public int TotalLineSearchIterations { get; private set; } - public int IterationsWithNoLineSearch { get; private set; } + public int IterationsWithNonTrivialLineSearch { get; private set; } - public MinimizationWithLineSearchOutput(IEvaluation function_info, int iterations, int total_line_search_iterations, int iterations_with_no_line_search) + public MinimizationWithLineSearchOutput(IEvaluation function_info, int iterations, int total_line_search_iterations, int iterations_with_non_trivial_line_search) : base(function_info, iterations) { this.TotalLineSearchIterations = total_line_search_iterations; - this.IterationsWithNoLineSearch = iterations_with_no_line_search; + 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..901b5ecc --- /dev/null +++ b/src/Numerics/Optimization/NewtonMinimizer.cs @@ -0,0 +1,128 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using MathNet.Numerics.LinearAlgebra; +using LU = MathNet.Numerics.LinearAlgebra.Factorization.LU; + +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.Evaluate(initial_guess); + + // Check that we're not already done + if (this.ExitCriteriaSatisfied(initial_guess, initial_eval.Gradient)) + return new MinimizationOutput(initial_eval, 0); + + // Set up line search algorithm + var line_searcher = new WeakWolfeLineSearch(1e-4, 0.9, 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 + { + candidate_point = objective.Evaluate(candidate_point.Point + search_direction); + } + + 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, 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(Vector gradient, Vector input) + { + foreach (var x in gradient) + { + if (Double.IsNaN(x) || Double.IsInfinity(x)) + throw new EvaluationException("Non-finite gradient returned."); + } + } + + private void ValidateObjective(double objective, Vector input) + { + if (Double.IsNaN(objective) || Double.IsInfinity(objective)) + throw new EvaluationException("Non-finite objective function returned."); + } + + private void ValidateHessian(Matrix hessian, Vector input) + { + for (int ii = 0; ii < hessian.RowCount; ++ii) + { + for (int jj = 0; jj < hessian.ColumnCount; ++jj) + { + if (Double.IsNaN(hessian[ii,jj]) || Double.IsInfinity(hessian[ii,jj])) + throw new EvaluationException("Non-finite Hessian returned."); + } + } + } + } +} diff --git a/src/Numerics/Optimization/WeakWolfeLineSearch.cs b/src/Numerics/Optimization/WeakWolfeLineSearch.cs index 0690a100..809a3068 100644 --- a/src/Numerics/Optimization/WeakWolfeLineSearch.cs +++ b/src/Numerics/Optimization/WeakWolfeLineSearch.cs @@ -23,6 +23,10 @@ namespace MathNet.Numerics.Optimization // 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; @@ -71,5 +75,24 @@ namespace MathNet.Numerics.Optimization return step > 0 && sufficient_decrease && not_too_steep; } + private void ValidateValue(double value, Vector input) + { + if (!this.IsFinite(value)) + throw new EvaluationException(String.Format("Non-finite value returned by objective function: {0}", value)); + } + + private void ValidateGradient(Vector gradient, Vector input) + { + foreach (double x in gradient) + if (!this.IsFinite(x)) + { + throw new EvaluationException(String.Format("Non-finite value returned by gradient: {0}", x)); + } + } + + private 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 1e82c22a..7dea0f51 100644 --- a/src/UnitTests/OptimizationTests/RosenbrockFunction.cs +++ b/src/UnitTests/OptimizationTests/RosenbrockFunction.cs @@ -17,9 +17,20 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests 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[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; + } } } diff --git a/src/UnitTests/OptimizationTests/TestBfgsMinimizer.cs b/src/UnitTests/OptimizationTests/TestBfgsMinimizer.cs new file mode 100644 index 00000000..742f7b59 --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestBfgsMinimizer.cs @@ -0,0 +1,50 @@ +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 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 })); + + 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 BfgsMinimizer(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 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 })); + + 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 new file mode 100644 index 00000000..e97ffd86 --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestNewtonMinimizer.cs @@ -0,0 +1,82 @@ +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 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 })); + + 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, 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 })); + + 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, 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 })); + + 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 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 })); + + 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 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 })); + + 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 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 })); + + 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 new file mode 100644 index 00000000..28098b39 --- /dev/null +++ b/src/UnitTests/OptimizationTests/TestRosenbrockFunction.cs @@ -0,0 +1,59 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using NUnit.Framework; +namespace MathNet.Numerics.UnitTests.OptimizationTests +{ + [TestFixture] + class TestRosenbrockFunction + { + [Test] + public void TestGradient() + { + var input = new LinearAlgebra.Double.DenseVector(new double[]{ -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; + + 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); + } + + [Test] + public void TestHessian() + { + var input = new LinearAlgebra.Double.DenseVector(new double[] { -0.9, -0.5 }); + + var v1 = RosenbrockFunction.Value(input); + var h = RosenbrockFunction.Hessian(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; + + 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); + + + 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); + } + } +} diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 2da443ab..ffca80d1 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -361,8 +361,11 @@ + + +