Optimization: cosmetics

Add IMinimizer interface so that same curve fitting approach can be used
for any function minimization algorithm.

src/NativeProviders/MKL/optimization.cpp

@@ -0,0 +1,227 @@
+#include <algorithm>
+#include <complex>
+#define MKL_Complex8 std::complex<float>
+#define MKL_Complex16 std::complex<double>
+
+#include "mkl_lapack.h"
+#include "mkl_cblas.h"
+#include "lapack_common.h"
+#include "wrapper_common.h"
+#include "mkl_lapacke.h"
+#include "mkl.h"
+#include "mkl_trans.h"
+
+template<typename T>
+inline MKL_INT lu_factor_2(MKL_INT m, T a[], MKL_INT ipiv[],
+                         void (*getrf)(const MKL_INT*, const MKL_INT*, T*, const MKL_INT*, MKL_INT*, MKL_INT*))
+{
+    std::complex<double> x = 5;
+	MKL_INT info = 0;
+    getrf(&m, &m, a, &m, ipiv, &info);
+    shift_ipiv_down(m, ipiv);
+    return info;
+};
+
+extern "C" {
+
+//_Mkl_Api(MKL_INT,dtrnlsp_init,(_TRNSP_HANDLE_t*, MKL_INT*, MKL_INT*, double*, double*, MKL_INT*, MKL_INT*, double*))
+//_Mkl_Api(MKL_INT,dtrnlsp_check,(_TRNSP_HANDLE_t*, MKL_INT*, MKL_INT*, double*, double*, double*, MKL_INT*))
+//_Mkl_Api(MKL_INT,dtrnlsp_solve,(_TRNSP_HANDLE_t*, double*, double*, MKL_INT*))
+//_Mkl_Api(MKL_INT,dtrnlsp_get,(_TRNSP_HANDLE_t*, MKL_INT*, MKL_INT*, double*, double*))
+//_Mkl_Api(MKL_INT,dtrnlsp_delete,(_TRNSP_HANDLE_t*))
+//
+//_Mkl_Api(MKL_INT,dtrnlspbc_init,(_TRNSPBC_HANDLE_t*, MKL_INT*, MKL_INT*, double*, double*, double*, double*, MKL_INT*, MKL_INT*, double*))
+//_Mkl_Api(MKL_INT,dtrnlspbc_check,(_TRNSPBC_HANDLE_t*, MKL_INT*, MKL_INT*, double*, double*, double*, double*, double*, MKL_INT*))
+//_Mkl_Api(MKL_INT,dtrnlspbc_solve,(_TRNSPBC_HANDLE_t*, double*, double*, MKL_INT*))
+//_Mkl_Api(MKL_INT,dtrnlspbc_get,(_TRNSPBC_HANDLE_t*, MKL_INT*, MKL_INT*, double*, double*))
+//_Mkl_Api(MKL_INT,dtrnlspbc_delete,(_TRNSPBC_HANDLE_t*))
+//
+//_Mkl_Api(MKL_INT,djacobi_init,(_JACOBIMATRIX_HANDLE_t*, MKL_INT*, MKL_INT*, double*, double*, double*))
+//_Mkl_Api(MKL_INT,djacobi_solve,(_JACOBIMATRIX_HANDLE_t*, double*, double*, MKL_INT*))
+//_Mkl_Api(MKL_INT,djacobi_delete,(_JACOBIMATRIX_HANDLE_t*))
+//_Mkl_Api(MKL_INT,djacobi,(USRFCND fcn, MKL_INT*, MKL_INT*, double*, double*, double*))
+//_Mkl_Api(MKL_INT,djacobix,(USRFCNXD fcn, MKL_INT*, MKL_INT*, double*, double*, double*,void*))
+
+	typedef double (__stdcall *function)();
+
+	DLLEXPORT double Test(function func, double input[])
+	{
+		input[0] = 5;
+		return func();
+	}
+
+    DLLEXPORT MKL_INT unbound_nonlinearleastsq_init(_TRNSPBC_HANDLE_t* handle, MKL_INT n, MKL_INT m, double x[], double eps[], MKL_INT iter1, MKL_INT iter2, double rs)
+    {
+		return dtrnlsp_init(handle, &n, &m, x, eps, &iter1, &iter2, &rs);
+	}
+
+    DLLEXPORT MKL_INT unbound_nonlinearleastsq_check(_TRNSPBC_HANDLE_t* handle, MKL_INT n, MKL_INT m, double fjac[], double fvec[], double eps[], int info[])
+    {
+		return dtrnlsp_check(handle, &n, &m, fjac, fvec, eps, info);
+	}
+
+	DLLEXPORT MKL_INT unbound_nonlinearleastsq_solve(_TRNSPBC_HANDLE_t* handle, double fvec[], double fjac[], int* RCI_Request)
+    {
+		return dtrnlsp_solve(handle, fvec, fjac, RCI_Request);
+    }
+
+	DLLEXPORT MKL_INT unbound_nonlinearleastsq_get(_TRNSPBC_HANDLE_t* handle, MKL_INT* iter, MKL_INT* st_cr, double* r1, double* r2)
+    {
+		return dtrnlsp_get(handle, iter, st_cr, r1, r2);
+    }
+
+	DLLEXPORT MKL_INT unbound_nonlinearleastsq_delete(_TRNSPBC_HANDLE_t* handle)
+    {
+		return 	dtrnlsp_delete(handle);
+    }
+
+	DLLEXPORT MKL_INT jacobi_init(_JACOBIMATRIX_HANDLE_t* handle, MKL_INT n, MKL_INT m, double x[], double fjac[], double eps)
+	{
+		return djacobi_init(handle, &n, &m, x, fjac, &eps);
+	}
+
+	DLLEXPORT MKL_INT jacobi_solve(_JACOBIMATRIX_HANDLE_t* handle, double f1[], double f2[], MKL_INT* RCI_Request)
+	{
+		return djacobi_solve(handle, f1, f2, RCI_Request);
+	}
+
+	DLLEXPORT MKL_INT jacobi_delete(_TRNSPBC_HANDLE_t* handle)
+    {
+		return djacobi_delete(handle);
+    }
+
+	DLLEXPORT void FreeBuffers()
+	{
+		MKL_FreeBuffers();
+	}
+
+	DLLEXPORT MKL_INT unbound_nonlinearleastsq(double parameters[], double parametersInitialGuess[], MKL_INT parametersLength, 
+		double residuals[], MKL_INT residualsLength, 
+		double jacobian[], // size parametersLength * residualsLength
+		double residualsMinus[], function updateResidualsMinus,
+		double residualsPlus[], function updateResidualsPlus,
+		function updateResiduals)
+	{
+		double eps[6]; // stop criteria
+		MKL_INT i;
+		for (i = 0; i < 6; i++)
+			eps[i] = 0.00001;
+
+		MKL_INT RCI_Request;     
+		MKL_INT successful;
+
+		MKL_INT maxIterations = 1000, maxTrialStepIterations = 100;
+
+		double rs = 0.0;
+		_TRNSP_HANDLE_t solverHandle;  
+		_JACOBIMATRIX_HANDLE_t jacobianHandle;
+
+		MKL_INT info[6]; // for parameter checking
+
+		double initialStepBound = 0.0;
+
+		double jacobianPrecision = 0.000001;
+
+		MKL_INT iteration;
+
+		double initialResidual, finalResidual;
+
+		// zero initial values:
+		for (i = 0; i < residualsLength; i++)
+			residuals[i] = 0.0;
+		for (i = 0; i < residualsLength * parametersLength; i++)
+			jacobian[i] = 0.0;
+
+		if (dtrnlsp_init(&solverHandle, &parametersLength, &residualsLength, parameters, eps, &maxIterations, &maxTrialStepIterations, &initialStepBound) !=
+			TR_SUCCESS)
+		{                             
+			MKL_FreeBuffers();
+			return 1;
+		}
+	 
+		if (dtrnlsp_check(&solverHandle, &parametersLength, &residualsLength, jacobian, residuals, eps, info) != TR_SUCCESS)
+		{                               
+			MKL_FreeBuffers();
+			return 1;
+		}
+		else
+		{
+			if (info[0] != 0 || // Handle invalid
+				info[1] != 0 || // Jacobian array not valid
+				info[2] != 0 || // Parameters array not valid
+				info[3] != 0)   // Eps array not valid
+			{                                  
+				MKL_FreeBuffers();
+				return 1;
+			}
+		}
+
+		if (djacobi_init(&jacobianHandle, &parametersLength, &residualsLength, parameters, jacobian, &jacobianPrecision) != TR_SUCCESS)
+		{
+			MKL_FreeBuffers ();
+			return 1;
+		}
+
+		RCI_Request = 0;
+		successful = 0;
+		while (successful == 0)
+		{
+			if (dtrnlsp_solve(&solverHandle, residuals, jacobian, &RCI_Request) != TR_SUCCESS)
+			{                             
+				MKL_FreeBuffers();
+				return 1;
+			}
+			if (RCI_Request == -1 || RCI_Request == -2 || RCI_Request == -3 ||
+				RCI_Request == -4 || RCI_Request == -5 || RCI_Request == -6)
+				successful = 1;
+			if (RCI_Request == 1) // recalculate function to update parameters
+			{
+				updateResiduals();
+			}
+			if (RCI_Request == 2)
+			{
+				MKL_INT rci_request = 0;
+				MKL_INT jacobianSuccessful = 0;
+				
+				// update Jacobian matrix:
+				while (jacobianSuccessful == 0)
+				{
+					if (djacobi_solve (&jacobianHandle, residualsPlus, residualsMinus, &rci_request) != TR_SUCCESS)
+					{
+						MKL_FreeBuffers ();
+						return 1;
+					}
+					if (rci_request == 1)
+						updateResidualsPlus();
+					else if (rci_request == 2)
+						updateResidualsMinus();
+					else if (rci_request == 0)
+						jacobianSuccessful = 1;
+				}
+			}
+		}
+
+		MKL_INT stopCriterionNumber, iterations;
+		if (dtrnlsp_get (&solverHandle, &iterations, &stopCriterionNumber, &initialResidual, &finalResidual) != TR_SUCCESS)
+		{                                  
+			MKL_FreeBuffers();
+			return 1;
+		}
+
+		if (dtrnlsp_delete (&solverHandle) != TR_SUCCESS)
+		{                                  
+			MKL_FreeBuffers();
+			return 1;
+		}
+
+		if (djacobi_delete (&jacobianHandle) != TR_SUCCESS)
+		{
+			MKL_FreeBuffers ();
+			return 1;
+		}
+                          
+		MKL_FreeBuffers ();
+		return 0;
+	}
+  
+}

