diff --git a/src/Numerics/Optimization/NonLinearLeastSquaresMinimizer.cs b/src/Numerics/Optimization/NonLinearLeastSquaresMinimizer.cs
index 533869b5..5611477c 100644
--- a/src/Numerics/Optimization/NonLinearLeastSquaresMinimizer.cs
+++ b/src/Numerics/Optimization/NonLinearLeastSquaresMinimizer.cs
@@ -13,14 +13,51 @@ namespace MathNet.Numerics.Optimization
     /// </summary>
     public class NonLinearLeastSquaresMinimizer
     {
+        public class Options
+        {
+            public int MaximumIterations = 1000;
+
+            public int MaximumTrialStepIterations = 100;
+
+            public ConvergenceType 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>
-        /// Criterion0: Δ &lt; eps(0) (trust region solvers only)
-        /// Criterion1: ||F(x)||2 &lt; eps(1)
-        /// Criterion2: The Jacobian matrix is singular.||J(x)(1:m,j)||2 &lt; eps(2), j = 1, ..., n
-        /// Criterion3: ||s||2 &lt; eps(3)
-        /// Criterion4: ||F(x)||2 - ||F(x) - J(x)s||2 &lt; eps(4)
+        /// For details of convergence criteria, see Options.
         /// </summary>
-        public enum ConvergenceType { NoneMaxIterationExceeded, Criterion0, Criterion1, Criterion2, Criterion3, Criterion4, SingularJacobian, Error };
+        public enum ConvergenceType { NoneMaxIterationExceeded, Criterion0, Criterion1, Criterion2, Criterion3, Criterion4, Error };
         
         public class Result
         {
@@ -36,7 +73,7 @@ namespace MathNet.Numerics.Optimization
         /// <param name="x"></param>
         /// <param name="y"></param>
         /// <param name="f"></param>
-        /// <param name="pStart"></param>
+        /// <param name="pStart">Initial guess of parameters, p.</param>
         /// <param name="jacobian">jac_j(x, p) = df / dp_j</param>
         /// <returns></returns>
         public static double[] CurveFit(double[] x, double[] y, Func<double, double[], double> f,
diff --git a/src/Numerics/Providers/Optimization/IOptimizationProvider.cs b/src/Numerics/Providers/Optimization/IOptimizationProvider.cs
index b5463849..3e2400ba 100644
--- a/src/Numerics/Providers/Optimization/IOptimizationProvider.cs
+++ b/src/Numerics/Providers/Optimization/IOptimizationProvider.cs
@@ -57,7 +57,7 @@ namespace MathNet.Numerics.Providers.Optimization
         where T : struct
     {
         NonLinearLeastSquaresMinimizer.Result NonLinearLeastSquaresUnboundedMinimize(
-            int residualsLength, T[] initialGuess, LeastSquaresForwardModel function, 
-            out T[] parameters, Jacobian jacobianFunction = null);
+            int residualsLength, T[] initialGuess, LeastSquaresForwardModel function,
+            out T[] parameters, Jacobian jacobianFunction = null, NonLinearLeastSquaresMinimizer.Options options = null);
     }
 }
diff --git a/src/Numerics/Providers/Optimization/Mkl/MklOptimizationProvider.cs b/src/Numerics/Providers/Optimization/Mkl/MklOptimizationProvider.cs
index 571ca076..eeac095e 100644
--- a/src/Numerics/Providers/Optimization/Mkl/MklOptimizationProvider.cs
+++ b/src/Numerics/Providers/Optimization/Mkl/MklOptimizationProvider.cs
@@ -42,8 +42,9 @@ namespace MathNet.Numerics.Providers.Optimization.Mkl
     {
         const int TR_SUCCESS = 1501;
 
-        public NonLinearLeastSquaresMinimizer.Result NonLinearLeastSquaresUnboundedMinimize(int residualsLength, double[] initialGuess, LeastSquaresForwardModel function, out double[] parameters, Jacobian jacobianFunction = null)
+        public NonLinearLeastSquaresMinimizer.Result NonLinearLeastSquaresUnboundedMinimize(int residualsLength, double[] initialGuess, LeastSquaresForwardModel function, out double[] parameters, Jacobian jacobianFunction = null, NonLinearLeastSquaresMinimizer.Options options = null)
         {
+            if (options == null) options = new NonLinearLeastSquaresMinimizer.Options();
             bool analyticJacobian = jacobianFunction != null;
             double[] residuals = new double[residualsLength];
             double[] residualsMinus = new double[residualsLength];
@@ -53,15 +54,15 @@ namespace MathNet.Numerics.Providers.Optimization.Mkl
 
             double[] eps = new double[6]; // stop criteria
             int i;
-            for (i = 0; i < 6; i++)
-                eps[i] = 1e-8;
+            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 = 1000, maxTrialStepIterations = 100;
+            int maxIterations = options.MaximumIterations, maxTrialStepIterations = options.MaximumTrialStepIterations;
 
             IntPtr solverHandle = IntPtr.Zero;
             IntPtr jacobianHandle = IntPtr.Zero;
@@ -70,7 +71,7 @@ namespace MathNet.Numerics.Providers.Optimization.Mkl
 
             double initialStepBound = 0.0;
 
-            double jacobianPrecision = 1e-8;
+            double jacobianPrecision = options.JacobianPrecision;
 
             // zero initial values:
             for (i = 0; i < residuals.Length; i++)
@@ -185,7 +186,7 @@ namespace MathNet.Numerics.Providers.Optimization.Mkl
                 case -3:
                     convergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.Criterion1; break;
                 case -4:
-                    convergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.SingularJacobian; break;
+                    convergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.Criterion2; break;
                 case -5:
                     convergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.Criterion3; break;
                 case -6:
@@ -193,7 +194,7 @@ namespace MathNet.Numerics.Providers.Optimization.Mkl
             }
 
             // no errors, find reason for stopping;
-            return new NonLinearLeastSquaresMinimizer.Result() { ConvergenceType = convergenceType };
+            return new NonLinearLeastSquaresMinimizer.Result() { ConvergenceType = convergenceType, NumberOfIterations = iterations };
         }
 
         public static NonLinearLeastSquaresMinimizer.Result ErrorResult()