Merge pull request #435 from scottstephens/unified_optimization_better_tests

Better tests for optimization

src/Numerics/Optimization/ObjectiveFunctions/LazyObjectiveFunctionBase.cs

@@ -6,14 +6,14 @@ namespace MathNet.Numerics.Optimization.ObjectiveFunctions
     {
         Vector<double> _point;
 
-        bool _hasFunctionValue;
-        double _functionValue;
+        protected bool _hasFunctionValue;
+        protected double _functionValue;
 
-        bool _hasGradientValue;
-        Vector<double> _gradientValue;
+        protected bool _hasGradientValue;
+        protected Vector<double> _gradientValue;
 
-        bool _hasHessianValue;
-        Matrix<double> _hessianValue;
+        protected bool _hasHessianValue;
+        protected Matrix<double> _hessianValue;
 
         protected LazyObjectiveFunctionBase(bool gradientSupported, bool hessianSupported)
         {

src/UnitTests/OptimizationTests/TestBfgsBMinimizer.cs → src/UnitTests/OptimizationTests/BfgsBMinimizerTests.cs

@@ -32,11 +32,16 @@ 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;
 
 namespace MathNet.Numerics.UnitTests.OptimizationTests
 {
 	[TestFixture]
-	public class TestBfgsBMinimizer
+	public class BfgsBMinimizerTests
 	{
 	    [Test]
 		public void FindMinimum_Rosenbrock_Easy()
@@ -149,6 +154,54 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
             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.");
+        }
+
+        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.IsBounded)
+                    .Select(x => new TestCaseData(x)
+                        .SetName(x.FullName)
+                    )
+                    .GetEnumerator();
+            }
+
+            IEnumerator IEnumerable.GetEnumerator()
+            {
+                return this.GetEnumerator();
+            }
+        }
     }
 }
 

src/UnitTests/OptimizationTests/TestBfgsMinimizer.cs → src/UnitTests/OptimizationTests/BfgsMinimizerTests.cs

@@ -1,12 +1,17 @@
 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;
 
 namespace MathNet.Numerics.UnitTests.OptimizationTests
 {
     [TestFixture]
-    public class TestBfgsMinimizer
+    public class BfgsMinimizerTests
     {
         [Test]
         public void FindMinimum_Rosenbrock_Easy()
@@ -73,5 +78,53 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
 			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/ConjugateGradientMinimizerTests.cs

@@ -0,0 +1,100 @@
+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;
+
+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/TestGoldenSectionMinimizer.cs → src/UnitTests/OptimizationTests/GoldenSectionMinimizerTests.cs

@@ -5,7 +5,7 @@ using NUnit.Framework;
 namespace MathNet.Numerics.UnitTests.OptimizationTests
 {
     [TestFixture]
-    public class TestGoldenSectionMinimizer
+    public class GoldenSectionMinimizerTests
     {
         [Test]
         public void Test_Works()

src/UnitTests/OptimizationTests/NelderMeadSimplexTests.cs

@@ -30,8 +30,12 @@
 
 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;
 
 namespace MathNet.Numerics.UnitTests.OptimizationTests
 {
@@ -65,5 +69,64 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
             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/TestNewtonMinimizer.cs → src/UnitTests/OptimizationTests/NewtonMinimizerTests.cs

@@ -3,6 +3,10 @@ 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;
 
 namespace MathNet.Numerics.UnitTests.OptimizationTests
 {
@@ -51,7 +55,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
     }
 
     [TestFixture]
-    public class TestNewtonMinimizer
+    public class NewtonMinimizerTests
     {
         [Test]
         public void FindMinimum_Rosenbrock_Easy()
@@ -118,5 +122,66 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
             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/TestRosenbrockFunction.cs → src/UnitTests/OptimizationTests/RosenbrockFunctionTests.cs

@@ -5,7 +5,7 @@ using NUnit.Framework;
 namespace MathNet.Numerics.UnitTests.OptimizationTests
 {
     [TestFixture]
-    class TestRosenbrockFunction
+    class RosenbrockFunctionTests
     {
         [Test]
         public void TestGradient()

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

@@ -1,33 +0,0 @@
-using System;
-using MathNet.Numerics.LinearAlgebra.Double;
-using MathNet.Numerics.Optimization;
-using NUnit.Framework;
-
-namespace MathNet.Numerics.UnitTests.OptimizationTests
-{
-    [TestFixture]
-    public class TestConjugateGradientMinimizer
-    {
-        [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));
-        }
-    }
-}

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

@@ -345,16 +345,32 @@
     <Compile Include="EuclidTests\IntegerTheoryTest.cs" />
     <Compile Include="LinearAlgebraTests\MatrixStorageCombinatorsTests.cs" />
     <Compile Include="LinearAlgebraTests\VectorStorageCombinatorsTests.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\TestGoldenSectionMinimizer.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\TestBfgsMinimizer.cs" />
-    <Compile Include="OptimizationTests\TestConjugateGradientMinimizer.cs" />
-    <Compile Include="OptimizationTests\TestNewtonMinimizer.cs" />
-    <Compile Include="OptimizationTests\TestRosenbrockFunction.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" />
@@ -401,7 +417,7 @@
     <Compile Include="StatisticsTests\StatTestData.cs" />
     <Compile Include="TrigonometryTest.cs" />
     <Compile Include="UseLinearAlgebraProvider.cs" />
-    <Compile Include="OptimizationTests\TestBfgsBMinimizer.cs" />
+    <Compile Include="OptimizationTests\BfgsBMinimizerTests.cs" />
   </ItemGroup>
   <ItemGroup>
     <None Include="..\..\data\Codeplex-5667.csv">