us1zcljz
/
mathnet-numerics
forked from tsai/mathnet-numerics
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

NonLinearLeastSquaresMinimizer change API

optimization-1
joemoorhouse 13 years ago
parent
59cebe9f01
commit
73c61e7e4f
4 changed files with 77 additions and 63 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 105
      src/Numerics/Optimization/NonLinearLeastSquaresMinimizer.cs
  2. 4
      src/Numerics/Providers/Optimization/IOptimizationProvider.cs
  3. 24
      src/Numerics/Providers/Optimization/Mkl/MklOptimizationProvider.cs
  4. 7
      src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs

105
src/Numerics/Optimization/NonLinearLeastSquaresMinimizer.cs
View File

@ -6,66 +6,77 @@ using MathNet.Numerics.Providers.Optimization;
using MathNet.Numerics.Providers.Optimization.Mkl; using MathNet.Numerics.Providers.Optimization.Mkl;
namespace MathNet.Numerics.Optimization namespace MathNet.Numerics.Optimization
{ {
/// <summary> /// <summary>
/// This class is a special function minimizer that minimizes functions of the form /// Options for Non-Linear Least Squares Minimization.
/// f(p) = |r(p)|^2 where r is a vector of residuals and p is a vector of model parameters.
/// </summary> /// </summary>
public class NonLinearLeastSquaresMinimizer public class NonLinearLeastSquaresOptions
{ {
public class Options public int MaximumIterations = 1000;
{
public int MaximumIterations = 1000;
public int MaximumTrialStepIterations = 100; public int MaximumTrialStepIterations = 100;
public ConvergenceType ConvergenceType; public NonLinearLeastSquaresConvergenceType ConvergenceType;
/// <summary> /// <summary>
/// Convergence if Δ &lt; Criterion0, Δ is trust region size. /// Convergence if Δ &lt; Criterion0, Δ is trust region size.
/// </summary> /// </summary>
public double Criterion0 = 1e-7; public double Criterion0 = 1e-7;
/// <summary> /// <summary>
/// Convergence if |r|2 &lt; Criterion1, r is residuals vector. /// Convergence if |r|2 &lt; Criterion1, r is residuals vector.
/// </summary> /// </summary>
public double Criterion1 = 1e-7; public double Criterion1 = 1e-7;
/// <summary> /// <summary>
/// Jacobian considered singular if |J(:,j)|2 &lt; Criterion2 for any j. /// Jacobian considered singular if |J(:,j)|2 &lt; Criterion2 for any j.
/// </summary> /// </summary>
public double Criterion2 = 1e-7; public double Criterion2 = 1e-7;
/// <summary> /// <summary>
/// Jacobian considered singular if |s|2 &lt; Criterion3. s is trial step. /// Jacobian considered singular if |s|2 &lt; Criterion3. s is trial step.
/// </summary> /// </summary>
public double Criterion3 = 1e-7; public double Criterion3 = 1e-7;
/// <summary> /// <summary>
/// |r|2 - |r - Js|2 &lt; Criterion4 /// |r|2 - |r - Js|2 &lt; Criterion4
/// </summary> /// </summary>
public double Criterion4 = 1e-7; public double Criterion4 = 1e-7;
public double TrialStepPrecision = 1e-10; public double TrialStepPrecision = 1e-10;
/// <summary>
/// Only if Jacobian calculated by central differences.
/// </summary>
public double JacobianPrecision = 1e-8;
}
/// <summary> /// <summary>
/// For details of convergence criteria, see Options. /// Only if Jacobian calculated by central differences.
/// </summary> /// </summary>
public enum ConvergenceType { NoneMaxIterationExceeded, Criterion0, Criterion1, Criterion2, Criterion3, Criterion4, Error }; public double JacobianPrecision = 1e-8;
}
public class Result
{ /// <summary>
public int NumberOfIterations; /// For details of convergence criteria, see Options.
/// </summary>
public enum NonLinearLeastSquaresConvergenceType { NoneMaxIterationExceeded, Criterion0, Criterion1, Criterion2, Criterion3, Criterion4, Error };
/// <summary>
/// Result of Non-Linear Least Squares Minimization.
/// </summary>
public class NonLinearLeastSquaresResult
{
public int NumberOfIterations;
public NonLinearLeastSquaresConvergenceType ConvergenceType;
}
/// <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
{
public NonLinearLeastSquaresResult Result { get; private set; }
public readonly NonLinearLeastSquaresOptions Options = new NonLinearLeastSquaresOptions();
public ConvergenceType ConvergenceType;
}
/// <summary> /// <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. /// 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. /// returning its best fitting parameters p.
@ -76,7 +87,7 @@ namespace MathNet.Numerics.Optimization
/// <param name="pStart">Initial guess of parameters, p.</param> /// <param name="pStart">Initial guess of parameters, p.</param>
/// <param name="jacobian">jac_j(x, p) = df / dp_j</param> /// <param name="jacobian">jac_j(x, p) = df / dp_j</param>
/// <returns></returns> /// <returns></returns>
public static double[] CurveFit(double[] x, double[] y, Func<double, double[], double> f, public double[] CurveFit(double[] x, double[] y, Func<double, double[], double> f,
double[] pStart, Func<double, double[], double[]> jacobian = null) double[] pStart, Func<double, double[], double[]> jacobian = null)
{ {
if (x.Length != y.Length) throw new ArgumentException("x and y lengths different"); if (x.Length != y.Length) throw new ArgumentException("x and y lengths different");
@ -101,7 +112,7 @@ namespace MathNet.Numerics.Optimization
}; };
double[] parameters; double[] parameters;
Result result = provider.NonLinearLeastSquaresUnboundedMinimize(y.Length, pStart, function, out parameters, jacobianFunction); Result = provider.NonLinearLeastSquaresUnboundedMinimize(y.Length, pStart, function, out parameters, jacobianFunction);
return parameters; return parameters;
} }
} }