src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj

@@ -292,6 +292,7 @@
     <ClCompile Include="..\..\Common\WindowsDLL.cpp" />
     <ClCompile Include="..\..\MKL\blas.c" />
     <ClCompile Include="..\..\MKL\lapack.cpp" />
+    <ClCompile Include="..\..\MKL\optimization.cpp" />
     <ClCompile Include="..\..\MKL\vector_functions.c" />
   </ItemGroup>
   <ItemGroup>

src/Numerics/Numerics.csproj

@@ -87,6 +87,10 @@
     <Reference Include="System.Xml" />
   </ItemGroup>
   <ItemGroup>
+    <Compile Include="Optimization\BrentMinimizer.cs" />
+    <Compile Include="Optimization\IMinimizer.cs" />
+    <Compile Include="Optimization\NonLinearLeastSquaresMinimizer.cs" />
+    <Compile Include="Optimization\PowellMinimizer.cs" />
     <Compile Include="Distributions\Triangular.cs" />
     <Compile Include="Euclid.cs" />
     <Compile Include="Generate.cs" />
@@ -186,6 +190,9 @@
     <Compile Include="LinearAlgebra\Vector.BCL.cs" />
     <Compile Include="LinearAlgebra\Vector.Operators.cs" />
     <Compile Include="NonConvergenceException.cs" />
+    <Compile Include="Providers\Optimization\IOptimizationProvider.cs" />
+    <Compile Include="Providers\Optimization\Mkl\MklOptimizationProvider.cs" />
+    <Compile Include="Providers\Optimization\Mkl\SafeNativeMethods.cs" />
     <Compile Include="Random\SystemRandomSource.cs" />
     <Compile Include="Random\RandomSeed.cs" />
     <Compile Include="RootFinding\Broyden.cs" />

