Optimization: initial commit of candidate API and Newton implementation

src/Numerics/Numerics.csproj

@@ -115,6 +115,19 @@
     <Compile Include="LinearRegression\Options.cs" />
     <Compile Include="OdeSolvers\AdamsBashforth.cs" />
     <Compile Include="OdeSolvers\RungeKutta.cs" />
+    <Compile Include="Optimization\BaseEvaluation.cs" />
+    <Compile Include="Optimization\Exceptions.cs" />
+    <Compile Include="Optimization\IEvaluation.cs" />
+    <Compile Include="Optimization\Implementation\LineSearchOutput.cs" />
+    <Compile Include="Optimization\Implementation\NullEvaluation.cs" />
+    <Compile Include="Optimization\Implementation\ObjectiveChecker.cs" />
+    <Compile Include="Optimization\Implementation\WeakWolfeLineSearch.cs" />
+    <Compile Include="Optimization\IObjectiveFunction.cs" />
+    <Compile Include="Optimization\BaseObjectiveFunction.cs" />
+    <Compile Include="Optimization\IUnconstrainedMinimizer.cs" />
+    <Compile Include="Optimization\MinimizationOutput.cs" />
+    <Compile Include="Optimization\MinimizationWithLineSearchOutput.cs" />
+    <Compile Include="Optimization\NewtonMinimizer.cs" />
     <Compile Include="Precision.Comparison.cs" />
     <Compile Include="Precision.Equality.cs" />
     <Compile Include="Distributions\Bernoulli.cs" />
@@ -487,5 +500,6 @@
       <LastGenOutput>Resources.Designer.cs</LastGenOutput>
     </EmbeddedResource>
   </ItemGroup>
+  <ItemGroup />
   <Import Project="$(MSBuildToolsPath)\Microsoft.CSharp.targets" />
 </Project>

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<double> Point { get; set; }
+        public double ValueRaw { get; set; }
+        public Vector<double> GradientRaw { get; set; }
+        public Matrix<double> 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<double> Gradient
+        {
+            get
+            {
+                if (!Status.HasFlag(EvaluationStatus.Gradient))
+                {
+                    setGradient();
+                    Status |= EvaluationStatus.Gradient;
+                }
+                return GradientRaw;
+            }
+        }
+        public Matrix<double> Hessian
+        {
+            get
+            {
+                if (!Status.HasFlag(EvaluationStatus.Hessian))
+                {
+                    setHessian();
+                    Status |= EvaluationStatus.Hessian;
+                }
+                return HessianRaw;
+            }
+        }
+
+        public void Reset(Vector<double> new_point)
+        {
+            this.Point = new_point;
+            this.Status = EvaluationStatus.None;
+        }
+
+        protected abstract void setValue();
+        protected abstract void setGradient();
+        protected abstract void setHessian();
+    }
+}

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<T> : 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<double> point, IEvaluation output)
+        {
+            output.Reset(point);
+        }
+
+        public virtual IEvaluation CreateEvaluationObject()
+        {
+            return default(T);
+        }
+    }
+}

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) { }
+    }
+}

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<double> Point { get; set; }
+        EvaluationStatus Status { get; set; }
+
+        // Used by algorithm
+        double Value { get; }
+        Vector<double> Gradient { get; }
+        Matrix<double> Hessian { get; }
+
+        // Used by ObjectiveFunction
+        void Reset(Vector<double> new_point);
+        double ValueRaw { get; set; }
+        Vector<double> GradientRaw { get; set; }
+        Matrix<double> HessianRaw { get; set; }
+    }
+}

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<double> point, IEvaluation output);   
+    }
+}

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<double> initial_guess);
+    }
+
+}

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;
+        }
+    }
+}

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<double> 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();
+        }
+    }
+}

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<double> 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<double> GradientRaw
+        {
+            get { return this.InnerEvaluation.GradientRaw;  }
+            set { this.InnerEvaluation.GradientRaw = value; }
+        }
+
+        public Matrix<double> 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<double> Gradient
+        {
+            get
+            {
+
+                if (!this.GradientChecked)
+                {
+                    Vector<double> 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<double> Hessian
+        {
+            get
+            {
+
+                if (!this.HessianChecked)
+                {
+                    Matrix<double> 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<double> new_point)
+        {
+            this.InnerEvaluation.Reset(new_point);
+        }
+    }
+
+    public class ObjectiveChecker : IObjectiveFunction
+    {
+        public IObjectiveFunction InnerObjective { get; private set; }
+        public Action<IEvaluation> ValueChecker { get; private set; }
+        public Action<IEvaluation> GradientChecker { get; private set; }
+        public Action<IEvaluation> HessianChecker { get; private set; }
+
+        public ObjectiveChecker(IObjectiveFunction objective, Action<IEvaluation> value_checker, Action<IEvaluation> gradient_checker, Action<IEvaluation> 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<double> 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();
+        }
+    }
+}

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<double> 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<double> 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<double> 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));
+        }
+    }
+}

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<double> 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;
+        }
+    }
+}

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;
+        }
+    }
+}

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<double>;
+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<double> 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<double> 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<double> candidate_point, Vector<double> 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);
+                }
+            }
+        }
+    }
+}

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<double> input)
+        {
+            return Math.Pow((1 - input[0]), 2) + 100 * Math.Pow((input[1] - input[0] * input[0]), 2);
+        }
+
+        public static Vector<double> Gradient(Vector<double> input)
+        {
+            Vector<double> 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<double> Hessian(Vector<double> input)
+        {
+
+            Matrix<double> 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<double> input)
+        {
+            return 1000.0 + 100.0 * RosenbrockFunction.Value(input / 100.0);
+        }
+
+        public static Vector<double> Gradient(Vector<double> input)
+        {
+            return 100.0 * RosenbrockFunction.Gradient(input / 100.0);
+        }
+
+        public static Matrix<double> Hessian(Vector<double> input)
+        {
+            return 100.0 * RosenbrockFunction.Hessian(input / 100.0);
+        }
+
+    }
+}

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>(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>(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>(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>(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>(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>(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));
+        }
+    }
+}

src/UnitTests/UnitTests.csproj

@@ -366,6 +366,8 @@
     <Compile Include="LinearAlgebraTests\VectorToStringTests.cs" />
     <Compile Include="OdeSolvers\OdeSolverTest.cs" />
     <Compile Include="Random\RandomSerializationTests.cs" />
+    <Compile Include="OptimizationTests\RosenbrockFunction.cs" />
+    <Compile Include="OptimizationTests\TestNewtonMinimizer.cs" />
     <Compile Include="Random\SystemRandomSourceTests.cs" />
     <Compile Include="RootFindingTests\BfgsTest.cs" />
     <Compile Include="RootFindingTests\BisectionTest.cs" />