4
src/Numerics/Providers/Optimization/IOptimizationProvider.cs
View File

@ -56,8 +56,8 @@ namespace MathNet.Numerics.Providers.Optimization
public interface IOptimizationProvider<T> public interface IOptimizationProvider<T>
where T : struct where T : struct
{ {
NonLinearLeastSquaresMinimizer.Result NonLinearLeastSquaresUnboundedMinimize( NonLinearLeastSquaresResult NonLinearLeastSquaresUnboundedMinimize(
int residualsLength, T[] initialGuess, LeastSquaresForwardModel function, int residualsLength, T[] initialGuess, LeastSquaresForwardModel function,
out T[] parameters, Jacobian jacobianFunction = null, NonLinearLeastSquaresMinimizer.Options options = null); out T[] parameters, Jacobian jacobianFunction = null, NonLinearLeastSquaresOptions options = null);
} }
} }

24
src/Numerics/Providers/Optimization/Mkl/MklOptimizationProvider.cs
View File

@ -42,9 +42,9 @@ namespace MathNet.Numerics.Providers.Optimization.Mkl
{ {
const int TR_SUCCESS = 1501; const int TR_SUCCESS = 1501;
public NonLinearLeastSquaresMinimizer.Result NonLinearLeastSquaresUnboundedMinimize(int residualsLength, double[] initialGuess, LeastSquaresForwardModel function, out double[] parameters, Jacobian jacobianFunction = null, NonLinearLeastSquaresMinimizer.Options options = null) public NonLinearLeastSquaresResult NonLinearLeastSquaresUnboundedMinimize(int residualsLength, double[] initialGuess, LeastSquaresForwardModel function, out double[] parameters, Jacobian jacobianFunction = null, NonLinearLeastSquaresOptions options = null)
{ {
if (options == null) options = new NonLinearLeastSquaresMinimizer.Options(); if (options == null) options = new NonLinearLeastSquaresOptions();
bool analyticJacobian = jacobianFunction != null; bool analyticJacobian = jacobianFunction != null;
double[] residuals = new double[residualsLength]; double[] residuals = new double[residualsLength];
double[] residualsMinus = new double[residualsLength]; double[] residualsMinus = new double[residualsLength];
@ -176,30 +176,30 @@ namespace MathNet.Numerics.Providers.Optimization.Mkl
SafeNativeMethods.FreeBuffers(); SafeNativeMethods.FreeBuffers();
NonLinearLeastSquaresMinimizer.ConvergenceType convergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.Error; NonLinearLeastSquaresConvergenceType convergenceType = NonLinearLeastSquaresConvergenceType.Error;
switch (rciRequest) switch (rciRequest)
{ {
case -1: case -1:
convergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.NoneMaxIterationExceeded; break; convergenceType = NonLinearLeastSquaresConvergenceType.NoneMaxIterationExceeded; break;
case -2: case -2:
convergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.Criterion0; break; convergenceType = NonLinearLeastSquaresConvergenceType.Criterion0; break;
case -3: case -3:
convergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.Criterion1; break; convergenceType = NonLinearLeastSquaresConvergenceType.Criterion1; break;
case -4: case -4:
convergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.Criterion2; break; convergenceType = NonLinearLeastSquaresConvergenceType.Criterion2; break;
case -5: case -5:
convergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.Criterion3; break; convergenceType = NonLinearLeastSquaresConvergenceType.Criterion3; break;
case -6: case -6:
convergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.Criterion4; break; convergenceType = NonLinearLeastSquaresConvergenceType.Criterion4; break;
} }
// no errors, find reason for stopping; // no errors, find reason for stopping;
return new NonLinearLeastSquaresMinimizer.Result() { ConvergenceType = convergenceType, NumberOfIterations = iterations }; return new NonLinearLeastSquaresResult() { ConvergenceType = convergenceType, NumberOfIterations = iterations };
} }
public static NonLinearLeastSquaresMinimizer.Result ErrorResult() public static NonLinearLeastSquaresResult ErrorResult()
{ {
return new NonLinearLeastSquaresMinimizer.Result() { ConvergenceType = NonLinearLeastSquaresMinimizer.ConvergenceType.Error }; return new NonLinearLeastSquaresResult() { ConvergenceType = NonLinearLeastSquaresConvergenceType.Error };
} }
} }
} }