src/Numerics/Optimization/BrentMinimizer.cs

@@ -0,0 +1,321 @@
+// <copyright file="BrentMinimizer.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2013 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using System;
+
+namespace MathNet.Numerics.Optimization
+{
+    /// <summary>
+    /// Options for Brent Minimization.
+    /// </summary>
+    public class BrentOptions
+    {
+        public int MaximumIterations = 1000;
+        public double Tolerance = 1e-4;
+    }
+
+    /// <summary>
+    /// Result of Brent Minimization.
+    /// </summary>
+    public class BrentResult
+    {
+        public int NumberOfIterations;
+        public double MinimumPoint;
+        public double MinimumFunctionValue;
+    }
+
+    /// <summary>
+    /// Minimizes f(p) where p is a model parameter scalar, i.e. a line-search.
+    /// Inspired by the SciPy implementation.
+    /// </summary>
+    public class BrentMinimizer
+    {
+        public struct Bracket
+        {
+            public double PointA;
+            public double PointB;
+            public double PointC;
+            public double FunctionA;
+            public double FunctionB;
+            public double FunctionC;
+        }
+
+        public BrentResult Result { get; private set; }
+
+        public readonly BrentOptions Options = new BrentOptions();
+
+        const double VerySmallNumber = 1e-21, GoldenRatio = 1.618034, MinimumTolerance = 1.0e-11;
+        const double GrowLimit = 110.0, ConjugateGradient = 0.3819660;
+
+        /// <summary>
+        /// Find the minimum of the supplied function using the Brent method.
+        /// </summary>
+        /// <param name="function"></param>
+        /// <returns></returns>
+        public double Minimize(Func<double, double> function)
+        {
+            Bracket bracket;
+            UpdateBracketInterval(function, new Bracket { PointA = 0, PointB = 1 }, out bracket);
+            int iterations = 0;
+            double x, w, v, fx, fw, fv;
+            double a, b, deltax, rat;
+            double cg = ConjugateGradient;
+            x = w = v = bracket.PointB; // x is the point with lowest function value encountered
+            fw = fv = fx = function(x);
+            if (bracket.PointA < bracket.PointC)
+            {
+                a = bracket.PointA;
+                b = bracket.PointC;
+            }
+            else
+            {
+                a = bracket.PointC;
+                b = bracket.PointA;
+            }
+            deltax = 0.0;
+            rat = 0;
+            while (iterations < Options.MaximumIterations)
+            {
+                double tol1, tol2, xmid;
+                double temp1, temp2, p;
+                double u, fu, dx_temp;
+                tol1 = Options.Tolerance*Math.Abs(x) + MinimumTolerance;
+                tol2 = 2.0*tol1;
+                xmid = 0.5*(a + b);
+                if (Math.Abs(x - xmid) < (tol2 - 0.5*(b - a))) // check for convergence
+                    break;
+                if (Math.Abs(deltax) <= tol1)
+                {
+                    if (x >= xmid) deltax = a - x; // do a golden section step
+                    else deltax = b - x;
+                    rat = cg*deltax;
+                }
+                else // do a parabolic step
+                {
+                    temp1 = (x - w)*(fx - fv);
+                    temp2 = (x - v)*(fx - fw);
+                    p = (x - v)*temp2 - (x - w)*temp1;
+                    temp2 = 2.0*(temp2 - temp1);
+                    if (temp2 > 0.0) p = -p;
+                    temp2 = Math.Abs(temp2);
+                    dx_temp = deltax;
+                    deltax = rat;
+                    // check parabolic fit
+                    if ((p > temp2*(a - x)) && (p < temp2*(b - x))
+                        && (Math.Abs(p) < Math.Abs(0.5*temp2*dx_temp)))
+                    {
+                        rat = p*1.0/temp2; // if parabolic step is useful.
+                        u = x + rat;
+                        if (((u - a) < tol2) || ((b - u) < tol2))
+                        {
+                            if ((xmid - x) >= 0) rat = tol1;
+                            else rat = -tol1;
+                        }
+                    }
+                    else
+                    {
+                        if (x >= xmid) deltax = a - x; // if it's not do a golden section step
+                        else deltax = b - x;
+                        rat = cg*deltax;
+                    }
+                }
+
+                if (Math.Abs(rat) < tol1) // update by at least tol1
+                {
+                    if (rat >= 0) u = x + tol1;
+                    else u = x - tol1;
+                }
+                else
+                {
+                    u = x + rat;
+                }
+                fu = function(u);
+                if (fu > fx) // if it's bigger than current
+                {
+                    if (u < x) a = u;
+                    else b = u;
+                    if ((fu <= fw) || (w == x))
+                    {
+                        v = w;
+                        w = u;
+                        fv = fw;
+                        fw = fu;
+                    }
+                    else if ((fu <= fv) || (v == x) || (v == w))
+                    {
+                        v = u;
+                        fv = fu;
+                    }
+                }
+                else
+                {
+                    if (u >= x) a = x;
+                    else b = x;
+                    v = w;
+                    w = x;
+                    x = u;
+                    fv = fw;
+                    fw = fx;
+                    fx = fu;
+                }
+                iterations++;
+            }
+            Result = new BrentResult { NumberOfIterations = iterations, MinimumPoint = x, MinimumFunctionValue = fx };
+            return x;
+        }
+
+        /// <summary>
+        /// Find the minimum of the supplied function along a specified line, using the Brent method.
+        /// </summary>
+        /// <param name="function"></param>
+        /// <param name="direction">Direction of line.</param>
+        /// <param name="startingPoint">Starting point of line.</param>
+        /// <returns></returns>
+        public double Minimize(Func<double[], double> function, double[] direction, double[] startingPoint, out double[] minimumPoint)
+        {
+            double[] point = new double[direction.Length];
+            Func<double, double> functionAlongLine = p =>
+            {
+                for (int i = 0; i < point.Length; ++i)
+                    point[i] = startingPoint[i] + direction[i]*p;
+                return function(point);
+            };
+            double result = Minimize(functionAlongLine);
+            minimumPoint = point;
+            return result;
+        }
+
+        /// <summary>
+        /// Updates the bracket.
+        /// </summary>
+        /// <param name="function"></param>
+        /// <param name="bracketInitial"></param>
+        /// <param name="newBracket"></param>
+        /// <returns></returns>
+        static bool UpdateBracketInterval(Func<double, double> function, Bracket bracketInitial, out Bracket newBracket)
+        {
+            double pointA = bracketInitial.PointA;
+            double pointB = bracketInitial.PointB;
+            const int maxIterations = 1000;
+            double functionA = function(pointA);
+            double functionB = function(pointB);
+            if (functionA < functionB) // Swap points over
+            {
+                double temp = functionA;
+                functionA = functionB;
+                functionB = temp;
+                temp = pointA;
+                pointA = pointB;
+                pointB = temp;
+            }
+            double pointC = pointB + GoldenRatio*(pointB - pointA);
+            double functionC = function(pointC);
+            int iterations = 0;
+            double temp1, temp2, value, denom;
+            while (functionC < functionB)
+            {
+                double functionW;
+                temp1 = (pointB - pointA)*(functionB - functionC);
+                temp2 = (pointB - pointC)*(functionB - functionA);
+                value = temp2 - temp1;
+                if (Math.Abs(value) < VerySmallNumber) denom = 2.0*VerySmallNumber;
+                else denom = 2.0*value;
+                double pointW = pointB - ((pointB - pointC)*temp2 - (pointB - pointA)*temp1)/denom;
+                double wlim = pointB + GrowLimit*(pointC - pointB);
+                if (iterations > maxIterations)
+                {
+                    newBracket = bracketInitial;
+                    return false;
+                }
+                iterations++;
+                if ((pointW - pointC)*(pointB - pointW) > 0.0)
+                {
+                    functionW = function(pointW);
+                    if (functionW < functionC)
+                    {
+                        pointA = pointB;
+                        pointB = pointW;
+                        functionA = functionB;
+                        functionB = functionW;
+                        break;
+                    }
+                    if (functionW > functionB)
+                    {
+                        pointC = pointW;
+                        functionC = functionW;
+                        break;
+                    }
+                    pointW = pointC + GoldenRatio*(pointC - pointB);
+                    functionW = function(pointW);
+                }
+                else if ((pointW - wlim)*(wlim - pointC) >= 0.0)
+                {
+                    pointW = wlim;
+                    functionW = function(pointW);
+                }
+                else if ((pointW - wlim)*(pointC - pointW) > 0.0)
+                {
+                    functionW = function(pointW);
+                    if (functionW < functionC)
+                    {
+                        pointB = pointC;
+                        pointC = pointW;
+                        pointW = pointC + GoldenRatio*(pointC - pointB);
+                        functionB = functionC;
+                        functionC = functionW;
+                        functionW = function(pointW);
+                    }
+                }
+                else
+                {
+                    pointW = pointC + GoldenRatio*(pointC - pointB);
+                    functionW = function(pointW);
+                }
+                pointA = pointB;
+                pointB = pointC;
+                pointC = pointW;
+                functionA = functionB;
+                functionB = functionC;
+                functionC = functionW;
+            }
+            newBracket = new Bracket
+            {
+                PointA = pointA,
+                PointB = pointB,
+                PointC = pointC,
+                FunctionA = functionA,
+                FunctionB = functionB,
+                FunctionC = functionC
+            };
+            return true;
+        }
+
+    }
+}

