Merge pull request #489 from eriove/unified_optimization

Unified optimization

src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs

@@ -1720,7 +1720,8 @@ namespace MathNet.Numerics.LinearAlgebra
         /// matrix and a given other matrix being the 'x' of atan2 and the
         /// 'this' matrix being the 'y'
         /// </summary>
-        /// <param name="other"></param>
+        /// <param name="other">The other matrix 'y'</param>
+        /// <param name="result">The matrix with the result and 'x'</param>
         /// <returns></returns>
         public void PointwiseAtan2(Matrix<T> other, Matrix<T> result)
         {

src/Numerics/LinearAlgebra/Vector.Arithmetic.cs

@@ -1028,6 +1028,7 @@ namespace MathNet.Numerics.LinearAlgebra
         /// </summary>
         /// <param name="f">Function which takes a scalar and a vector, modifies the vector in place and returns void</param>
         /// <param name="x">The scalar to be passed to the function</param>
+        /// <param name="result">The vector where the result will be placed</param>
         /// <exception cref="ArgumentException">If this vector and <paramref name="result"/> are not the same size.</exception>
         protected void PointwiseBinary(Action<T, Vector<T>> f, T x, Vector<T> result)
         {

src/Numerics/Numerics.csproj

@@ -44,7 +44,7 @@
     <WarningLevel>4</WarningLevel>
     <CodeAnalysisRuleSet>AllRules.ruleset</CodeAnalysisRuleSet>
     <NoWarn>1591</NoWarn>
-    <LangVersion>5</LangVersion>
+    <LangVersion>6</LangVersion>
   </PropertyGroup>
   <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Debug|AnyCPU' ">
     <OutputPath>..\..\out\lib-debug\Net40\</OutputPath>
@@ -58,7 +58,7 @@
     <ErrorReport>prompt</ErrorReport>
     <WarningLevel>4</WarningLevel>
     <CodeAnalysisRuleSet>AllRules.ruleset</CodeAnalysisRuleSet>
-    <LangVersion>5</LangVersion>
+    <LangVersion>6</LangVersion>
   </PropertyGroup>
   <PropertyGroup Condition="'$(Configuration)|$(Platform)' == 'Release-Signed|AnyCPU'">
     <OutputPath>..\..\out\lib-signed\Net40\</OutputPath>
@@ -75,7 +75,7 @@
     <WarningLevel>4</WarningLevel>
     <CodeAnalysisRuleSet>AllRules.ruleset</CodeAnalysisRuleSet>
     <NoWarn>1591</NoWarn>
-    <LangVersion>5</LangVersion>
+    <LangVersion>6</LangVersion>
   </PropertyGroup>
   <ItemGroup>
     <Reference Include="System" />
@@ -115,6 +115,24 @@
     <Compile Include="LinearRegression\Options.cs" />
     <Compile Include="OdeSolvers\AdamsBashforth.cs" />
     <Compile Include="OdeSolvers\RungeKutta.cs" />
+    <Compile Include="Optimization\BfgsMinimizerBase.cs" />
+    <Compile Include="Optimization\LineSearch\WolfeLineSearch.cs" />
+    <Compile Include="Optimization\NelderMeadSimplex.cs" />
+    <Compile Include="Optimization\ObjectiveFunctions\ForwardDifferenceGradientObjectiveFunction.cs" />
+    <Compile Include="Optimization\ObjectiveFunctions\LazyObjectiveFunctionBase.cs" />
+    <Compile Include="Optimization\Exceptions.cs" />
+    <Compile Include="Optimization\ObjectiveFunctions\ObjectiveFunctionBase.cs" />
+    <Compile Include="Optimization\ObjectiveFunctions\LazyObjectiveFunction.cs" />
+    <Compile Include="Optimization\ObjectiveFunctions\ValueObjectiveFunction.cs" />
+    <Compile Include="Optimization\ObjectiveFunctions\HessianObjectiveFunction.cs" />
+    <Compile Include="Optimization\ObjectiveFunctions\GradientObjectiveFunction.cs" />
+    <Compile Include="Optimization\ObjectiveFunctions\GradientHessianObjectiveFunction.cs" />
+    <Compile Include="Optimization\IObjectiveFunction.cs" />
+    <Compile Include="Optimization\LineSearch\LineSearchResult.cs" />
+    <Compile Include="Optimization\LineSearch\WeakWolfeLineSearch.cs" />
+    <Compile Include="Optimization\IUnconstrainedMinimizer.cs" />
+    <Compile Include="Optimization\MinimizationResult.cs" />
+    <Compile Include="Optimization\MinimizationWithLineSearchResult.cs" />
     <Compile Include="Precision.Comparison.cs" />
     <Compile Include="Precision.Equality.cs" />
     <Compile Include="Distributions\Bernoulli.cs" />
@@ -233,6 +251,7 @@
     <Compile Include="Providers\Common\NativeProviderLoader.cs" />
     <Compile Include="Random\SystemRandomSource.cs" />
     <Compile Include="Random\RandomSeed.cs" />
+    <Compile Include="Optimization\BfgsSolver.cs" />
     <Compile Include="RootFinding\Broyden.cs" />
     <Compile Include="RootFinding\Cubic.cs" />
     <Compile Include="RootFinding\NewtonRaphson.cs" />
@@ -242,6 +261,17 @@
     <Compile Include="RootFinding\Brent.cs" />
     <Compile Include="FindRoots.cs" />
     <Compile Include="RootFinding\Bisection.cs" />
+    <Compile Include="Optimization\BfgsBMinimizer.cs" />
+    <Compile Include="Optimization\BfgsMinimizer.cs" />
+    <Compile Include="Optimization\ConjugateGradientMinimizer.cs" />
+    <Compile Include="Optimization\GoldenSectionMinimizer.cs" />
+    <Compile Include="Optimization\MinimizationResult1D.cs" />
+    <Compile Include="Optimization\NewtonMinimizer.cs" />
+    <Compile Include="Optimization\ObjectiveFunction.cs" />
+    <Compile Include="Optimization\ObjectiveFunction1D.cs" />
+    <Compile Include="Optimization\OptimizationResult.cs" />
+    <Compile Include="Optimization\LineSearch\StrongWolfeLineSearch.cs" />
+    <Compile Include="Optimization\QuadraticGradientProjectionSearch.cs" />
     <Compile Include="SpecialFunctions\Evaluate.cs" />
     <Compile Include="ExcelFunctions.cs" />
     <Compile Include="SpecialFunctions\ExponentialIntegral.cs" />

src/Numerics/Optimization/BfgsBMinimizer.cs

@@ -0,0 +1,305 @@
+// <copyright file="BfgsTest.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2016 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using System;
+using System.Collections.Generic;
+using MathNet.Numerics.LinearAlgebra;
+using MathNet.Numerics.LinearAlgebra.Double;
+using MathNet.Numerics.Optimization.LineSearch;
+
+namespace MathNet.Numerics.Optimization
+{
+    /// <summary>
+    /// Broyden–Fletcher–Goldfarb–Shanno Bounded (BFGS-B) algorithm is an iterative method for solving box-constrained nonlinear optimization problems
+    /// http://www.ece.northwestern.edu/~nocedal/PSfiles/limited.ps.gz
+    /// </summary>
+    public class BfgsBMinimizer : BfgsMinimizerBase
+    {
+        public BfgsBMinimizer(double gradientTolerance, double parameterTolerance, double functionProgressTolerance, int maximumIterations = 1000)
+            : base(gradientTolerance,parameterTolerance,functionProgressTolerance,maximumIterations)
+        {
+        }
+
+        /// <summary>
+        /// Find the minimum of the objective function given lower and upper bounds
+        /// </summary>
+        /// <param name="objective">The objective function, must support a gradient</param>
+        /// <param name="lowerBound">The lower bound</param>
+        /// <param name="upperBound">The upper bound</param>
+        /// <param name="initialGuess">The initial guess</param>
+        /// <returns>The MinimizationResult which contains the minimum and the ExitCondition</returns>
+        public MinimizationResult FindMinimum(IObjectiveFunction objective, Vector<double> lowerBound, Vector<double> upperBound, Vector<double> initialGuess)
+        {
+            _lowerBound = lowerBound;
+            _upperBound = upperBound;
+            if (!objective.IsGradientSupported)
+                throw new IncompatibleObjectiveException("Gradient not supported in objective function, but required for BFGS minimization.");
+
+            // Check that dimensions match
+            if (lowerBound.Count != upperBound.Count || lowerBound.Count != initialGuess.Count)
+                throw new ArgumentException("Dimensions of bounds and/or initial guess do not match.");
+
+            // Check that initial guess is feasible
+            for (int ii = 0; ii < initialGuess.Count; ++ii)
+                if (initialGuess[ii] < lowerBound[ii] || initialGuess[ii] > upperBound[ii])
+                    throw new ArgumentException("Initial guess is not in the feasible region");
+
+            objective.EvaluateAt(initialGuess);
+            ValidateGradientAndObjective(objective);
+
+            // Check that we're not already done
+            MinimizationResult.ExitCondition currentExitCondition = ExitCriteriaSatisfied(objective, null, 0);
+            if (currentExitCondition != MinimizationResult.ExitCondition.None)
+                return new MinimizationResult(objective, 0, currentExitCondition);
+
+            // Set up line search algorithm
+            var lineSearcher = new StrongWolfeLineSearch(1e-4, 0.9, Math.Max(ParameterTolerance, 1e-5), maxIterations: 1000);
+
+            // Declare state variables
+            Vector<double> reducedSolution1, reducedGradient, reducedInitialPoint, reducedCauchyPoint, solution1;
+            Matrix<double> reducedHessian;
+            List<int> reducedMap;
+
+            // First step
+            var pseudoHessian = CreateMatrix.DiagonalIdentity<double>(initialGuess.Count);
+
+            // Determine active set
+            var gradientProjectionResult = QuadraticGradientProjectionSearch.Search(objective.Point, objective.Gradient, pseudoHessian, lowerBound, upperBound);
+            var cauchyPoint = gradientProjectionResult.CauchyPoint;
+            var fixedCount = gradientProjectionResult.FixedCount;
+            var isFixed = gradientProjectionResult.IsFixed;
+            var freeCount = lowerBound.Count - fixedCount;
+
+            if (freeCount > 0)
+            {
+                reducedGradient = new DenseVector(freeCount);
+                reducedHessian = new DenseMatrix(freeCount, freeCount);
+                reducedMap = new List<int>(freeCount);
+                reducedInitialPoint = new DenseVector(freeCount);
+                reducedCauchyPoint = new DenseVector(freeCount);
+
+                CreateReducedData(objective.Point, cauchyPoint, isFixed, lowerBound, upperBound, objective.Gradient, pseudoHessian, reducedInitialPoint, reducedCauchyPoint, reducedGradient, reducedHessian, reducedMap);
+
+                // Determine search direction and maximum step size
+                reducedSolution1 = reducedInitialPoint + reducedHessian.Cholesky().Solve(-reducedGradient);
+
+                solution1 = ReducedToFull(reducedMap, reducedSolution1, cauchyPoint);
+            }
+            else
+            {
+                solution1 = cauchyPoint;
+            }
+
+            var directionFromCauchy = solution1 - cauchyPoint;
+            var maxStepFromCauchyPoint = FindMaxStep(cauchyPoint, directionFromCauchy, lowerBound, upperBound);
+
+            var solution2 = cauchyPoint + Math.Min(maxStepFromCauchyPoint, 1.0)*directionFromCauchy;
+
+            var lineSearchDirection = solution2 - objective.Point;
+            var maxLineSearchStep = FindMaxStep(objective.Point, lineSearchDirection, lowerBound, upperBound);
+            var estStepSize = -objective.Gradient*lineSearchDirection/(lineSearchDirection*pseudoHessian*lineSearchDirection);
+
+            var startingStepSize = Math.Min(Math.Max(estStepSize, 1.0), maxLineSearchStep);
+
+            // Line search
+            LineSearchResult lineSearchResult;
+            try
+            {
+                lineSearchResult = lineSearcher.FindConformingStep(objective, lineSearchDirection, startingStepSize, upperBound: maxLineSearchStep);
+            }
+            catch (Exception e)
+            {
+                throw new InnerOptimizationException("Line search failed.", e);
+            }
+
+            var previousPoint = objective.Fork();
+            var candidatePoint = lineSearchResult.FunctionInfoAtMinimum;
+            var gradient = candidatePoint.Gradient;
+            var step = candidatePoint.Point - initialGuess;
+
+            // Subsequent steps
+            int totalLineSearchSteps = lineSearchResult.Iterations;
+            int iterationsWithNontrivialLineSearch = lineSearchResult.Iterations > 0 ? 0 : 1;
+
+            int iterations = DoBfgsUpdate(ref currentExitCondition, lineSearcher, ref pseudoHessian, ref lineSearchDirection, ref previousPoint, ref lineSearchResult, ref candidatePoint, ref step, ref totalLineSearchSteps, ref iterationsWithNontrivialLineSearch);
+
+            if (iterations == MaximumIterations && currentExitCondition == MinimizationResult.ExitCondition.None)
+                throw new MaximumIterationsException(string.Format("Maximum iterations ({0}) reached.", MaximumIterations));
+
+            return new MinimizationWithLineSearchResult(candidatePoint, iterations, currentExitCondition, totalLineSearchSteps, iterationsWithNontrivialLineSearch);
+        }
+
+        protected override Vector<double> CalculateSearchDirection(ref Matrix<double> pseudoHessian, 
+            out double maxLineSearchStep,
+            out double startingStepSize, 
+            IObjectiveFunction previousPoint, 
+            IObjectiveFunction candidatePoint, 
+            Vector<double> step)
+        {
+            Vector<double> lineSearchDirection;
+            var y = candidatePoint.Gradient - previousPoint.Gradient;
+
+            double sy = step * y;
+            if (sy > 0.0) // only do update if it will create a positive definite matrix
+            {
+                double sts = step * step;
+                
+                var Hs = pseudoHessian * step;
+                var sHs = step * pseudoHessian * step;
+                pseudoHessian = pseudoHessian + y.OuterProduct(y) * (1.0 / sy) - Hs.OuterProduct(Hs) * (1.0 / sHs);
+            }
+            else
+            {
+                //pseudo_hessian = LinearAlgebra.Double.DiagonalMatrix.Identity(initial_guess.Count);
+            }
+
+            // Determine active set
+            var gradientProjectionResult = QuadraticGradientProjectionSearch.Search(candidatePoint.Point, candidatePoint.Gradient, pseudoHessian, _lowerBound, _upperBound);
+            var cauchyPoint = gradientProjectionResult.CauchyPoint;
+            var fixedCount = gradientProjectionResult.FixedCount;
+            var isFixed = gradientProjectionResult.IsFixed;
+            var freeCount = _lowerBound.Count - fixedCount;
+            Vector<double> solution1;
+            if (freeCount > 0)
+            {
+                var reducedGradient = new DenseVector(freeCount);
+                var reducedHessian = new DenseMatrix(freeCount, freeCount);
+                var reducedMap = new List<int>(freeCount);
+                var reducedInitialPoint = new DenseVector(freeCount);
+                var reducedCauchyPoint = new DenseVector(freeCount);
+
+                CreateReducedData(candidatePoint.Point, cauchyPoint, isFixed, _lowerBound, _upperBound, candidatePoint.Gradient, pseudoHessian, reducedInitialPoint, reducedCauchyPoint, reducedGradient, reducedHessian, reducedMap);
+
+                // Determine search direction and maximum step size
+                Vector<double> reducedSolution1 = reducedInitialPoint + reducedHessian.Cholesky().Solve(-reducedGradient);
+
+                solution1 = ReducedToFull(reducedMap, reducedSolution1, cauchyPoint);
+            }
+            else
+            {
+                solution1 = cauchyPoint;
+            }
+
+            var directionFromCauchy = solution1 - cauchyPoint;
+            var maxStepFromCauchyPoint = FindMaxStep(cauchyPoint, directionFromCauchy, _lowerBound, _upperBound);
+
+            var solution2 = cauchyPoint + Math.Min(maxStepFromCauchyPoint, 1.0) * directionFromCauchy;
+
+            lineSearchDirection = solution2 - candidatePoint.Point;
+            maxLineSearchStep = FindMaxStep(candidatePoint.Point, lineSearchDirection, _lowerBound, _upperBound);
+
+            if (maxLineSearchStep == 0.0)
+            {
+                lineSearchDirection = cauchyPoint - candidatePoint.Point;
+                maxLineSearchStep = FindMaxStep(candidatePoint.Point, lineSearchDirection, _lowerBound, _upperBound);
+            }
+
+            double estStepSize = -candidatePoint.Gradient * lineSearchDirection / (lineSearchDirection * pseudoHessian * lineSearchDirection);
+
+            startingStepSize = Math.Min(Math.Max(estStepSize, 1.0), maxLineSearchStep);
+            return lineSearchDirection;
+        }
+
+        private static Vector<double> ReducedToFull(List<int> reducedMap, Vector<double> reducedVector, Vector<double> fullVector)
+        {
+            var output = fullVector.Clone();
+            for (int ii = 0; ii < reducedMap.Count; ++ii)
+                output[reducedMap[ii]] = reducedVector[ii];
+            return output;
+        }
+
+        private Vector<double> _lowerBound;
+        private Vector<double> _upperBound;
+
+        private static double FindMaxStep(Vector<double> startingPoint, Vector<double> searchDirection, Vector<double> lowerBound, Vector<double> upperBound)
+        {
+            double maxStep = Double.PositiveInfinity;
+            for (int ii = 0; ii < startingPoint.Count; ++ii)
+            {
+                double paramMaxStep;
+                if (searchDirection[ii] > 0)
+                    paramMaxStep = (upperBound[ii] - startingPoint[ii])/searchDirection[ii];
+                else if (searchDirection[ii] < 0)
+                    paramMaxStep = (startingPoint[ii] - lowerBound[ii])/-searchDirection[ii];
+                else
+                    paramMaxStep = Double.PositiveInfinity;
+
+                if (paramMaxStep < maxStep)
+                    maxStep = paramMaxStep;
+            }
+            return maxStep;
+        }
+
+        private static void CreateReducedData(Vector<double> initialPoint, Vector<double> cauchyPoint, List<bool> isFixed, Vector<double> lowerBound, Vector<double> upperBound, Vector<double> gradient, Matrix<double> pseudoHessian, Vector<double> reducedInitialPoint, Vector<double> reducedCauchyPoint, Vector<double> reducedGradient, Matrix<double> reducedHessian, List<int> reducedMap)
+        {
+            int ll = 0;
+            for (int ii = 0; ii < lowerBound.Count; ++ii)
+            {
+                if (!isFixed[ii])
+                {
+                    // hessian
+                    int mm = 0;
+                    for (int jj = 0; jj < lowerBound.Count; ++jj)
+                    {
+                        if (!isFixed[jj])
+                        {
+                            reducedHessian[ll, mm++] = pseudoHessian[ii, jj];
+                        }
+                    }
+
+                    // gradient
+                    reducedInitialPoint[ll] = initialPoint[ii];
+                    reducedCauchyPoint[ll] = cauchyPoint[ii];
+                    reducedGradient[ll] = gradient[ii];
+                    ll += 1;
+                    reducedMap.Add(ii);
+
+                }
+            }
+        }
+
+        protected override double GetProjectedGradient(IObjectiveFunction candidatePoint, int ii)
+        {
+            double projectedGradient;
+            bool atLowerBound = candidatePoint.Point[ii] - _lowerBound[ii] < VerySmall;
+            bool atUpperBound = _upperBound[ii] - candidatePoint.Point[ii] < VerySmall;
+
+            if (atLowerBound && atUpperBound)
+                projectedGradient = 0.0;
+            else if (atLowerBound)
+                projectedGradient = Math.Min(candidatePoint.Gradient[ii], 0.0);
+            else if (atUpperBound)
+                projectedGradient = Math.Max(candidatePoint.Gradient[ii], 0.0);
+            else
+                projectedGradient = base.GetProjectedGradient(candidatePoint, ii);
+            return projectedGradient;
+        }
+    }
+}

src/Numerics/Optimization/BfgsMinimizer.cs

@@ -0,0 +1,142 @@
+// <copyright file="BfgsTest.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2016 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using System;
+using MathNet.Numerics.LinearAlgebra;
+using MathNet.Numerics.Optimization.LineSearch;
+
+namespace MathNet.Numerics.Optimization
+{
+    /// <summary>
+    /// Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm is an iterative method for solving unconstrained nonlinear optimization problems
+    /// </summary>
+    public class BfgsMinimizer : BfgsMinimizerBase
+    {
+        /// <summary>
+        /// Creates BFGS minimizer
+        /// </summary>
+        /// <param name="gradientTolerance">The gradient tolerance</param>
+        /// <param name="parameterTolerance">The parameter tolerance</param>
+        /// <param name="functionProgressTolerance">The funciton progress tolerance</param>
+        /// <param name="maximumIterations">The maximum number of iterations</param>
+        public BfgsMinimizer(double gradientTolerance, double parameterTolerance, double functionProgressTolerance, int maximumIterations=1000)
+            :base(gradientTolerance,parameterTolerance,functionProgressTolerance,maximumIterations)
+        {
+        }
+
+        /// <summary>
+        /// Find the minimum of the objective function given lower and upper bounds
+        /// </summary>
+        /// <param name="objective">The objective function, must support a gradient</param>
+        /// <param name="initialGuess">The initial guess</param>
+        /// <returns>The MinimizationResult which contains the minimum and the ExitCondition</returns>
+        public MinimizationResult FindMinimum(IObjectiveFunction objective, Vector<double> initialGuess)
+        {
+            if (!objective.IsGradientSupported)
+                throw new IncompatibleObjectiveException("Gradient not supported in objective function, but required for BFGS minimization.");
+
+            objective.EvaluateAt(initialGuess);
+            ValidateGradientAndObjective(objective);
+
+            // Check that we're not already done
+            MinimizationResult.ExitCondition currentExitCondition = ExitCriteriaSatisfied(objective, null, 0);
+            if (currentExitCondition != MinimizationResult.ExitCondition.None)
+                return new MinimizationResult(objective, 0, currentExitCondition);
+
+            // Set up line search algorithm
+            var lineSearcher = new WeakWolfeLineSearch(1e-4, 0.9, Math.Max(ParameterTolerance, 1e-10), 1000);
+
+            // First step
+            var inversePseudoHessian = CreateMatrix.DenseIdentity<double>(initialGuess.Count);
+            var lineSearchDirection = -objective.Gradient;
+            var stepSize = 100 * GradientTolerance / (lineSearchDirection * lineSearchDirection);
+
+            var previousPoint = objective;
+
+            LineSearchResult lineSearchResult;
+            try
+            {
+                lineSearchResult = lineSearcher.FindConformingStep(objective, lineSearchDirection, stepSize);
+            }
+            catch (OptimizationException e)
+            {
+                throw new InnerOptimizationException("Line search failed.", e);
+            }
+            catch (ArgumentException e)
+            {
+                throw new InnerOptimizationException("Line search failed.", e);
+            }
+
+            var candidate = lineSearchResult.FunctionInfoAtMinimum;
+            ValidateGradientAndObjective(candidate);
+
+            var gradient = candidate.Gradient;
+            var step = candidate.Point - initialGuess;
+
+            // Subsequent steps
+            Matrix<double> I = CreateMatrix.DiagonalIdentity<double>(initialGuess.Count);
+            int iterations;
+            int totalLineSearchSteps = lineSearchResult.Iterations;
+            int iterationsWithNontrivialLineSearch = lineSearchResult.Iterations > 0 ? 0 : 1;
+            iterations = DoBfgsUpdate(ref currentExitCondition, lineSearcher, ref inversePseudoHessian, ref lineSearchDirection, ref previousPoint, ref lineSearchResult, ref candidate, ref step, ref totalLineSearchSteps, ref iterationsWithNontrivialLineSearch);
+
+            if (iterations == MaximumIterations && currentExitCondition == MinimizationResult.ExitCondition.None)
+                throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", MaximumIterations));
+
+            return new MinimizationWithLineSearchResult(candidate, iterations, MinimizationResult.ExitCondition.AbsoluteGradient, totalLineSearchSteps, iterationsWithNontrivialLineSearch);
+        }
+
+        protected override Vector<double> CalculateSearchDirection(ref Matrix<double> inversePseudoHessian,
+            out double maxLineSearchStep,
+            out double startingStepSize,
+            IObjectiveFunction previousPoint,
+            IObjectiveFunction candidate, 
+            Vector<double> step)
+        {
+            startingStepSize = 1.0;
+            maxLineSearchStep = double.PositiveInfinity;
+
+            Vector<double> lineSearchDirection;
+            var y = candidate.Gradient - previousPoint.Gradient;
+
+            double sy = step * y;
+            inversePseudoHessian = inversePseudoHessian + ((sy + y * inversePseudoHessian * y) / Math.Pow(sy, 2.0)) * step.OuterProduct(step) - ((inversePseudoHessian * y.ToColumnMatrix()) * step.ToRowMatrix() + step.ToColumnMatrix() * (y.ToRowMatrix() * inversePseudoHessian)) * (1.0 / sy);
+            lineSearchDirection = -inversePseudoHessian * candidate.Gradient;
+
+            if (lineSearchDirection * candidate.Gradient >= 0.0)
+            {
+                lineSearchDirection = -candidate.Gradient;
+                inversePseudoHessian = CreateMatrix.DenseIdentity<double>(candidate.Point.Count);
+            }
+
+            return lineSearchDirection;
+        }
+    }
+}

src/Numerics/Optimization/BfgsMinimizerBase.cs

@@ -0,0 +1,158 @@
+// <copyright file="BfgsTest.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2016 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using MathNet.Numerics.LinearAlgebra;
+using MathNet.Numerics.LinearAlgebra.Double;
+using MathNet.Numerics.Optimization.LineSearch;
+using System;
+
+namespace MathNet.Numerics.Optimization
+{
+    public abstract class BfgsMinimizerBase
+    {
+        public double GradientTolerance { get; set; }
+        public double ParameterTolerance { get; set; }
+        public double FunctionProgressTolerance { get; set; }
+        public int MaximumIterations { get; set; }
+
+        protected const double VerySmall = 1e-15;
+
+        /// <summary>
+        /// Creates a base class for BFGS minimization
+        /// </summary>
+        /// <param name="gradientTolerance">The gradient tolerance</param>
+        /// <param name="parameterTolerance">The parameter tolerance</param>
+        /// <param name="functionProgressTolerance">The funciton progress tolerance</param>
+        /// <param name="maximumIterations">The maximum number of iterations</param>
+        public BfgsMinimizerBase(double gradientTolerance, double parameterTolerance, double functionProgressTolerance, int maximumIterations)
+        {
+            GradientTolerance = gradientTolerance;
+            ParameterTolerance = parameterTolerance;
+            FunctionProgressTolerance = functionProgressTolerance;
+            MaximumIterations = maximumIterations;
+        }
+
+        protected MinimizationResult.ExitCondition ExitCriteriaSatisfied(IObjectiveFunction candidatePoint, IObjectiveFunction lastPoint, int iterations)
+        {
+            Vector<double> relGrad = new DenseVector(candidatePoint.Point.Count);
+            double relativeGradient = 0.0;
+            double normalizer = Math.Max(Math.Abs(candidatePoint.Value), 1.0);
+            for (int ii = 0; ii < relGrad.Count; ++ii)
+            {
+                double projectedGradient = GetProjectedGradient(candidatePoint, ii);
+
+                double tmp = projectedGradient * 
+                    Math.Max(Math.Abs(candidatePoint.Point[ii]), 1.0) / normalizer;
+                relativeGradient = Math.Max(relativeGradient, Math.Abs(tmp));
+            }
+            if (relativeGradient < GradientTolerance)
+            {
+                return MinimizationResult.ExitCondition.RelativeGradient;
+            }
+
+            if (lastPoint != null)
+            {
+                double mostProgress = 0.0;
+                for (int ii = 0; ii < candidatePoint.Point.Count; ++ii)
+                {
+                    var tmp = Math.Abs(candidatePoint.Point[ii] - lastPoint.Point[ii]) / 
+                        Math.Max(Math.Abs(lastPoint.Point[ii]), 1.0);
+                    mostProgress = Math.Max(mostProgress, tmp);
+                }
+                if (mostProgress < ParameterTolerance)
+                {
+                    return MinimizationResult.ExitCondition.LackOfProgress;
+                }
+
+                double functionChange = candidatePoint.Value - lastPoint.Value;
+                if (iterations > 500 && functionChange < 0 && Math.Abs(functionChange) < FunctionProgressTolerance)
+                    return MinimizationResult.ExitCondition.LackOfProgress;
+            }
+
+            return MinimizationResult.ExitCondition.None;
+        }
+
+        protected virtual double GetProjectedGradient(IObjectiveFunction candidatePoint, int ii)
+        {
+            return candidatePoint.Gradient[ii];
+        }
+
+        protected void ValidateGradientAndObjective(IObjectiveFunction eval)
+        {
+            foreach (var x in eval.Gradient)
+            {
+                if (Double.IsNaN(x) || Double.IsInfinity(x))
+                    throw new EvaluationException("Non-finite gradient returned.", eval);
+            }
+            if (Double.IsNaN(eval.Value) || Double.IsInfinity(eval.Value))
+                throw new EvaluationException("Non-finite objective function returned.", eval);
+        }
+
+        protected int DoBfgsUpdate(ref MinimizationResult.ExitCondition currentExitCondition, WolfeLineSearch lineSearcher, ref Matrix<double> inversePseudoHessian, ref Vector<double> lineSearchDirection, ref IObjectiveFunction previousPoint, ref LineSearchResult lineSearchResult, ref IObjectiveFunction candidate, ref Vector<double> step, ref int totalLineSearchSteps, ref int iterationsWithNontrivialLineSearch)
+        {
+            int iterations;
+            for (iterations = 1; iterations < MaximumIterations; ++iterations)
+            {
+                double startingStepSize;
+                double maxLineSearchStep;
+                lineSearchDirection = CalculateSearchDirection(ref inversePseudoHessian, out maxLineSearchStep, out startingStepSize, previousPoint, candidate, step);
+
+                try
+                {
+                    lineSearchResult = lineSearcher.FindConformingStep(candidate, lineSearchDirection, startingStepSize, maxLineSearchStep);
+                }
+                catch (Exception e)
+                {
+                    throw new InnerOptimizationException("Line search failed.", e);
+                }
+
+                iterationsWithNontrivialLineSearch += lineSearchResult.Iterations > 0 ? 1 : 0;
+                totalLineSearchSteps += lineSearchResult.Iterations;
+
+                step = lineSearchResult.FunctionInfoAtMinimum.Point - candidate.Point;
+                previousPoint = candidate;
+                candidate = lineSearchResult.FunctionInfoAtMinimum;
+
+                currentExitCondition = ExitCriteriaSatisfied(candidate, previousPoint, iterations);
+                if (currentExitCondition != MinimizationResult.ExitCondition.None)
+                    break;
+            }
+
+            return iterations;
+        }
+
+        protected abstract Vector<double> CalculateSearchDirection(ref Matrix<double> inversePseudoHessian,
+            out double maxLineSearchStep,
+            out double startingStepSize,
+            IObjectiveFunction previousPoint,
+            IObjectiveFunction candidate,
+            Vector<double> step);
+    }
+}

src/Numerics/Optimization/BfgsSolver.cs

@@ -0,0 +1,110 @@
+// <copyright file="BfgsSolver.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2015 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using System;
+using MathNet.Numerics.LinearAlgebra;
+using MathNet.Numerics.LinearAlgebra.Double;
+using MathNet.Numerics.Optimization.LineSearch;
+
+namespace MathNet.Numerics.Optimization
+{
+    /// <summary>
+    /// Broyden-Fletcher-Goldfarb-Shanno solver for finding function minima
+    /// See http://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm
+    /// Inspired by implementation: https://github.com/PatWie/CppNumericalSolvers/blob/master/src/BfgsSolver.cpp
+    /// </summary>
+    public static class BfgsSolver
+    {
+        private const double GradientTolerance = 1e-5;
+        private const int MaxIterations = 100000;
+
+        /// <summary>
+        /// Finds a minimum of a function by the BFGS quasi-Newton method
+        /// This uses the function and it's gradient (partial derivatives in each direction) and approximates the Hessian
+        /// </summary>
+        /// <param name="initialGuess">An initial guess</param>
+        /// <param name="functionValue">Evaluates the function at a point</param>
+        /// <param name="functionGradient">Evaluates the gradient of the function at a point</param>
+        /// <returns>The minimum found</returns>
+        public static Vector<double> Solve(Vector initialGuess, Func<Vector<double>, double> functionValue, Func<Vector<double>, Vector<double>> functionGradient)
+        {
+            var objectiveFunction = ObjectiveFunction.Gradient(functionValue, functionGradient);
+            objectiveFunction.EvaluateAt(initialGuess);
+
+            int dim = initialGuess.Count;
+            int iter = 0;
+            // H represents the approximation of the inverse hessian matrix
+            // it is updated via the Sherman–Morrison formula (http://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula)
+            Matrix<double> H = DenseMatrix.CreateIdentity(dim);
+
+            Vector<double> x = initialGuess;
+            Vector<double> x_old = x;
+            Vector<double> grad;
+            WolfeLineSearch wolfeLineSearch = new WeakWolfeLineSearch(1e-4, 0.9, 1e-5, 200);
+            do
+            {
+                // search along the direction of the gradient
+                grad = objectiveFunction.Gradient;
+                Vector<double> p = -1 * H * grad;
+                var lineSearchResult = wolfeLineSearch.FindConformingStep(objectiveFunction, p, 1.0);
+                double rate = lineSearchResult.FinalStep;
+                x = x + rate * p;
+                Vector<double> grad_old = grad;
+
+                // update the gradient
+                objectiveFunction.EvaluateAt(x);
+                grad = objectiveFunction.Gradient;// functionGradient(x);
+
+                Vector<double> s = x - x_old;
+                Vector<double> y = grad - grad_old;
+
+                double rho = 1.0 / (y * s);
+                if (iter == 0)
+                {
+                    // set up an initial hessian
+                    H = (y * s) / (y * y) * DenseMatrix.CreateIdentity(dim);
+                }
+
+                var sM = s.ToColumnMatrix();
+                var yM = y.ToColumnMatrix();
+
+                // Update the estimate of the hessian
+                H = H
+                    - rho * (sM * (yM.TransposeThisAndMultiply(H)) + (H * yM).TransposeAndMultiply(sM))
+                    + rho * rho * (y.DotProduct(H * y) + 1.0 / rho) * (sM.TransposeAndMultiply(sM));
+                x_old = x;
+                iter++;
+            }
+            while ((grad.InfinityNorm() > GradientTolerance) && (iter < MaxIterations));
+
+            return x;
+        }
+    }
+}

src/Numerics/Optimization/ConjugateGradientMinimizer.cs

@@ -0,0 +1,115 @@
+using System;
+using MathNet.Numerics.LinearAlgebra;
+using MathNet.Numerics.Optimization.LineSearch;
+
+namespace MathNet.Numerics.Optimization
+{
+    public class ConjugateGradientMinimizer
+    {
+        public double GradientTolerance { get; set; }
+        public int MaximumIterations { get; set; }
+
+        public ConjugateGradientMinimizer(double gradientTolerance, int maximumIterations)
+        {
+            GradientTolerance = gradientTolerance;
+            MaximumIterations = maximumIterations;
+        }
+
+        public MinimizationResult FindMinimum(IObjectiveFunction objective, Vector<double> initialGuess)
+        {
+            if (!objective.IsGradientSupported)
+                throw new IncompatibleObjectiveException("Gradient not supported in objective function, but required for ConjugateGradient minimization.");
+
+            objective.EvaluateAt(initialGuess);
+            var gradient = objective.Gradient;
+            ValidateGradient(objective);
+
+            // Check that we're not already done
+            if (ExitCriteriaSatisfied(initialGuess, gradient))
+                return new MinimizationResult(objective, 0, MinimizationResult.ExitCondition.AbsoluteGradient);
+
+            // Set up line search algorithm
+            var lineSearcher = new WeakWolfeLineSearch(1e-4, 0.1, 1e-4, 1000);
+
+            // First step
+            var steepestDirection = -gradient;
+            var searchDirection = steepestDirection;
+            double initialStepSize = 100 * GradientTolerance / (gradient * gradient);
+
+            LineSearchResult result;
+            try
+            {
+                result = lineSearcher.FindConformingStep(objective, searchDirection, initialStepSize);
+            }
+            catch (Exception e)
+            {
+                throw new InnerOptimizationException("Line search failed.", e);
+            }
+
+            objective = result.FunctionInfoAtMinimum;
+            ValidateGradient(objective);
+
+            double stepSize = result.FinalStep;
+
+            // Subsequent steps
+            int iterations = 1;
+            int totalLineSearchSteps = result.Iterations;
+            int iterationsWithNontrivialLineSearch = result.Iterations > 0 ? 0 : 1;
+            int steepestDescentResets = 0;
+            while (!ExitCriteriaSatisfied(objective.Point, objective.Gradient) && iterations < MaximumIterations)
+            {
+                var previousSteepestDirection = steepestDirection;
+                steepestDirection = -objective.Gradient;
+                var searchDirectionAdjuster = Math.Max(0, steepestDirection*(steepestDirection - previousSteepestDirection)/(previousSteepestDirection*previousSteepestDirection));
+                searchDirection = steepestDirection + searchDirectionAdjuster * searchDirection;
+                if (searchDirection * objective.Gradient >= 0)
+                {
+                    searchDirection = steepestDirection;
+                    steepestDescentResets += 1;
+                }
+
+                try
+                {
+                    result = lineSearcher.FindConformingStep(objective, searchDirection, stepSize);
+                }
+                catch (Exception e)
+                {
+                    throw new InnerOptimizationException("Line search failed.", e);
+                }
+
+                iterationsWithNontrivialLineSearch += result.Iterations == 0 ? 1 : 0;
+                totalLineSearchSteps += result.Iterations;
+				stepSize = result.FinalStep;
+                objective = result.FunctionInfoAtMinimum;
+                iterations += 1;
+            }
+
+            if (iterations == MaximumIterations)
+            {
+                throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", MaximumIterations));
+            }
+
+            return new MinimizationWithLineSearchResult(objective, iterations, MinimizationResult.ExitCondition.AbsoluteGradient, totalLineSearchSteps, iterationsWithNontrivialLineSearch);
+        }
+
+        bool ExitCriteriaSatisfied(Vector<double> candidatePoint, Vector<double> gradient)
+        {
+            return gradient.Norm(2.0) < GradientTolerance;
+        }
+
+        void ValidateGradient(IObjectiveFunction objective)
+        {
+            foreach (var x in objective.Gradient)
+            {
+                if (Double.IsNaN(x) || Double.IsInfinity(x))
+                    throw new EvaluationException("Non-finite gradient returned.", objective);
+            }
+        }
+
+        void ValidateObjective(IObjectiveFunction objective)
+        {
+            if (Double.IsNaN(objective.Value) || Double.IsInfinity(objective.Value))
+                throw new EvaluationException("Non-finite objective function returned.", objective);
+        }
+    }
+}

src/Numerics/Optimization/Exceptions.cs

@@ -0,0 +1,52 @@
+using System;
+
+namespace MathNet.Numerics.Optimization
+{
+    public class OptimizationException : Exception
+    {
+        public OptimizationException(string message)
+            : base(message) { }
+
+        public OptimizationException(string message, Exception innerException)
+            : base(message, innerException) { }
+    }
+
+    public class MaximumIterationsException : OptimizationException
+    {
+        public MaximumIterationsException(string message)
+            : base(message) { }
+    }
+
+    public class EvaluationException : OptimizationException
+    {
+        public IObjectiveFunction ObjectiveFunction { get; private set; }
+
+        public EvaluationException(string message, IObjectiveFunction eval)
+            : base(message)
+        {
+            ObjectiveFunction = eval;
+        }
+
+        public EvaluationException(string message, IObjectiveFunction eval, Exception innerException)
+            : base(message, innerException)
+        {
+            ObjectiveFunction = 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/GoldenSectionMinimizer.cs

@@ -0,0 +1,107 @@
+using System;
+
+namespace MathNet.Numerics.Optimization
+{
+    public class GoldenSectionMinimizer
+    {
+        public double XTolerance { get; set; }
+        public int MaximumIterations { get; set; }
+        public int MaximumExpansionSteps { get; set; }
+        public double LowerExpansionFactor { get; set; }
+        public double UpperExpansionFactor { get; set; }
+
+        public GoldenSectionMinimizer(double xTolerance = 1e-5, int maxIterations = 1000, int maxExpansionSteps = 10, double lowerExpansionFactor = 2.0, double upperExpansionFactor = 2.0)
+        {
+            XTolerance = xTolerance;
+            MaximumIterations = maxIterations;
+            MaximumExpansionSteps = maxExpansionSteps;
+            LowerExpansionFactor = lowerExpansionFactor;
+            UpperExpansionFactor = upperExpansionFactor;
+        }
+
+        public MinimizationResult1D FindMinimum(IObjectiveFunction1D objective, double lowerBound, double upperBound)
+        {
+            if (upperBound <= lowerBound)
+                throw new OptimizationException("Lower bound must be lower than upper bound.");
+
+            double middlePointX = lowerBound + (upperBound - lowerBound)/(1 + Constants.GoldenRatio);
+            IEvaluation1D lower = objective.Evaluate(lowerBound);
+            IEvaluation1D middle = objective.Evaluate(middlePointX);
+            IEvaluation1D upper = objective.Evaluate(upperBound);
+
+            ValueChecker(lower.Value, lowerBound);
+            ValueChecker(middle.Value, middlePointX);
+            ValueChecker(upper.Value, upperBound);
+
+            int expansion_steps = 0;
+            while ((expansion_steps < this.MaximumExpansionSteps) && (upper.Value < middle.Value || lower.Value < middle.Value))
+            {
+                if (lower.Value < middle.Value)
+                {
+                    lowerBound = 0.5*(upperBound + lowerBound) - this.LowerExpansionFactor*0.5*(upperBound - lowerBound);
+                    lower = objective.Evaluate(lowerBound);
+                }
+
+                if (upper.Value < middle.Value)
+                {
+                    upperBound = 0.5*(upperBound + lowerBound) + this.UpperExpansionFactor*0.5*(upperBound - lowerBound);
+                    upper = objective.Evaluate(upperBound);
+                }
+
+                middlePointX = lowerBound + (upperBound - lowerBound)/(1 + Constants.GoldenRatio);
+                middle = objective.Evaluate(middlePointX);
+
+                expansion_steps += 1;
+            }
+
+            if (upper.Value < middle.Value || lower.Value < middle.Value)
+                throw new OptimizationException("Lower and upper bounds do not necessarily bound a minimum.");
+
+            int iterations = 0;
+            while (Math.Abs(upper.Point - lower.Point) > XTolerance && iterations < MaximumIterations)
+            {
+                double testX = lower.Point + (upper.Point - middle.Point);
+                var test = objective.Evaluate(testX);
+                ValueChecker(test.Value, testX);
+
+                if (test.Point < middle.Point)
+                {
+                    if (test.Value > middle.Value)
+                    {
+                        lower = test;
+                    }
+                    else
+                    {
+                        upper = middle;
+                        middle = test;
+                    }
+                }
+                else
+                {
+                    if (test.Value > middle.Value)
+                    {
+                        upper = test;
+                    }
+                    else
+                    {
+                        lower = middle;
+                        middle = test;
+                    }
+                }
+
+                iterations += 1;
+            }
+
+            if (iterations == MaximumIterations)
+                throw new MaximumIterationsException("Max iterations reached.");
+
+            return new MinimizationResult1D(middle, iterations, MinimizationResult.ExitCondition.BoundTolerance);
+        }
+
+        void ValueChecker(double value, double point)
+        {
+            if (Double.IsNaN(value) || Double.IsInfinity(value))
+                throw new Exception("Objective function returned non-finite value.");
+        }
+    }
+}

src/Numerics/Optimization/IObjectiveFunction.cs

@@ -0,0 +1,33 @@
+using MathNet.Numerics.LinearAlgebra;
+
+namespace MathNet.Numerics.Optimization
+{
+    /// <summary>
+    /// Objective function with a frozen evaluation that must not be changed from the outside.
+    /// </summary>
+    public interface IObjectiveFunctionEvaluation
+    {
+        /// <summary>Create a new unevaluated and independent copy of this objective function</summary>
+        IObjectiveFunction CreateNew();
+
+        /// <summary>Create a new independent copy of this objective function, evaluated at the same point.</summary>
+        IObjectiveFunction Fork();
+
+        Vector<double> Point { get; }
+        double Value { get; }
+
+        bool IsGradientSupported { get; }
+        Vector<double> Gradient { get; }
+
+        bool IsHessianSupported { get; }
+        Matrix<double> Hessian { get; }
+    }
+
+    /// <summary>
+    /// Objective function with a mutable evaluation.
+    /// </summary>
+    public interface IObjectiveFunction : IObjectiveFunctionEvaluation
+    {
+        void EvaluateAt(Vector<double> point);
+    }
+}

src/Numerics/Optimization/IUnconstrainedMinimizer.cs

@@ -0,0 +1,9 @@
+using MathNet.Numerics.LinearAlgebra;
+
+namespace MathNet.Numerics.Optimization
+{
+    public interface IUnconstrainedMinimizer
+    {
+        MinimizationResult FindMinimum(IObjectiveFunction objective, Vector<double> initialGuess);
+    }
+}

src/Numerics/Optimization/LineSearch/LineSearchResult.cs

@@ -0,0 +1,43 @@
+// <copyright file="BfgsTest.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2016 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.Optimization.LineSearch
+{
+    public class LineSearchResult : MinimizationResult
+    {
+        public double FinalStep { get; private set; }
+
+        public LineSearchResult(IObjectiveFunction functionInfo, int iterations, double finalStep, ExitCondition reasonForExit)
+            : base(functionInfo, iterations, reasonForExit)
+        {
+            FinalStep = finalStep;
+        }
+    }
+}

src/Numerics/Optimization/LineSearch/StrongWolfeLineSearch.cs

@@ -0,0 +1,50 @@
+// <copyright file="BfgsTest.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2016 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using System;
+
+namespace MathNet.Numerics.Optimization.LineSearch
+{
+    public class StrongWolfeLineSearch : WolfeLineSearch
+    {
+        public StrongWolfeLineSearch(double c1, double c2, double parameterTolerance, int maxIterations = 10)
+            : base(c1, c2, parameterTolerance, maxIterations)
+        {
+            // Argument validation in base class
+        }
+
+        protected override MinimizationResult.ExitCondition WolfeExitCondition { get { return MinimizationResult.ExitCondition.StrongWolfeCriteria; } }
+
+        protected override bool WolfeCondition(double stepDd, double initialDd)
+        {
+            return Math.Abs(stepDd) > C2 * Math.Abs(initialDd);
+        }
+    }
+}

src/Numerics/Optimization/LineSearch/WeakWolfeLineSearch.cs

@@ -0,0 +1,97 @@
+// <copyright file="WeakWolfeLineSearch.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2016 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using System;
+using MathNet.Numerics.LinearAlgebra;
+
+namespace MathNet.Numerics.Optimization.LineSearch
+{
+    /// <summary>
+    /// Search for a step size alpha that satisfies the weak wolfe conditions. The weak Wolfe
+    /// Conditions are
+    /// i)  Armijo Rule:         f(x_k + alpha_k p_k) &lt;= f(x_k) + c1 alpha_k p_k^T g(x_k)
+    /// ii) Curvature Condition: p_k^T g(x_k + alpha_k p_k) &gt;= c2 p_k^T g(x_k)
+    /// where g(x) is the gradient of f(x), 0 &lt; c1 &lt; c2 &lt; 1.
+    /// 
+    /// Implementation is based on http://www.math.washington.edu/~burke/crs/408/lectures/L9-weak-Wolfe.pdf
+    /// 
+    /// references:
+    /// http://en.wikipedia.org/wiki/Wolfe_conditions
+    /// http://www.math.washington.edu/~burke/crs/408/lectures/L9-weak-Wolfe.pdf
+    /// </summary>
+    public class WeakWolfeLineSearch : WolfeLineSearch
+    {
+        public WeakWolfeLineSearch(double c1, double c2, double parameterTolerance, int maxIterations = 10)
+            : base(c1,c2,parameterTolerance,maxIterations)
+        {
+            // Validation in base class
+        }
+
+        protected override MinimizationResult.ExitCondition WolfeExitCondition
+        {
+            get { return MinimizationResult.ExitCondition.WeakWolfeCriteria; }
+        }
+
+        protected override bool WolfeCondition(double stepDd, double initialDd)
+        {
+            return stepDd < C2 * initialDd;
+        }
+
+        protected override void ValidateValue(IObjectiveFunction eval)
+        {
+            if (!IsFinite(eval.Value))
+            {
+                throw new EvaluationException(String.Format("Non-finite value returned by objective function: {0}", eval.Value), eval);
+            }
+        }
+
+        protected override void ValidateInputArguments(IObjectiveFunctionEvaluation startingPoint, Vector<double> searchDirection, double initialStep, double upperBound)
+        {
+            if (!startingPoint.IsGradientSupported)
+                throw new ArgumentException("objective function does not support gradient");
+        }
+
+        protected override void ValidateGradient(IObjectiveFunction eval)
+        {
+            foreach (double x in eval.Gradient)
+            {
+                if (!IsFinite(x))
+                {
+                    throw new EvaluationException(string.Format("Non-finite value returned by gradient: {0}", x), eval);
+                }
+            }
+        }
+
+        static bool IsFinite(double x)
+        {
+            return !(double.IsNaN(x) || double.IsInfinity(x));
+        }
+    }
+}

src/Numerics/Optimization/LineSearch/WolfeLineSearch.cs

@@ -0,0 +1,158 @@
+// <copyright file="BfgsTest.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2016 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using MathNet.Numerics.LinearAlgebra;
+using System;
+
+namespace MathNet.Numerics.Optimization.LineSearch
+{
+    public abstract class WolfeLineSearch
+    {
+        protected double C1 { get; }
+        protected double C2 { get; }
+        protected double ParameterTolerance { get; }
+        protected int MaximumIterations { get; }
+
+        public WolfeLineSearch(double c1, double c2, double parameterTolerance, int maxIterations = 10)
+        {
+            if (c1 <= 0)
+                throw new ArgumentException(string.Format("c1 {0} should be greater than 0", c1));
+            if (c2 <= c1)
+                throw new ArgumentException(string.Format("c1 {0} should be less than c2 {1}", c1, c2));
+            if (c2 >= 1)
+                throw new ArgumentException(string.Format("c2 {0} should be less than 1", c2));
+
+            C1 = c1;
+            C2 = c2;
+            ParameterTolerance = parameterTolerance;
+            MaximumIterations = maxIterations;
+        }
+
+        /// <summary>Implemented following http://www.math.washington.edu/~burke/crs/408/lectures/L9-weak-Wolfe.pdf</summary>
+        /// <param name="startingPoint">The objective function being optimized, evaluated at the starting point of the search</param>
+        /// <param name="searchDirection">Search direction</param>
+        /// <param name="initialStep">Initial size of the step in the search direction</param>
+        public LineSearchResult FindConformingStep(IObjectiveFunctionEvaluation startingPoint, Vector<double> searchDirection, double initialStep)
+        {
+            return FindConformingStep(startingPoint, searchDirection, initialStep, double.PositiveInfinity);
+        }
+
+        /// <summary></summary>
+        /// <param name="startingPoint">The objective function being optimized, evaluated at the starting point of the search</param>
+        /// <param name="searchDirection">Search direction</param>
+        /// <param name="initialStep">Initial size of the step in the search direction</param>
+        /// <param name="upperBound">The upper bound</param>
+        public LineSearchResult FindConformingStep(IObjectiveFunctionEvaluation startingPoint, Vector<double> searchDirection, double initialStep, double upperBound)
+        {
+            ValidateInputArguments(startingPoint, searchDirection, initialStep, upperBound);
+
+            double lowerBound = 0.0;
+            double step = initialStep;
+
+            double initialValue = startingPoint.Value;
+            Vector<double> initialGradient = startingPoint.Gradient;
+
+            double initialDd = searchDirection * initialGradient;
+
+            IObjectiveFunction objective = startingPoint.CreateNew();
+            int ii;
+            MinimizationResult.ExitCondition reasonForExit = MinimizationResult.ExitCondition.None;
+            for (ii = 0; ii < MaximumIterations; ++ii)
+            {
+                objective.EvaluateAt(startingPoint.Point + searchDirection * step);
+                ValidateGradient(objective);  
+                ValidateValue(objective);       
+
+                double stepDd = searchDirection * objective.Gradient;
+
+                if (objective.Value > initialValue + C1 * step * initialDd)
+                {
+                    upperBound = step;
+                    step = 0.5 * (lowerBound + upperBound);
+                }
+                else if (WolfeCondition(stepDd,initialDd))
+                {
+                    lowerBound = step;
+                    step = double.IsPositiveInfinity(upperBound) ? 2 * lowerBound : 0.5 * (lowerBound + upperBound);
+                }
+                else
+                {
+                    reasonForExit = WolfeExitCondition; 
+                    break;
+                }
+
+                if (!double.IsInfinity(upperBound))
+                {
+                    double maxRelChange = 0.0;
+                    for (int jj = 0; jj < objective.Point.Count; ++jj)
+                    {
+                        double tmp = Math.Abs(searchDirection[jj] * (upperBound - lowerBound)) / Math.Max(Math.Abs(objective.Point[jj]), 1.0);
+                        maxRelChange = Math.Max(maxRelChange, tmp);
+                    }
+                    if (maxRelChange < ParameterTolerance)
+                    {
+                        reasonForExit = MinimizationResult.ExitCondition.LackOfProgress;
+                        break;
+                    }
+                }
+            }
+
+            if (ii == MaximumIterations && Double.IsPositiveInfinity(upperBound))
+            {
+                throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached. Function appears to be unbounded in search direction.", MaximumIterations));
+            }
+
+            if (ii == MaximumIterations)
+            {
+                throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", MaximumIterations));
+            }
+
+            return new LineSearchResult(objective, ii, step, reasonForExit);
+        }
+        protected abstract MinimizationResult.ExitCondition WolfeExitCondition { get; }
+
+        protected abstract bool WolfeCondition(double stepDd, double initialDd);
+
+        protected virtual void ValidateGradient(IObjectiveFunction objective)
+        {
+        }
+        protected virtual void ValidateValue(IObjectiveFunction objective)
+        {
+        }
+
+        protected virtual void ValidateInputArguments(IObjectiveFunctionEvaluation startingPoint, Vector<double> searchDirection, double initialStep, double upperBound)
+        {
+
+        }
+    }
+
+        
+    
+}

src/Numerics/Optimization/MinimizationResult.cs

@@ -0,0 +1,32 @@
+using MathNet.Numerics.LinearAlgebra;
+
+namespace MathNet.Numerics.Optimization
+{
+    public class MinimizationResult
+    {
+        public enum ExitCondition
+        {
+            None,
+            RelativeGradient,
+            LackOfProgress,
+            AbsoluteGradient,
+            WeakWolfeCriteria,
+            BoundTolerance,
+            StrongWolfeCriteria,
+            LackOfFunctionImprovement,
+            Converged
+        }
+
+        public Vector<double> MinimizingPoint { get { return FunctionInfoAtMinimum.Point; } }
+        public IObjectiveFunction FunctionInfoAtMinimum { get; private set; }
+        public int Iterations { get; private set; }
+        public ExitCondition ReasonForExit { get; private set; }
+
+        public MinimizationResult(IObjectiveFunction functionInfo, int iterations, ExitCondition reasonForExit)
+        {
+            FunctionInfoAtMinimum = functionInfo;
+            Iterations = iterations;
+            ReasonForExit = reasonForExit;
+        }
+    }
+}

src/Numerics/Optimization/MinimizationResult1D.cs

@@ -0,0 +1,17 @@
+namespace MathNet.Numerics.Optimization
+{
+    public class MinimizationResult1D
+    {
+        public double MinimizingPoint { get { return FunctionInfoAtMinimum.Point; } }
+        public IEvaluation1D FunctionInfoAtMinimum { get; private set; }
+        public int Iterations { get; private set; }
+        public MinimizationResult.ExitCondition ReasonForExit { get; private set; }
+
+        public MinimizationResult1D(IEvaluation1D functionInfo, int iterations, MinimizationResult.ExitCondition reasonForExit)
+        {
+            FunctionInfoAtMinimum = functionInfo;
+            Iterations = iterations;
+            ReasonForExit = reasonForExit;
+        }
+    }
+}

src/Numerics/Optimization/MinimizationWithLineSearchResult.cs

@@ -0,0 +1,15 @@
+namespace MathNet.Numerics.Optimization
+{
+    public class MinimizationWithLineSearchResult : MinimizationResult
+    {
+        public int TotalLineSearchIterations { get; private set; }
+        public int IterationsWithNonTrivialLineSearch { get; private set; }
+
+        public MinimizationWithLineSearchResult(IObjectiveFunction functionInfo, int iterations, ExitCondition reasonForExit, int totalLineSearchIterations, int iterationsWithNonTrivialLineSearch)
+            : base(functionInfo, iterations, reasonForExit)
+        {
+            TotalLineSearchIterations = totalLineSearchIterations;
+            IterationsWithNonTrivialLineSearch = iterationsWithNonTrivialLineSearch;
+        }
+    }
+}

src/Numerics/Optimization/NelderMeadSimplex.cs

@@ -0,0 +1,415 @@
+// <copyright file="NelderMeadSimplex.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+// 
+// Copyright (c) 2009-2015 Math.NET
+// 
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+// 
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+// 
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+// Converted from code relased with a MIT liscense available at https://code.google.com/p/nelder-mead-simplex/
+
+using MathNet.Numerics.LinearAlgebra;
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+
+namespace MathNet.Numerics.Optimization
+{
+    /// <summary>
+    /// Class implementing the Nelder-Mead simplex algorithm, used to find a minima when no gradient is available. 
+    /// Called fminsearch() in Matlab. A description of the algorithm can be found at 
+    /// http://se.mathworks.com/help/matlab/math/optimizing-nonlinear-functions.html#bsgpq6p-11
+    /// or
+    /// https://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method
+    /// </summary>
+    public sealed class NelderMeadSimplex
+    {
+        private static readonly double JITTER = 1e-10d;           // a small value used to protect against floating point noise
+
+        public double ConvergenceTolerance { get; set; }
+        public int MaximumIterations { get; set; }
+
+        public NelderMeadSimplex(double convergenceTolerance, int maximumIterations)
+        {
+            ConvergenceTolerance = convergenceTolerance;
+            MaximumIterations = maximumIterations;
+        }
+
+        /// <summary>
+        /// Finds the minimum of the objective function without an intial pertubation, the default values used
+        /// by fminsearch() in Matlab are used instead
+        /// http://se.mathworks.com/help/matlab/math/optimizing-nonlinear-functions.html#bsgpq6p-11
+        /// </summary>
+        /// <param name="objectiveFunction">The objective function, no gradient or hessian needed</param>
+        /// <param name="initialGuess">The intial guess</param>
+        /// <returns>The minimum point</returns>
+        public MinimizationResult FindMinimum(IObjectiveFunction objectiveFunction, Vector<double> initialGuess)
+        {
+            var initalPertubation = new MathNet.Numerics.LinearAlgebra.Double.DenseVector(initialGuess.Count);
+            for (int i = 0; i < initialGuess.Count; i++)
+            {
+                initalPertubation[i] = initialGuess[i] == 0.0 ? 0.00025 : initialGuess[i] * 0.05;
+            }
+            return FindMinimum(objectiveFunction, initialGuess, initalPertubation);
+        }
+
+        /// <summary>
+        /// Finds the minimum of the objective function with an intial pertubation
+        /// </summary>
+        /// <param name="objectiveFunction">The objective function, no gradient or hessian needed</param>
+        /// <param name="initialGuess">The intial guess</param>
+        /// <param name="initalPertubation">The inital pertubation</param>
+        /// <returns>The minimum point</returns>
+        public MinimizationResult FindMinimum(IObjectiveFunction objectiveFunction, Vector<double> initialGuess, Vector<double> initalPertubation)
+        {
+            // confirm that we are in a position to commence
+            if (objectiveFunction == null)
+                throw new ArgumentNullException("objectiveFunction","ObjectiveFunction must be set to a valid ObjectiveFunctionDelegate");
+
+            if (initialGuess == null)
+                throw new ArgumentNullException("initialGuess", "initialGuess must be initialized");
+
+            if (initialGuess == null)
+                throw new ArgumentNullException("initalPertubation", "initalPertubation must be initialized, if unknown use overloaded version of FindMinimum()");
+
+            SimplexConstant[] simplexConstants = SimplexConstant.CreateSimplexConstantsFromVectors(initialGuess,initalPertubation);
+
+            // create the initial simplex
+            int numDimensions = simplexConstants.Length;
+            int numVertices = numDimensions + 1;
+            Vector<double>[] vertices = InitializeVertices(simplexConstants);
+            double[] errorValues = new double[numVertices];
+
+            int evaluationCount = 0;
+            MinimizationResult.ExitCondition exitCondition = MinimizationResult.ExitCondition.None;
+            ErrorProfile errorProfile;
+
+            errorValues = InitializeErrorValues(vertices, objectiveFunction);
+
+            // iterate until we converge, or complete our permitted number of iterations
+            while (true)
+            {
+                errorProfile = EvaluateSimplex(errorValues);
+
+                // see if the range in point heights is small enough to exit
+                if (HasConverged(ConvergenceTolerance, errorProfile, errorValues))
+                {
+                    exitCondition = MinimizationResult.ExitCondition.Converged;
+                    break;
+                }
+
+                // attempt a reflection of the simplex
+                double reflectionPointValue = TryToScaleSimplex(-1.0, ref errorProfile, vertices, errorValues, objectiveFunction);
+                ++evaluationCount;
+                if (reflectionPointValue <= errorValues[errorProfile.LowestIndex])
+                {
+                    // it's better than the best point, so attempt an expansion of the simplex
+                    double expansionPointValue = TryToScaleSimplex(2.0, ref errorProfile, vertices, errorValues, objectiveFunction);
+                    ++evaluationCount;
+                }
+                else if (reflectionPointValue >= errorValues[errorProfile.NextHighestIndex])
+                {
+                    // it would be worse than the second best point, so attempt a contraction to look
+                    // for an intermediate point
+                    double currentWorst = errorValues[errorProfile.HighestIndex];
+                    double contractionPointValue = TryToScaleSimplex(0.5, ref errorProfile, vertices, errorValues, objectiveFunction);
+                    ++evaluationCount;
+                    if (contractionPointValue >= currentWorst)
+                    {
+                        // that would be even worse, so let's try to contract uniformly towards the low point; 
+                        // don't bother to update the error profile, we'll do it at the start of the
+                        // next iteration
+                        ShrinkSimplex(errorProfile, vertices, errorValues, objectiveFunction);
+                        evaluationCount += numVertices; // that required one function evaluation for each vertex; keep track
+                    }
+                }
+                // check to see if we have exceeded our alloted number of evaluations
+                if (evaluationCount >= MaximumIterations)
+                {
+                    throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", MaximumIterations));
+                }
+            }
+            var regressionResult = new MinimizationResult(objectiveFunction, evaluationCount, exitCondition);
+            return regressionResult;
+        }
+
+        /// <summary>
+        /// Evaluate the objective function at each vertex to create a corresponding
+        /// list of error values for each vertex
+        /// </summary>
+        /// <param name="vertices"></param>
+        /// <param name="objectiveFunction"></param>
+        /// <returns></returns>
+        private static double[] InitializeErrorValues(Vector<double>[] vertices, IObjectiveFunction objectiveFunction)
+        {
+            double[] errorValues = new double[vertices.Length];
+            for (int i = 0; i < vertices.Length; i++)
+            {
+                objectiveFunction.EvaluateAt(vertices[i]);
+                errorValues[i] = objectiveFunction.Value;
+            }
+            return errorValues;
+        }
+
+        /// <summary>
+        /// Check whether the points in the error profile have so little range that we
+        /// consider ourselves to have converged
+        /// </summary>
+        /// <param name="convergenceTolerance"></param>
+        /// <param name="errorProfile"></param>
+        /// <param name="errorValues"></param>
+        /// <returns></returns>
+        private static bool HasConverged(double convergenceTolerance, ErrorProfile errorProfile, double[] errorValues)
+        {
+            double range = 2 * Math.Abs(errorValues[errorProfile.HighestIndex] - errorValues[errorProfile.LowestIndex]) /
+                (Math.Abs(errorValues[errorProfile.HighestIndex]) + Math.Abs(errorValues[errorProfile.LowestIndex]) + JITTER);
+
+            if (range < convergenceTolerance)
+            {
+                return true;
+            }
+            else
+            {
+                return false;
+            }
+        }
+
+        /// <summary>
+        /// Examine all error values to determine the ErrorProfile
+        /// </summary>
+        /// <param name="errorValues"></param>
+        /// <returns></returns>
+        private static ErrorProfile EvaluateSimplex(double[] errorValues)
+        {
+            ErrorProfile errorProfile = new ErrorProfile();
+            if (errorValues[0] > errorValues[1])
+            {
+                errorProfile.HighestIndex = 0;
+                errorProfile.NextHighestIndex = 1;
+            }
+            else
+            {
+                errorProfile.HighestIndex = 1;
+                errorProfile.NextHighestIndex = 0;
+            }
+
+            for (int index = 0; index < errorValues.Length; index++)
+            {
+                double errorValue = errorValues[index];
+                if (errorValue <= errorValues[errorProfile.LowestIndex])
+                {
+                    errorProfile.LowestIndex = index;
+                }
+                if (errorValue > errorValues[errorProfile.HighestIndex])
+                {
+                    errorProfile.NextHighestIndex = errorProfile.HighestIndex; // downgrade the current highest to next highest
+                    errorProfile.HighestIndex = index;
+                }
+                else if (errorValue > errorValues[errorProfile.NextHighestIndex] && index != errorProfile.HighestIndex)
+                {
+                    errorProfile.NextHighestIndex = index;
+                }
+            }
+
+            return errorProfile;
+        }
+
+        /// <summary>
+        /// Construct an initial simplex, given starting guesses for the constants, and
+        /// initial step sizes for each dimension
+        /// </summary>
+        /// <param name="simplexConstants"></param>
+        /// <returns></returns>
+        private static Vector<double>[] InitializeVertices(SimplexConstant[] simplexConstants)
+        {
+            int numDimensions = simplexConstants.Length;
+            Vector<double>[] vertices = new Vector<double>[numDimensions + 1];
+
+            // define one point of the simplex as the given initial guesses
+            var p0 = new MathNet.Numerics.LinearAlgebra.Double.DenseVector(numDimensions);
+            for (int i = 0; i < numDimensions; i++)
+            {
+                p0[i] = simplexConstants[i].Value;
+            }
+
+            // now fill in the vertices, creating the additional points as:
+            // P(i) = P(0) + Scale(i) * UnitVector(i)
+            vertices[0] = p0;
+            for (int i = 0; i < numDimensions; i++)
+            {
+                double scale = simplexConstants[i].InitialPerturbation;
+                Vector<double> unitVector = new MathNet.Numerics.LinearAlgebra.Double.DenseVector(numDimensions);
+                unitVector[i] = 1;
+                vertices[i + 1] = p0.Add(unitVector.Multiply(scale));
+            }
+            return vertices;
+        }
+
+        /// <summary>
+        /// Test a scaling operation of the high point, and replace it if it is an improvement
+        /// </summary>
+        /// <param name="scaleFactor"></param>
+        /// <param name="errorProfile"></param>
+        /// <param name="vertices"></param>
+        /// <param name="errorValues"></param>
+        /// <param name="objectiveFunction"></param>
+        /// <returns></returns>
+        private static double TryToScaleSimplex(double scaleFactor, ref ErrorProfile errorProfile, Vector<double>[] vertices,
+                                          double[] errorValues, IObjectiveFunction objectiveFunction)
+        {
+            // find the centroid through which we will reflect
+            Vector<double> centroid = ComputeCentroid(vertices, errorProfile);
+
+            // define the vector from the centroid to the high point
+            Vector<double> centroidToHighPoint = vertices[errorProfile.HighestIndex].Subtract(centroid);
+
+            // scale and position the vector to determine the new trial point
+            Vector<double> newPoint = centroidToHighPoint.Multiply(scaleFactor).Add(centroid);
+
+            // evaluate the new point
+            objectiveFunction.EvaluateAt(newPoint);
+            double newErrorValue = objectiveFunction.Value;
+
+            // if it's better, replace the old high point
+            if (newErrorValue < errorValues[errorProfile.HighestIndex])
+            {
+                vertices[errorProfile.HighestIndex] = newPoint;
+                errorValues[errorProfile.HighestIndex] = newErrorValue;
+            }
+
+            return newErrorValue;
+        }
+
+        /// <summary>
+        /// Contract the simplex uniformly around the lowest point
+        /// </summary>
+        /// <param name="errorProfile"></param>
+        /// <param name="vertices"></param>
+        /// <param name="errorValues"></param>
+        /// <param name="objectiveFunction"></param>
+        private static void ShrinkSimplex(ErrorProfile errorProfile, Vector<double>[] vertices, double[] errorValues,
+                                      IObjectiveFunction objectiveFunction)
+        {
+            Vector<double> lowestVertex = vertices[errorProfile.LowestIndex];
+            for (int i = 0; i < vertices.Length; i++)
+            {
+                if (i != errorProfile.LowestIndex)
+                {
+                    vertices[i] = (vertices[i].Add(lowestVertex)).Multiply(0.5);
+                    objectiveFunction.EvaluateAt(vertices[i]);
+                    errorValues[i] = objectiveFunction.Value;
+                }
+            }
+        }
+
+        /// <summary>
+        /// Compute the centroid of all points except the worst
+        /// </summary>
+        /// <param name="vertices"></param>
+        /// <param name="errorProfile"></param>
+        /// <returns></returns>
+        private static Vector<double> ComputeCentroid(Vector<double>[] vertices, ErrorProfile errorProfile)
+        {
+            int numVertices = vertices.Length;
+            // find the centroid of all points except the worst one
+            Vector<double> centroid = new MathNet.Numerics.LinearAlgebra.Double.DenseVector(numVertices - 1);
+            for (int i = 0; i < numVertices; i++)
+            {
+                if (i != errorProfile.HighestIndex)
+                {
+                    centroid = centroid.Add(vertices[i]);
+                }
+            }
+            return centroid.Multiply(1.0d / (numVertices - 1));
+        }
+
+        private sealed class SimplexConstant
+        {
+            private double _value;
+            private double _initialPerturbation;
+
+            public SimplexConstant(double value, double initialPerturbation)
+            {
+                _value = value;
+                _initialPerturbation = initialPerturbation;
+            }
+
+            /// <summary>
+            /// The value of the constant
+            /// </summary>
+            public double Value
+            {
+                get { return _value; }
+                set { _value = value; }
+            }
+
+            // The size of the initial perturbation
+            public double InitialPerturbation
+            {
+                get { return _initialPerturbation; }
+                set { _initialPerturbation = value; }
+            }
+
+            public static SimplexConstant[] CreateSimplexConstantsFromVectors(Vector<double> initialGuess, Vector<double> initialPertubation)
+            {
+                var constants = new SimplexConstant[initialGuess.Count];
+                for (int i = 0; i < constants.Length;i++ )
+                {
+                    constants[i] = new SimplexConstant(initialGuess[i], initialPertubation[i]);
+                }
+                return constants;
+            }
+        }
+
+        private sealed class ErrorProfile
+        {
+            private int _highestIndex;
+            private int _nextHighestIndex;
+            private int _lowestIndex;
+
+            public int HighestIndex
+            {
+                get { return _highestIndex; }
+                set { _highestIndex = value; }
+            }
+
+            public int NextHighestIndex
+            {
+                get { return _nextHighestIndex; }
+                set { _nextHighestIndex = value; }
+            }
+
+            public int LowestIndex
+            {
+                get { return _lowestIndex; }
+                set { _lowestIndex = value; }
+            }
+        }
+    }
+
+    
+
+}

src/Numerics/Optimization/NewtonMinimizer.cs

@@ -0,0 +1,130 @@
+using System;
+using MathNet.Numerics.LinearAlgebra;
+using MathNet.Numerics.Optimization.LineSearch;
+
+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 gradientTolerance, int maximumIterations, bool useLineSearch = false)
+        {
+            GradientTolerance = gradientTolerance;
+            MaximumIterations = maximumIterations;
+            UseLineSearch = useLineSearch;
+        }
+
+        public MinimizationResult FindMinimum(IObjectiveFunction objective, Vector<double> initialGuess)
+        {
+            if (!objective.IsGradientSupported)
+            {
+                throw new IncompatibleObjectiveException("Gradient not supported in objective function, but required for Newton minimization.");
+            }
+
+            if (!objective.IsHessianSupported)
+            {
+                throw new IncompatibleObjectiveException("Hessian not supported in objective function, but required for Newton minimization.");
+            }
+
+            // Check that we're not already done
+            objective.EvaluateAt(initialGuess);
+            ValidateGradient(objective);
+            if (ExitCriteriaSatisfied(objective.Gradient))
+            {
+                return new MinimizationResult(objective, 0, MinimizationResult.ExitCondition.AbsoluteGradient);
+            }
+
+            // Set up line search algorithm
+            var lineSearcher = new WeakWolfeLineSearch(1e-4, 0.9, 1e-4, maxIterations: 1000);
+
+            // Subsequent steps
+            int iterations = 0;
+            int totalLineSearchSteps = 0;
+            int iterationsWithNontrivialLineSearch = 0;
+            bool tmpLineSearch = false;
+            while (!ExitCriteriaSatisfied(objective.Gradient) && iterations < MaximumIterations)
+            {
+                ValidateHessian(objective);
+
+                var searchDirection = objective.Hessian.LU().Solve(-objective.Gradient);
+                if (searchDirection * objective.Gradient >= 0)
+                {
+                    searchDirection = -objective.Gradient;
+                    tmpLineSearch = true;
+                }
+
+                if (UseLineSearch || tmpLineSearch)
+                {
+                    LineSearchResult result;
+                    try
+                    {
+                        result = lineSearcher.FindConformingStep(objective, searchDirection, 1.0);
+                    }
+                    catch (Exception e)
+                    {
+                        throw new InnerOptimizationException("Line search failed.", e);
+                    }
+
+                    iterationsWithNontrivialLineSearch += result.Iterations > 0 ? 1 : 0;
+                    totalLineSearchSteps += result.Iterations;
+                    objective = result.FunctionInfoAtMinimum;
+                }
+                else
+                {
+                    objective.EvaluateAt(objective.Point + searchDirection);
+                }
+
+                ValidateGradient(objective);
+
+                tmpLineSearch = false;
+                iterations += 1;
+            }
+
+            if (iterations == MaximumIterations)
+            {
+                throw new MaximumIterationsException(String.Format("Maximum iterations ({0}) reached.", MaximumIterations));
+            }
+
+            return new MinimizationWithLineSearchResult(objective, iterations, MinimizationResult.ExitCondition.AbsoluteGradient, totalLineSearchSteps, iterationsWithNontrivialLineSearch);
+        }
+
+        bool ExitCriteriaSatisfied(Vector<double> gradient)
+        {
+            return gradient.Norm(2.0) < GradientTolerance;
+        }
+
+        static void ValidateGradient(IObjectiveFunction 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(IObjectiveFunction eval)
+        {
+            if (Double.IsNaN(eval.Value) || Double.IsInfinity(eval.Value))
+                throw new EvaluationException("Non-finite objective function returned.", eval);
+        }
+
+        private void ValidateHessian(IObjectiveFunction 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/Numerics/Optimization/ObjectiveFunction.cs

@@ -0,0 +1,65 @@
+using System;
+using MathNet.Numerics.LinearAlgebra;
+using MathNet.Numerics.Optimization.ObjectiveFunctions;
+
+namespace MathNet.Numerics.Optimization
+{
+    public static class ObjectiveFunction
+    {
+        /// <summary>
+        /// Objective function where neither Gradient nor Hessian is available.
+        /// </summary>
+        public static IObjectiveFunction Value(Func<Vector<double>, double> function)
+        {
+            return new ValueObjectiveFunction(function);
+        }
+
+        /// <summary>
+        /// Objective function where the Gradient is available. Greedy evaluation.
+        /// </summary>
+        public static IObjectiveFunction Gradient(Func<Vector<double>, Tuple<double, Vector<double>>> function)
+        {
+            return new GradientObjectiveFunction(function);
+        }
+
+        /// <summary>
+        /// Objective function where the Gradient is available. Lazy evaluation.
+        /// </summary>
+        public static IObjectiveFunction Gradient(Func<Vector<double>, double> function, Func<Vector<double>, Vector<double>> gradient)
+        {
+            return new LazyObjectiveFunction(function, gradient: gradient);
+        }
+
+        /// <summary>
+        /// Objective function where the Hessian is available. Greedy evaluation.
+        /// </summary>
+        public static IObjectiveFunction Hessian(Func<Vector<double>, Tuple<double, Matrix<double>>> function)
+        {
+            return new HessianObjectiveFunction(function);
+        }
+
+        /// <summary>
+        /// Objective function where the Hessian is available. Lazy evaluation.
+        /// </summary>
+        public static IObjectiveFunction Hessian(Func<Vector<double>, double> function, Func<Vector<double>, Matrix<double>> hessian)
+        {
+            return new LazyObjectiveFunction(function, hessian: hessian);
+        }
+
+        /// <summary>
+        /// Objective function where both Gradient and Hessian are available. Greedy evaluation.
+        /// </summary>
+        public static IObjectiveFunction GradientHessian(Func<Vector<double>, Tuple<double, Vector<double>, Matrix<double>>> function)
+        {
+            return new GradientHessianObjectiveFunction(function);
+        }
+
+        /// <summary>
+        /// Objective function where both Gradient and Hessian are available. Lazy evaluation.
+        /// </summary>
+        public static IObjectiveFunction GradientHessian(Func<Vector<double>, double> function, Func<Vector<double>, Vector<double>> gradient, Func<Vector<double>, Matrix<double>> hessian)
+        {
+            return new LazyObjectiveFunction(function, gradient: gradient, hessian: hessian);
+        }
+    }
+}

src/Numerics/Optimization/ObjectiveFunction1D.cs

@@ -0,0 +1,98 @@
+using System;
+
+namespace MathNet.Numerics.Optimization
+{
+    public interface IEvaluation1D
+    {
+        double Point { get; }
+        double Value { get; }
+        double Derivative { get; }
+        double SecondDerivative { get; }
+    }
+
+    public interface IObjectiveFunction1D
+    {
+        bool DerivativeSupported { get; }
+        bool SecondDerivativeSupported { get; }
+        IEvaluation1D Evaluate(double point);
+    }
+
+    public class CachedEvaluation1D : IEvaluation1D
+    {
+        private double? _value;
+        private double? _derivative;
+        private double? _secondDerivative;
+        private readonly SimpleObjectiveFunction1D _objectiveObject;
+        private readonly double _point;
+
+        public CachedEvaluation1D(SimpleObjectiveFunction1D f, double point)
+        {
+            _objectiveObject = f;
+            _point = point;
+        }
+        private double SetValue()
+        {
+            _value = _objectiveObject.Objective(_point);
+            return _value.Value;
+        }
+        private double SetDerivative()
+        {
+            _derivative = _objectiveObject.Derivative(_point);
+            return _derivative.Value;
+        }
+        private double SetSecondDerivative()
+        {
+            _secondDerivative = _objectiveObject.SecondDerivative(_point);
+            return _secondDerivative.Value;
+        }
+
+        public double Point { get { return _point; } }
+        public double Value { get { return _value ?? SetValue(); } }
+        public double Derivative { get { return _derivative ?? SetDerivative(); } }
+        public double SecondDerivative { get { return _secondDerivative ?? SetSecondDerivative(); } }
+
+    }
+
+    public class SimpleObjectiveFunction1D : IObjectiveFunction1D
+    {
+        public Func<double, double> Objective { get; private set; }
+        public Func<double, double> Derivative { get; private set; }
+        public Func<double, double> SecondDerivative { get; private set; }
+
+        public SimpleObjectiveFunction1D(Func<double, double> objective)
+        {
+            Objective = objective;
+            Derivative = null;
+            SecondDerivative = null;
+        }
+
+        public SimpleObjectiveFunction1D(Func<double, double> objective, Func<double, double> derivative)
+        {
+            Objective = objective;
+            Derivative = derivative;
+            SecondDerivative = null;
+        }
+
+        public SimpleObjectiveFunction1D(Func<double, double> objective, Func<double, double> derivative, Func<double,double> secondDerivative)
+        {
+            Objective = objective;
+            Derivative = derivative;
+            SecondDerivative = secondDerivative;
+        }
+
+        public bool DerivativeSupported
+        {
+            get { return Derivative != null; }
+        }
+
+        public bool SecondDerivativeSupported
+        {
+            get { return SecondDerivative != null; }
+        }
+
+        public IEvaluation1D Evaluate(double point)
+        {
+            return new CachedEvaluation1D(this, point);
+        }
+    }
+}

src/Numerics/Optimization/ObjectiveFunctions/ForwardDifferenceGradientObjectiveFunction.cs

@@ -0,0 +1,145 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using MathNet.Numerics.LinearAlgebra;
+
+namespace MathNet.Numerics.Optimization.ObjectiveFunctions
+{
+    /// <summary>
+    /// Adapts an objective function with only value implemented 
+    /// to provide a gradient as well. Gradient calculation is 
+    /// done using the finite difference method, specifically 
+    /// forward differences.
+    /// 
+    /// For each gradient computed, the algorithm requires an 
+    /// additional number of function evaluations equal to the 
+    /// functions's number of input parameters.
+    /// </summary>
+    public class ForwardDifferenceGradientObjectiveFunction : IObjectiveFunction
+    {
+        public IObjectiveFunction InnerObjectiveFunction { get; protected set; }
+        protected Vector<double> LowerBound { get; set; }
+        protected Vector<double> UpperBound { get; set; }
+
+        protected bool ValueEvaluated { get; set; } = false;
+        protected bool GradientEvaluated { get; set; } = false;
+        private Vector<double> _gradient;
+
+        public double MinimumIncrement { get; set; }
+        public double RelativeIncrement { get; set; }
+
+        public ForwardDifferenceGradientObjectiveFunction(IObjectiveFunction valueOnlyObj, Vector<double> lowerBound, Vector<double> upperBound, double relativeIncrement=1e-5, double minimumIncrement=1e-8)
+        {
+            InnerObjectiveFunction = valueOnlyObj;
+            LowerBound = lowerBound;
+            UpperBound = upperBound;
+            _gradient = new LinearAlgebra.Double.DenseVector(LowerBound.Count);
+            RelativeIncrement = relativeIncrement;
+            MinimumIncrement = minimumIncrement;
+        }
+
+        protected void EvaluateValue()
+        {
+            ValueEvaluated = true;
+        }
+
+        protected void EvaluateGradient()
+        {
+            if (!ValueEvaluated)
+                EvaluateValue();
+
+            var tmp_point = Point.Clone();
+            var tmp_obj = InnerObjectiveFunction.CreateNew();
+            for (int ii = 0; ii < _gradient.Count; ++ii)
+            {
+                var orig_point = tmp_point[ii];
+                var rel_incr = orig_point * RelativeIncrement;
+                var h = Math.Max(rel_incr, MinimumIncrement);
+                var mult = 1;
+                if (orig_point + h > UpperBound[ii])
+                    mult = -1;
+
+                tmp_point[ii] = orig_point + mult*h;
+                tmp_obj.EvaluateAt(tmp_point);
+                double bumped_value = tmp_obj.Value;
+                _gradient[ii] = (mult * bumped_value - mult * InnerObjectiveFunction.Value) / h;
+
+                tmp_point[ii] = orig_point;
+            }
+            GradientEvaluated = true;
+        }
+
+        public Vector<double> Gradient
+        {
+            get
+            {
+                if (!GradientEvaluated)
+                    EvaluateGradient();
+                return _gradient;
+            }
+            protected set { _gradient = value; }
+        }
+
+        public Matrix<double> Hessian
+        {
+            get
+            {
+                throw new NotImplementedException();
+            }
+        }
+
+        public bool IsGradientSupported
+        {
+            get
+            {
+                return true;
+            }
+        }
+
+        public bool IsHessianSupported
+        {
+            get
+            {
+                return false;
+            }
+        }
+
+        public Vector<double> Point { get; protected set; }
+
+        public double Value
+        {
+            get
+            {
+                if (!ValueEvaluated)
+                    EvaluateValue();
+                return this.InnerObjectiveFunction.Value;
+            }
+        }
+
+        public IObjectiveFunction CreateNew()
+        {
+            var tmp = new ForwardDifferenceGradientObjectiveFunction(this.InnerObjectiveFunction.CreateNew(), LowerBound, UpperBound, this.RelativeIncrement, this.MinimumIncrement);
+            return tmp;
+        }
+
+        public void EvaluateAt(Vector<double> point)
+        {
+            Point = point;
+            ValueEvaluated = false;
+            GradientEvaluated = false;
+            InnerObjectiveFunction.EvaluateAt(point);
+        }
+
+        public IObjectiveFunction Fork()
+        {
+            return new ForwardDifferenceGradientObjectiveFunction(this.InnerObjectiveFunction.Fork(), LowerBound, UpperBound, this.RelativeIncrement, this.MinimumIncrement)
+            {
+                Point = Point?.Clone(),
+                GradientEvaluated = GradientEvaluated,
+                ValueEvaluated = ValueEvaluated,
+                _gradient = _gradient?.Clone()
+            };
+        }
+    }
+}

src/Numerics/Optimization/ObjectiveFunctions/GradientHessianObjectiveFunction.cs

@@ -0,0 +1,57 @@
+using System;
+using MathNet.Numerics.LinearAlgebra;
+
+namespace MathNet.Numerics.Optimization.ObjectiveFunctions
+{
+    internal class GradientHessianObjectiveFunction : IObjectiveFunction
+    {
+        readonly Func<Vector<double>, Tuple<double, Vector<double>, Matrix<double>>> _function;
+
+        public GradientHessianObjectiveFunction(Func<Vector<double>, Tuple<double, Vector<double>, Matrix<double>>> function)
+        {
+            _function = function;
+        }
+
+        public IObjectiveFunction CreateNew()
+        {
+            return new GradientHessianObjectiveFunction(_function);
+        }
+
+        public IObjectiveFunction Fork()
+        {
+            // no need to deep-clone values since they are replaced on evaluation
+            return new GradientHessianObjectiveFunction(_function)
+            {
+                Point = Point,
+                Value = Value,
+                Gradient = Gradient,
+                Hessian = Hessian
+            };
+        }
+
+        public bool IsGradientSupported
+        {
+            get { return true; }
+        }
+
+        public bool IsHessianSupported
+        {
+            get { return true; }
+        }
+
+        public void EvaluateAt(Vector<double> point)
+        {
+            Point = point;
+
+            var result = _function(point);
+            Value = result.Item1;
+            Gradient = result.Item2;
+            Hessian = result.Item3;
+        }
+
+        public Vector<double> Point { get; private set; }
+        public double Value { get; private set; }
+        public Vector<double> Gradient { get; private set; }
+        public Matrix<double> Hessian { get; private set; }
+    }
+}

src/Numerics/Optimization/ObjectiveFunctions/GradientObjectiveFunction.cs

@@ -0,0 +1,59 @@
+using System;
+using MathNet.Numerics.LinearAlgebra;
+
+namespace MathNet.Numerics.Optimization.ObjectiveFunctions
+{
+    internal class GradientObjectiveFunction : IObjectiveFunction
+    {
+        readonly Func<Vector<double>, Tuple<double, Vector<double>>> _function;
+
+        public GradientObjectiveFunction(Func<Vector<double>, Tuple<double, Vector<double>>> function)
+        {
+            _function = function;
+        }
+
+        public IObjectiveFunction CreateNew()
+        {
+            return new GradientObjectiveFunction(_function);
+        }
+
+        public IObjectiveFunction Fork()
+        {
+            // no need to deep-clone values since they are replaced on evaluation
+            return new GradientObjectiveFunction(_function)
+            {
+                Point = Point,
+                Value = Value,
+                Gradient = Gradient
+            };
+        }
+
+        public bool IsGradientSupported
+        {
+            get { return true; }
+        }
+
+        public bool IsHessianSupported
+        {
+            get { return false; }
+        }
+
+        public void EvaluateAt(Vector<double> point)
+        {
+            Point = point;
+
+            var result = _function(point);
+            Value = result.Item1;
+            Gradient = result.Item2;
+        }
+
+        public Vector<double> Point { get; private set; }
+        public double Value { get; private set; }
+        public Vector<double> Gradient { get; private set; }
+
+        public Matrix<double> Hessian
+        {
+            get { throw new NotSupportedException(); }
+        }
+    }
+}

src/Numerics/Optimization/ObjectiveFunctions/HessianObjectiveFunction.cs

@@ -0,0 +1,59 @@
+using System;
+using MathNet.Numerics.LinearAlgebra;
+
+namespace MathNet.Numerics.Optimization.ObjectiveFunctions
+{
+    internal class HessianObjectiveFunction : IObjectiveFunction
+    {
+        readonly Func<Vector<double>, Tuple<double, Matrix<double>>> _function;
+
+        public HessianObjectiveFunction(Func<Vector<double>, Tuple<double, Matrix<double>>> function)
+        {
+            _function = function;
+        }
+
+        public IObjectiveFunction CreateNew()
+        {
+            return new HessianObjectiveFunction(_function);
+        }
+
+        public IObjectiveFunction Fork()
+        {
+            // no need to deep-clone values since they are replaced on evaluation
+            return new HessianObjectiveFunction(_function)
+            {
+                Point = Point,
+                Value = Value,
+                Hessian = Hessian
+            };
+        }
+
+        public bool IsGradientSupported
+        {
+            get { return false; }
+        }
+
+        public bool IsHessianSupported
+        {
+            get { return true; }
+        }
+
+        public void EvaluateAt(Vector<double> point)
+        {
+            Point = point;
+
+            var result = _function(point);
+            Value = result.Item1;
+            Hessian = result.Item2;
+        }
+
+        public Vector<double> Point { get; private set; }
+        public double Value { get; private set; }
+        public Matrix<double> Hessian { get; private set; }
+
+        public Vector<double> Gradient
+        {
+            get { throw new NotSupportedException(); }
+        }
+    }
+}

src/Numerics/Optimization/ObjectiveFunctions/LazyObjectiveFunction.cs

@@ -0,0 +1,142 @@
+// <copyright file="LazyObjectiveFunction.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2016 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using System;
+using MathNet.Numerics.LinearAlgebra;
+
+namespace MathNet.Numerics.Optimization.ObjectiveFunctions
+{
+    internal class LazyObjectiveFunction : IObjectiveFunction
+    {
+        readonly Func<Vector<double>, double> _function;
+        readonly Func<Vector<double>, Vector<double>> _gradient;
+        readonly Func<Vector<double>, Matrix<double>> _hessian;
+
+        Vector<double> _point;
+
+        bool _hasFunctionValue;
+        double _functionValue;
+
+        bool _hasGradientValue;
+        Vector<double> _gradientValue;
+
+        bool _hasHessianValue;
+        Matrix<double> _hessianValue;
+
+        public LazyObjectiveFunction(Func<Vector<double>, double> function, Func<Vector<double>, Vector<double>> gradient = null, Func<Vector<double>, Matrix<double>> hessian = null)
+        {
+            _function = function;
+            _gradient = gradient;
+            _hessian = hessian;
+
+            IsGradientSupported = gradient != null;
+            IsHessianSupported = hessian != null;
+        }
+
+        public IObjectiveFunction CreateNew()
+        {
+            return new LazyObjectiveFunction(_function, _gradient, _hessian);
+        }
+
+        public IObjectiveFunction Fork()
+        {
+            // no need to deep-clone values since they are replaced on evaluation
+            return new LazyObjectiveFunction(_function, _gradient, _hessian)
+            {
+                _point = _point,
+                _hasFunctionValue = _hasFunctionValue,
+                _functionValue = _functionValue,
+                _hasGradientValue = _hasGradientValue,
+                _gradientValue = _gradientValue,
+                _hasHessianValue = _hasHessianValue,
+                _hessianValue = _hessianValue
+            };
+        }
+
+        public bool IsGradientSupported { get; private set; }
+        public bool IsHessianSupported { get; private set; }
+
+        public void EvaluateAt(Vector<double> point)
+        {
+            _point = point;
+            _hasFunctionValue = false;
+            _hasGradientValue = false;
+            _hasHessianValue = false;
+
+            // don't keep references unnecessarily
+            _gradientValue = null;
+            _hessianValue = null;
+        }
+
+        public Vector<double> Point
+        {
+            get { return _point; }
+        }
+
+        public double Value
+        {
+            get
+            {
+                if (!_hasFunctionValue)
+                {
+                    _functionValue = _function(_point);
+                    _hasFunctionValue = true;
+                }
+                return _functionValue;
+            }
+        }
+
+        public Vector<double> Gradient
+        {
+            get
+            {
+                if (!_hasGradientValue)
+                {
+                    _gradientValue = _gradient(_point);
+                    _hasGradientValue = true;
+                }
+                return _gradientValue;
+            }
+        }
+
+        public Matrix<double> Hessian
+        {
+            get
+            {
+                if (!_hasHessianValue)
+                {
+                    _hessianValue = _hessian(_point);
+                    _hasHessianValue = true;
+                }
+                return _hessianValue;
+            }
+        }
+    }
+}

src/Numerics/Optimization/ObjectiveFunctions/LazyObjectiveFunctionBase.cs

@@ -0,0 +1,149 @@
+// <copyright file="LazyObjectiveFunctionBase.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2016 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using MathNet.Numerics.LinearAlgebra;
+
+namespace MathNet.Numerics.Optimization.ObjectiveFunctions
+{
+    public abstract class LazyObjectiveFunctionBase : IObjectiveFunction
+    {
+        Vector<double> _point;
+
+        protected bool HasFunctionValue { get; set; }
+        protected double FunctionValue { get; set; }
+
+        protected bool HasGradientValue { get; set; }
+        protected Vector<double> GradientValue { get; set; }
+
+        protected bool HasHessianValue { get; set; }
+        protected Matrix<double> HessianValue { get; set; }
+
+        protected LazyObjectiveFunctionBase(bool gradientSupported, bool hessianSupported)
+        {
+            IsGradientSupported = gradientSupported;
+            IsHessianSupported = hessianSupported;
+        }
+
+        public abstract IObjectiveFunction CreateNew();
+
+        public virtual IObjectiveFunction Fork()
+        {
+            // we need to deep-clone values since they may be updated inplace on evaluation
+            LazyObjectiveFunctionBase fork = (LazyObjectiveFunctionBase)CreateNew();
+            fork._point = _point?.Clone();
+            fork.HasFunctionValue = HasFunctionValue;
+            fork.FunctionValue = FunctionValue;
+            fork.HasGradientValue = HasGradientValue;
+            fork.GradientValue = GradientValue?.Clone();
+            fork.HasHessianValue = HasHessianValue;
+            fork.HessianValue = HessianValue?.Clone();
+            return fork;
+        }
+
+        public bool IsGradientSupported { get; private set; }
+        public bool IsHessianSupported { get; private set; }
+
+        public void EvaluateAt(Vector<double> point)
+        {
+            _point = point;
+            HasFunctionValue = false;
+            HasGradientValue = false;
+            HasHessianValue = false;
+        }
+
+        protected abstract void EvaluateValue();
+
+        protected virtual void EvaluateGradient()
+        {
+            Gradient = null;
+        }
+
+        protected virtual void EvaluateHessian()
+        {
+            Hessian = null;
+        }
+
+        public Vector<double> Point
+        {
+            get { return _point; }
+        }
+
+        public double Value
+        {
+            get
+            {
+                if (!HasFunctionValue)
+                {
+                    EvaluateValue();
+                }
+                return FunctionValue;
+            }
+            protected set
+            {
+                FunctionValue = value;
+                HasFunctionValue = true;
+            }
+        }
+
+        public Vector<double> Gradient
+        {
+            get
+            {
+                if (!HasGradientValue)
+                {
+                    EvaluateGradient();
+                }
+                return GradientValue;
+            }
+            protected set
+            {
+                GradientValue = value;
+                HasGradientValue = true;
+            }
+        }
+
+        public Matrix<double> Hessian
+        {
+            get
+            {
+                if (!HasHessianValue)
+                {
+                    EvaluateHessian();
+                }
+                return HessianValue;
+            }
+            protected set
+            {
+                HessianValue = value;
+                HasHessianValue = true;
+            }
+        }
+    }
+}

src/Numerics/Optimization/ObjectiveFunctions/ObjectiveFunctionBase.cs

@@ -0,0 +1,42 @@
+using MathNet.Numerics.LinearAlgebra;
+
+namespace MathNet.Numerics.Optimization.ObjectiveFunctions
+{
+    public abstract class ObjectiveFunctionBase : IObjectiveFunction
+    {
+        protected ObjectiveFunctionBase(bool isGradientSupported, bool isHessianSupported)
+        {
+            IsGradientSupported = isGradientSupported;
+            IsHessianSupported = isHessianSupported;
+        }
+
+        public abstract IObjectiveFunction CreateNew();
+
+        public virtual IObjectiveFunction Fork()
+        {
+            // we need to deep-clone values since they may be updated inplace on evaluation
+            ObjectiveFunctionBase objective = (ObjectiveFunctionBase)CreateNew();
+            objective.Point = Point == null ? null : Point.Clone();
+            objective.Value = Value;
+            objective.Gradient = Gradient == null ? null : Gradient.Clone();
+            objective.Hessian = Hessian == null ? null : Hessian.Clone();
+            return objective;
+        }
+
+        public bool IsGradientSupported { get; private set; }
+        public bool IsHessianSupported { get; private set; }
+
+        public void EvaluateAt(Vector<double> point)
+        {
+            Point = point;
+            Evaluate();
+        }
+
+        protected abstract void Evaluate();
+
+        public Vector<double> Point { get; private set; }
+        public double Value { get; protected set; }
+        public Vector<double> Gradient { get; protected set; }
+        public Matrix<double> Hessian { get; protected set; }
+    }
+}

src/Numerics/Optimization/ObjectiveFunctions/ValueObjectiveFunction.cs

@@ -0,0 +1,59 @@
+using System;
+using MathNet.Numerics.LinearAlgebra;
+
+namespace MathNet.Numerics.Optimization.ObjectiveFunctions
+{
+    internal class ValueObjectiveFunction : IObjectiveFunction
+    {
+        readonly Func<Vector<double>, double> _function;
+
+        public ValueObjectiveFunction(Func<Vector<double>, double> function)
+        {
+            _function = function;
+        }
+
+        public IObjectiveFunction CreateNew()
+        {
+            return new ValueObjectiveFunction(_function);
+        }
+
+        public IObjectiveFunction Fork()
+        {
+            // no need to deep-clone values since they are replaced on evaluation
+            return new ValueObjectiveFunction(_function)
+            {
+                Point = Point,
+                Value = Value,
+            };
+        }
+
+        public bool IsGradientSupported
+        {
+            get { return false; }
+        }
+
+        public bool IsHessianSupported
+        {
+            get { return false; }
+        }
+
+        public void EvaluateAt(Vector<double> point)
+        {
+            Point = point;
+            Value = _function(point);
+        }
+
+        public Vector<double> Point { get; private set; }
+        public double Value { get; private set; }
+
+        public Matrix<double> Hessian
+        {
+            get { throw new NotSupportedException(); }
+        }
+
+        public Vector<double> Gradient
+        {
+            get { throw new NotSupportedException(); }
+        }
+    }
+}

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/QuadraticGradientProjectionSearch.cs

@@ -0,0 +1,99 @@
+using System;
+using System.Collections.Generic;
+using MathNet.Numerics.LinearAlgebra;
+
+namespace MathNet.Numerics.Optimization
+{
+    public static class QuadraticGradientProjectionSearch
+    {
+        public static GradientProjectionResult Search(Vector<double> x0, Vector<double> gradient, Matrix<double> hessian, Vector<double> lowerBound, Vector<double> upperBound)
+        {
+            List<bool> isFixed = new List<bool>(x0.Count);
+            List<double> breakpoint = new List<double>(x0.Count);
+            for (int ii = 0; ii < x0.Count; ++ii)
+            {
+                breakpoint.Add(0.0);
+                isFixed.Add(false);
+                if (gradient[ii] < 0)
+                    breakpoint[ii] = (x0[ii] - upperBound[ii]) / gradient[ii];
+                else if (gradient[ii] > 0)
+                    breakpoint[ii] = (x0[ii] - lowerBound[ii]) / gradient[ii];
+                else
+                {
+                    if (Math.Abs(x0[ii] - upperBound[ii]) < 100 * Double.Epsilon || Math.Abs(x0[ii] - lowerBound[ii]) < 100 * Double.Epsilon)
+                        breakpoint[ii] = 0.0;
+                    else
+                        breakpoint[ii] = Double.PositiveInfinity;
+                }
+            }
+
+            var orderedBreakpoint = new List<double>(x0.Count);
+            orderedBreakpoint.AddRange(breakpoint);
+            orderedBreakpoint.Sort();
+
+            // Compute initial state variables
+            var d = -gradient;
+            for (int ii = 0; ii < d.Count; ++ii)
+                if (breakpoint[ii] <= 0.0)
+                    d[ii] *= 0.0;
+
+
+            int jj = -1;
+            var x = x0;
+            var f1 = gradient * d;
+            var f2 = 0.5 * d * hessian * d;
+            var sMin = -f1 / f2;
+            var maxS = orderedBreakpoint[0];
+
+            if (sMin < maxS)
+                return new GradientProjectionResult(x + sMin * d, 0,isFixed);
+
+            // while minimum of the last quadratic piece observed is beyond the interval searched
+            while (true)
+            {
+                // update data to the beginning of the interval we're searching
+                jj += 1;
+                x = x + d * maxS;
+                maxS = orderedBreakpoint[jj+1] - orderedBreakpoint[jj];
+
+                int fixedCount = 0;
+                for (int ii = 0; ii < d.Count; ++ii)
+                    if (orderedBreakpoint[jj] >= breakpoint[ii])
+                    {
+                        d[ii] *= 0.0;
+                        isFixed[ii] = true;
+                        fixedCount += 1;
+                    }
+
+                if (Double.IsPositiveInfinity(orderedBreakpoint[jj + 1]))
+                    return new GradientProjectionResult(x, fixedCount, isFixed);
+
+                f1 = gradient * d + (x - x0) * hessian * d;
+                f2 = d * hessian * d;
+
+                sMin = -f1 / f2;
+
+                if (sMin < maxS)
+                    return new GradientProjectionResult(x + sMin * d, fixedCount, isFixed);
+                else if (jj + 1 >= orderedBreakpoint.Count - 1)
+                {
+                    isFixed[isFixed.Count - 1] = true;
+                    return new GradientProjectionResult(x + maxS * d, lowerBound.Count, isFixed);
+                }
+            }
+        }
+
+        public struct GradientProjectionResult
+        {
+            public GradientProjectionResult(Vector<double> cauchyPoint, int fixedCount, List<bool> isFixed)
+            {
+                CauchyPoint = cauchyPoint;
+                FixedCount = fixedCount;
+                IsFixed = isFixed;
+            }
+            public Vector<double> CauchyPoint { get; }
+            public int FixedCount { get; }
+            public List<bool> IsFixed { get; }
+        }
+    }
+}

src/TestData/TestData.csproj

@@ -26,7 +26,7 @@
     <CodeAnalysisRuleSet>AllRules.ruleset</CodeAnalysisRuleSet>
     <NoWarn>1591</NoWarn>
     <Prefer32Bit>false</Prefer32Bit>
-    <LangVersion>5</LangVersion>
+    <LangVersion>6</LangVersion>
   </PropertyGroup>
   <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Debug|AnyCPU' ">
     <DefineConstants>DEBUG;TRACE</DefineConstants>
@@ -41,7 +41,7 @@
     <PlatformTarget>AnyCPU</PlatformTarget>
     <NoWarn>1591</NoWarn>
     <Prefer32Bit>false</Prefer32Bit>
-    <LangVersion>5</LangVersion>
+    <LangVersion>6</LangVersion>
   </PropertyGroup>
   <PropertyGroup Condition="'$(Configuration)|$(Platform)' == 'Release-Signed|AnyCPU'">
     <OutputPath>..\..\out\test-signed\Net35\</OutputPath>
@@ -57,7 +57,7 @@
     <CodeAnalysisRuleSet>AllRules.ruleset</CodeAnalysisRuleSet>
     <NoWarn>1591</NoWarn>
     <Prefer32Bit>false</Prefer32Bit>
-    <LangVersion>5</LangVersion>
+    <LangVersion>6</LangVersion>
   </PropertyGroup>
   <ItemGroup>
     <Reference Include="System" />

src/UnitTests/OptimizationTests/BfgsBMinimizerTests.cs

@@ -0,0 +1,253 @@
+// <copyright file="BfgsTest.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2016 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using System;
+using MathNet.Numerics.LinearAlgebra.Double;
+using MathNet.Numerics.Optimization;
+using NUnit.Framework;
+using System.Linq;
+using System.Text;
+using System.Collections.Generic;
+using MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions;
+using System.Collections;
+using MathNet.Numerics.Optimization.ObjectiveFunctions;
+using NUnit.Framework.Interfaces;
+
+namespace MathNet.Numerics.UnitTests.OptimizationTests
+{
+	[TestFixture]
+	public class BfgsBMinimizerTests
+	{
+	    [Test]
+		public void FindMinimum_Rosenbrock_Easy()
+		{
+            var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
+			var solver = new BfgsBMinimizer (1e-5, 1e-5, 1e-5, maximumIterations: 1000);
+			var lowerBound = new DenseVector(new[]{ -5.0, -5.0 });
+			var upperBound = new DenseVector(new[]{ 5.0, 5.0 });
+			var initialGuess = new DenseVector(new[] { 1.2, 1.2 });
+
+			var result = solver.FindMinimum(obj, lowerBound, upperBound, initialGuess);
+
+			Assert.That(Math.Abs(result.MinimizingPoint[0] - RosenbrockFunction.Minimum[0]), Is.LessThan(1e-3));
+			Assert.That(Math.Abs(result.MinimizingPoint[1] - RosenbrockFunction.Minimum[1]), Is.LessThan(1e-3));
+		}
+
+		[Test]
+		public void FindMinimum_Rosenbrock_Hard()
+		{
+            var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
+			var solver = new BfgsBMinimizer (1e-5, 1e-5, 1e-5, maximumIterations: 1000);
+
+			var lowerBound = new DenseVector(new[]{ -5.0, -5.0 });
+			var upperBound = new DenseVector(new[]{ 5.0, 5.0 });
+
+			var initialGuess = new DenseVector (new[]{ -1.2, 1.0 });
+
+			var result = solver.FindMinimum(obj, lowerBound, upperBound, initialGuess);
+
+			Assert.That(Math.Abs(result.MinimizingPoint[0] - RosenbrockFunction.Minimum[0]), Is.LessThan(1e-3));
+			Assert.That(Math.Abs(result.MinimizingPoint[1] - RosenbrockFunction.Minimum[1]), Is.LessThan(1e-3));
+		}
+
+		[Test]
+		public void FindMinimum_Rosenbrock_Overton()
+		{
+            var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
+			var solver = new BfgsBMinimizer (1e-5, 1e-5, 1e-5, maximumIterations: 1000);
+
+			var lowerBound = new DenseVector(new[]{ -5.0, -5.0 });
+			var upperBound = new DenseVector(new[]{ 5.0, 5.0 });
+			var initialGuess = new DenseVector (new[]{ -0.9, -0.5 });
+
+			var result = solver.FindMinimum (obj, lowerBound, upperBound, initialGuess);
+
+			Assert.That(Math.Abs(result.MinimizingPoint[0] - RosenbrockFunction.Minimum[0]), Is.LessThan(1e-3));
+			Assert.That(Math.Abs(result.MinimizingPoint[1] - RosenbrockFunction.Minimum[1]), Is.LessThan(1e-3));
+		}
+
+		[Test]
+		public void FindMinimum_Rosenbrock_Easy_OneBoundary()
+		{
+            var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
+			var solver = new BfgsBMinimizer (1e-5, 1e-5, 1e-5, maximumIterations: 1000);
+			var lowerBound = new DenseVector(new[]{ 1.0, -5.0 });
+			var upperBound = new DenseVector(new[]{ 5.0, 5.0 });
+			var initialGuess = new DenseVector(new[] { 1.2, 1.2 });
+
+			var result = solver.FindMinimum(obj, lowerBound, upperBound, initialGuess);
+
+			Assert.That(Math.Abs(result.MinimizingPoint[0] - RosenbrockFunction.Minimum[0]), Is.LessThan(1e-3));
+			Assert.That(Math.Abs(result.MinimizingPoint[1] - RosenbrockFunction.Minimum[1]), Is.LessThan(1e-3));
+		}
+
+		[Test]
+		public void FindMinimum_Rosenbrock_Easy_TwoBoundaries()
+		{
+            var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
+			var solver = new BfgsBMinimizer (1e-5, 1e-5, 1e-5, maximumIterations: 1000);
+			var lowerBound = new DenseVector(new[]{ 1.0, 1.0 });
+			var upperBound = new DenseVector(new[]{ 5.0, 5.0 });
+			var initialGuess = new DenseVector(new[] { 1.2, 1.2 });
+
+			var result = solver.FindMinimum(obj, lowerBound, upperBound, initialGuess);
+
+			Assert.That(Math.Abs(result.MinimizingPoint[0] - RosenbrockFunction.Minimum[0]), Is.LessThan(1e-3));
+			Assert.That(Math.Abs(result.MinimizingPoint[1] - RosenbrockFunction.Minimum[1]), Is.LessThan(1e-3));
+		}
+
+        [Test]
+        public void FindMinimum_Rosenbrock_MinimumGreateerOrEqualToLowerBoundary()
+        {
+            var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
+            var solver = new BfgsBMinimizer(1e-5, 1e-5, 1e-5, maximumIterations: 1000);
+
+            var lowerBound = new DenseVector(new[] { 2, 2.0 });
+            var upperBound = new DenseVector(new[] { 5.0, 5.0 });
+
+            var initialGuess = new DenseVector(new[] { 2.5, 2.5 });
+
+            var result = solver.FindMinimum(obj, lowerBound, upperBound, initialGuess);
+
+            Assert.GreaterOrEqual(result.MinimizingPoint[0],lowerBound[0]);
+            Assert.GreaterOrEqual(result.MinimizingPoint[1], lowerBound[1]);
+        }
+
+        [Test]
+        public void FindMinimum_Rosenbrock_MinimumLesserOrEqualToUpperBoundary()
+        {
+            var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
+            var solver = new BfgsBMinimizer(1e-5, 1e-5, 1e-5, maximumIterations: 1000);
+
+            var lowerBound = new DenseVector(new[] { -2.0, -2.0 });
+            var upperBound = new DenseVector(new[] { 0.5, 0.5 });
+
+            var initialGuess = new DenseVector(new[] { -0.9, -0.5 });
+
+            var result = solver.FindMinimum(obj, lowerBound, upperBound, initialGuess);
+
+            Assert.LessOrEqual(result.MinimizingPoint[0],upperBound[0]);
+            Assert.LessOrEqual(result.MinimizingPoint[1],upperBound[1]);
+        }
+
+        [Test]
+        [TestCaseSource(typeof(MghTestCaseEnumerator))]
+        public void Mgh_Tests(TestFunctions.TestCase test_case)
+        {
+            var obj = new MghObjectiveFunction(test_case.Function, true, true);
+            var solver = new BfgsBMinimizer(1e-8, 1e-8, 1e-8, 1000);
+
+            var result = solver.FindMinimum(obj, test_case.LowerBound, test_case.UpperBound, test_case.InitialGuess);
+
+            if (test_case.MinimizingPoint != null)
+            {
+                Assert.That((result.MinimizingPoint - test_case.MinimizingPoint).L2Norm(), Is.LessThan(1e-3));
+            }
+
+            var val1 = result.FunctionInfoAtMinimum.Value;
+            var val2 = test_case.MinimalValue;
+            var abs_min = Math.Min(Math.Abs(val1), Math.Abs(val2));
+            var abs_err = Math.Abs(val1 - val2);
+            var rel_err = abs_err / abs_min;
+            var success = (abs_min <= 1 && abs_err < 1e-3) || (abs_min > 1 && rel_err < 1e-3);
+            Assert.That(success, "Minimal function value is not as expected.");
+        }
+
+        [Test]
+        [TestCaseSource(typeof(FdMghTestCaseEnumerator))]
+        public void Mgh_FiniteDifference_Tests(TestFunctions.TestCase test_case)
+        {
+            var obj1 = new MghObjectiveFunction(test_case.Function, true, true);
+            var obj = new ForwardDifferenceGradientObjectiveFunction(obj1, test_case.LowerBound, test_case.UpperBound, 1e-10, 1e-10);
+            var solver = new BfgsBMinimizer(1e-8, 1e-8, 1e-8, 1000);
+
+            var result = solver.FindMinimum(obj, test_case.LowerBound, test_case.UpperBound, test_case.InitialGuess);
+
+            if (test_case.MinimizingPoint != null)
+            {
+                Assert.That((result.MinimizingPoint - test_case.MinimizingPoint).L2Norm(), Is.LessThan(1e-3));
+            }
+
+            var val1 = result.FunctionInfoAtMinimum.Value;
+            var val2 = test_case.MinimalValue;
+            var abs_min = Math.Min(Math.Abs(val1), Math.Abs(val2));
+            var abs_err = Math.Abs(val1 - val2);
+            var rel_err = abs_err / abs_min;
+            var success = (abs_min <= 1 && abs_err < 1e-3) || (abs_min > 1 && rel_err < 1e-3);
+            Assert.That(success, "Minimal function value is not as expected.");
+        }
+
+
+        private class BaseMghTestCaseEnumerator : IEnumerable<ITestCaseData>
+        {
+            private string _prefix = "";
+
+            public BaseMghTestCaseEnumerator(string prefix)
+            {
+                if (prefix.EndsWith(" "))
+                    _prefix = prefix;
+                else
+                    _prefix = prefix + " ";
+            }
+            public IEnumerator<ITestCaseData> GetEnumerator()
+            {
+                return 
+                    RosenbrockFunction2.TestCases
+                        .Concat(BealeFunction.TestCases)
+                        .Concat(HelicalValleyFunction.TestCases)
+                        .Concat(MeyerFunction.TestCases)
+                        .Concat(PowellSingularFunction.TestCases)
+                        .Concat(WoodFunction.TestCases)
+                        .Concat(BrownAndDennisFunction.TestCases)
+                    .Where(x => x.IsBounded)
+                    .Select(x => new TestCaseData(x)
+                        .SetName(_prefix + x.FullName)
+                    )
+                    .GetEnumerator();
+            }
+
+            IEnumerator IEnumerable.GetEnumerator()
+            {
+                return this.GetEnumerator();
+            }
+        }
+
+        private class MghTestCaseEnumerator : BaseMghTestCaseEnumerator
+        {
+            public MghTestCaseEnumerator() : base("") { }
+        }
+
+        private class FdMghTestCaseEnumerator : BaseMghTestCaseEnumerator
+        {
+            public FdMghTestCaseEnumerator() : base("FD") { }
+        }
+    }
+}
+

src/UnitTests/OptimizationTests/BfgsMinimizerTests.cs

@@ -0,0 +1,131 @@
+using System;
+using System.Linq;
+using MathNet.Numerics.LinearAlgebra.Double;
+using MathNet.Numerics.Optimization;
+using NUnit.Framework;
+using System.Text;
+using MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions;
+using System.Collections.Generic;
+using System.Collections;
+using NUnit.Framework.Interfaces;
+
+namespace MathNet.Numerics.UnitTests.OptimizationTests
+{
+    [TestFixture]
+    public class BfgsMinimizerTests
+    {
+        [Test]
+        public void FindMinimum_Rosenbrock_Easy()
+        {
+            var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
+            var solver = new BfgsMinimizer(1e-5, 1e-5, 1e-5, 1000);
+            var result = solver.FindMinimum(obj, new DenseVector(new[] { 1.2, 1.2 }));
+
+            Assert.That(Math.Abs(result.MinimizingPoint[0] - RosenbrockFunction.Minimum[0]), Is.LessThan(1e-3));
+            Assert.That(Math.Abs(result.MinimizingPoint[1] - RosenbrockFunction.Minimum[1]), Is.LessThan(1e-3));
+        }
+
+        [Test]
+        public void FindMinimum_Rosenbrock_Hard()
+        {
+            var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
+            var solver = new BfgsMinimizer(1e-5, 1e-5, 1e-5, 1000);
+            var result = solver.FindMinimum(obj, new DenseVector(new[] { -1.2, 1.0 }));
+
+            Assert.That(Math.Abs(result.MinimizingPoint[0] - RosenbrockFunction.Minimum[0]), Is.LessThan(1e-3));
+            Assert.That(Math.Abs(result.MinimizingPoint[1] - RosenbrockFunction.Minimum[1]), Is.LessThan(1e-3));
+        }
+
+        [Test]
+        public void FindMinimum_Rosenbrock_Overton()
+        {
+            var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
+            var solver = new BfgsMinimizer(1e-5, 1e-5, 1e-5, 1000);
+            var result = solver.FindMinimum(obj, new DenseVector(new[] { -0.9, -0.5 }));
+
+            Assert.That(Math.Abs(result.MinimizingPoint[0] - RosenbrockFunction.Minimum[0]), Is.LessThan(1e-3));
+            Assert.That(Math.Abs(result.MinimizingPoint[1] - RosenbrockFunction.Minimum[1]), Is.LessThan(1e-3));
+        }
+
+		[Test]
+		public void FindMinimum_BigRosenbrock_Easy()
+		{
+            var obj = ObjectiveFunction.Gradient(BigRosenbrockFunction.Value, BigRosenbrockFunction.Gradient);
+			var solver = new BfgsMinimizer(1e-10, 1e-5, 1e-5, 1000);
+			var result = solver.FindMinimum(obj, new DenseVector(new[] { 1.2*100.0, 1.2*100.0 }));
+
+			Assert.That(Math.Abs(result.MinimizingPoint[0] - BigRosenbrockFunction.Minimum[0]), Is.LessThan(1e-3));
+			Assert.That(Math.Abs(result.MinimizingPoint[1] - BigRosenbrockFunction.Minimum[1]), Is.LessThan(1e-3));
+		}
+
+		[Test]
+		public void FindMinimum_BigRosenbrock_Hard()
+		{
+            var obj = ObjectiveFunction.Gradient(BigRosenbrockFunction.Value, BigRosenbrockFunction.Gradient);
+			var solver = new BfgsMinimizer(1e-5, 1e-5, 1e-5, 1000);
+			var result = solver.FindMinimum(obj, new DenseVector(new[] { -1.2*100.0, 1.0*100.0 }));
+
+			Assert.That(Math.Abs(result.MinimizingPoint[0] - BigRosenbrockFunction.Minimum[0]), Is.LessThan(1e-3));
+			Assert.That(Math.Abs(result.MinimizingPoint[1] - BigRosenbrockFunction.Minimum[1]), Is.LessThan(1e-3));
+		}
+
+		[Test]
+		public void FindMinimum_BigRosenbrock_Overton()
+		{
+            var obj = ObjectiveFunction.Gradient(BigRosenbrockFunction.Value, BigRosenbrockFunction.Gradient);
+			var solver = new BfgsMinimizer(1e-5, 1e-5, 1e-5, 1000);
+			var result = solver.FindMinimum(obj, new DenseVector(new[] { -0.9*100.0, -0.5*100.0 }));
+
+			Assert.That(Math.Abs(result.MinimizingPoint[0] - BigRosenbrockFunction.Minimum[0]), Is.LessThan(1e-3));
+			Assert.That(Math.Abs(result.MinimizingPoint[1] - BigRosenbrockFunction.Minimum[1]), Is.LessThan(1e-3));
+		}
+
+        private class MghTestCaseEnumerator : IEnumerable<ITestCaseData>
+        {
+            public IEnumerator<ITestCaseData> GetEnumerator()
+            {
+                return
+                    RosenbrockFunction2.TestCases
+                        .Concat(BealeFunction.TestCases)
+                        .Concat(HelicalValleyFunction.TestCases)
+                        .Concat(MeyerFunction.TestCases)
+                        .Concat(PowellSingularFunction.TestCases)
+                        .Concat(WoodFunction.TestCases)
+                        .Concat(BrownAndDennisFunction.TestCases)
+                    .Where(x => x.IsUnbounded)
+                    .Select(x => new TestCaseData(x)
+                        .SetName(x.FullName)
+                    )
+                    .GetEnumerator();
+            }
+
+            IEnumerator IEnumerable.GetEnumerator()
+            {
+                return this.GetEnumerator();
+            }
+        }
+
+        [Test]
+        [TestCaseSource(typeof(MghTestCaseEnumerator))]
+        public void Mgh_Tests(TestFunctions.TestCase test_case)
+        {
+            var obj = new MghObjectiveFunction(test_case.Function, true, true);
+            var solver = new BfgsMinimizer(1e-8, 1e-8, 1e-8, 1000);
+
+            var result = solver.FindMinimum(obj, test_case.InitialGuess);
+
+            if (test_case.MinimizingPoint != null)
+            {
+                Assert.That((result.MinimizingPoint - test_case.MinimizingPoint).L2Norm(), Is.LessThan(1e-3));
+            }
+
+            var val1 = result.FunctionInfoAtMinimum.Value;
+            var val2 = test_case.MinimalValue;
+            var abs_min = Math.Min(Math.Abs(val1), Math.Abs(val2));
+            var abs_err = Math.Abs(val1 - val2);
+            var rel_err = abs_err / abs_min;
+            var success = (abs_min <= 1 && abs_err < 1e-3) || (abs_min > 1 && rel_err < 1e-3);
+            Assert.That(success, "Minimal function value is not as expected.");
+        }
+    }
+}

src/UnitTests/OptimizationTests/BfgsTest.cs

@@ -0,0 +1,56 @@
+// <copyright file="BfgsTest.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2016 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using MathNet.Numerics.LinearAlgebra.Double;
+using MathNet.Numerics.Optimization;
+using NUnit.Framework;
+
+namespace MathNet.Numerics.UnitTests.OptimizationTests
+{
+    [TestFixture, Category("RootFinding")]
+    internal class BfgsTest
+    {
+        private const double Precision = 1e-4;
+
+        [Test]
+        public void MinimizeRosenbrock()
+        {
+            CheckRosenbrock(15.0, 8.0, expectedMin: 0.0);
+            CheckRosenbrock(-1.2, 1.0, expectedMin: 0.0);
+            CheckRosenbrock(-1.2, 100.0, expectedMin: 0.0);
+        }
+
+        private static void CheckRosenbrock(double a, double b, double expectedMin)
+        {
+            var x = BfgsSolver.Solve(new DenseVector(new[] { a, b }), RosenbrockFunction.Value, RosenbrockFunction.Gradient);
+            Numerics.Precision.AlmostEqual(expectedMin, RosenbrockFunction.Value(x), Precision);
+        }
+    }
+}

src/UnitTests/OptimizationTests/ConjugateGradientMinimizerTests.cs

@@ -0,0 +1,101 @@
+using System;
+using MathNet.Numerics.LinearAlgebra.Double;
+using MathNet.Numerics.Optimization;
+using NUnit.Framework;
+using MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions;
+using System.Collections;
+using System.Collections.Generic;
+using System.Linq;
+using NUnit.Framework.Interfaces;
+
+namespace MathNet.Numerics.UnitTests.OptimizationTests
+{
+    [TestFixture]
+    public class ConjugateGradientMinimizerTests
+    {
+        [Test]
+        public void FindMinimum_Rosenbrock_Easy()
+        {
+            var obj = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
+            var solver = new ConjugateGradientMinimizer(1e-5, 1000);
+            var result = solver.FindMinimum(obj, new DenseVector(new[]{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 = ObjectiveFunction.Gradient(RosenbrockFunction.Value, RosenbrockFunction.Gradient);
+            var solver = new ConjugateGradientMinimizer(1e-5, 1000);
+            var result = solver.FindMinimum(obj, new DenseVector(new[] { -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));
+        }
+
+        private class MghTestCaseEnumerator : IEnumerable<ITestCaseData>
+        {
+            private static readonly string[] _ignore_list =
+            {
+                "Beale fun (MGH #5) unbounded",
+                "Meyer fun (MGH #10) unbounded",
+                "Powell singular fun (MGH #13) unbounded",
+                "Rosenbrock fun (MGH #1) hard start",
+                "Rosenbrock fun (MGH #1) Overton start",
+            };
+
+            private static bool in_ignore_list(string test_name)
+            {
+                return _ignore_list.Contains(test_name);
+            }
+
+            public IEnumerator<ITestCaseData> GetEnumerator()
+            {
+                return
+                    RosenbrockFunction2.TestCases
+                        .Concat(BealeFunction.TestCases)
+                        .Concat(HelicalValleyFunction.TestCases)
+                        .Concat(MeyerFunction.TestCases)
+                        .Concat(PowellSingularFunction.TestCases)
+                        .Concat(WoodFunction.TestCases)
+                        .Concat(BrownAndDennisFunction.TestCases)
+                    .Where(x => x.IsUnbounded)
+                    .Select(x => new TestCaseData(x)
+                        .SetName(x.FullName)
+                        .IgnoreIf(in_ignore_list(x.FullName),"Algo error, not implementation error.")
+                    )
+                    .GetEnumerator();
+            }
+
+            IEnumerator IEnumerable.GetEnumerator()
+            {
+                return this.GetEnumerator();
+            }
+        }
+
+        [Test]
+        [TestCaseSource(typeof(MghTestCaseEnumerator))]
+        public void Mgh_Tests(TestFunctions.TestCase test_case)
+        {
+            var obj = new MghObjectiveFunction(test_case.Function, true, true);
+            var solver = new ConjugateGradientMinimizer(1e-8, 1000);
+
+            var result = solver.FindMinimum(obj, test_case.InitialGuess);
+
+            if (test_case.MinimizingPoint != null)
+            {
+                Assert.That((result.MinimizingPoint - test_case.MinimizingPoint).L2Norm(), Is.LessThan(1e-3));
+            }
+
+            var val1 = result.FunctionInfoAtMinimum.Value;
+            var val2 = test_case.MinimalValue;
+            var abs_min = Math.Min(Math.Abs(val1), Math.Abs(val2));
+            var abs_err = Math.Abs(val1 - val2);
+            var rel_err = abs_err / abs_min;
+            var success = (abs_min <= 1 && abs_err < 1e-3) || (abs_min > 1 && rel_err < 1e-3);
+            Assert.That(success, "Minimal function value is not as expected.");
+        }
+    }
+}

src/UnitTests/OptimizationTests/GoldenSectionMinimizerTests.cs

@@ -0,0 +1,32 @@
+using System;
+using MathNet.Numerics.Optimization;
+using NUnit.Framework;
+
+namespace MathNet.Numerics.UnitTests.OptimizationTests
+{
+    [TestFixture]
+    public class GoldenSectionMinimizerTests
+    {
+        [Test]
+        public void Test_Works()
+        {
+            var algorithm = new GoldenSectionMinimizer(1e-5, 1000);
+            var f1 = new Func<double, double>(x => (x - 3)*(x - 3));
+            var obj = new SimpleObjectiveFunction1D(f1);
+            var r1 = algorithm.FindMinimum(obj, -100, 100);
+
+            Assert.That(Math.Abs(r1.MinimizingPoint - 3.0), Is.LessThan(1e-4));
+        }
+
+        [Test]
+        public void Test_ExpansionWorks()
+        {
+            var algorithm = new GoldenSectionMinimizer(1e-5, 1000);
+            var f1 = new Func<double, double>(x => (x - 3)*(x - 3));
+            var obj = new SimpleObjectiveFunction1D(f1);
+            var r1 = algorithm.FindMinimum(obj, -5, 5);
+
+            Assert.That(Math.Abs(r1.MinimizingPoint - 3.0), Is.LessThan(1e-4));
+        }
+    }
+}

src/UnitTests/OptimizationTests/NelderMeadSimplexTests.cs

@@ -0,0 +1,133 @@
+// <copyright file="NelderMeadSimplexTests.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2016 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using MathNet.Numerics.LinearAlgebra.Double;
+using MathNet.Numerics.Optimization;
+using MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions;
+using NUnit.Framework;
+using System;
+using System.Collections;
+using System.Collections.Generic;
+using System.Linq;
+using NUnit.Framework.Interfaces;
+
+namespace MathNet.Numerics.UnitTests.OptimizationTests
+{
+    [TestFixture]
+    public class NelderMeadSimplexTests
+    {
+        [Test]
+        public void NMS_FindMinimum_Rosenbrock_Easy()
+        {
+            var obj = ObjectiveFunction.Value(RosenbrockFunction.Value);
+            var solver = new NelderMeadSimplex(1e-5, maximumIterations: 1000);
+            var initialGuess = new DenseVector(new[] { 1.2, 1.2 });
+
+            var result = solver.FindMinimum(obj, initialGuess);
+
+            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 NMS_FindMinimum_Rosenbrock_Hard()
+        {
+            var obj = ObjectiveFunction.Value(RosenbrockFunction.Value);
+            var solver = new NelderMeadSimplex(1e-5, maximumIterations: 1000);
+
+            var initialGuess = new DenseVector(new[] { -1.2, 1.0 });
+
+            var result = solver.FindMinimum(obj,initialGuess);
+
+            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));
+        }
+
+        private class MghTestCaseEnumerator : IEnumerable<ITestCaseData>
+        {
+            private static readonly string[] _ignore_list =
+            {
+                "Meyer fun (MGH #10) unbounded",
+            };
+
+            private static bool in_ignore_list(string test_name)
+            {
+                return _ignore_list.Contains(test_name);
+            }
+
+            public IEnumerator<ITestCaseData> GetEnumerator()
+            {
+                return
+                    RosenbrockFunction2.TestCases
+                        .Concat(BealeFunction.TestCases)
+                        .Concat(HelicalValleyFunction.TestCases)
+                        .Concat(MeyerFunction.TestCases)
+                        .Concat(PowellSingularFunction.TestCases)
+                        .Concat(WoodFunction.TestCases)
+                        .Concat(BrownAndDennisFunction.TestCases)
+                    .Where(x => x.IsUnbounded)
+                    .Select(x => new TestCaseData(x)
+                        .SetName(x.FullName)
+                        .IgnoreIf(in_ignore_list(x.FullName), "Algo error, not implementation error")
+                    )
+                    .GetEnumerator();
+            }
+
+            IEnumerator IEnumerable.GetEnumerator()
+            {
+                return this.GetEnumerator();
+            }
+        }
+
+        [Test]
+        [TestCaseSource(typeof(MghTestCaseEnumerator))]
+        public void Mgh_Tests(TestFunctions.TestCase test_case)
+        {
+            var obj = new MghObjectiveFunction(test_case.Function, true, true);
+            var solver = new NelderMeadSimplex(1e-8, 1000);
+
+            var result = solver.FindMinimum(obj, test_case.InitialGuess);
+
+            if (test_case.MinimizingPoint != null)
+            {
+                Assert.That((result.MinimizingPoint - test_case.MinimizingPoint).L2Norm(), Is.LessThan(1e-3));
+            }
+
+            var val1 = result.FunctionInfoAtMinimum.Value;
+            var val2 = test_case.MinimalValue;
+            var abs_min = Math.Min(Math.Abs(val1), Math.Abs(val2));
+            var abs_err = Math.Abs(val1 - val2);
+            var rel_err = abs_err / abs_min;
+            var success = (abs_min <= 1 && abs_err < 1e-3) || (abs_min > 1 && rel_err < 1e-3);
+            Assert.That(success, "Minimal function value is not as expected.");
+        }
+    }
+}

src/UnitTests/OptimizationTests/NewtonMinimizerTests.cs

@@ -0,0 +1,188 @@
+using System;
+using MathNet.Numerics.LinearAlgebra.Double;
+using MathNet.Numerics.Optimization;
+using MathNet.Numerics.Optimization.ObjectiveFunctions;
+using NUnit.Framework;
+using MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions;
+using System.Collections.Generic;
+using System.Collections;
+using System.Linq;
+using NUnit.Framework.Interfaces;
+
+namespace MathNet.Numerics.UnitTests.OptimizationTests
+{
+    public class LazyRosenbrockObjectiveFunction : LazyObjectiveFunctionBase
+    {
+        public LazyRosenbrockObjectiveFunction() : base(true, true) { }
+
+        public override IObjectiveFunction CreateNew()
+        {
+            return new LazyRosenbrockObjectiveFunction();
+        }
+
+        protected override void EvaluateValue()
+        {
+            Value = RosenbrockFunction.Value(Point);
+        }
+
+        protected override void EvaluateGradient()
+        {
+            Gradient = RosenbrockFunction.Gradient(Point);
+        }
+
+        protected override void EvaluateHessian()
+        {
+            Hessian = RosenbrockFunction.Hessian(Point);
+        }
+    }
+
+    public class RosenbrockObjectiveFunction : ObjectiveFunctionBase
+    {
+        public RosenbrockObjectiveFunction() : base(true, true) { }
+
+        public override IObjectiveFunction CreateNew()
+        {
+            return new RosenbrockObjectiveFunction();
+        }
+
+        protected override void Evaluate()
+        {
+            // here we could directly overwrite the existing matrix cells instead.
+            // note: values must then be initialized manually first, if null.
+            Value = RosenbrockFunction.Value(Point);
+            Gradient = RosenbrockFunction.Gradient(Point);
+            Hessian = RosenbrockFunction.Hessian(Point);
+        }
+    }
+
+    [TestFixture]
+    public class NewtonMinimizerTests
+    {
+        [Test]
+        public void FindMinimum_Rosenbrock_Easy()
+        {
+            var obj = ObjectiveFunction.GradientHessian(RosenbrockFunction.Value, RosenbrockFunction.Gradient, RosenbrockFunction.Hessian);
+            var solver = new NewtonMinimizer(1e-5, 1000);
+            var result = solver.FindMinimum(obj, new DenseVector(new[] { 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 = ObjectiveFunction.GradientHessian(point => Tuple.Create(RosenbrockFunction.Value(point), RosenbrockFunction.Gradient(point), RosenbrockFunction.Hessian(point)));
+            var solver = new NewtonMinimizer(1e-5, 1000);
+            var result = solver.FindMinimum(obj, new DenseVector(new[] { -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 LazyRosenbrockObjectiveFunction();
+            var solver = new NewtonMinimizer(1e-5, 1000);
+            var result = solver.FindMinimum(obj, new DenseVector(new[] { -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 RosenbrockObjectiveFunction();
+            var solver = new NewtonMinimizer(1e-5, 1000, true);
+            var result = solver.FindMinimum(obj, new DenseVector(new[] { 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 LazyRosenbrockObjectiveFunction();
+            var solver = new NewtonMinimizer(1e-5, 1000, true);
+            var result = solver.FindMinimum(obj, new DenseVector(new[] { -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 LazyRosenbrockObjectiveFunction();
+            var solver = new NewtonMinimizer(1e-5, 1000, true);
+            var result = solver.FindMinimum(obj, new DenseVector(new[] { -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));
+        }
+
+        private class MghTestCaseEnumerator : IEnumerable<ITestCaseData>
+        {
+            private static readonly string[] _ignore_list =
+            {
+                "Beale fun (MGH #5) unbounded",
+                "Meyer fun (MGH #10) unbounded",
+                "Wood fun (MGH #14) unbounded",
+            };
+
+            private static bool in_ignore_list(string test_name)
+            {
+                return _ignore_list.Contains(test_name);
+            }
+
+            public IEnumerator<ITestCaseData> GetEnumerator()
+            {
+                return
+                    RosenbrockFunction2.TestCases
+                        .Concat(BealeFunction.TestCases)
+                        .Concat(HelicalValleyFunction.TestCases)
+                        .Concat(MeyerFunction.TestCases)
+                        .Concat(PowellSingularFunction.TestCases)
+                        .Concat(WoodFunction.TestCases)
+                        .Concat(BrownAndDennisFunction.TestCases)
+                    .Where(x => x.IsUnbounded)
+                    .Select(x => new TestCaseData(x)
+                        .SetName(x.FullName)
+                        .IgnoreIf(in_ignore_list(x.FullName),"Algo error, not implementation error")
+                    )
+                    .GetEnumerator();
+            }
+
+            IEnumerator IEnumerable.GetEnumerator()
+            {
+                return this.GetEnumerator();
+            }
+        }
+
+        [Test]
+        [TestCaseSource(typeof(MghTestCaseEnumerator))]
+        public void Mgh_Tests(TestFunctions.TestCase test_case)
+        {
+            var obj = new MghObjectiveFunction(test_case.Function, true, true);
+            var solver = new NewtonMinimizer(1e-8, 1000, useLineSearch: false);
+
+            var result = solver.FindMinimum(obj, test_case.InitialGuess);
+
+            if (test_case.MinimizingPoint != null)
+            {
+                Assert.That((result.MinimizingPoint - test_case.MinimizingPoint).L2Norm(), Is.LessThan(1e-3));
+            }
+
+            var val1 = result.FunctionInfoAtMinimum.Value;
+            var val2 = test_case.MinimalValue;
+            var abs_min = Math.Min(Math.Abs(val1), Math.Abs(val2));
+            var abs_err = Math.Abs(val1 - val2);
+            var rel_err = abs_err / abs_min;
+            var success = (abs_min <= 1 && abs_err < 1e-3) || (abs_min > 1 && rel_err < 1e-3);
+            Assert.That(success, "Minimal function value is not as expected.");
+        }
+    }
+}

src/UnitTests/OptimizationTests/RosenbrockFunction.cs

@@ -0,0 +1,66 @@
+using System;
+using MathNet.Numerics.LinearAlgebra;
+using MathNet.Numerics.LinearAlgebra.Double;
+
+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 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 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 Vector<double> Minimum
+        {
+            get
+            {
+                return new DenseVector(new double[] { 1, 1 });
+            }
+        }
+    }
+
+    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);
+        }
+
+        public static Vector<double> Minimum
+        {
+            get
+            {
+                return new DenseVector(new double[] { 100, 100 });
+            }
+        }
+    }
+}

src/UnitTests/OptimizationTests/RosenbrockFunctionTests.cs

@@ -0,0 +1,55 @@
+using System;
+using MathNet.Numerics.LinearAlgebra.Double;
+using NUnit.Framework;
+
+namespace MathNet.Numerics.UnitTests.OptimizationTests
+{
+    [TestFixture]
+    class RosenbrockFunctionTests
+    {
+        [Test]
+        public void TestGradient()
+        {
+            var input = new DenseVector(new[]{ -0.9, -0.5 } );
+
+            var v1 = RosenbrockFunction.Value(input);
+            var g = RosenbrockFunction.Gradient(input);
+
+            var eps = 1e-5;
+            var eps0 = (new DenseVector(new[] { 1.0, 0.0 })) * eps;
+            var eps1 = (new DenseVector(new[] { 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 DenseVector(new[] { -0.9, -0.5 });
+
+            var v1 = RosenbrockFunction.Value(input);
+            var h = RosenbrockFunction.Hessian(input);
+
+            var eps = 1e-5;
+
+            var eps0 = (new DenseVector(new[] { 1.0, 0.0 })) * eps;
+            var eps1 = (new DenseVector(new[] { 0.0, 1.0 })) * eps;
+
+            var epsuu = (new DenseVector(new[] { 1.0, 1.0 })) * eps;
+            var epsud = (new DenseVector(new[] { 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);
+        }
+    }
+}

src/UnitTests/OptimizationTests/TestCaseDataExtensions.cs

@@ -0,0 +1,20 @@
+using NUnit.Framework;
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace MathNet.Numerics.UnitTests.OptimizationTests
+{
+    internal static class TestCaseDataExtensions
+    {
+        public static TestCaseData IgnoreIf(this TestCaseData input, bool do_ignore, string reason)
+        {
+            if (do_ignore)
+                return input.Ignore(reason);
+            else
+                return input;
+        }
+    }
+}

src/UnitTests/OptimizationTests/TestFunctionAdapters.cs

@@ -0,0 +1,54 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+using MathNet.Numerics.LinearAlgebra;
+using MathNet.Numerics.Optimization;
+using MathNet.Numerics.Optimization.ObjectiveFunctions;
+using MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions;
+using MathNet.Numerics.LinearAlgebra.Double;
+
+namespace MathNet.Numerics.UnitTests.OptimizationTests
+{
+    public class MghObjectiveFunction : LazyObjectiveFunctionBase
+    {
+        private ITestFunction TestFunction;
+
+        public MghObjectiveFunction(ITestFunction testFunction, bool use_gradient, bool use_hessian)
+            : base(use_gradient, use_hessian)
+        {
+            this.TestFunction = testFunction;
+        }
+
+        public override IObjectiveFunction CreateNew()
+        {
+            return new MghObjectiveFunction(this.TestFunction, this.IsGradientSupported, this.IsHessianSupported);
+        }
+
+        protected override void EvaluateValue()
+        {
+            this.Value = this.TestFunction.SsqValue(this.Point);
+        }
+
+        protected override void EvaluateGradient()
+        {
+            if (this.IsGradientSupported)
+            {
+                if (this.GradientValue == null)
+                    this.Gradient = new DenseVector(this.TestFunction.ParameterDimension);
+                this.TestFunction.SsqGradientByRef(this.Point, GradientValue);
+            }
+        }
+
+        protected override void EvaluateHessian()
+        {
+            if (this.IsHessianSupported)
+            {
+                if (this.HessianValue == null)
+                    this.Hessian = new DenseMatrix(this.TestFunction.ParameterDimension, this.TestFunction.ParameterDimension);
+                this.TestFunction.SsqHessianByRef(this.Point, HessianValue);
+            }
+        }
+    }
+}

src/UnitTests/OptimizationTests/TestFunctionTests.cs

@@ -0,0 +1,309 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+using MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions;
+using NUnit.Framework;
+using MathNet.Numerics.LinearAlgebra;
+using System.Collections;
+
+namespace MathNet.Numerics.UnitTests.OptimizationTests
+{
+    [TestFixture]
+    public class TestFunctionTests
+    {
+        private static IEnumerable<TestFunctions.TestCase> MghCases
+        {
+            get
+            {
+                return Enumerable.Empty<TestFunctions.TestCase>()
+                    .Concat(RosenbrockFunction2.TestCases)
+                    .Concat(BealeFunction.TestCases)
+                    .Concat(HelicalValleyFunction.TestCases)
+                    .Concat(MeyerFunction.TestCases)
+                    .Concat(PowellSingularFunction.TestCases)
+                    .Concat(WoodFunction.TestCases)
+                    .Concat(BrownAndDennisFunction.TestCases);
+            }
+        }
+
+        private class MghCaseEnumerator : IEnumerable<TestCaseData>
+        {
+            public string CategoryName { get; protected set; }
+
+            public MghCaseEnumerator(string category_name)
+            {
+                this.CategoryName = category_name;
+            }
+
+            public virtual IEnumerator<TestCaseData> GetEnumerator()
+            {
+                return MghCases
+                    .Select(x =>
+                        new TestCaseData(x)
+                            .SetName($"{x.FullName} {this.CategoryName}")
+                    ).GetEnumerator();
+            }
+
+            IEnumerator IEnumerable.GetEnumerator()
+            {
+                return this.GetEnumerator();
+            }
+        }
+
+        [Test]
+        public void Smoke_Construction()
+        {
+            var c = new TestCase()
+            {
+                InitialGuess = new double[] { 1, 2, 3 },
+                MinimizingPoint = new double[] { 1, 1, 1 },
+                MinimalValue = 0
+            };
+        }
+
+        private class ValueAtMinimumSource : MghCaseEnumerator
+        {
+            public ValueAtMinimumSource() : base("ValueAtMinimum") { }
+
+            public override IEnumerator<TestCaseData> GetEnumerator()
+            {
+                return MghCases
+                    .Where(x => x.MinimizingPoint != null)
+                    .Select(x =>
+                        new TestCaseData(x)
+                            .SetName($"{x.FullName} {this.CategoryName}")
+                    )
+                    .GetEnumerator();
+            }
+        }
+
+        [Test]
+        [TestCaseSource(typeof(ValueAtMinimumSource))]
+        public void ValueAtMinimum(TestFunctions.TestCase test_case)
+        {
+            if (test_case.MinimizingPoint != null)
+            {
+                var value_at_minimum = test_case.Function.SsqValue(test_case.MinimizingPoint);
+                Assert.That(
+                    Math.Abs(value_at_minimum - test_case.MinimalValue) < 1e-3,
+                    $"Function value at minimum not as expected."
+                );
+            }
+        }
+
+        private class GradientAtStartSource : MghCaseEnumerator
+        {
+            public GradientAtStartSource() : base("GradientAtStart") { }
+        }
+
+        [Test]
+        [TestCaseSource(typeof(GradientAtStartSource))]
+        public void GradientAtStart(TestFunctions.TestCase test_case)
+        {
+            var a_grad = test_case.Function.SsqGradient(test_case.InitialGuess);
+            var fd_grad = Vector<double>.Build.Dense(test_case.Function.ParameterDimension, 0.0);
+
+            for (int ii = 0; ii < test_case.Function.ParameterDimension; ++ii)
+            {
+                var h = 1e-6;
+
+                var bump_up = test_case.InitialGuess.Clone();
+                bump_up[ii] += h;
+                var bump_down = test_case.InitialGuess.Clone();
+                bump_down[ii] -= h;
+
+                var up_val = test_case.Function.SsqValue(bump_up);
+                var down_val = test_case.Function.SsqValue(bump_down);
+
+                fd_grad[ii] = 0.5 * (up_val - down_val) / h;
+            }
+
+            for (int ii = 0; ii < test_case.Function.ParameterDimension; ++ii)
+            {
+                var val1 = a_grad[ii];
+                var val2 = fd_grad[ii];
+                var min_abs_val = Math.Min(Math.Abs(val1), Math.Abs(val2));
+                if (min_abs_val <= 1)
+                    Assert.That(Math.Abs(val1 - val2) < 1e-3, $"Problem with gradient value at start point.");
+                else
+                    Assert.That(Math.Abs(val1 - val2) / min_abs_val < 1e-3, $"Problem with gradient value at start point.");
+            }
+        }
+
+        private class HessianAtStartSource : MghCaseEnumerator
+        {
+            public HessianAtStartSource() : base("HessianAtStart") { }
+
+            public override IEnumerator<TestCaseData> GetEnumerator()
+            {
+                return MghCases
+                    .Where(x => x.MinimizingPoint != null)
+                    .Select(x =>
+                        new TestCaseData(x)
+                            .SetName($"{x.FullName} {this.CategoryName}")
+                    )
+                    .GetEnumerator();
+            }
+        }
+        [Test]
+        [TestCaseSource(typeof(HessianAtStartSource))]
+        public void HessianAtStart(TestFunctions.TestCase test_case)
+        {
+            var a_hess = test_case.Function.SsqHessian(test_case.InitialGuess);
+            var fd_hess = Matrix<double>.Build.Dense(test_case.Function.ParameterDimension, test_case.Function.ParameterDimension);
+            for (int ii = 0; ii < test_case.Function.ParameterDimension; ++ii)
+            {
+                for (int jj = 0; jj < test_case.Function.ParameterDimension; ++jj)
+                {
+                    var h1 = 1e-3 * Math.Max(1.0, Math.Abs(test_case.InitialGuess[ii]));
+                    var h2 = 1e-3 * Math.Max(1.0, Math.Abs(test_case.InitialGuess[jj]));
+
+                    var bump_uu = test_case.InitialGuess.Clone();
+                    bump_uu[ii] += h1;
+                    bump_uu[jj] += h2;
+
+                    var bump_dd = test_case.InitialGuess.Clone();
+                    bump_dd[ii] -= h1;
+                    bump_dd[jj] -= h2;
+
+                    var bump_ud = test_case.InitialGuess.Clone();
+                    bump_ud[ii] += h1;
+                    bump_ud[jj] -= h2;
+
+                    var bump_du = test_case.InitialGuess.Clone();
+                    bump_du[ii] -= h1;
+                    bump_du[jj] += h2;
+
+                    var val_uu = test_case.Function.SsqValue(bump_uu);
+                    var val_dd = test_case.Function.SsqValue(bump_dd);
+                    var val_ud = test_case.Function.SsqValue(bump_ud);
+                    var val_du = test_case.Function.SsqValue(bump_du);
+
+                    fd_hess[ii, jj] = (val_uu - val_ud + val_dd - val_du) / (4 * h1 * h2);
+                }
+            }
+
+            for (int ii = 0; ii < test_case.Function.ParameterDimension; ++ii)
+            {
+                for (int jj = 0; jj < test_case.Function.ParameterDimension; ++jj)
+                {
+                    var val1 = fd_hess[ii, jj];
+                    var val2 = a_hess[ii, jj];
+
+                    var abs_min = Math.Min(Math.Abs(val1), Math.Abs(val2));
+                    if (abs_min <= 1)
+                    {
+                        Assert.That(Math.Abs(val1 - val2) < 1e-3, $"Problem with hessian at start point.");
+                    }
+                    else
+                    {
+                        Assert.That(Math.Abs(val1 - val2) / abs_min < 0.05, $"Problem with hessian at start point.");
+                    }
+                }
+            }
+        }
+
+        private class ItemGradientAtStartSource : MghCaseEnumerator
+        {
+            public ItemGradientAtStartSource() : base("ItemGradientAtStart") { }
+        }
+
+        [Test]
+        [TestCaseSource(typeof(ItemGradientAtStartSource))]
+        public void ItemGradientAtStart(TestFunctions.TestCase test_case)
+        {
+            for (var item_index = 0; item_index < test_case.Function.ItemDimension; ++item_index)
+            {
+
+                var a_grad = test_case.Function.ItemGradient(test_case.InitialGuess, item_index);
+                var h = 1e-4;
+                var fd_grad = Vector<double>.Build.Dense(test_case.Function.ParameterDimension, 0.0);
+
+                for (int ii = 0; ii < test_case.Function.ParameterDimension; ++ii)
+                {
+                    var bump_up = test_case.InitialGuess.Clone();
+                    bump_up[ii] += h;
+                    var bump_down = test_case.InitialGuess.Clone();
+                    bump_down[ii] -= h;
+
+                    var up_val = test_case.Function.ItemValue(bump_up, item_index);
+                    var down_val = test_case.Function.ItemValue(bump_down, item_index);
+
+                    fd_grad[ii] = 0.5 * (up_val - down_val) / h;
+                }
+
+                for (int ii = 0; ii < test_case.Function.ParameterDimension; ++ii)
+                {
+                    Assert.That(Math.Abs(fd_grad[ii] - a_grad[ii]) < 1e-3, $"Failed for parameter {ii}");
+                }
+            }
+        }
+
+        private class ItemHessianAtStartSource : MghCaseEnumerator
+        {
+            public ItemHessianAtStartSource() : base("ItemHessianAtStart") { }
+        }
+
+        [Test]
+        [TestCaseSource(typeof(ItemHessianAtStartSource))]
+        public void ItemHessianAtStart(TestFunctions.TestCase test_case)
+        {
+            for (var item_index = 0; item_index < test_case.Function.ItemDimension; ++item_index)
+            {
+                var a_hess = test_case.Function.ItemHessian(test_case.InitialGuess, item_index);
+                var h = 1e-4;
+                var fd_hess = Matrix<double>.Build.Dense(test_case.Function.ParameterDimension, test_case.Function.ParameterDimension);
+                for (int ii = 0; ii < test_case.Function.ParameterDimension; ++ii)
+                {
+                    for (int jj = 0; jj < test_case.Function.ParameterDimension; ++jj)
+                    {
+                        var bump_uu = test_case.InitialGuess.Clone();
+                        bump_uu[ii] += h;
+                        bump_uu[jj] += h;
+
+                        var bump_dd = test_case.InitialGuess.Clone();
+                        bump_dd[ii] -= h;
+                        bump_dd[jj] -= h;
+
+                        var bump_ud = test_case.InitialGuess.Clone();
+                        bump_ud[ii] += h;
+                        bump_ud[jj] -= h;
+
+                        var bump_du = test_case.InitialGuess.Clone();
+                        bump_du[ii] -= h;
+                        bump_du[jj] += h;
+
+                        var val_uu = test_case.Function.ItemValue(bump_uu, item_index);
+                        var val_dd = test_case.Function.ItemValue(bump_dd, item_index);
+                        var val_ud = test_case.Function.ItemValue(bump_ud, item_index);
+                        var val_du = test_case.Function.ItemValue(bump_du, item_index);
+
+                        fd_hess[ii, jj] = (val_uu - val_ud + val_dd - val_du) / (4 * h * h);
+                    }
+                }
+
+                for (int ii = 0; ii < test_case.Function.ParameterDimension; ++ii)
+                {
+                    for (int jj = 0; jj < test_case.Function.ParameterDimension; ++jj)
+                    {
+                        var val1 = fd_hess[ii, jj];
+                        var val2 = a_hess[ii, jj];
+
+                        var abs_min = Math.Min(Math.Abs(val1), Math.Abs(val2));
+                        if (abs_min <= 1)
+                        {
+                            Assert.That(Math.Abs(val1 - val2) < 1e-3, $"Problem with hessian at start point.");
+                        }
+                        else
+                        {
+                            Assert.That(Math.Abs(val1 - val2) / abs_min < 0.05, $"Problem with hessian at start point.");
+                        }
+                    }
+                }
+            }
+        }
+
+    }
+}

src/UnitTests/OptimizationTests/TestFunctions/BaseTestFunction.cs

@@ -0,0 +1,125 @@
+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.TestFunctions
+{
+    public abstract class BaseTestFunction : ITestFunction
+    {
+        public abstract string Description { get; }
+        public abstract int ParameterDimension { get; }
+        public abstract int ItemDimension { get; }
+
+        public abstract double ItemValue(Vector<double> x, int itemIndex);
+        public abstract void ItemGradientByRef(Vector<double> x, int itemIndex, Vector<double> output);
+        public abstract void ItemHessianByRef(Vector<double> x, int itemIndex, Matrix<double> output);
+
+        public virtual Vector<double> ItemGradient(Vector<double> x, int itemIndex)
+        {
+            var output = new LinearAlgebra.Double.DenseVector(this.ParameterDimension);
+            this.ItemGradientByRef(x, itemIndex, output);
+            return output;
+        }
+        
+        public virtual Matrix<double> ItemHessian(Vector<double> x, int itemIndex)
+        {
+            var output = new LinearAlgebra.Double.DenseMatrix(this.ParameterDimension, this.ParameterDimension);
+            this.ItemHessianByRef(x, itemIndex, output);
+            return output;
+        }
+
+
+        public virtual void JacobianbyRef(Vector<double> x, Matrix<double> output)
+        {
+            for (int ii = 0; ii < this.ItemDimension; ++ii)
+            {
+                var grad = this.ItemGradient(x, ii);
+                output.SetRow(ii, grad);
+            }
+        }
+
+        public virtual Matrix<double> Jacobian(Vector<double> x)
+        {
+            var output = new LinearAlgebra.Double.DenseMatrix(this.ItemDimension, this.ParameterDimension);
+            this.JacobianbyRef(x, output);
+            return output;
+        }
+
+        public virtual void SsqGradientByRef(Vector<double> x, Vector<double> output)
+        {
+            if (output.Count != this.ParameterDimension)
+                throw new ArgumentException($"Output vector must match parameter dimension of function; expected {this.ParameterDimension}, got {output.Count}.");
+
+            for (int jj = 0; jj < this.ParameterDimension; ++jj)
+                output[jj] = 0.0;
+            var tmp_grad = new LinearAlgebra.Double.DenseVector(this.ParameterDimension);
+            double tmp_value = 0.0;
+
+            for (int ii = 0; ii < this.ItemDimension; ++ii)
+            {
+                tmp_value = this.ItemValue(x, ii);
+                this.ItemGradientByRef(x, ii, tmp_grad);
+                for (int jj = 0; jj < this.ParameterDimension; ++jj)
+                    output[jj] += 2 * tmp_value * tmp_grad[jj];
+            }
+        }
+
+        public virtual Vector<double> SsqGradient(Vector<double> x)
+        {
+            var output = new LinearAlgebra.Double.DenseVector(this.ParameterDimension);
+            this.SsqGradientByRef(x, output);
+            return output;
+        }
+
+        public virtual void SsqHessianByRef(Vector<double> x, Matrix<double> output)
+        {
+            if (output.RowCount != this.ParameterDimension || output.ColumnCount != this.ParameterDimension)
+                throw new ArgumentException($"Output matrix must match parameter dimension of function; expected {this.ParameterDimension}x{this.ParameterDimension}, got {output.RowCount}x{output.ColumnCount}.");
+
+            for (int ii = 0; ii < this.ParameterDimension; ++ii)
+                for (int jj = 0; jj < this.ParameterDimension; ++jj)
+                    output[ii,jj] = 0.0;
+
+            var tmp_grad = new LinearAlgebra.Double.DenseVector(this.ParameterDimension);
+            var tmp_hess = new LinearAlgebra.Double.DenseMatrix(this.ParameterDimension, this.ParameterDimension);
+            double tmp_value = 0.0;
+
+            for (int ii = 0; ii < this.ItemDimension; ++ii)
+            {
+                tmp_value = this.ItemValue(x, ii);
+                this.ItemGradientByRef(x, ii, tmp_grad);
+                this.ItemHessianByRef(x, ii, tmp_hess);
+                for (int jj = 0; jj < this.ParameterDimension; ++jj)
+                {
+                    for (int kk = 0; kk < this.ParameterDimension; ++kk)
+                    {
+                        var increment = 2 * (tmp_value * tmp_hess[jj, kk] + tmp_grad[jj] * tmp_grad[kk]);
+                        output[jj, kk] += increment;
+                    }
+                }
+            }
+        }
+
+        public virtual Matrix<double> SsqHessian(Vector<double> x)
+        {
+            var output = new LinearAlgebra.Double.DenseMatrix(this.ParameterDimension, this.ParameterDimension);
+            this.SsqHessianByRef(x, output);
+            return output;
+        }
+
+        public virtual double SsqValue(Vector<double> x)
+        {
+            double ssq = 0.0;
+            for (int ii = 0; ii < this.ItemDimension; ++ii)
+            {
+                var tmp = this.ItemValue(x, ii);
+                ssq += tmp * tmp;
+            }
+            return ssq;
+        }
+
+    }
+}

src/UnitTests/OptimizationTests/TestFunctions/BealeFunction.cs

@@ -0,0 +1,97 @@
+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.TestFunctions
+{
+    public class BealeFunction : BaseTestFunction
+    {
+        public static IEnumerable<TestCase> TestCases
+        {
+            get
+            {
+                yield return new TestCase()
+                {
+                    Function = new BealeFunction(),
+                    InitialGuess = new double[] { 1, 1 },
+                    MinimalValue = 0,
+                    MinimizingPoint = new double[] { 3, 0.5 }, 
+                    CaseName = "unbounded"
+                };
+                yield return new TestCase()
+                {
+                    Function = new BealeFunction(),
+                    InitialGuess = new double[] { 1, 1 },
+                    MinimalValue = 0,
+                    MinimizingPoint = new double[] { 3, 0.5 },
+                    LowerBound = new double[] { -1000, -1000},
+                    UpperBound = new double[] { 1000, 1000},
+                    CaseName = "loose bounds"
+                };
+                yield return new TestCase()
+                {
+                    Function = new BealeFunction(),
+                    InitialGuess = new double[] { 1, 1 },
+                    MinimalValue = 0,
+                    MinimizingPoint = new double[] { 3, 0.5 },
+                    LowerBound = new double[] { 0.6, 0.5 },
+                    UpperBound = new double[] { 10, 100 },
+                    CaseName = "tight bounds"
+                };
+            }
+        }
+
+        public BealeFunction() { }
+
+        public override string Description
+        {
+            get
+            {
+                return "Beale fun (MGH #5)";
+            }
+        }
+
+        public override int ItemDimension
+        {
+            get
+            {
+                return 3;
+            }
+        }
+
+        public override int ParameterDimension
+        {
+            get
+            {
+                return 2;
+            }
+        }
+
+        public override void ItemGradientByRef(Vector<double> x, int itemIndex, Vector<double> output)
+        {
+            int ii = itemIndex + 1;
+            output[0] = -1 + Math.Pow(x[1], ii);
+            output[1] = ii * x[0] * Math.Pow(x[1], ii - 1);
+        }
+
+        public override void ItemHessianByRef(Vector<double> x, int itemIndex, Matrix<double> output)
+        {
+            int ii = itemIndex + 1;
+            output[0, 0] = 0;
+            output[0, 1] = ii * Math.Pow(x[1], ii - 1);
+            output[1, 0] = ii * Math.Pow(x[1], ii - 1);
+            output[1, 1] = (ii - 1) * ii * x[0] * Math.Pow(x[1], ii - 2);
+        }
+
+        private static readonly double[] y = { 1.5, 2.25, 2.625};
+
+        public override double ItemValue(Vector<double> x, int itemIndex)
+        {
+            int ii = itemIndex + 1;
+            return y[itemIndex] - x[0] * (1 - Math.Pow(x[1], ii));
+        }
+    }
+}

src/UnitTests/OptimizationTests/TestFunctions/BrownAndDennisFunction.cs

@@ -0,0 +1,114 @@
+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.TestFunctions
+{
+    public class BrownAndDennisFunction : BaseTestFunction
+    {
+        public static IEnumerable<TestCase> TestCases
+        {
+            get
+            {
+                yield return new TestCase()
+                {
+                    Function = new BrownAndDennisFunction(20),
+                    InitialGuess = new double[] { 25, 5, -5, -1 },
+                    MinimalValue = 85822.2,
+                    MinimizingPoint = null, 
+                    CaseName = "unbounded"
+                };
+                yield return new TestCase()
+                {
+                    Function = new BrownAndDennisFunction(20),
+                    InitialGuess = new double[] { 25, 5, -5, -1 },
+                    MinimalValue = 85822.2,
+                    MinimizingPoint = null,
+                    LowerBound = new double[] { -1000, -1000, -1000, -1000 },
+                    UpperBound = new double[] {1000, 1000, 1000, 1000 },
+                    CaseName = "loose bounds"                    
+                };
+                yield return new TestCase()
+                {
+                    Function = new BrownAndDennisFunction(20),
+                    InitialGuess = new double[] { 25, 5, -5, -1 },
+                    MinimalValue = 0.88860479e5,
+                    MinimizingPoint = null,
+                    LowerBound = new double[] { -10, 0, -100, -20 },
+                    UpperBound = new double[] { 100, 15, 0, 0.2 },
+                    CaseName = "tight bounds"
+                };
+            }
+        }
+        
+        private readonly int _items;
+
+        public BrownAndDennisFunction(int items)
+        {
+            if (items < 4)
+                throw new ArgumentException("items must be >= 4");
+            _items = items;
+        }
+
+        public override string Description
+        {
+            get
+            {
+                return "Brown & Dennis fun (MGH #16)";
+            }
+        }
+
+        public override int ItemDimension
+        {
+            get
+            {
+                return _items;
+            }
+        }
+
+        public override int ParameterDimension
+        {
+            get
+            {
+                return 4;
+            }
+        }
+
+        public override void ItemGradientByRef(Vector<double> x, int itemIndex, Vector<double> output)
+        {
+            var ii = itemIndex + 1;
+            var t = ii / 5.0;
+            output[0] = 2 * (x[0] + t * x[1] - Math.Exp(t));
+            output[1] = (2*ii/25.0) * (5 * x[0] + ii * x[1] - 5 * Math.Exp(t));
+            output[2] = 2 * (x[2] + x[3] * Math.Sin(t) - Math.Cos(t));
+            output[3] = 2 * Math.Sin(t) * (x[2] + Math.Sin(t) * x[3] - Math.Cos(t));
+        }
+
+        public override void ItemHessianByRef(Vector<double> x, int itemIndex, Matrix<double> output)
+        {
+            for (int ii = 0; ii < 4; ++ii)
+                for (int jj = 0; jj < 4; ++jj)
+                    output[ii, jj] = 0;
+            var i = itemIndex + 1;
+            var t = i / 5.0;
+            output[0, 0] = 2;
+            output[0, 1] = 2 * t;
+            output[1, 0] = 2 * t;
+            output[1, 1] = 2 * t * t;
+            output[2, 2] = 2;
+            output[2, 3] = 2 * Math.Sin(t);
+            output[3, 2] = 2 * Math.Sin(t);
+            output[3, 3] = 2 * Math.Pow(Math.Sin(t), 2);
+        }
+
+        public override double ItemValue(Vector<double> x, int itemIndex)
+        {
+            var ii = itemIndex + 1;
+            var t = ii / 5.0;
+            return Math.Pow(x[0] + t * x[1] - Math.Exp(t), 2.0) + Math.Pow(x[2] + x[3] * Math.Sin(t) - Math.Cos(t), 2);
+        }
+    }
+}

src/UnitTests/OptimizationTests/TestFunctions/BrownBadlyScaledFunction.cs

@@ -0,0 +1,131 @@
+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.TestFunctions
+{
+    public class BrownBadlyScaledFunction : BaseTestFunction
+    {
+        public static IEnumerable<TestCase> TestCases
+        {
+            get
+            {
+                yield return new TestCase()
+                {
+                    Function = new BrownBadlyScaledFunction(),
+                    InitialGuess = new double[] { 1, 1 },
+                    MinimalValue = 0,
+                    MinimizingPoint = new double[] { 1e6, 2e-6 },
+                    CaseName = "unbounded"
+                };
+                yield return new TestCase()
+                {
+                    Function = new BrownBadlyScaledFunction(),
+                    InitialGuess = new double[] { 1, 1 },
+                    MinimalValue = 0,
+                    MinimizingPoint = new double[] { 1e6, 2e-6 },
+                    LowerBound = new double[] { -1e8, -1e8 },
+                    UpperBound = new double[] { 1e8, 1e8 },
+                    CaseName = "loose bounds"
+                };
+                yield return new TestCase()
+                {
+                    Function = new BrownBadlyScaledFunction(),
+                    InitialGuess = new double[] { 1, 1 },
+                    MinimalValue = 0.784e3,
+                    MinimizingPoint = new double[] { 1e6, 2e-6 },
+                    LowerBound = new double[] { 0, 3e-5 },
+                    UpperBound = new double[] { 1e6, 100 },
+                    CaseName = "tight bounds"                 
+                };
+            }
+        }
+
+        public BrownBadlyScaledFunction() { }
+
+        public override string Description
+        {
+            get
+            {
+                return "Brown badly scaled fun (MGH #4)";
+            }
+        }
+
+        public override int ItemDimension
+        {
+            get
+            {
+                return 3;
+            }
+        }
+
+        public override int ParameterDimension
+        {
+            get
+            {
+                return 2;
+            }
+        }
+
+        public override void ItemGradientByRef(Vector<double> x, int itemIndex, Vector<double> output)
+        {
+            switch (itemIndex)
+            {
+                case 0:
+                    output[0] = 1;
+                    output[1] = 0;
+                    break;
+                case 1:
+                    output[0] = 0;
+                    output[1] = 1;
+                    break;
+                case 2:
+                    output[0] = x[1];
+                    output[1] = x[0];
+                    break;
+                default:
+                    throw new ArgumentException("itemIndex must be <= 2");
+            }
+        }
+
+        public override void ItemHessianByRef(Vector<double> x, int itemIndex, Matrix<double> output)
+        {
+            switch (itemIndex)
+            {
+                case 0:
+                case 1:
+                    output[0, 0] = 0;
+                    output[0, 1] = 0;
+                    output[1, 0] = 0;
+                    output[1, 1] = 0;
+                    break;
+                case 2:
+                    output[0, 0] = 0;
+                    output[0, 1] = 1;
+                    output[1, 0] = 1;
+                    output[1, 1] = 0;
+                    break;
+                default:
+                    throw new ArgumentException("itemIndex must be <= 2");
+            }
+        }
+
+        public override double ItemValue(Vector<double> x, int itemIndex)
+        {
+            switch (itemIndex)
+            {
+                case 0:
+                    return x[0] - 1e6;
+                case 1:
+                    return x[1] - 2e-6;
+                case 2:
+                    return x[0] * x[1] - 2;
+                default:
+                    throw new ArgumentException("itemIndex must be <= 2");
+            }
+        }
+    }
+}

src/UnitTests/OptimizationTests/TestFunctions/FreudensteinAndRothFunction.cs

@@ -0,0 +1,98 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+using MathNet.Numerics.LinearAlgebra;
+using DenseVector = MathNet.Numerics.LinearAlgebra.Double.DenseVector;
+using DenseMatrix = MathNet.Numerics.LinearAlgebra.Double.DenseMatrix;
+
+namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions
+{
+    public class FreudensteinAndRothFunction : BaseTestFunction
+    {
+
+        public static IEnumerable<TestCase> TestCases
+        {
+            get
+            {
+                yield return new TestCase()
+                {
+                    Function = new FreudensteinAndRothFunction(),
+                    InitialGuess = new double[] { 0.5, -2 },
+                    MinimizingPoint = new double[] { 5, 4 },
+                    MinimalValue = 0,
+                    CaseName = "unbounded"
+                };
+                yield return new TestCase()
+                {
+                    Function = new FreudensteinAndRothFunction(),
+                    InitialGuess = new double[] { 0.5, -2 },
+                    MinimizingPoint = new double[] {5, 4},
+                    MinimalValue = 0,
+                    LowerBound = new double[] { -1000, -1000 },
+                    UpperBound = new double[] { 1000, 1000},
+                    CaseName = "loose bounds"
+                };
+            }
+        }
+
+        public override string Description {  get { return "Freudenstein & Roth fun (MGH #2)"; } }
+
+        public override int ParameterDimension
+        {
+            get
+            {
+                return 2;
+            }
+        }
+
+        public override int ItemDimension
+        {
+            get
+            {
+                return 2;
+            }
+        }
+
+        public override double ItemValue(Vector<double> x, int itemIndex)
+        {
+            if (itemIndex == 0)
+                return -13 + x[0] + ((5 - x[1]) * x[1] - 2) * x[1];
+            else 
+                return -29 + x[0] + ((x[1] + 1) * x[1] - 14) * x[1];
+        }
+
+        public override void ItemGradientByRef(Vector<double> x, int itemIndex, Vector<double> output)
+        {
+            if (itemIndex == 0)
+            {
+                output[0] = 1;
+                output[1] = -2 + (5 - 2 * x[1]) * x[1] + (5 - x[1]) * x[1];
+            }
+            else
+            {
+                output[0] = 1;
+                output[1] = -14 + x[1] * (1 + x[1]) + x[1] * (1 + 2 * x[1]);
+            }
+        }
+
+        public override void ItemHessianByRef(Vector<double> x, int itemIndex, Matrix<double> output)
+        {
+            if (itemIndex == 0)
+            {
+                output[0, 0] = 0;
+                output[0, 1] = 0;
+                output[1, 0] = 0;
+                output[1, 1] = 10 - 6 * x[1];
+            }
+            else
+            {
+                output[0, 0] = 0;
+                output[0, 1] = 0;
+                output[1, 0] = 0;
+                output[1, 1] = 2 + 6 * x[1];
+            }
+        }
+    }
+}

src/UnitTests/OptimizationTests/TestFunctions/HelicalValleyFunction.cs

@@ -0,0 +1,178 @@
+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.TestFunctions
+{
+    public class HelicalValleyFunction : BaseTestFunction
+    {
+        public static IEnumerable<TestCase> TestCases
+        {
+            get
+            {
+                yield return new TestCase()
+                {
+                    Function = new HelicalValleyFunction(),
+                    InitialGuess = new double[] { -1, 0, 0 },
+                    MinimalValue = 0,
+                    MinimizingPoint = new double[] { 1, 0, 0 },
+                    CaseName = "unbounded"
+                };
+                yield return new TestCase()
+                {
+                    Function = new HelicalValleyFunction(),
+                    InitialGuess = new double[] { -1, 0, 0 },
+                    MinimalValue = 0,
+                    MinimizingPoint = new double[] { 1, 0, 0 },
+                    LowerBound = new double[] { -1000, -1000, -1000 },
+                    UpperBound = new double[] { 1000, 1000, 1000 },
+                    CaseName = "loose bounds"
+                };
+                yield return new TestCase()
+                {
+                    Function = new HelicalValleyFunction(),
+                    InitialGuess = new double[] { -1, 0, 0 },
+                    MinimalValue = 0.99042212,
+                    LowerBound = new double[] { -100, -1, -1 },
+                    UpperBound = new double[] { 0.8, 1, 1 },
+                    CaseName = "tight bounds"
+                };
+            }
+        }
+
+        public override string Description
+        {
+            get
+            {
+                return "Helical valley fun (MGH #7)";
+            }
+        }
+
+        public override int ItemDimension
+        {
+            get
+            {
+                return 3;
+            }
+        }
+
+        public override int ParameterDimension
+        {
+            get
+            {
+                return 3;
+            }
+        }
+
+        public override void ItemGradientByRef(Vector<double> x, int itemIndex, Vector<double> output)
+        {
+            switch (itemIndex)
+            {
+                case 0:
+                    output[0] = -100 * theta10(x[0], x[1]);
+                    output[1] = -100 * theta01(x[0], x[1]);
+                    output[2] = 10;
+                    break;
+                case 1:
+                    output[0] = (10 * x[0]) / Math.Sqrt(x[0]*x[0] + x[1]*x[1]);
+                    output[1] = (10 * x[1]) / Math.Sqrt(x[0]*x[0] + x[1]*x[1]);
+                    output[2] = 0;
+                    break;
+                case 2:
+                    output[0] = 0;
+                    output[1] = 0;
+                    output[2] = 1;
+                    break;
+                default:
+                    throw new ArgumentException("itemIndex must be <= 2");
+            }
+        }
+
+        public override void ItemHessianByRef(Vector<double> x, int itemIndex, Matrix<double> output)
+        {
+            switch (itemIndex)
+            {
+                case 0:
+                    output[0, 0] = -100 * theta20(x[0], x[1]);
+                    output[0, 1] = -100 * theta11(x[0], x[1]);
+                    output[0, 2] = 0;
+                    output[1, 0] = -100 * theta11(x[0], x[1]);
+                    output[1, 1] = -100 * theta02(x[0], x[1]);
+                    output[1, 2] = 0;
+                    output[2, 0] = 0;
+                    output[2, 1] = 0;
+                    output[2, 2] = 0;
+                    break;
+                case 1:
+                    output[0, 0] = (10 * x[1]*x[1]) / Math.Pow(x[0]*x[0] + x[1]*x[1],1.5);
+                    output[0, 1] = (-10 * x[0] * x[1]) / Math.Pow(x[0]*x[0] + x[1]*x[1],1.5);
+                    output[0, 2] = 0;
+                    output[1, 0] = (-10 * x[0] * x[1]) / Math.Pow(x[0] * x[0] + x[1] * x[1], 1.5);
+                    output[1, 1] = (10 * x[0]*x[0]) / Math.Pow(x[0] * x[0] + x[1] * x[1], 1.5);
+                    output[1, 2] = 0;
+                    output[2, 0] = 0;
+                    output[2, 1] = 0;
+                    output[2, 2] = 0;
+                    break;
+                case 2:
+                    for (int ii = 0; ii < 2; ++ii)
+                        for (int jj = 0; jj < 2; ++jj)
+                            output[ii, jj] = 0;
+                    break;
+                default:
+                    throw new ArgumentException("itemIndex must be <= 2");
+            }
+        }
+
+        private static double theta(double x1, double x2)
+        {
+            if (x1 >= 0)
+                return 0.5 * Math.Atan(x2 / x1) / Math.PI;
+            else
+                return 0.5 * Math.Atan(x2 / x1) / Math.PI + 0.5;
+        }
+
+        private static double theta10(double x1, double x2)
+        {
+            return -(x2 / (2 * Math.PI * Math.Pow(x1,2) + 2 * Math.PI * Math.Pow(x2,2)));
+        }
+
+        private static double theta01(double x1, double x2)
+        {
+            return x1 / (2 * Math.PI * x1*x1 + 2 * Math.PI * x2*x2);
+        }
+
+        private static double theta20(double x1,double x2)
+        {
+            return (x1 * x2) / (Math.PI * Math.Pow(x1 * x1 + x2 * x2, 2));
+        }
+
+        private static double theta11(double x1, double x2)
+        {
+            return (-x1 * x1 + x2 * x2) / (2 * Math.PI * Math.Pow(x1 * x1 + x2 * x2, 2));
+        }
+
+        private static double theta02(double x1, double x2)
+        {
+            return -((x1 * x2) / (Math.PI * Math.Pow(x1*x1 + x2*x2, 2)));
+        }
+        
+        public override double ItemValue(Vector<double> x, int itemIndex)
+        {
+            switch (itemIndex)
+            {
+                case 0:
+                    return 10 * (x[2] - 10 * theta(x[0], x[1]));
+                case 1:
+                    return 10 * (Math.Sqrt(x[0] * x[0] + x[1] * x[1]) - 1);
+                case 2:
+                    return x[2];
+                default:
+                    throw new ArgumentException("itemIndex must be <= 2");
+            }
+        }
+    }
+}

src/UnitTests/OptimizationTests/TestFunctions/ITestFunction.cs

@@ -0,0 +1,69 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+using MathNet.Numerics.LinearAlgebra;
+using DenseVector = MathNet.Numerics.LinearAlgebra.Double.DenseVector;
+
+namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions
+{
+    public class TestCase
+    {
+        public string CaseName;
+        public ITestFunction Function;
+        public DenseVector InitialGuess;
+        public DenseVector LowerBound;
+        public DenseVector UpperBound;
+        public double MinimalValue;
+        public DenseVector MinimizingPoint;
+
+        public bool IsBounded
+        {
+            get
+            {
+                return this.LowerBound != null && this.UpperBound != null;
+            }
+        }
+
+        public bool IsUnbounded
+        {
+            get
+            {
+                return this.IsUnboundedOverride ?? this.LowerBound == null || this.UpperBound == null;
+            }
+        }
+
+        public bool? IsUnboundedOverride;
+
+        public string FullName
+        {
+            get
+            {
+                return $"{this.Function.Description} {this.CaseName}";
+            }
+        }
+    }
+
+    public interface ITestFunction
+    {
+        string Description { get; }
+        int ParameterDimension { get; }
+        int ItemDimension { get; }
+
+        double ItemValue(Vector<double> x, int itemIndex);
+        Vector<double> ItemGradient(Vector<double> x, int itemIndex);
+        void ItemGradientByRef(Vector<double> x, int itemIndex, Vector<double> output);
+        Matrix<double> ItemHessian(Vector<double> x, int itemIndex);
+        void ItemHessianByRef(Vector<double> x, int itemIndex, Matrix<double> output);
+
+        Matrix<double> Jacobian(Vector<double> x);
+        void JacobianbyRef(Vector<double> x, Matrix<double> output);
+
+        double SsqValue(Vector<double> x);
+        Vector<double> SsqGradient(Vector<double> x);
+        void SsqGradientByRef(Vector<double> x, Vector<double> output);
+        Matrix<double> SsqHessian(Vector<double> x);
+        void SsqHessianByRef(Vector<double> x, Matrix<double> output);
+    }
+}

src/UnitTests/OptimizationTests/TestFunctions/JennrichAndSampsonFunction.cs

@@ -0,0 +1,102 @@
+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.TestFunctions
+{
+    public class JennrichAndSampsonFunction : BaseTestFunction
+    {
+        private readonly int _m;
+
+        public JennrichAndSampsonFunction(int itemDimension)
+        {
+            if (itemDimension < 2)
+                throw new ArgumentException("itemDimension must be at least 2.");
+            _m = itemDimension;
+        }
+
+        public static IEnumerable<TestCase> TestCases
+        {
+            get
+            {
+                yield return new TestCase()
+                {
+                    Function = new JennrichAndSampsonFunction(10),
+                    InitialGuess = new double[] { 0.3, 0.4 },
+                    MinimalValue = 124.362,
+                    MinimizingPoint = new double[] { 0.2578, 0.2578 },
+                    CaseName = "unbounded"
+                };
+                //yield return new TestCase()
+                //{
+                //    Function = new JennrichAndSampsonFunction(10),
+                //    LowerBound = new double[] { 0.6, 0.5 },
+                //    UpperBound = new double[] { 10, 50 },
+                //    StartPoint = new double[] { 1.0, 1.0 },
+                //    MinimizingInput = null,
+                //    MinimizingValue = 0,
+                //    CaseName = "tight bounds"
+                //};
+                yield return new TestCase()
+                {
+                    Function = new JennrichAndSampsonFunction(10),
+                    LowerBound = new double[] { -50, -50 },
+                    UpperBound = new double[] { 50, 50 },
+                    InitialGuess = new double[] { 0.3, 0.4 },
+                    MinimizingPoint = null,
+                    MinimalValue = 0,
+                    CaseName = "loose bounds"
+                };
+            }
+        }
+
+        public override string Description
+        {
+            get
+            {
+                return $"Jennrich & Sampson fun (MGH #6) (n={this.ItemDimension})";
+            }
+        }
+
+        public override int ItemDimension
+        {
+            get
+            {
+                return _m;
+            }
+        }
+
+        public override int ParameterDimension
+        {
+            get
+            {
+                return 2;
+            }
+        }
+
+        public override void ItemGradientByRef(Vector<double> x, int itemIndex, Vector<double> output)
+        {
+            int ii = itemIndex + 1;
+            output[0] = -(Math.Exp(ii * x[0]) * ii);
+            output[1] = -(Math.Exp(ii * x[1]) * ii);
+        }
+
+        public override void ItemHessianByRef(Vector<double> x, int itemIndex, Matrix<double> output)
+        {
+            int ii = itemIndex + 1;
+            output[0, 0] = -(Math.Exp(ii * x[0]) * ii*ii);
+            output[0, 1] = 0;
+            output[1, 0] = 0;
+            output[1, 1] = -(Math.Exp(ii * x[1]) * ii*ii);
+        }
+
+        public override double ItemValue(Vector<double> x, int itemIndex)
+        {
+            int ii = itemIndex + 1;
+            return 2 + 2 * ii - (Math.Exp(ii * x[0]) + Math.Exp(ii * x[1]));
+        }
+    }
+}

src/UnitTests/OptimizationTests/TestFunctions/MeyerFunction.cs

@@ -0,0 +1,99 @@
+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.TestFunctions
+{
+    public class MeyerFunction : BaseTestFunction
+    {
+        public static IEnumerable<TestCase> TestCases
+        {
+            get
+            {
+                yield return new TestCase()
+                {
+                    Function = new MeyerFunction(),
+                    InitialGuess = new double[] { 0.02, 4000, 250 },
+                    MinimalValue = 87.9458,
+                    MinimizingPoint = null,
+                    CaseName = "unbounded"
+                };
+                yield return new TestCase()
+                {
+                    Function = new MeyerFunction(),
+                    InitialGuess = new double[] { 0.02, 4000, 250 },
+                    MinimalValue = 87.9458,
+                    MinimizingPoint = null,
+                    LowerBound = new double[] { -1e6, -1e6, -1e6 },
+                    UpperBound = new double[] { 1e6, 1e6, 1e6 },
+                    CaseName = "loose bounds"
+                };
+            }
+        }
+
+        public MeyerFunction() { }
+
+        public override string Description
+        {
+            get
+            {
+                return "Meyer fun (MGH #10)";
+            }
+        }
+
+        public override int ItemDimension
+        {
+            get
+            {
+                return 16;
+            }
+        }
+
+        public override int ParameterDimension
+        {
+            get
+            {
+                return 3;
+            }
+        }
+
+        public override void ItemGradientByRef(Vector<double> x, int itemIndex, Vector<double> output)
+        {
+            int ii = itemIndex + 1;
+            output[0] = Math.Exp(x[1] / (45.0 + 5 * ii + x[2]));
+            output[1] = (Math.Exp(x[1] / (45.0 + 5 * ii + x[2])) * x[0]) / (45 + 5 * ii + x[2]);
+            output[2] = -(Math.Exp(x[1] / (45.0 + 5 * ii + x[2])) * x[0] * x[1]) / Math.Pow(45 + 5 * ii + x[2], 2);
+
+        }
+
+        public override void ItemHessianByRef(Vector<double> x, int itemIndex, Matrix<double> output)
+        {
+            var ii = itemIndex + 1;
+
+            var t0 = (45.0 + 5 * ii + x[2]);
+            var t1 = Math.Exp(x[1] / t0);
+
+            output[0, 0] = 0;
+            output[0, 1] = t1 / t0;
+            output[0, 2] = -t1 * x[1] / Math.Pow(t0, 2);
+            output[1, 0] = t1 / t0;
+            output[1, 1] = t1 * x[0] / Math.Pow(t0, 2);
+            output[1, 2] = -t1 * x[0] * (t0 + x[1]) / Math.Pow(t0, 3);
+            output[2, 0] = -t1 * x[1] / Math.Pow(t0, 2);
+            output[2, 1] = -t1 * x[0] * (t0 + x[1]) / Math.Pow(t0, 3);
+            output[2, 2] = t1 * x[0] * x[1] * (2*t0 + x[1]) / Math.Pow(t0, 4);
+        }
+
+        private static readonly double[] y = { 34780, 28610, 23650, 19630, 16370, 13720, 11540, 9744, 8261, 7030, 6005, 5147, 4427, 3820, 3307, 2872 };
+
+        public override double ItemValue(Vector<double> x, int itemIndex)
+        {
+            var ii = itemIndex + 1;
+            var t = 45.0 + 5 * ii;
+            return x[0] * Math.Exp(x[1] / (t + x[2])) - y[itemIndex];
+        }
+    }
+}

src/UnitTests/OptimizationTests/TestFunctions/PowellBadlyScaledFunction.cs

@@ -0,0 +1,111 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+using MathNet.Numerics.LinearAlgebra;
+using MathNet.Numerics.LinearAlgebra.Double;
+
+namespace MathNet.Numerics.UnitTests.OptimizationTests.TestFunctions
+{
+    public class PowellBadlyScaledFunction : BaseTestFunction
+    {
+        public static IEnumerable<TestCase> TestCases
+        {
+            get
+            {
+                yield return new TestCase()
+                {
+                    Function = new PowellBadlyScaledFunction(),
+                    InitialGuess = new double[] { 0, 1 },
+                    MinimizingPoint = new double[] { 1.098e-5, 9.106 },
+                    MinimalValue = 0,
+                    CaseName = "unbounded"
+                };
+                yield return new TestCase()
+                {
+                    Function = new PowellBadlyScaledFunction(),
+                    InitialGuess = new double[] { 0, 1 },
+                    MinimizingPoint = new double[] { 1.098e-5, 9.106 },
+                    MinimalValue = 0,
+                    LowerBound = new double[] { -1000, -1000 },
+                    UpperBound = new double[] { 1000, 1000 },
+                    CaseName = "loose bounds"
+                };
+                yield return new TestCase()
+                {
+                    Function = new PowellBadlyScaledFunction(),
+                    LowerBound = new double[] { 0, 1 },
+                    UpperBound = new double[] { 1, 9 },
+                    InitialGuess = new double[] { 0, 1 },
+                    MinimalValue = 0.15125900e-9,
+                    CaseName = "tight bounds"
+                };
+            }
+        }
+
+        public override string Description
+        {
+            get
+            {
+                return "Powell badly scaled fun (MGH #3)";
+            }
+        }
+
+        public override int ItemDimension
+        {
+            get
+            {
+                return 2;
+            }
+        }
+
+        public override int ParameterDimension
+        {
+            get
+            {
+                return 2;
+            }
+        }
+
+        public override void ItemGradientByRef(Vector<double> x, int itemIndex, Vector<double> output)
+        {
+            if (itemIndex == 0)
+            {
+                output[0] = 10000 * x[1];
+                output[1] = 10000 * x[0];
+            }
+            else if (itemIndex == 1)
+            {
+                output[0] = -Math.Exp(-x[0]);
+                output[1] = -Math.Exp(-x[1]);
+            }
+        }
+
+        public override void ItemHessianByRef(Vector<double> x, int itemIndex, Matrix<double> output)
+        {
+            if (itemIndex == 0)
+            {
+                output[0, 0] = 0;
+                output[0, 1] = 10000;
+                output[1, 0] = 10000;
+                output[1, 1] = 0;
+            }
+            else
+            {
+                output[0, 0] = Math.Exp(-x[0]);
+                output[0, 1] = 0;
+                output[1, 0] = 0;
+                output[1,1] = Math.Exp(-x[1]);
+            }
+        }
+
+        public override double ItemValue(Vector<double> x, int itemIndex)
+        {
+            if (itemIndex == 0)
+                return 10000.0 * x[0] * x[1] - 1;
+            else
+                return Math.Exp(-x[0]) + Math.Exp(-x[1]) - 1.0001;
+        }
+    }
+}

src/UnitTests/OptimizationTests/TestFunctions/PowellSingularFunction.cs

@@ -0,0 +1,141 @@
+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.TestFunctions
+{
+    public class PowellSingularFunction : BaseTestFunction
+    {
+        public static IEnumerable<TestCase> TestCases
+        {
+            get
+            {
+                yield return new TestCase()
+                {
+                    Function = new PowellSingularFunction(),
+                    InitialGuess = new double[] { 3, -1, 0, 1 },
+                    MinimalValue = 0,
+                    MinimizingPoint = new double[] {0,0,0,0},
+                    CaseName = "unbounded"
+                };
+                yield return new TestCase()
+                {
+                    Function = new PowellSingularFunction(),
+                    InitialGuess = new double[] { 3, -1, 0, 1 },
+                    MinimalValue = 0,
+                    MinimizingPoint = new double[] { 0, 0, 0, 0 },
+                    LowerBound = new double[] {-1000, -1000, -1000, -1000},
+                    UpperBound = new double[] { 1000, 1000, 1000, 1000 },
+                    CaseName = "loose bounds"
+                };
+            }
+        }
+
+
+        public PowellSingularFunction() { }
+
+        public override string Description
+        {
+            get
+            {
+                return "Powell singular fun (MGH #13)";
+            }
+        }
+
+        public override int ItemDimension
+        {
+            get
+            {
+                return 4;
+            }
+        }
+
+        public override int ParameterDimension
+        {
+            get
+            {
+                return 4;
+            }
+        }
+
+        public override void ItemGradientByRef(Vector<double> x, int itemIndex, Vector<double> output)
+        {
+            switch (itemIndex)
+            {
+                case 0:
+                    output[0] = 1;
+                    output[1] = 10;
+                    output[2] = 0;
+                    output[3] = 0;
+                    break;
+                case 1:
+                    output[0] = 0;
+                    output[1] = 0;
+                    output[2] = Math.Sqrt(5);
+                    output[3] = -Math.Sqrt(5);
+                    break;
+                case 2:
+                    output[0] = 0;
+                    output[1] = 2*(x[1]-2*x[2]);
+                    output[2] = -4*x[1] + 8*x[2];
+                    output[3] = 0;
+                    break;
+                case 3:
+                    output[0] = 2*Math.Sqrt(10)*(x[0] - x[3]);
+                    output[1] = 0;
+                    output[2] = 0;
+                    output[3] = -2*Math.Sqrt(10)*(x[0] - x[3]);
+                    break;
+                default:
+                    throw new ArgumentException("itemIndex must be <= 3");
+            }
+        }
+
+        public override void ItemHessianByRef(Vector<double> x, int itemIndex, Matrix<double> output)
+        {
+            for (int ii = 0; ii < 4; ++ii)
+                for (int jj = 0; jj < 4; ++jj)
+                    output[ii, jj] = 0;
+            switch(itemIndex)
+            {
+                case 0:
+                case 1:
+                    break;
+                case 2:
+                    output[1, 1] = 2;
+                    output[1, 2] = -4;
+                    output[2, 1] = -4;
+                    output[2, 2] = 8;
+                    break;
+                case 3:
+                    output[0, 0] = 2 * Math.Sqrt(10);
+                    output[0, 3] = -2 * Math.Sqrt(10);
+                    output[3, 0] = -2 * Math.Sqrt(10);
+                    output[3, 3] = 2 * Math.Sqrt(10);
+                    break;
+                default:
+                    throw new ArgumentException("itemIndex must be <= 3");
+            }
+        }
+
+        public override double ItemValue(Vector<double> x, int itemIndex)
+        {
+            switch (itemIndex)
+            {
+                case 0:
+                    return x[0] + 10 * x[1];
+                case 1:
+                    return Math.Sqrt(5) * (x[2] - x[3]);
+                case 2:
+                    return Math.Pow(x[1] - 2 * x[2], 2);
+                case 3:
+                    return Math.Sqrt(10.0) * Math.Pow(x[0] - x[3], 2);
+                default:
+                    throw new ArgumentException("itemIndex must be <= 3");
+            }
+        }
+    }
+}

src/UnitTests/OptimizationTests/TestFunctions/RosenbrockFunction2.cs

@@ -0,0 +1,156 @@
+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.TestFunctions
+{
+    public class RosenbrockFunction2 : BaseTestFunction
+    {
+        public static IEnumerable<TestCase> TestCases
+        {
+            get
+            {
+                yield return new TestCase()
+                {
+                    Function = new RosenbrockFunction2(),
+                    InitialGuess = new double[] { -1.2, 1 },
+                    MinimizingPoint = new double[] { 1, 1 },
+                    MinimalValue = 0,
+                    LowerBound = new double[] { -1000, -1000 },
+                    UpperBound = new double[] { 1000, 1000 },
+                    CaseName = "hard start",
+                    IsUnboundedOverride = true
+                };
+                yield return new TestCase()
+                {
+                    Function = new RosenbrockFunction2(),
+                    InitialGuess = new double[] { 1.2, 1.2 },
+                    MinimizingPoint = new double[] { 1, 1 },
+                    MinimalValue = 0,
+                    LowerBound = new double[] { -5, -5 },
+                    UpperBound = new double[] { 5, 5 },
+                    CaseName = "easy start"
+                };
+                yield return new TestCase()
+                {
+                    Function = new RosenbrockFunction2(),
+                    InitialGuess = new double[] { -0.9, -0.5 },
+                    MinimizingPoint = new double[] { 1, 1 },
+                    MinimalValue = 0,
+                    LowerBound = new double[] { -5, -5 },
+                    UpperBound = new double[] { 5, 5 },
+                    CaseName = "Overton start",
+                    IsUnboundedOverride = true
+                };
+                yield return new TestCase()
+                {
+                    Function = new RosenbrockFunction2(),
+                    InitialGuess = new double[] { 1.2, 1.2 },
+                    MinimizingPoint = new double[] { 1, 1 },
+                    MinimalValue = 0,
+                    LowerBound = new double[] { 1, -5 },
+                    UpperBound = new double[] { 5, 5 },
+                    CaseName = "easy one active bound"
+                };
+                yield return new TestCase()
+                {
+                    Function = new RosenbrockFunction2(),
+                    InitialGuess = new double[] { 1.2, 1.2 },
+                    MinimizingPoint = new double[] { 1, 1 },
+                    MinimalValue = 0,
+                    LowerBound = new double[] { 1, 1 },
+                    UpperBound = new double[] { 5, 5 },
+                    CaseName = "easy two active bounds"
+                };
+                yield return new TestCase()
+                {
+                    Function = new RosenbrockFunction2(),
+                    InitialGuess = new double[] { 2.5, 2.5 },
+                    MinimizingPoint = new double[] { 2, 4 },
+                    MinimalValue = 1,
+                    LowerBound = new double[] { 2, 2 },
+                    UpperBound = new double[] { 5, 5 },
+                    CaseName = "min on lower bound, not local"
+                };
+                yield return new TestCase()
+                {
+                    Function = new RosenbrockFunction2(),
+                    InitialGuess = new double[] { -0.9, -0.5 },
+                    MinimizingPoint = new double[] { 0.5, 0.25 },
+                    MinimalValue = 0.25,
+                    LowerBound = new double[] { -2, -2 },
+                    UpperBound = new double[] { 0.5, 0.5 },
+                    CaseName = "min on upper bound, not local"
+                };
+
+
+            }
+        }
+        public RosenbrockFunction2() { }
+
+        public override string Description
+        {
+            get
+            {
+                return "Rosenbrock fun (MGH #1)";
+            }
+        }
+
+        public override int ItemDimension
+        {
+            get
+            {
+                return 2;
+            }
+        }
+
+        public override int ParameterDimension
+        {
+            get
+            {
+                return 2;
+            }
+        }
+
+        public override void ItemGradientByRef(Vector<double> x, int itemIndex, Vector<double> output)
+        {
+            if (itemIndex == 0)
+            {
+                output[0] = -20 * x[0];
+                output[1] = 10;
+            } else
+            {
+                output[0] = -1;
+                output[1] = 0;
+            }
+        }
+
+        public override void ItemHessianByRef(Vector<double> x, int itemIndex, Matrix<double> output)
+        {
+            if (itemIndex == 0)
+            {
+                output[0, 0] = -20;
+                output[0, 1] = 0;
+                output[1, 0] = 0;
+                output[1, 1] = 0;
+            } else
+            {
+                output[0, 0] = 0;
+                output[0, 1] = 0;
+                output[1, 0] = 0;
+                output[1, 1] = 0;
+            }
+        }
+
+        public override double ItemValue(Vector<double> x, int itemIndex)
+        {
+            if (itemIndex == 0)
+                return 10 * (x[1] - x[0] * x[0]);
+            else
+                return 1 - x[0];
+        }
+    }
+}

src/UnitTests/OptimizationTests/TestFunctions/WoodFunction.cs

@@ -0,0 +1,164 @@
+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.TestFunctions
+{
+    public class WoodFunction : BaseTestFunction
+    {
+        public static IEnumerable<TestCase> TestCases
+        {
+            get
+            {
+                yield return new TestCase()
+                {
+                    Function = new WoodFunction(),
+                    InitialGuess = new double[] { -3, -1, -3, -1 },
+                    MinimalValue = 0,
+                    MinimizingPoint = new double[] { 1, 1, 1, 1 },
+                    CaseName = "unbounded"
+                };
+                yield return new TestCase()
+                {
+                    Function = new WoodFunction(),
+                    InitialGuess = new double[] { -3, -1, -3, -1 },
+                    MinimalValue = 0,
+                    MinimizingPoint = new double[] { 1, 1, 1, 1 },
+                    LowerBound = new double[] { -1000, -1000, -1000, -1000 },
+                    UpperBound = new double[] { 1000, 1000, 1000, 1000 },
+                    CaseName = "loose bounds"
+                };
+                yield return new TestCase()
+                {
+                    Function = new WoodFunction(),
+                    InitialGuess = new double[] { -3, -1, -3, -1 },
+                    MinimalValue = 1.5567008,
+                    MinimizingPoint = null,
+                    LowerBound = new double[] { -100, -100, -100, -100 },
+                    UpperBound = new double[] { 0, 10, 100, 100 },
+                    CaseName = "tight bounds"
+                };
+            }
+        }
+
+        public WoodFunction() { }
+
+        public override string Description
+        {
+            get
+            {
+                return "Wood fun (MGH #14)";
+            }
+        }
+
+        public override int ItemDimension
+        {
+            get
+            {
+                return 6;
+            }
+        }
+
+        public override int ParameterDimension
+        {
+            get
+            {
+                return 4;
+            }
+        }
+
+        public override void ItemGradientByRef(Vector<double> x, int itemIndex, Vector<double> output)
+        {
+            switch (itemIndex)
+            {
+                case 0:
+                    output[0] = -20 * x[0];
+                    output[1] = 10;
+                    output[2] = 0;
+                    output[3] = 0;
+                    break;
+                case 1:
+                    output[0] = -1;
+                    output[1] = 0;
+                    output[2] = 0;
+                    output[3] = 0;
+                    break;
+                case 2:
+                    output[0] = 0;
+                    output[1] = 0;
+                    output[2] = -6 * Math.Sqrt(10) * x[2];
+                    output[3] = 3 * Math.Sqrt(10);
+                    break;
+                case 3:
+                    output[0] = 0;
+                    output[1] = 0;
+                    output[2] = -1;
+                    output[3] = 0;
+                    break;
+                case 4:
+                    output[0] = 0;
+                    output[1] = Math.Sqrt(10);
+                    output[2] = 0;
+                    output[3] = Math.Sqrt(10);
+                    break;
+                case 5:
+                    output[0] = 0;
+                    output[1] = 1.0 / Math.Sqrt(10);
+                    output[2] = 0;
+                    output[3] = -1.0 / Math.Sqrt(10);
+                    break;
+                default:
+                    throw new ArgumentException("itemIndex must be <= 5");
+            }
+        }
+
+        public override void ItemHessianByRef(Vector<double> x, int itemIndex, Matrix<double> output)
+        {
+            for (int ii = 0; ii < 4; ++ii)
+                for (int jj = 0; jj < 4; ++jj)
+                    output[ii, jj] = 0;
+            switch (itemIndex)
+            {
+                case 0:
+                    output[0, 0] = -20;
+                    break;
+                case 1:
+                    break;
+                case 2:
+                    output[2, 2] = -6 * Math.Sqrt(10);
+                    break;
+                case 3:
+                case 4:
+                case 5:
+                    break;
+                default:
+                    throw new ArgumentException("itemIndex must be <= 5");
+
+            }
+        }
+
+        public override double ItemValue(Vector<double> x, int itemIndex)
+        {
+            switch (itemIndex)
+            {
+                case 0:
+                    return 10 * (x[1] - x[0] * x[0]);
+                case 1:
+                    return 1 - x[0];
+                case 2:
+                    return Math.Sqrt(90) * (x[3] - x[2] * x[2]);
+                case 3:
+                    return 1 - x[2];
+                case 4:
+                    return Math.Sqrt(10) * (x[1] + x[3] - 2);
+                case 5:
+                    return (x[1] - x[3]) / Math.Sqrt(10);
+                default:
+                    throw new ArgumentException("itemIndex must be <= 5");
+            }
+        }
+    }
+}

src/UnitTests/UnitTests.csproj

@@ -30,7 +30,7 @@
     <CodeAnalysisRuleSet>AllRules.ruleset</CodeAnalysisRuleSet>
     <NoWarn>1591</NoWarn>
     <Prefer32Bit>false</Prefer32Bit>
-    <LangVersion>5</LangVersion>
+    <LangVersion>6</LangVersion>
   </PropertyGroup>
   <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Debug|AnyCPU' ">
     <DefineConstants>DEBUG;TRACE</DefineConstants>
@@ -46,7 +46,7 @@
     <PlatformTarget>AnyCPU</PlatformTarget>
     <NoWarn>1591</NoWarn>
     <Prefer32Bit>false</Prefer32Bit>
-    <LangVersion>5</LangVersion>
+    <LangVersion>6</LangVersion>
   </PropertyGroup>
   <PropertyGroup Condition="'$(Configuration)|$(Platform)' == 'Release-Signed|AnyCPU'">
     <OutputPath>..\..\out\test-signed\Net40\</OutputPath>
@@ -62,7 +62,7 @@
     <CodeAnalysisRuleSet>AllRules.ruleset</CodeAnalysisRuleSet>
     <NoWarn>1591</NoWarn>
     <Prefer32Bit>false</Prefer32Bit>
-    <LangVersion>5</LangVersion>
+    <LangVersion>6</LangVersion>
   </PropertyGroup>
   <ItemGroup>
     <Reference Include="System" />
@@ -366,8 +366,32 @@
     <Compile Include="LinearAlgebraTests\VectorToStringTests.cs" />
     <Compile Include="OdeSolvers\OdeSolverTest.cs" />
     <Compile Include="Random\RandomSerializationTests.cs" />
+    <Compile Include="OptimizationTests\TestCaseDataExtensions.cs" />
+    <Compile Include="OptimizationTests\TestFunctions\BrownAndDennisFunction.cs" />
+    <Compile Include="OptimizationTests\TestFunctions\HelicalValleyFunction.cs" />
+    <Compile Include="OptimizationTests\NelderMeadSimplexTests.cs" />
+    <Compile Include="OptimizationTests\TestFunctionAdapters.cs" />
+    <Compile Include="OptimizationTests\TestFunctions\BaseTestFunction.cs" />
+    <Compile Include="OptimizationTests\TestFunctions\BealeFunction.cs" />
+    <Compile Include="OptimizationTests\TestFunctions\BrownBadlyScaledFunction.cs" />
+    <Compile Include="OptimizationTests\TestFunctions\FreudensteinAndRothFunction.cs" />
+    <Compile Include="OptimizationTests\TestFunctions\ITestFunction.cs" />
+    <Compile Include="OptimizationTests\TestFunctions\JennrichAndSampsonFunction.cs" />
+    <Compile Include="OptimizationTests\TestFunctions\MeyerFunction.cs" />
+    <Compile Include="OptimizationTests\TestFunctions\PowellBadlyScaledFunction.cs" />
+    <Compile Include="OptimizationTests\TestFunctions\PowellSingularFunction.cs" />
+    <Compile Include="OptimizationTests\TestFunctions\RosenbrockFunction2.cs" />
+    <Compile Include="OptimizationTests\TestFunctions\WoodFunction.cs" />
+    <Compile Include="OptimizationTests\GoldenSectionMinimizerTests.cs" />
+    <Compile Include="OptimizationTests\TestFunctionTests.cs" />
     <Compile Include="Random\SystemRandomSourceTests.cs" />
+    <Compile Include="OptimizationTests\BfgsTest.cs" />
     <Compile Include="RootFindingTests\BisectionTest.cs" />
+    <Compile Include="OptimizationTests\RosenbrockFunction.cs" />
+    <Compile Include="OptimizationTests\BfgsMinimizerTests.cs" />
+    <Compile Include="OptimizationTests\ConjugateGradientMinimizerTests.cs" />
+    <Compile Include="OptimizationTests\NewtonMinimizerTests.cs" />
+    <Compile Include="OptimizationTests\RosenbrockFunctionTests.cs" />
     <Compile Include="PermutationTest.cs" />
     <Compile Include="PrecisionTest.cs" />
     <Compile Include="Properties\AssemblyInfo.cs" />
@@ -416,10 +440,16 @@
     <Compile Include="StatisticsTests\StatTestData.cs" />
     <Compile Include="TrigonometryTest.cs" />
     <Compile Include="UseLinearAlgebraProvider.cs" />
+    <Compile Include="OptimizationTests\BfgsBMinimizerTests.cs" />
   </ItemGroup>
   <ItemGroup>
     <None Include="paket.references" />
   </ItemGroup>
+  <ItemGroup>
+    <Service Include="{508349B6-6B84-4DF5-91F0-309BEEBAD82D}" />
+    <Service Include="{82A7F48D-3B50-4B1E-B82E-3ADA8210C358}" />
+  </ItemGroup>
+  <ItemGroup />
   <Import Project="$(MSBuildToolsPath)\Microsoft.CSharp.targets" />
   <Choose>
     <When Condition="$(TargetFrameworkIdentifier) == '.NETFramework' And ($(TargetFrameworkVersion) == 'v2.0' Or $(TargetFrameworkVersion) == 'v3.0')">