7
src/UnitTests/OptimizationTests/NonLinearLeastSquaresTest.cs
View File

@ -40,17 +40,20 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests
[Test] [Test]
public void CurveFit() public void CurveFit()
{ {
var minimizer = new NonLinearLeastSquaresMinimizer();
minimizer.Options.MaximumIterations = 1010;
minimizer.Options.Criterion4 = 1e-8;
// y = b1*(1-exp[-b2*x]) + e // y = b1*(1-exp[-b2*x]) + e
var xin = new double[] { 1, 2, 3, 5, 7, 10 }; var xin = new double[] { 1, 2, 3, 5, 7, 10 };
var yin = new double[] { 109, 149, 149, 191, 213, 224 }; var yin = new double[] { 109, 149, 149, 191, 213, 224 };
var popt = NonLinearLeastSquaresMinimizer.CurveFit(xin, yin, (x, p) => p[0] * (1 - Math.Exp(-p[1] * x)), new double[] { 1, 1 }); 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> function = (x, p) => p[0] * (1 - Math.Exp(-p[1] * x));
Func<double, double[], double[]> jacobian = (x, p) => new double[] { Func<double, double[], double[]> jacobian = (x, p) => new double[] {
1 - Math.Exp(-p[1] * x), 1 - Math.Exp(-p[1] * x),
p[0] * x * Math.Exp(-p[1] * x) }; p[0] * x * Math.Exp(-p[1] * x) };
popt = NonLinearLeastSquaresMinimizer.CurveFit(xin, yin, function, new double[] { 1, 1 }, jacobian); // 100, 0.75 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[] expected = new double[] { 2.1380940889E+02, 5.4723748542E-01 };

Write Preview
Loading…
Cancel
Save