src/Numerics/Optimization/IMinimizer.cs

@@ -0,0 +1,72 @@
+// <copyright file="IMinimizer.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2013 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using System;
+
+namespace MathNet.Numerics.Optimization
+{
+    /// <summary>
+    /// Interface implemented by a class that minimizes f(p) where p is a vector of model parameters.
+    /// A class implenting this interface can then be used to solve curve fitting problems.
+    /// </summary>
+    public interface IMinimizer
+    {
+        double[] Minimize(Func<double[], double> function, double[] pInitialGuess);
+    }
+
+    public static class NonLinearCurveFitExtensions
+    {
+        /// <summary>
+        /// Non-Linear Least-Squares fitting the points (x,y) to a specified function of y : x -> f(x, p), p being a vector of parameters.
+        /// returning its best fitting parameters p.
+        /// </summary>
+        /// <param name="x"></param>
+        /// <param name="y"></param>
+        /// <param name="f"></param>
+        /// <param name="pStart">Initial guess of parameters, p.</param>
+        /// <returns></returns>
+        public static double[] CurveFit(this IMinimizer minimizer, double[] x, double[] y, Func<double, double[], double> f,
+            double[] pStart)
+        {
+            // Need to minimize sum of squares of residuals; create this function:
+            Func<double[], double> function = p =>
+            {
+                double sum = 0;
+                for (int i = 0; i < x.Length; ++i)
+                {
+                    double temp = y[i] - f(x[i], p);
+                    sum += temp*temp;
+                }
+                return sum;
+            };
+            return minimizer.Minimize(function, pStart);
+        }
+    }
+}

src/Numerics/Optimization/NonLinearLeastSquaresMinimizer.cs

