First commit working on optimization tools

(Updated by cdrnet)

src/Numerics/Numerics.csproj

@@ -195,6 +195,15 @@
     <Compile Include="RootFinding\Brent.cs" />
     <Compile Include="FindRoots.cs" />
     <Compile Include="RootFinding\Bisection.cs" />
+    <Compile Include="Optimization\BFGS.cs" />
+    <Compile Include="Optimization\BisectionRootFinder.cs" />
+    <Compile Include="Optimization\ConjugateGradientMinimizer.cs" />
+    <Compile Include="Optimization\GoldenSectionMinimizer.cs" />
+    <Compile Include="Optimization\MinimizationOutput.cs" />
+    <Compile Include="Optimization\ObjectiveFunction.cs" />
+    <Compile Include="Optimization\OptimizationResult.cs" />
+    <Compile Include="Optimization\StrongWolfeLineSearch.cs" />
+    <Compile Include="Optimization\WeakWolfeLineSearch.cs" />
     <Compile Include="SpecialFunctions\Evaluate.cs" />
     <Compile Include="ExcelFunctions.cs" />
     <Compile Include="SpecialFunctions\ModifiedStruve.cs" />

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

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

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<double> 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<double> 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<double> candidate_point, Vector<double> gradient)
+        {
+            return gradient.Norm(2.0) < this.GradientTolerance;
+        }
+
+        private void ValidateGradient(Vector<double> gradient, Vector<double> 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<double> input)
+        {
+            if (Double.IsNaN(objective) || Double.IsInfinity(objective))
+                throw new Exception("Non-finite objective function returned.");
+        }
+    }
+}

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

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

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<double> Point { get; }
+        double Value { get; }
+        Vector<double> Gradient { get; }
+        Matrix<double> Hessian { get; }
+    }
+
+    public interface IObjectiveFunction
+    {
+        bool GradientSupported { get; }
+        bool HessianSupported { get; }
+        IEvaluation Evaluate(Vector<double> point);
+    }
+
+    public class CachedEvaluation : IEvaluation
+    {
+        private double? _value;
+        private Vector<double> _gradient;
+        private Matrix<double> _hessian;
+        private SimpleObjectiveFunction _objective_object;
+        private Vector<double> _point;
+
+        public CachedEvaluation(SimpleObjectiveFunction f, Vector<double> point)
+        {
+            _objective_object = f;
+			_point = point;
+        }
+        private double setValue()
+        {
+            _value = _objective_object.Objective(_point);
+            return _value.Value;
+        }
+        private Vector<double> setGradient()
+        {
+            _gradient = _objective_object.Gradient(_point);
+            return _gradient;
+        }
+        private Matrix<double> setHessian()
+        {
+            _hessian = _objective_object.Hessian(_point);
+            return _hessian;
+        }
+
+        public Vector<double> Point { get { return _point; } }
+        public double Value { get { return _value ?? setValue(); } }
+        public Vector<double> Gradient { get { return _gradient ?? setGradient(); } }
+        public Matrix<double> Hessian { get { return _hessian ?? setHessian(); } }
+
+    }
+
+    public class SimpleObjectiveFunction : IObjectiveFunction
+    {
+        public Func<Vector<double>,double> Objective { get; private set; }
+        public Func<Vector<double>,Vector<double>> Gradient { get; private set; }
+        public Func<Vector<double>, Matrix<double>> Hessian { get; private set; }
+
+        public SimpleObjectiveFunction(Func<Vector<double>,double> objective)
+        {
+            this.Objective = objective;
+            this.Gradient = null;
+            this.Hessian = null;
+        }
+
+        public SimpleObjectiveFunction(Func<Vector<double>,double> objective, Func<Vector<double>,Vector<double>> gradient)
+        {
+            this.Objective = objective;
+            this.Gradient = gradient;
+            this.Hessian = null;
+        }
+
+        public SimpleObjectiveFunction(Func<Vector<double>,double> objective, Func<Vector<double>,Vector<double>> gradient, Func<Vector<double>, Matrix<double>> 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<double> point)
+        {
+            return new CachedEvaluation(this, point);
+        }
+    }
+}

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

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

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

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<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]) + 2 * 100 * (input[1] - input[0] * input[0]) * (-2 * input[0]);
+            output[1] = 2 * 100 * (input[1] - input[0] * input[0]);
+            return output;
+        }
+    }
+}

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

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

src/UnitTests/UnitTests.csproj

@@ -360,6 +360,9 @@
     <Compile Include="EuclidTests\IntegerTheoryTest.cs" />
     <Compile Include="Random\SystemRandomSourceTests.cs" />
     <Compile Include="RootFindingTests\BisectionTest.cs" />
+    <Compile Include="OptimizationTests\RosenbrockFunction.cs" />
+    <Compile Include="OptimizationTests\TestBisectionRootFinder.cs" />
+    <Compile Include="OptimizationTests\TestConjugateGradientMinimizer.cs" />
     <Compile Include="PermutationTest.cs" />
     <Compile Include="PrecisionTest.cs" />
     <Compile Include="Properties\AssemblyInfo.cs" />
@@ -682,5 +685,6 @@
   <ItemGroup>
     <Service Include="{508349B6-6B84-4DF5-91F0-309BEEBAD82D}" />
   </ItemGroup>
+  <ItemGroup />
   <Import Project="$(MSBuildToolsPath)\Microsoft.CSharp.targets" />
 </Project>