@@ -0,0 +1,160 @@
+// <copyright file="NonLinearLeastSquaresMinimizer.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2013 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using System;
+using MathNet.Numerics.Providers.Optimization;
+
+#if NATIVEMKL
+using MathNet.Numerics.Providers.Optimization.Mkl;
+#endif
+
+namespace MathNet.Numerics.Optimization
+{
+    /// <summary>
+    /// Options for Non-Linear Least Squares Minimization.
+    /// </summary>
+    public class NonLinearLeastSquaresOptions
+    {
+        public int MaximumIterations = 1000;
+        public int MaximumTrialStepIterations = 100;
+        public NonLinearLeastSquaresConvergenceType ConvergenceType;
+
+        /// <summary>
+        ///  Convergence if Δ &lt; Criterion0, Δ is trust region size.
+        /// </summary>
+        public double Criterion0 = 1e-7;
+
+        /// <summary>
+        /// Convergence if |r|2 &lt; Criterion1, r is residuals vector.
+        /// </summary>
+        public double Criterion1 = 1e-7;
+
+        /// <summary>
+        /// Jacobian considered singular if |J(:,j)|2 &lt; Criterion2 for any j.
+        /// </summary>
+        public double Criterion2 = 1e-7;
+
+        /// <summary>
+        /// Jacobian considered singular if |s|2 &lt; Criterion3. s is trial step.
+        /// </summary>
+        public double Criterion3 = 1e-7;
+
+        /// <summary>
+        /// |r|2 - |r - Js|2 &lt; Criterion4
+        /// </summary>
+        public double Criterion4 = 1e-7;
+
+        public double TrialStepPrecision = 1e-10;
+
+        /// <summary>
+        /// Only if Jacobian calculated by central differences.
+        /// </summary>
+        public double JacobianPrecision = 1e-8;
+    }
+
+    /// <summary>
+    /// For details of convergence criteria, see Options.
+    /// </summary>
+    public enum NonLinearLeastSquaresConvergenceType
+    {
+        MaxIterationsExceeded,
+        Criterion0,
+        Criterion1,
+        Criterion2,
+        Criterion3,
+        Criterion4,
+        Error
+    };
+
+    /// <summary>
+    /// Result of Non-Linear Least Squares Minimization.
+    /// </summary>
+    public class NonLinearLeastSquaresResult
+    {
+        public int NumberOfIterations;
+        public NonLinearLeastSquaresConvergenceType ConvergenceType;
+    }
+
+#if NATIVEMKL
+
+    /// <summary>
+    /// This class is a special function minimizer that minimizes functions of the form
+    /// f(p) = |r(p)|^2 where r is a vector of residuals and p is a vector of model parameters.
+    /// </summary>
+    public class NonLinearLeastSquaresMinimizer // Note does not implement IMinimizer, since it curve fitting problems can be solved more efficiently.
+    {
+        public NonLinearLeastSquaresResult Result { get; private set; }
+
+        public readonly NonLinearLeastSquaresOptions Options = new NonLinearLeastSquaresOptions();
+
+        /// <summary>
+        /// Non-Linear Least-Squares fitting the points (x,y) to a specified function of y : x -> f(x, p), p being a vector of parameters.
+        /// returning its best fitting parameters p.
+        /// </summary>
+        /// <param name="x"></param>
+        /// <param name="y"></param>
+        /// <param name="f"></param>
+        /// <param name="pStart">Initial guess of parameters, p.</param>
+        /// <param name="jacobian">jac_j(x, p) = df / dp_j</param>
+        /// <returns></returns>
+        public double[] CurveFit(double[] x, double[] y, Func<double, double[], double> f,
+            double[] pStart, Func<double, double[], double[]> jacobian = null)
+        {
+            if (x.Length != y.Length) throw new ArgumentException("x and y lengths different");
+            var provider = new MklOptimizationProvider();
+            LeastSquaresForwardModel function = (p, r) =>
+            {
+                for (int i = 0; i < r.Length; ++i)
+                    r[i] = y[i] - f(x[i], p);
+            };
+
+            // jac is df_i / dp_j
+
+            Jacobian jacobianFunction = null;
+            if (jacobian != null)
+                jacobianFunction = (p, jac) =>
+                {
+                    for (int i = 0; i < y.Length; ++i)
+                    {
+                        double[] values = jacobian(x[i], p);
+                        for (int j = 0; j < values.Length; ++j)
+                            jac[j*y.Length + i] = -values[j];
+                    }
+                };
+
+            double[] parameters;
+            Result = provider.NonLinearLeastSquaresUnboundedMinimize(y.Length, pStart, function, out parameters, jacobianFunction);
+            return parameters;
+        }
+    }
+
+#endif
+
+}

src/Numerics/Optimization/PowellMinimizer.cs

@@ -0,0 +1,204 @@
+// <copyright file="PowellMinimizer.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2013 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using System;
+
+namespace MathNet.Numerics.Optimization
+{
+    /// <summary>
+    /// Options for Powell Minimization.
+    /// </summary>
+    public class PowellOptions
+    {
+        public int? MaximumIterations = null;
+        public int? MaximumFunctionCalls = null;
+        public double PointTolerance = 1e-4; // absolute
+        public double FunctionTolerance = 1e-4; // relative
+    }
+
+    public enum PowellConvergenceType
+    {
+        Success,
+        MaxIterationsExceeded,
+        MaxFunctionCallsExceeded
+    }
+
+    /// <summary>
+    /// Result of Powell Minimization.
+    /// </summary>
+    public class PowellResult
+    {
+        public int NumberOfIterations;
+        public int NumberOfFunctionCalls;
+        public double[] MinimumPoint;
+        public double MinimumFunctionValue;
+        public PowellConvergenceType ConvergenceType;
+    }
+
+    /// <summary>
+    /// Minimizes f(p) where p is a vector of model parameters using the Powell method.
+    /// </summary>
+    public class PowellMinimizer : IMinimizer
+    {
+        public PowellResult Result { get; private set; }
+
+        public readonly PowellOptions Options = new PowellOptions();
+
+        /// <summary>
+        /// Find the minimum of the supplied function using the Powell method.
+        /// </summary>
+        /// <param name="function"></param>
+        /// <param name="pInitialGuess"></param>
+        /// <returns></returns>
+        public double[] Minimize(Func<double[], double> function, double[] pInitialGuess)
+        {
+            BrentMinimizer brentMinimizer = new BrentMinimizer();
+            int n = pInitialGuess.Length; // number of dimensions
+            // used in closure:
+            double[] point = new double[n];
+            double[] startingPoint = new double[n];
+            double[] direction = new double[n];
+            double lineMiniumum = 0;
+            int functionCalls = 0;
+            Func<double, double> functionAlongLine = p =>
+            {
+                for (int i = 0; i < point.Length; ++i)
+                    point[i] = startingPoint[i] + direction[i]*p;
+                lineMiniumum = function(point);
+                functionCalls++;
+                return lineMiniumum;
+            };
+
+            double fval;
+
+            int iterations = 0;
+            int maxIterations = (Options.MaximumIterations == null) ? n*1000 : (int)Options.MaximumIterations;
+            int maxFunctionCalls = (Options.MaximumFunctionCalls == null) ? n*1000 : (int)Options.MaximumFunctionCalls;
+
+            // An array of n directions:
+            double[][] directionSet = new double[n][];
+            for (int i = 0; i < n; ++i)
+            {
+                directionSet[i] = new double[n];
+                directionSet[i][i] = 1.0;
+            }
+            double[] x = pInitialGuess;
+            double[] x1 = (double[])x.Clone();
+            double[] x2 = new double[n];
+            double[] direction1 = new double[n];
+
+            brentMinimizer.Options.Tolerance = Options.PointTolerance*100;
+
+            fval = function(x);
+
+            double fx;
+
+            while (true)
+            {
+                fx = fval;
+                int bigIndex = 0;
+                double delta = 0.0;
+                double fx2;
+                // loop over all directions:
+                for (int i = 0; i < n; ++i)
+                {
+                    fx2 = fval;
+
+                    // Do a linesearch along direction
+                    direction = directionSet[i];
+                    startingPoint = x;
+                    double u = brentMinimizer.Minimize(functionAlongLine);
+                    fval = functionAlongLine(u); // updates point
+
+                    for (int j = 0; j < n; ++j) x[j] = point[j];
+
+                    if ((fx2 - fval) > delta)
+                    {
+                        delta = fx2 - fval;
+                        bigIndex = i; // record direction index of the largest decrease seen
+                    }
+                }
+                iterations++;
+                if (2.0*(fx - fval) <= Options.FunctionTolerance*((Math.Abs(fx) + Math.Abs(fval)) + 1e-20)) break;
+                if (functionCalls >= maxFunctionCalls) break;
+                if (iterations >= maxIterations) break;
+
+                // Construct the extrapolated point
+                for (int i = 0; i < n; ++i)
+                {
+                    direction1[i] = x[i] - x1[i];
+                    x2[i] = 2.0*x[i] - x1[i];
+                    x1[i] = x[i];
+                }
+
+                fx2 = function(x2);
+                if (fx > fx2)
+                {
+                    double t = 2.0*(fx + fx2 - 2.0*fval);
+                    double temp = (fx - fval - delta);
+                    t *= temp*temp;
+                    temp = fx - fx2;
+                    t -= delta*temp*temp;
+                    if (t < 0.0)
+                    {
+                        // Do a linesearch along direction
+                        direction = direction1;
+                        startingPoint = x;
+                        double u = brentMinimizer.Minimize(functionAlongLine);
+                        fval = functionAlongLine(u); // updates point
+
+                        for (int i = 0; i < n; ++i)
+                        {
+                            directionSet[bigIndex][i] = directionSet[n - 1][i];
+                            directionSet[n - 1][i] = point[i] - x[i];
+                            x[i] = point[i];
+                        }
+                    }
+                }
+            }
+            var convergenceType = PowellConvergenceType.Success;
+            if (functionCalls >= maxFunctionCalls)
+                convergenceType = PowellConvergenceType.MaxFunctionCallsExceeded;
+            else if (iterations > maxIterations)
+                convergenceType = PowellConvergenceType.MaxFunctionCallsExceeded;
+
+            Result = new PowellResult
+            {
+                MinimumPoint = (double[])x.Clone(),
+                MinimumFunctionValue = fx,
+                ConvergenceType = convergenceType,
+                NumberOfIterations = iterations,
+                NumberOfFunctionCalls = functionCalls
+            };
+
+            return Result.MinimumPoint;
+        }
+    }
+}

src/Numerics/Providers/Optimization/IOptimizationProvider.cs

@@ -0,0 +1,63 @@
+// <copyright file="IOptimizationProvider.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2013 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using MathNet.Numerics.Optimization;
+
+namespace MathNet.Numerics.Providers.Optimization
+{
+    /// <summary>
+    /// Function specifying the model. This takes in model parameters, calculates residuals
+    /// and updates the residuals array with these.
+    /// </summary>
+    /// <param name="parameters">The model parameters. The function must not change these.</param>
+    /// <param name="r">The residuals to be updated. The existing array should be updated.</param>
+    /// <returns></returns>
+    public delegate void LeastSquaresForwardModel(double[] p, double[] r);
+
+    /// <summary>
+    /// Function providing the Jacobian matrix in column-major format for a set of model parameter values.
+    /// Jacobian is dr_i / dp_i, r being the residuals vector and p the vector of parameters
+    /// </summary>
+    /// <param name="p">THe model parameters. The function must not change these.</param>
+    /// <param name="jacobian">THe Jacobian matrix in column-major format.</param>
+    public delegate void Jacobian(double[] p, double[] jacobian);
+
+    /// <summary>
+    /// Interface to linear algebra algorithms that work off 1-D arrays.
+    /// </summary>
+    /// <typeparam name="T">Supported data types are Double, Single, Complex, and Complex32.</typeparam>
+    public interface IOptimizationProvider<T>
+        where T : struct
+    {
+        NonLinearLeastSquaresResult NonLinearLeastSquaresUnboundedMinimize(
+            int residualsLength, T[] initialGuess, LeastSquaresForwardModel function,
+            out T[] parameters, Jacobian jacobianFunction = null, NonLinearLeastSquaresOptions options = null);
+    }
+}

src/Numerics/Providers/Optimization/Mkl/MklOptimizationProvider.cs

@@ -0,0 +1,204 @@
+// <copyright file="MklOptimizationProvider.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2013 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+#if NATIVEMKL
+
+using System;
+using MathNet.Numerics.Optimization;
+
+namespace MathNet.Numerics.Providers.Optimization.Mkl
+{
+    public class MklOptimizationProvider : IOptimizationProvider<double>
+    {
+        const int TR_SUCCESS = 1501;
+
+        public NonLinearLeastSquaresResult NonLinearLeastSquaresUnboundedMinimize(int residualsLength, double[] initialGuess, LeastSquaresForwardModel function, out double[] parameters, Jacobian jacobianFunction = null, NonLinearLeastSquaresOptions options = null)
+        {
+            if (options == null) options = new NonLinearLeastSquaresOptions();
+            bool analyticJacobian = jacobianFunction != null;
+            double[] residuals = new double[residualsLength];
+            double[] residualsMinus = new double[residualsLength];
+            double[] residualsPlus = new double[residualsLength];
+            double[] jacobian = new double[residualsLength * initialGuess.Length];
+            parameters = new double[initialGuess.Length];
+
+            double[] eps = new double[6]; // stop criteria
+            int i;
+            eps[0] = options.Criterion0; eps[1] = options.Criterion1; eps[2] = options.Criterion2;
+            eps[3] = options.Criterion3; eps[4] = options.Criterion4; eps[5] = options.TrialStepPrecision;
+
+            for (i = 0; i < initialGuess.Length; i++)
+                parameters[i] = initialGuess[i];
+
+            int successful;
+
+            int maxIterations = options.MaximumIterations, maxTrialStepIterations = options.MaximumTrialStepIterations;
+
+            IntPtr solverHandle = IntPtr.Zero;
+            IntPtr jacobianHandle = IntPtr.Zero;
+
+            int[] info = new int[6]; // for parameter checking
+
+            double initialStepBound = 0.0;
+
+            double jacobianPrecision = options.JacobianPrecision;
+
+            // zero initial values:
+            for (i = 0; i < residuals.Length; i++)
+                residuals[i] = 0.0;
+            for (i = 0; i < residuals.Length * parameters.Length; i++)
+                jacobian[i] = 0.0;
+
+            if (SafeNativeMethods.unbound_nonlinearleastsq_init(ref solverHandle, parameters.Length, residualsLength, parameters, eps, maxIterations, maxTrialStepIterations, initialStepBound) !=
+                TR_SUCCESS)
+            {
+                SafeNativeMethods.FreeBuffers();
+                return ErrorResult();
+            }
+
+            if (SafeNativeMethods.unbound_nonlinearleastsq_check(ref solverHandle, parameters.Length, residualsLength, jacobian, residuals, eps, info) != TR_SUCCESS)
+            {
+                SafeNativeMethods.FreeBuffers();
+                return ErrorResult();
+            }
+            else
+            {
+                if (info[0] != 0 || // Handle invalid
+                    info[1] != 0 || // Jacobian array not valid
+                    info[2] != 0 || // Parameters array not valid
+                    info[3] != 0)   // Eps array not valid
+                {
+                    SafeNativeMethods.FreeBuffers();
+                    return ErrorResult();
+                }
+            }
+
+            if (SafeNativeMethods.jacobi_init(ref jacobianHandle, parameters.Length, residuals.Length, parameters, jacobian, jacobianPrecision) != TR_SUCCESS)
+            {
+                SafeNativeMethods.FreeBuffers();
+                return ErrorResult();
+            }
+
+            int rciRequest = 0;
+            successful = 0;
+            while (successful == 0)
+            {
+                if (SafeNativeMethods.unbound_nonlinearleastsq_solve(ref solverHandle, residuals, jacobian, ref rciRequest) != TR_SUCCESS)
+                {
+                    SafeNativeMethods.FreeBuffers();
+                    return ErrorResult();
+                }
+                if (rciRequest == -1 || rciRequest == -2 || rciRequest == -3 ||
+                    rciRequest == -4 || rciRequest == -5 || rciRequest == -6)
+                    successful = 1;
+                if (rciRequest == 1) // recalculate function to update parameters
+                {
+                    function(parameters, residuals);
+                }
+                if (rciRequest == 2)
+                {
+                    if (analyticJacobian)
+                        jacobianFunction(parameters, jacobian);
+                    else
+                    {
+                        // calculate by central differences:
+                        int rciRequestJacobian = 0;
+                        int jacobianSuccessful = 0;
+
+                        // update Jacobian matrix:
+                        while (jacobianSuccessful == 0)
+                        {
+                            if (SafeNativeMethods.jacobi_solve(ref jacobianHandle, residualsPlus, residualsMinus, ref rciRequestJacobian) != TR_SUCCESS)
+                            {
+                                SafeNativeMethods.FreeBuffers();
+                                return ErrorResult();
+                            }
+                            if (rciRequestJacobian == 1)
+                                function(parameters, residualsPlus);
+                            else if (rciRequestJacobian == 2)
+                                function(parameters, residualsMinus);
+                            else if (rciRequestJacobian == 0)
+                                jacobianSuccessful = 1;
+                        }
+                    }
+                }
+            }
+
+            int stopCriterionNumber = 0, iterations = 0;
+            double initialResidual = 0, finalResidual = 0;
+            if (SafeNativeMethods.unbound_nonlinearleastsq_get(ref solverHandle, ref iterations, ref stopCriterionNumber, ref initialResidual, ref finalResidual) != TR_SUCCESS)
+            {
+                SafeNativeMethods.FreeBuffers();
+                return ErrorResult();
+            }
+
+            if (SafeNativeMethods.unbound_nonlinearleastsq_delete(ref solverHandle) != TR_SUCCESS)
+            {
+                SafeNativeMethods.FreeBuffers();
+                return ErrorResult();
+            }
+
+            if (SafeNativeMethods.jacobi_delete(ref jacobianHandle) != TR_SUCCESS)
+            {
+                SafeNativeMethods.FreeBuffers();
+                return ErrorResult();
+            }
+
+            SafeNativeMethods.FreeBuffers();
+
+            NonLinearLeastSquaresConvergenceType convergenceType = NonLinearLeastSquaresConvergenceType.Error;
+            switch (rciRequest)
+            {
+                case -1:
+                    convergenceType = NonLinearLeastSquaresConvergenceType.MaxIterationsExceeded; break;
+                case -2:
+                    convergenceType = NonLinearLeastSquaresConvergenceType.Criterion0; break;
+                case -3:
+                    convergenceType = NonLinearLeastSquaresConvergenceType.Criterion1; break;
+                case -4:
+                    convergenceType = NonLinearLeastSquaresConvergenceType.Criterion2; break;
+                case -5:
+                    convergenceType = NonLinearLeastSquaresConvergenceType.Criterion3; break;
+                case -6:
+                    convergenceType = NonLinearLeastSquaresConvergenceType.Criterion4; break;
+            }
+
+            // no errors, find reason for stopping;
+            return new NonLinearLeastSquaresResult() { ConvergenceType = convergenceType, NumberOfIterations = iterations };
+        }
+
+        public static NonLinearLeastSquaresResult ErrorResult()
+        {
+            return new NonLinearLeastSquaresResult() { ConvergenceType = NonLinearLeastSquaresConvergenceType.Error };
+        }
+    }
+}
+
+#endif

src/Numerics/Providers/Optimization/Mkl/SafeNativeMethods.cs

@@ -0,0 +1,83 @@
+// <copyright file="SafeNativeMethods.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009-2013 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+#if NATIVEMKL
+
+using System.Runtime.InteropServices;
+using System.Security;
+using System;
+
+namespace MathNet.Numerics.Providers.Optimization.Mkl
+{
+    /// <summary>
+    /// P/Invoke methods to the native math libraries.
+    /// </summary>
+    [SuppressUnmanagedCodeSecurity]
+    [SecurityCritical]
+    internal static class SafeNativeMethods
+    {
+        /// <summary>
+        /// Name of the native DLL.
+        /// </summary>
+        const string DllName = "MathNet.Numerics.MKL.dll";
+
+        #region Non-Linear Least Squares Unbounded
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern int unbound_nonlinearleastsq_init(ref IntPtr handle, int n, int m, double[] x, double[] eps, int iter1, int iter2, double rs);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern int unbound_nonlinearleastsq_check(ref IntPtr handle, int n, int m, double[] fjac, double[] fvec, double[] eps, int[] info);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern int unbound_nonlinearleastsq_solve(ref IntPtr handle, double[] fvec, double[] fjac, ref int RCI_Request);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern int unbound_nonlinearleastsq_get(ref IntPtr handle, ref int iter, ref int st_cr, ref double r1, ref double r2);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern int unbound_nonlinearleastsq_delete(ref IntPtr handle);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern int jacobi_init(ref IntPtr handle, int n, int m, double[] x, double[] fjac, double eps);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern int jacobi_solve(ref IntPtr handle, double[] f1, double[] f2, ref int RCI_Request);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern int jacobi_delete(ref IntPtr handle);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern int FreeBuffers();
+
+        #endregion
+
+    }
+}
+
+#endif

src/UnitTests/OptimizationTests/FunctionMinimizationTests.cs

@@ -0,0 +1,71 @@
+// <copyright file="FunctionMinimizationTests.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2013 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using System;
+using System.Diagnostics;
+using MathNet.Numerics.Optimization;
+using NUnit.Framework;
+
+namespace MathNet.Numerics.UnitTests.OptimizationTests
+{
+    [TestFixture]
+    public class FunctionMinimizationTests
+    {
+        [Test]
+        public void PowellCurveFit()
+        {
+            // y = b1*(1-exp[-b2*x])  +  e
+            var xin = new double[] { 1, 2, 3, 5, 7, 10 };
+            var yin = new double[] { 109, 149, 149, 191, 213, 224 };
+
+            Func<double, double[], double> function = (x, p) => p[0]*(1 - Math.Exp(-p[1]*x));
+
+            var minimizer = new PowellMinimizer();
+
+            var watch = Stopwatch.StartNew();
+            double[] popt = null;
+            minimizer.Options.PointTolerance = 1e-8;
+
+            for (int i = 0; i < 1000; ++i)
+            {
+                popt = minimizer.CurveFit(xin, yin, function, new double[] { 1, 1 }); // 100, 0.75
+            }
+            double elapsed = watch.ElapsedMilliseconds;
+
+            double[] expected = new double[] { 2.1380940889E+02, 5.4723748542E-01 };
+
+            double residual = 0;
+            for (int i = 0; i < yin.Length; ++i) residual += (yin[i] - function(xin[i], popt))*(yin[i] - function(xin[i], popt));
+
+            Assert.AreEqual(expected[0], popt[0], 1e-4);
+            Assert.AreEqual(expected[1], popt[1], 1e-4);
+        }
+    }
+}

src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs

@@ -0,0 +1,73 @@
+// <copyright file="NonLinearLeastSquaresTest.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2013 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+#if NATIVEMKL
+
+using System;
+using MathNet.Numerics.Optimization;
+using NUnit.Framework;
+
+namespace MathNet.Numerics.UnitTests.OptimizationTests
+{
+    [TestFixture]
+    public class NonLinearLeastSquaresTest
+    {
+        [Test]
+        public void CurveFit()
+        {
+            var minimizer = new NonLinearLeastSquaresMinimizer();
+            minimizer.Options.MaximumIterations = 1010;
+            minimizer.Options.Criterion4 = 1e-8;
+            // y = b1*(1-exp[-b2*x])  +  e
+            var xin = new double[] { 1, 2, 3, 5, 7, 10 };
+            var yin = new double[] { 109, 149, 149, 191, 213, 224 };
+
+            // estimated derivative method: does not find best solution.
+            //var popt = minimizer.CurveFit(xin, yin, (x, p) => p[0] * (1 - Math.Exp(-p[1] * x)), new double[] { 1, 1 });
+
+            Func<double, double[], double> function = (x, p) => p[0] * (1 - Math.Exp(-p[1] * x));
+            Func<double, double[], double[]> jacobian = (x, p) => new double[] {
+                1 - Math.Exp(-p[1] * x),
+                p[0] * x * Math.Exp(-p[1] * x) };
+
+            var popt = minimizer.CurveFit(xin, yin, function, new double[] { 1, 1 }, jacobian); // 100, 0.75
+
+            double[] expected = new double[] { 2.1380940889E+02, 5.4723748542E-01 };
+
+            double residual = 0;
+            for (int i = 0; i < yin.Length; ++i) residual += (yin[i] - function(xin[i], popt)) * (yin[i] - function(xin[i], popt));
+
+            Assert.AreEqual(expected[0], popt[0], 1e-3);
+            Assert.AreEqual(expected[1], popt[1], 1e-3);
+        }
+    }
+}
+
+#endif

src/UnitTests/UnitTests.csproj

@@ -355,6 +355,8 @@
     <Compile Include="LinearAlgebraTests\Double\VectorArithmeticTheory.cs" />
     <Compile Include="LinearAlgebraTests\VectorArithmeticTheory.cs" />
     <Compile Include="MatrixHelpers.cs" />
+    <Compile Include="OptimizationTests\FunctionMinimizationTests.cs" />
+    <Compile Include="OptimizationTests\NonLinearLeastSquaresTest.cs" />
     <Compile Include="EuclidTests\GcdRelatedTest.cs" />
     <Compile Include="EuclidTests\GcdRelatedTestBigInteger.cs" />
     <Compile Include="EuclidTests\IntegerTheoryTest.cs" />