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

Optimization: added Levenberg-Marquardt and trust region dogleg minimizer.

pull/614/head
diluculo 8 years ago
parent
7e45066368
commit
aca34e13b9
8 changed files with 1885 additions and 1 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 404
      src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs
  2. 6
      src/Numerics/Optimization/ExitCondition.cs
  3. 70
      src/Numerics/Optimization/IObjectiveModel.cs
  4. 263
      src/Numerics/Optimization/LevenbergMarquardtMinimizer.cs
  5. 52
      src/Numerics/Optimization/ModelMinimizationResult.cs
  6. 40
      src/Numerics/Optimization/ObjectiveModel.cs
  7. 729
      src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs
  8. 322
      src/Numerics/Optimization/TrustRegionDogLegMinimizer.cs

404
src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs
View File

@ -0,0 +1,404 @@
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearAlgebra.Double;
using MathNet.Numerics.Optimization;
using NUnit.Framework;
using System;
namespace MathNet.Numerics.UnitTests.OptimizationTests
{
[TestFixture]
public class NonLinearCurveFittingTests
{
// model: Rosenbrock
// f(x; a, b) = (1 - a)^2 + 100*(b - a^2)^2
// derivatives:
// df/da = 400*a^3 - 400*a*b + 2*a - 2
// df/db = 200*(b - a^2)
// best fitted parameters:
// a = 1
// b = 1
private double RosenbrockModel(Vector<double> p, double x)
{
var y = Math.Pow(1.0 - p[0], 2) + 100.0 * Math.Pow(p[1] - p[0] * p[0], 2);
return y;
}
private Vector<double> RosenbrockPrime(Vector<double> p, double x)
{
var prime = Vector<double>.Build.Dense(p.Count);
prime[0] = 400.0 * p[0] * p[0] * p[0] - 400.0 * p[0] * p[1] + 2.0 * p[0] - 2.0;
prime[1] = 200.0 * (p[1] - p[0] * p[0]);
return prime;
}
private Vector<double> RosenbrockX = Vector<double>.Build.Dense(2);
private Vector<double> RosenbrockY = Vector<double>.Build.Dense(2);
private Vector<double> RosenbrockPbest = new DenseVector(new double[] { 1.0, 1.0 });
private Vector<double> RosenbrockStart1 = new DenseVector(new double[] { -1.2, 1.0 });
private Vector<double> RosebbrockLowerBound = new DenseVector(new double[] { -5.0, -5.0 });
private Vector<double> RosenbrockUpperBound = new DenseVector(new double[] { 5.0, 5.0 });
[Test]
public void LMDER_FindMinimum_Rosenbrock_Unconstrained()
{
var obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockPrime, RosenbrockX, RosenbrockY);
var solver = new LevenbergMarquardtMinimizer(maximumIterations: 10000);
var result = solver.FindMinimum(obj, RosenbrockStart1);
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[0], result.BestFitParameters[0], 3);
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[1], result.BestFitParameters[1], 3);
}
[Test]
public void LMDIF_FindMinimum_Rosenbrock_Unconstrained()
{
var obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockX, RosenbrockY, accuracyOrder:2);
var solver = new LevenbergMarquardtMinimizer(maximumIterations: 10000);
var result = solver.FindMinimum(obj, RosenbrockStart1);
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[0], result.BestFitParameters[0], 3);
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[1], result.BestFitParameters[1], 3);
}
[Test]
public void LMDER_FindMinimum_Rosenbrock_BoxConstrained()
{
var obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockPrime, RosenbrockX, RosenbrockY,
lowerBound : RosebbrockLowerBound, upperBound : RosenbrockUpperBound);
var solver = new LevenbergMarquardtMinimizer(maximumIterations: 10000);
var result = solver.FindMinimum(obj, RosenbrockStart1);
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[0], result.BestFitParameters[0], 3);
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[1], result.BestFitParameters[1], 3);
}
[Test]
public void LMDIF_FindMinimum_Rosenbrock_BoxConstrained()
{
var obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockX, RosenbrockY,
lowerBound: RosebbrockLowerBound, upperBound: RosenbrockUpperBound,
accuracyOrder: 2);
var solver = new LevenbergMarquardtMinimizer(maximumIterations: 10000);
var result = solver.FindMinimum(obj, RosenbrockStart1);
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[0], result.BestFitParameters[0], 3);
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[1], result.BestFitParameters[1], 3);
}
[Test]
public void TRLMDER_FindMinimum_Rosenbrock_Unconstrained()
{
var obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockPrime, RosenbrockX, RosenbrockY);
var solver = new TrustRegionDogLegMinimizer(maximumIterations: 10000);
var result = solver.FindMinimum(obj, RosenbrockStart1);
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[0], result.BestFitParameters[0], 1);
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[1], result.BestFitParameters[1], 1);
}
[Test]
public void TRLMDIF_FindMinimum_Rosenbrock_Unconstrained()
{
var obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockPrime, RosenbrockX, RosenbrockY);
var solver = new TrustRegionDogLegMinimizer(maximumIterations: 10000);
var result = solver.FindMinimum(obj, RosenbrockStart1);
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[0], result.BestFitParameters[0], 1);
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[1], result.BestFitParameters[1], 1);
}
// model: Rat43 (https://www.itl.nist.gov/div898/strd/nls/data/ratkowsky3.shtml)
// f(x; a, b, c, d) = a / ((1 + exp(b - c * x))^(1 / d))
// best fitted parameters:
// a = 6.9964151270E+02 +/- 1.6302297817E+01
// b = 5.2771253025E+00 +/- 2.0828735829E+00
// c = 7.5962938329E-01 +/- 1.9566123451E-01
// d = 1.2792483859E+00 +/- 6.8761936385E-01
private double Rat43Model(Vector<double> p, double x)
{
var y = p[0] / Math.Pow(1.0 + Math.Exp(p[1] - p[2] * x), 1.0 / p[3]);
return y;
}
private Vector<double> Rat43X = new DenseVector(new double[] {
1.00, 2.00, 3.00, 4.00, 5.00, 6.00, 7.00, 8.00, 9.00, 10.00,
11.00, 12.00, 13.00, 14.00, 15.00
});
private Vector<double> Rat43Y = new DenseVector(new double[] {
16.08, 33.83, 65.80, 97.20, 191.55, 326.20, 386.87, 520.53, 590.03, 651.92,
724.93, 699.56, 689.96, 637.56, 717.41
});
private Vector<double> Rat43Pbest = new DenseVector(new double[] {
6.9964151270E+02, 5.2771253025E+00, 7.5962938329E-01, 1.2792483859E+00
});
private Vector<double> Rat43Pstd = new DenseVector(new double[]{
1.6302297817E+01, 2.0828735829E+00, 1.9566123451E-01, 6.8761936385E-01
});
private Vector<double> Rat43Start1 = new DenseVector(new double[] { 100, 10, 1, 1 });
private Vector<double> Rat43Start2 = new DenseVector(new double[] { 700, 5, 0.75, 1.3 });
[Test]
public void LMDIF_FindMinimum_Rat43_Unconstrained()
{
var obj = ObjectiveModel.FittingModel(Rat43Model, Rat43X, Rat43Y, accuracyOrder: 6);
var solver = new LevenbergMarquardtMinimizer();
var result = solver.FindMinimum(obj, Rat43Start1);
for (int i = 0; i < result.BestFitParameters.Count; i++)
{
AssertHelpers.AlmostEqualRelative(Rat43Pbest[i], result.BestFitParameters[i], 6);
AssertHelpers.AlmostEqualRelative(Rat43Pstd[i], result.StandardErrors[i], 6);
}
}
[Test]
public void TRLMDIF_FindMinimum_Rat43_Unconstrained()
{
var obj = ObjectiveModel.FittingModel(Rat43Model, Rat43X, Rat43Y, accuracyOrder: 6);
var solver = new TrustRegionDogLegMinimizer();
var result = solver.FindMinimum(obj, Rat43Start2);
for (int i = 0; i < result.BestFitParameters.Count; i++)
{
AssertHelpers.AlmostEqualRelative(Rat43Pbest[i], result.BestFitParameters[i], 2);
AssertHelpers.AlmostEqualRelative(Rat43Pstd[i], result.StandardErrors[i], 2);
}
}
// model: BoxBod (https://www.itl.nist.gov/div898/strd/nls/data/boxbod.shtml)
// f(x; a, b) = a*(1 - exp(-b*x))
// derivatives:
// df/da = 1 - exp(-b*x)
// df/db = a*x*exp(-b*x)
// best fitted parameters:
// a = 2.1380940889E+02 +/- 1.2354515176E+01
// b = 5.4723748542E-01 +/- 1.0455993237E-01
private double BoxBodModel(Vector<double> p, double x)
{
var y = p[0] * (1.0 - Math.Exp(-p[1] * x));
return y;
}
private Vector<double> BoxBodPrime(Vector<double> p, double x)
{
var prime = Vector<double>.Build.Dense(p.Count);
prime[0] = 1.0 - Math.Exp(-p[1] * x);
prime[1] = p[0] * x * Math.Exp(-p[1] * x);
return prime;
}
private Vector<double> BoxBodX = new DenseVector(new double[] { 1, 2, 3, 5, 7, 10 });
private Vector<double> BoxBodY = new DenseVector(new double[] { 109, 149, 149, 191, 213, 224 });
private Vector<double> BoxBodPbest = new DenseVector(new double[] { 2.1380940889E+02, 5.4723748542E-01 });
private Vector<double> BoxBodPstd = new DenseVector(new double[] { 1.2354515176E+01, 1.0455993237E-01 });
private Vector<double> BoxBodStart1 = new DenseVector(new double[] { 1.0, 1.0 });
private Vector<double> BoxBodStart2 = new DenseVector(new double[] { 100.0, 0.75 });
private Vector<double> BoxBodLowerBound = new DenseVector(new double[] { 0, 0 });
private Vector<double> BoxBodUpperBound = new DenseVector(new double[] { 500.0, 10 });
private Vector<double> BoxBodScales = new DenseVector(new double[] { 100.0, 1 });
[Test]
public void LMDER_FindMinimum_BoxBod_Unconstrained()
{
var obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY);
var solver = new LevenbergMarquardtMinimizer();
var result = solver.FindMinimum(obj, BoxBodStart1);
AssertHelpers.AlmostEqualRelative(BoxBodPbest[0], result.BestFitParameters[0], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPbest[1], result.BestFitParameters[1], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[0], result.StandardErrors[0], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[1], result.StandardErrors[1], 6);
}
[Test]
public void LMDIF_FindMinimum_BoxBod_Unconstrained()
{
var obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodX, BoxBodY, accuracyOrder:6);
var solver = new LevenbergMarquardtMinimizer();
var result = solver.FindMinimum(obj, BoxBodStart1);
AssertHelpers.AlmostEqualRelative(BoxBodPbest[0], result.BestFitParameters[0], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPbest[1], result.BestFitParameters[1], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[0], result.StandardErrors[0], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[1], result.StandardErrors[1], 6);
}
[Test]
public void LMDER_FindMinimum_BoxBod_BoxConstrained()
{
var obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY,
lowerBound: BoxBodLowerBound, upperBound: BoxBodUpperBound);
var solver = new LevenbergMarquardtMinimizer();
var result = solver.FindMinimum(obj, BoxBodStart1);
AssertHelpers.AlmostEqualRelative(BoxBodPbest[0], result.BestFitParameters[0], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPbest[1], result.BestFitParameters[1], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[0], result.StandardErrors[0], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[1], result.StandardErrors[1], 6);
}
[Test]
public void LMDIF_FindMinimum_BoxBod_BoxConstrained()
{
var obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodX, BoxBodY,
lowerBound: BoxBodLowerBound, upperBound: BoxBodUpperBound,
accuracyOrder: 6);
var solver = new LevenbergMarquardtMinimizer();
var result = solver.FindMinimum(obj, BoxBodStart1);
AssertHelpers.AlmostEqualRelative(BoxBodPbest[0], result.BestFitParameters[0], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPbest[1], result.BestFitParameters[1], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[0], result.StandardErrors[0], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[1], result.StandardErrors[1], 6);
}
[Test]
public void TRLMDIF_FindMinimum_BoxBod_Unconstrained()
{
var obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodX, BoxBodY, accuracyOrder: 6);
var solver = new TrustRegionDogLegMinimizer();
var result = solver.FindMinimum(obj, BoxBodStart1);
AssertHelpers.AlmostEqualRelative(BoxBodPbest[0], result.BestFitParameters[0], 3);
AssertHelpers.AlmostEqualRelative(BoxBodPbest[1], result.BestFitParameters[1], 3);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[0], result.StandardErrors[0], 3);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[1], result.StandardErrors[1], 3);
}
// model : Thurber (https://www.itl.nist.gov/div898/strd/nls/data/thurber.shtml)
// f(x; b1 ... b7) = (b1 + b2*x + b3*x^2 + b4*x^3) / (1 + b5*x + b6*x^2 + b7*x^3)
// derivatives:
// df/db1 = 1/(b5*x + b6*x^2 + b7*x^3 + 1)
// df/db2 = x/(b5*x + b6*x^2 + b7*x^3 + 1)
// df/db3 = x^2/(b5*x + b6*x^2 + b7*x^3 + 1)
// df/db4 = x^3/(b5*x + b6*x^2 + b7*x^3 + 1)
// df/db5 = -(x*(b1 + x*(b2 + x*(b3 + b4*x))))/(b5*x + b6*x^2 + b7*x^3 + 1)^2
// df/db6 = -(x^2*(b1 + x*(b2 + x*(b3 + b4*x))))/(b5*x + b6*x^2 + b7*x^3 + 1)^2
// df/db7 = -(x^3*(b1 + x*(b2 + x*(b3 + b4*x))))/(b5*x + b6*x^2 + b7*x^3 + 1)^2
// best fitted parameters:
// b1 = 1.2881396800E+03 +/- 4.6647963344E+00
// b2 = 1.4910792535E+03 +/- 3.9571156086E+01
// b3 = 5.8323836877E+02 +/- 2.8698696102E+01
// b4 = 7.5416644291E+01 +/- 5.5675370270E+00
// b5 = 9.6629502864E-01 +/- 3.1333340687E-02
// b6 = 3.9797285797E-01 +/- 1.4984928198E-02
// b7 = 4.9727297349E-02 +/- 6.5842344623E-03
private double ThurberModel(Vector<double> p, double x)
{
var xSq = x * x;
var xCb = xSq * x;
var y = (p[0] + p[1] * x + p[2] * xSq + p[3] * xCb)
/ (1 + p[4] * x + p[5] * xSq + p[6] * xCb);
return y;
}
private Vector<double> ThurberPrime(Vector<double> p, double x)
{
var prime = Vector<double>.Build.Dense(p.Count);
var xSq = x * x;
var xCb = xSq * x;
var num = (p[0] + x * (p[1] + x * (p[2] + p[3] * x)));
var den = (p[4] * x + p[5] * xSq + p[6] * xCb + 1.0);
var denSq = den * den;
prime[0] = 1.0 / den;
prime[1] = x / den;
prime[2] = xSq / den;
prime[3] = xCb / den;
prime[4] = -(x * num) / denSq;
prime[5] = -(xSq * num) / denSq;
prime[6] = -(xCb * num) / denSq;
return prime;
}
private Vector<double> ThurberX = new DenseVector(new double[] {
-3.067, -2.981, -2.921, -2.912, -2.84,
-2.797, -2.702, -2.699, -2.633, -2.481,
-2.363, -2.322, -1.501, -1.460, -1.274,
-1.212, -1.100, -1.046, -0.915, -0.714,
-0.566, -0.545, -0.400, -0.309, -0.109,
-0.103, 0.01, 0.119, 0.377, 0.79,
0.963, 1.006, 1.115, 1.572, 1.841,
2.047, 2.2});
private Vector<double> ThurberY = new DenseVector(new double[] {
80.574, 084.248, 087.264, 087.195, 089.076,
089.608, 089.868, 090.101, 092.405, 095.854,
100.696, 101.060, 401.672, 390.724, 567.534,
635.316, 733.054, 759.087, 894.206, 990.785,
1090.109, 1080.914, 1122.643, 1178.351, 1260.531,
1273.514, 1288.339, 1327.543, 1353.863, 1414.509,
1425.208, 1421.384, 1442.962, 1464.350, 1468.705,
1447.894, 1457.628});
private Vector<double> ThurberPbest = new DenseVector(new double[] {
1.2881396800E+03, 1.4910792535E+03, 5.8323836877E+02, 7.5416644291E+01, 9.6629502864E-01,
3.9797285797E-01, 4.9727297349E-02 });
private Vector<double> ThurberPstd = new DenseVector(new double[] {
4.6647963344E+00, 3.9571156086E+01, 2.8698696102E+01, 5.5675370270E+00, 3.1333340687E-02,
1.4984928198E-02, 6.5842344623E-03 });
private Vector<double> ThurberInitialGuess = new DenseVector(new double[] { 1000.0, 1000.0, 400.0, 40.0, 0.7, 0.3, 0.03 });
private Vector<double> ThurberScales = new DenseVector(new double[7] { 1000, 1000, 400, 40, 0.7, 0.3, 0.03 });
[Test]
public void LMDER_FindMinimum_Thurber_Unconstrained()
{
var obj = ObjectiveModel.FittingModel(ThurberModel, ThurberPrime, ThurberX, ThurberY);
var solver = new LevenbergMarquardtMinimizer();
var result = solver.FindMinimum(obj, ThurberInitialGuess);
for (int i = 0; i < result.BestFitParameters.Count; i++)
{
AssertHelpers.AlmostEqualRelative(ThurberPbest[i], result.BestFitParameters[i], 6);
AssertHelpers.AlmostEqualRelative(ThurberPstd[i], result.StandardErrors[i], 6);
}
}
[Test]
public void LMDIF_FindMinimum_Thurber_Unconstrained()
{
var obj = ObjectiveModel.FittingModel(ThurberModel, ThurberX, ThurberY, accuracyOrder: 6);
var solver = new LevenbergMarquardtMinimizer();
var result = solver.FindMinimum(obj, ThurberInitialGuess);
for (int i = 0; i < result.BestFitParameters.Count; i++)
{
AssertHelpers.AlmostEqualRelative(ThurberPbest[i], result.BestFitParameters[i], 6);
AssertHelpers.AlmostEqualRelative(ThurberPstd[i], result.StandardErrors[i], 6);
}
}
[Test]
public void TRLMDIF_FindMinimum_Thurber_Scaled()
{
var obj = ObjectiveModel.FittingModel(ThurberModel, ThurberX, ThurberY,
scales: ThurberScales,
accuracyOrder: 6);
var solver = new TrustRegionDogLegMinimizer();
var result = solver.FindMinimum(obj, ThurberInitialGuess);
for (int i = 0; i < result.BestFitParameters.Count; i++)
{
AssertHelpers.AlmostEqualRelative(ThurberPbest[i], result.BestFitParameters[i], 3);
AssertHelpers.AlmostEqualRelative(ThurberPstd[i], result.StandardErrors[i], 3);
}
}
}
}

6
src/Numerics/Optimization/ExitCondition.cs
View File

@ -32,12 +32,16 @@ namespace MathNet.Numerics.Optimization
public enum ExitCondition
{
None,
InvalidValues,
ExceedIterations,
RelativePoints,
RelativeGradient,
LackOfProgress,
AbsoluteGradient,
WeakWolfeCriteria,
BoundTolerance,
StrongWolfeCriteria,
Converged
Converged,
ManuallyStopped
}
}

70
src/Numerics/Optimization/IObjectiveModel.cs
View File

@ -0,0 +1,70 @@
using MathNet.Numerics.LinearAlgebra;
namespace MathNet.Numerics.Optimization
{
public interface IObjectiveModelEvaluation
{
IObjectiveModel CreateNew();
/// <summary>
/// Get the y-values of the fitted model that correspond to the independent values.
/// </summary>
Vector<double> Values { get; }
/// <summary>
/// Get the values of the parameters.
/// </summary>
Vector<double> Parameters { get; }
/// <summary>
/// Get the residual sum of squares.
/// </summary>
double Residue { get; }
/// <summary>
/// Get the Jacobian matrix, J(x; p) = df(x; p)/dp.
/// </summary>
Matrix<double> Jacobian { get; }
/// <summary>
/// Get the Gradient vector. G = J'(y - f(x; p))
/// </summary>
Vector<double> Gradient { get; }
/// <summary>
/// Get the approximated Hessian matrix. H = J'J
/// </summary>
Matrix<double> Hessian { get; }
/// <summary>
/// Get the covariance matrix.
/// </summary>
Matrix<double> Covariance { get; }
/// <summary>
/// Get the number of calls to function.
/// </summary>
int FunctionEvaluations { get; }
/// <summary>
/// Get the number of calls to jacobian.
/// </summary>
int JacobianEvaluations { get; }
/// <summary>
/// Get the degree of freedom.
/// </summary>
int DegreeOfFreedom { get; }
/// <summary>
/// Get whether or not the analytical jacobian is supported.
/// </summary>
bool IsJacobianSupported { get; }
}
public interface IObjectiveModel : IObjectiveModelEvaluation
{
void EvaluateFunction(Vector<double> parameters);
void EvaluateJacobian(Vector<double> parameters);
void EvaluateCovariance(Vector<double> parameters);
/// <summary>Create a new independent copy of this objective function, evaluated at the same point.</summary>
IObjectiveModel Fork();
}
}

263
src/Numerics/Optimization/LevenbergMarquardtMinimizer.cs
View File

@ -0,0 +1,263 @@
using MathNet.Numerics.LinearAlgebra;
using System;
using System.Linq;
namespace MathNet.Numerics.Optimization
{
public sealed class LevenbergMarquardtMinimizer
{
#region Tolerances and options
/// <summary>
/// The scale factor for initial mu
/// </summary>
public static double InitialMu { get; set; }
/// <summary>
/// The stopping threshold for infinity norm of the gradient.
/// </summary>
public static double GradientTolerance { get; set; }
/// <summary>
/// The stopping threshold for L2 norm of the change of the parameters.
/// </summary>
public static double StepTolerance { get; set; }
/// <summary>
/// The stopping threshold for the function value or L2 norm of the residuals.
/// </summary>
public static double FunctionTolerance { get; set; }
/// <summary>
/// The maximum number of iterations.
/// </summary>
public int MaximumIterations { get; set; }
#endregion Tolerances and options
public LevenbergMarquardtMinimizer(double initialMu = 1E-3, double gradientTolerance = 1E-18, double stepTolerance = 1E-18, double functionTolerance = 1E-18, int maximumIterations = -1)
{
InitialMu = initialMu;
GradientTolerance = gradientTolerance;
StepTolerance = stepTolerance;
FunctionTolerance = functionTolerance;
MaximumIterations = maximumIterations;
}
public ModelMinimizationResult FindMinimum(IObjectiveModel objective, Vector<double> initialGuess)
{
if (objective == null)
throw new ArgumentNullException("objective");
if (initialGuess == null)
throw new ArgumentNullException("initialGuess");
return Minimum(objective, initialGuess, InitialMu, FunctionTolerance, GradientTolerance, StepTolerance, MaximumIterations);
}
public ModelMinimizationResult FindMinimum(IObjectiveModel objective, double[] initialGuess)
{
if (objective == null)
throw new ArgumentNullException("objective");
if (initialGuess == null)
throw new ArgumentNullException("initialGuess");
return Minimum(objective, CreateVector.DenseOfArray<double>(initialGuess), InitialMu, GradientTolerance, StepTolerance, FunctionTolerance, MaximumIterations);
}
/// <summary>
/// Non-linear least square fitting by the Levenberg-Marduardt algorithm.
/// </summary>
/// <param name="objective">The objective function, including model, observations, and parameter bounds.</param>
/// <param name="initialGuess">The initial guess values.</param>
/// <param name="initialMu">The initial damping parameter of mu.</param>
/// <param name="gradientTolerance">The stopping threshold for infinity norm of the gradient vector.</param>
/// <param name="stepTolerance">The stopping threshold for L2 norm of the change of parameters.</param>
/// <param name="functionTolerance">The stopping threshold for L2 norm of the residuals.</param>
/// <param name="maximumIterations">The max iterations.</param>
/// <returns>The result of the Levenberg-Marquardt minimization</returns>
public static ModelMinimizationResult Minimum(IObjectiveModel objective, Vector<double> initialGuess, double initialMu = 1E-3, double gradientTolerance = 1E-18, double stepTolerance = 1E-18, double functionTolerance = 1E-18, int maximumIterations = -1)
{
// Non-linear least square fitting by the Levenberg-Marduardt algorithm.
//
// Levenberg-Marquardt is finding the minimum of a function F(p) that is a sum of squares of nonlinear functions.
//
// For given datum pair (x, y), uncertainties σ (or weighting W = 1 / σ^2) and model function f = f(x; p),
// let's find the parameters of the model so that the sum of the quares of the deviations is minimized.
//
// F(p) = 1/2 * ∑{ Wi * (yi - f(xi; p))^2 }
// pbest = argmin F(p)
//
// We will use the following terms:
// Weighting W is the diagonal matrix and can be decomposed as LL', so L = 1/σ
// Residuals, R = L(y - f(x; p))
// Residual sum of squares, RSS = ||R||^2 = R.DotProduct(R)
// Jacobian J = df(x; p)/dp
// Gradient g = J'W(y − f(x; p)) = J'LR
// Approximated Hessian H = J'WJ
//
// The Levenberg-Marquardt algorithm is summarized as follows:
// initially let μ = τ * max(diag(J'WJ)).
// repeat
// solve linear equations: (J'WJ + μI)ΔP = J'R
// let ρ = (||R||^2 - ||Rnew||^2) / (Δp'(μΔp + J'R)).
// if ρ > ε, P = P + ΔP; μ = μ * max(1/3, 1 - (2ρ - 1)^3); ν = 2;
// otherwise μ = μ*ν; ν = 2*ν;
//
// References:
// [1]. Madsen, K., H. B. Nielsen, and O. Tingleff.
// "Methods for Non-Linear Least Squares Problems. Technical University of Denmark, 2004. Lecture notes." (2004).
// Available Online from: http://orbit.dtu.dk/files/2721358/imm3215.pdf
// [2]. Gavin, Henri.
// "The Levenberg-Marquardt method for nonlinear least squares curve-fitting problems."
// Department of Civil and Environmental Engineering, Duke University (2017): 1-19.
// Availble Online from: http://people.duke.edu/~hpgavin/ce281/lm.pdf
if (objective == null)
throw new ArgumentNullException("objective");
if (initialGuess == null)
throw new ArgumentNullException("initialGuess");
ExitCondition exitCondition = ExitCondition.None;
// First, calculate function values and setup variables
objective.EvaluateFunction(initialGuess);
var P = objective.Parameters; // current parameters
var Pstep = Vector<double>.Build.Dense(P.Count); // the change of parameters
var RSS = objective.Residue; // Residual Sum of Squares = R'R
if (maximumIterations < 0)
{
maximumIterations = (objective.IsJacobianSupported)
? 100 * (initialGuess.Count + 1)
: 200 * (initialGuess.Count + 1);
}
// if RSS == NaN, stop
if (double.IsNaN(RSS))
{
exitCondition = ExitCondition.InvalidValues;
return new ModelMinimizationResult(objective, -1, exitCondition);
}
// When only function evaluation is needed, set maximumIterations to zero,
if (maximumIterations == 0)
{
exitCondition = ExitCondition.ManuallyStopped;
}
// if RSS <= fTol, stop
if (RSS <= functionTolerance)
{
exitCondition = ExitCondition.Converged; // SmallRSS
}
// Evaluate gradient and Hessian
objective.EvaluateJacobian(P);
var Gradient = objective.Gradient;
var Hessian = objective.Hessian;
var diagonalOfHessian = Hessian.Diagonal(); // diag(H)
// if ||g||oo <= gtol, found and stop
if (Gradient.InfinityNorm() <= gradientTolerance)
{
exitCondition = ExitCondition.RelativeGradient;
}
if (exitCondition != ExitCondition.None)
{
objective.EvaluateCovariance(P);
return new ModelMinimizationResult(objective, -1, exitCondition);
}
double mu = initialMu * diagonalOfHessian.Max(); // μ
double nu = 2; // ν
int iterations = 0;
while (iterations < maximumIterations && exitCondition == ExitCondition.None)
{
iterations++;
while (true)
{
Hessian.SetDiagonal(Hessian.Diagonal() + mu); // hessian[i, i] = hessian[i, i] + mu;
// solve normal equations
Pstep = Hessian.Solve(Gradient);
// if ||ΔP|| <= xTol * (||P|| + xTol), found and stop
if (Pstep.L2Norm() <= stepTolerance * (stepTolerance + P.DotProduct(P)))
{
exitCondition = ExitCondition.RelativePoints; // SmallRelativeParameters
break;
}
var Pnew = P + Pstep; // new parameters to test
objective.EvaluateFunction(Pnew);
var RSSnew = objective.Residue;
if (double.IsNaN(RSSnew))
{
exitCondition = ExitCondition.InvalidValues;
break;
}
// calculate the ratio of the actual to the predicted reduction.
// ρ = (RSS - RSSnew) / (Δp'(μΔp + g))
var predictedReduction = Pstep.DotProduct(mu * Pstep + Gradient);
var rho = (predictedReduction != 0)
? (RSS - RSSnew) / predictedReduction
: 0;
if (rho > 0.0)
{
// accepted
Pnew.CopyTo(P);
RSS = RSSnew;
// update gradient and Hessian
objective.EvaluateJacobian(P);
Gradient = objective.Gradient;
Hessian = objective.Hessian;
diagonalOfHessian = Hessian.Diagonal();
// if ||g||_oo <= gtol, found and stop
if (Gradient.InfinityNorm() <= gradientTolerance)
{
exitCondition = ExitCondition.RelativeGradient;
}
// if ||R||^2 < fTol, found and stop
if (RSS <= functionTolerance)
{
exitCondition = ExitCondition.Converged; // SmallRSS
}
mu = mu * Math.Max(1.0 / 3.0, 1.0 - Math.Pow(2.0 * rho - 1.0, 3));
nu = 2;
break;
}
else
{
// rejected, increased μ
mu = mu * nu;
nu = 2 * nu;
Hessian.SetDiagonal(diagonalOfHessian);
}
}
}
if (iterations >= maximumIterations)
{
exitCondition = ExitCondition.ExceedIterations;
}
// finalize
objective.EvaluateCovariance(P);
return new ModelMinimizationResult(objective, iterations, exitCondition);
}
}
}

52
src/Numerics/Optimization/ModelMinimizationResult.cs
View File

@ -0,0 +1,52 @@
using MathNet.Numerics.LinearAlgebra;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace MathNet.Numerics.Optimization
{
public class ModelMinimizationResult
{
public IObjectiveModel ModelInfoAtMinimum { get; private set; }
/// <summary>
/// Returns the best fit parameters.
/// </summary>
public Vector<double> BestFitParameters { get { return ModelInfoAtMinimum.Parameters; } }
/// <summary>
/// Returns the standard errors of the corresponding parameters
/// </summary>
public Vector<double> StandardErrors
{
get
{
if (ModelInfoAtMinimum.Covariance == null)
return null;
return ModelInfoAtMinimum.Covariance.Diagonal().PointwiseSqrt();
}
}
/// <summary>
/// Returns the y-values of the fitted model that correspond to the independent values.
/// </summary>
public Vector<double> BestFitValues { get { return ModelInfoAtMinimum.Values; } }
/// <summary>
/// Returns the residual sum of squares.
/// </summary>
public double Residue { get { return ModelInfoAtMinimum.Residue; } }
public double DegreeOfFreedom { get { return ModelInfoAtMinimum.DegreeOfFreedom; } }
public int Iterations { get; private set; }
public ExitCondition ReasonForExit { get; private set; }
public ModelMinimizationResult(IObjectiveModel modelInfo, int iterations, ExitCondition reasonForExit)
{
ModelInfoAtMinimum = modelInfo;
Iterations = iterations;
ReasonForExit = reasonForExit;
}
}
}

40
src/Numerics/Optimization/ObjectiveModel.cs
View File

@ -0,0 +1,40 @@
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.Optimization.ObjectiveModels;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace MathNet.Numerics.Optimization
{
public static class ObjectiveModel
{
/// <summary>
/// Fitting model with a user supplied jacobian for non-linear least squares regression.
/// </summary>
public static IObjectiveModel FittingModel(Func<Vector<double>, double, double> function, Func<Vector<double>, double, Vector<double>> derivatives,
Vector<double> observedX, Vector<double> observedY, Vector<double> weight = null,
Vector<double> lowerBound = null, Vector<double> upperBound = null, Vector<double> scales = null, List<bool> isFixed = null)
{
var objective = new FittingObjectiveModel(function, derivatives);
objective.SetObserved(observedX, observedY, weight);
objective.SetParameters(lowerBound, upperBound, scales, isFixed);
return objective;
}
/// <summary>
/// Fitting model for non-linear least squares regression.
/// The numerical jacobian with accuracy order is used.
/// </summary>
public static IObjectiveModel FittingModel(Func<Vector<double>, double, double> function,
Vector<double> observedX, Vector<double> observedY, Vector<double> weight = null,
Vector<double> lowerBound = null, Vector<double> upperBound = null, Vector<double> scales = null, List<bool> isFixed = null,
int accuracyOrder = 2)
{
var objective = new FittingObjectiveModel(function, null, accuracyOrder: accuracyOrder);
objective.SetObserved(observedX, observedY, weight);
objective.SetParameters(lowerBound, upperBound, scales, isFixed);
return objective;
}
}
}

729
src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs
View File

@ -0,0 +1,729 @@
using MathNet.Numerics.LinearAlgebra;
using System;
using System.Collections.Generic;
using System.Linq;
namespace MathNet.Numerics.Optimization.ObjectiveModels
{
internal class FittingObjectiveModel : IObjectiveModel
{
readonly Func<Vector<double>, double, double> userFunction; // (p, x) => f(x; p)
readonly Func<Vector<double>, double, Vector<double>> userDerivatives; // (p, x) => df(x; p)/dp
#region Public Variables
/// <summary>
/// Set or get the values of the independent variable.
/// </summary>
public Vector<double> ObservedX { get; private set; }
/// <summary>
/// Set or get the values of the observations.
/// </summary>
public Vector<double> ObservedY { get; private set; }
/// <summary>
/// Set or get the values of the weights for the observations.
/// inverse of the standard measurement errors
/// If null, unity weighting is used.
/// </summary>
public Matrix<double> Weights { get; private set; }
// W = LL'
private Vector<double> L;
/// <summary>
/// Set or get the values of the parameters.
/// </summary>
public Vector<double> Parameters { get; private set; }
/// <summary>
/// Set or get the values of the parameters.
/// </summary>
public List<bool> IsFixed { get; set; }
/// <summary>
/// Set or get the values of the parameters.
/// </summary>
public Vector<double> LowerBound { get; set; }
/// <summary>
/// Set or get the values of the parameters.
/// </summary>
public Vector<double> UpperBound { get; set; }
/// <summary>
/// Set or get the scale factor of the parameters.
/// </summary>
public Vector<double> Scales { get; set; }
/// <summary>
/// Set of get whether or not the parameters are bounded.
/// </summary>
public bool IsBounded { get; set; }
/// <summary>
/// Get the y-values of the fitted model that correspond to the independent values.
/// </summary>
public Vector<double> Values { get; private set; }
/// <summary>
/// Get the error values, R(x; p) = L * (y - f(x; p)) where L = sqrt(W)
/// </summary>
private Vector<double> Residuals;
/// <summary>
/// Get the residual sum of squares, R.DotProduct(R)
/// </summary>
public double Residue { get; private set; }
/// <summary>
/// Get the Jacobian matrix of x and p, J(x; p).
/// </summary>
public Matrix<double> Jacobian { get; private set; }
/// <summary>
/// Get the Gradient vector of x and p, J'WR
/// </summary>
public Vector<double> Gradient { get; private set; }
/// <summary>
/// Get the Hessian matrix of x and p, J'WJ
/// </summary>
public Matrix<double> Hessian { get; private set; }
/// <summary>
/// Get the number of observations.
/// </summary>
public int NumberOfObservations { get { return (ObservedY == null) ? 0 : ObservedY.Count; } }
/// <summary>
/// Get the number of unknown parameters.
/// </summary>
public int NumberOfParameters { get { return (Parameters == null) ? 0 : Parameters.Count; } }
/// <summary>
/// Get the degree of freedom
/// </summary>
public int DegreeOfFreedom
{
get
{
var dof = NumberOfObservations - NumberOfParameters;
if (IsFixed != null)
{
dof = dof + IsFixed.Count(p => p == true);
}
return dof;
}
}
/// <summary>
/// Get the covariance matrix.
/// </summary>
public Matrix<double> Covariance { get; private set; }
/// <summary>
/// Get the number of calls to function.
/// </summary>
public int FunctionEvaluations { get; private set; }
/// <summary>
/// Get the number of calls to jacobian.
/// </summary>
public int JacobianEvaluations { get; private set; }
/// <summary>
/// Set or get the desired accuracy order of the numerical jacobian.
/// </summary>
public int AccuracyOrder { get; set; }
/// <summary>
/// Get whether or not the analytical jacobian is supported.
/// </summary>
public bool IsJacobianSupported { get { return userDerivatives != null; } }
#endregion Public Variables
public FittingObjectiveModel(Func<Vector<double>, double, double>function, Func<Vector<double>, double, Vector<double>> derivatives, int accuracyOrder = 2)
{
userFunction = function;
userDerivatives = derivatives;
AccuracyOrder = Math.Min(6, Math.Max(1, accuracyOrder));
}
public IObjectiveModel Fork()
{
return new FittingObjectiveModel(userFunction, userDerivatives, AccuracyOrder)
{
ObservedX = ObservedX,
ObservedY = ObservedY,
Weights = Weights,
Parameters = Parameters,
LowerBound = LowerBound,
UpperBound = UpperBound,
IsFixed = IsFixed,
Scales = Scales,
IsBounded = IsBounded,
Residue = Residue,
Jacobian = Jacobian
};
}
public IObjectiveModel CreateNew()
{
return new FittingObjectiveModel(userFunction, userDerivatives);
}
/// <summary>
/// Set observed data to fit.
/// </summary>
public void SetObserved(Vector<double> observedX, Vector<double> observedY, Vector<double> weights = null)
{
if (observedX == null || observedY == null)
{
throw new ArgumentNullException("The data set can't be null.");
}
if (observedX.Count != observedY.Count)
{
throw new ArgumentException("The observed x data can't have different from observed y data.");
}
ObservedX = observedX;
ObservedY = observedY;
if (weights != null && weights.Count != observedY.Count)
{
throw new ArgumentException("The weightings can't have different from observations.");
}
if (weights != null && weights.Count(x => double.IsInfinity(x) || double.IsNaN(x)) > 0)
{
throw new ArgumentException("The weightings are not well-defined.");
}
if (weights != null && weights.Count(x => x == 0) == weights.Count)
{
throw new ArgumentException("All the weightings can't be zero.");
}
if (weights != null && weights.Count(x => x < 0) > 0)
{
weights = weights.PointwiseAbs();
}
Weights = (weights == null)
? null
: Matrix<double>.Build.DenseOfDiagonalVector(weights);
L = (weights == null)
? null
: Weights.Diagonal().PointwiseSqrt();
}
/// <summary>
/// Set observed data to fit.
/// </summary>
public void SetObserved(double[] observedX, double[] observedY, double[] weights = null)
{
if (observedX == null || observedY == null)
{
throw new ArgumentNullException("The data set can't be null.");
}
if (observedX.Length != observedY.Length)
{
throw new ArgumentException("The observed x data can't have different from observed y data.");
}
var wVector = (weights == null)
? null
: Vector<double>.Build.DenseOfArray(weights);
SetObserved(Vector<double>.Build.DenseOfArray(observedX), Vector<double>.Build.DenseOfArray(observedY), wVector);
}
/// <summary>
/// Set parameters.
/// <para/>
/// If bounded, the paramneters will be projected to unconstrained range by the mapping rule from the MINPACK.
/// If the projection is not needed, set IsBounded = false befre calling the Minimization method.
/// </summary>
/// <param name="lowerBound">The lower bounds of parameters.</param>
/// <param name="upperBound">The upper bounds of parameters.</param>
/// /// <param name="scales">The scaling constants of parameters</param>
/// <param name="isFixed">The list to the parameters fix or free.</param>
public void SetParameters(Vector<double> lowerBound = null, Vector<double> upperBound = null, Vector<double> scales = null, List<bool> isFixed = null)
{
if (lowerBound != null && lowerBound.Count(x => double.IsInfinity(x) || double.IsNaN(x)) > 0)
{
throw new ArgumentException("The lower bounds must be finite.");
}
LowerBound = lowerBound;
if (upperBound != null && upperBound.Count(x => double.IsInfinity(x) || double.IsNaN(x)) > 0)
{
throw new ArgumentException("The upper bounds must be finite.");
}
if (upperBound != null && lowerBound != null && upperBound.Count != lowerBound.Count)
{
throw new ArgumentException("The upper bounds can't have different elements from the lower bounds.");
}
UpperBound = upperBound;
if (scales != null && scales.Count(x => double.IsInfinity(x) || double.IsNaN(x) || x == 0) > 0)
{
throw new ArgumentException("The scales must be finite.");
}
if (scales != null && lowerBound != null && scales.Count != lowerBound.Count)
{
throw new ArgumentException("The upper bounds can't have different elements from the lower bounds.");
}
if (scales != null && upperBound != null && scales.Count != upperBound.Count)
{
throw new ArgumentException("The upper bounds can't have different elements from the upper bounds.");
}
if (scales != null && scales.Count(x => x < 0) > 0)
{
scales.PointwiseAbs();
}
Scales = scales;
IsBounded = (LowerBound != null || UpperBound != null || Scales != null);
if (isFixed != null && lowerBound != null && isFixed.Count != lowerBound.Count)
{
throw new ArgumentException("The initial guess can't have different elements from the lower bounds.");
}
if (isFixed != null && upperBound != null && isFixed.Count != upperBound.Count)
{
throw new ArgumentException("The initial guess can't have different elements from the upper bounds.");
}
if (isFixed != null && scales != null && isFixed.Count != scales.Count)
{
throw new ArgumentException("The initial guess can't have different elements from the scales.");
}
if (isFixed != null && isFixed.Count(p => p == true) == isFixed.Count)
{
throw new ArgumentException("All the parameters can't be fixed.");
}
IsFixed = isFixed;
}
/// <summary>
/// Set parameters.
/// <para/>
/// If bounded, the paramneters will be projected to unconstrained range by the mapping rule.
/// If the projection is not needed, set IsBounded = false befre calling the Minimization method.
/// </summary>
/// <param name="lowerBound">The lower bounds of parameters.</param>
/// <param name="upperBound">The upper bounds of parameters.</param>
/// <param name="scales">The scaling constants of parameters</param>
/// <param name="isFixed">The list to the parameters fix or free.</param>
public void SetParameters(double[] lowerBound = null, double[] upperBound = null, double[] scales = null, bool[] isFixed = null)
{
var lb = (lowerBound == null) ? null : Vector<double>.Build.DenseOfArray(lowerBound);
var ub = (upperBound == null) ? null : Vector<double>.Build.DenseOfArray(upperBound);
var sc = (scales == null) ? null : Vector<double>.Build.DenseOfArray(scales);
var fp = (isFixed == null) ? null : isFixed.ToList();
SetParameters(lb, ub, sc, fp);
}
public void EvaluateFunction(Vector<double> parameters)
{
ValidateParameters(parameters);
// To handle the box constrained minimization as the unconstrained minimization,
// parameters are mapping by the following rule.
//
// 1. lower < P < upper
// Pint = asin(2 * (Pext - lower) / (upper - lower) - 1)
// Pext = lower + (sin(Pint) + 1) * (upper - lower) / 2
// dPext/dPint = (upper - lower) / 2 * cos(Pint)
// 2. lower < P
// Pint = sqrt((Pext - lower + 1)^2 - 1)
// Pext = lower - 1 + sqrt(Pint^2 + 1)
// dPext/dPint = Pint / sqrt(Pint^2 + 1)
// 3. P < upper
// Pint = sqrt((upper - Pext + 1)^2 - 1)
// Pext = upper + 1 - sqrt(Pint^2 + 1)
// dPext/dPint = - Pint / sqrt(Pint^2 + 1)
// 4. no bounds, but scales
// Pint = Pext / scale
// Pext = Pint * scale
// dPext/dPint = scale
//
// see ProjectParametersToInternal(Pext), ProjectParametersToExternal(Pint), ScaleFactorsOfJacobian(Pint) methods.
//
// References:
// [1] https://lmfit.github.io/lmfit-py/bounds.html
//
//
// Except when it is initial guess, the parameters argument is always internal parameter.
// So, first map the parameters argument to the external parameters in order to calculate function values.
var Pext = (FunctionEvaluations > 0 && this.IsBounded)
? ProjectParametersToExternal(parameters)
: parameters.Clone();
// Project parameters, now this.Parameters are the internal parameters.
Parameters = (this.IsBounded)
? ProjectParametersToInternal(Pext)
: Pext;
// Calculates the residuals, (y[i] - f(x[i]; p)) * L[i]
if (Values == null)
{
Values = Vector<double>.Build.Dense(NumberOfObservations);
}
for (int i = 0; i < NumberOfObservations; i++)
{
Values[i] = userFunction(Pext, ObservedX[i]);
}
FunctionEvaluations++;
// calculate the weighted residuals
Residuals = (Weights == null)
? ObservedY - Values
: (ObservedY - Values).PointwiseMultiply(L);
// Calculate the residual sum of squares
Residue = Residuals.DotProduct(Residuals);
return;
}
public void EvaluateJacobian(Vector<double> parameters)
{
var Pext = (IsBounded)
? ProjectParametersToExternal(parameters)
: parameters.Clone();
// Calculates the jacobian of x and p.
if (userDerivatives != null)
{
// analytical jacobian
if (Jacobian == null)
{
Jacobian = Matrix<double>.Build.Dense(NumberOfObservations, NumberOfParameters);
}
for (int i = 0; i < NumberOfObservations; i++)
{
Jacobian.SetRow(i, userDerivatives(Pext, ObservedX[i]));
}
JacobianEvaluations++;
}
else
{
// numerical jacobian
Jacobian = NumericalJacobian(Pext, Values, AccuracyOrder);
FunctionEvaluations += AccuracyOrder;
}
var scaleFactors = (this.IsBounded)
? ScaleFactorsOfJacobian(Parameters)
: Vector<double>.Build.Dense(Parameters.Count, 1.0);
// Jint(x; Pint) = Jext(x; Pext) * scale where scale = dPext/dPint
for (int i = 0; i < NumberOfObservations; i++)
{
for (int j = 0; j < NumberOfParameters; j++)
{
if (IsFixed != null && IsFixed[j])
{
// if j-th parameter is fixed, set J[i, j] = 0
Jacobian[i, j] = 0.0;
}
else
{
Jacobian[i, j] = Jacobian[i, j] * scaleFactors[j];
}
}
}
// Gradient, g = J'W(y − f(x; p)) = J'L(L'E) = J'LR
Gradient = (Weights == null)
? Jacobian.Transpose() * (ObservedY - Values)
: Jacobian.Transpose() * Weights * (ObservedY - Values);
// approximated Hessian, H = J'WJ + ∑LRiHi ~ J'WJ near the minimum
Hessian = (Weights == null)
? Jacobian.Transpose() * Jacobian
: Jacobian.Transpose() * Weights * Jacobian;
}
public void EvaluateCovariance(Vector<double> parameters)
{
// convert to bounded(external) parameters
var Pext = (IsBounded)
? ProjectParametersToExternal(parameters)
: parameters.Clone();
// set IsBounded = false to get external Parameters and covariance matrix
this.IsBounded = false;
EvaluateFunction(Pext);
EvaluateJacobian(Pext);
if (Hessian == null || Residuals == null || DegreeOfFreedom < 1)
{
Covariance = null;
return;
}
var covariance = Hessian.PseudoInverse() * Residuals.DotProduct(Residuals) / DegreeOfFreedom;
Covariance = covariance;
// restore isBounded
this.IsBounded = (LowerBound != null || UpperBound != null);
return;
}
private void ValidateParameters(Vector<double> parameters)
{
if (parameters == null)
{
throw new ArgumentNullException("parameters");
}
else if (parameters.Count(p => double.IsNaN(p) || double.IsInfinity(p)) > 0)
{
throw new ArgumentException("the parameters must be finite.");
}
if (LowerBound != null && parameters.Count != LowerBound.Count)
{
throw new ArgumentException("The parameters can't have different size from the lower bounds.");
}
if (UpperBound != null && parameters.Count != UpperBound.Count)
{
throw new ArgumentException("The parameters can't have different size from the upper bounds.");
}
if (Scales != null && parameters.Count != Scales.Count)
{
throw new ArgumentException("The parameters can't have different size from the scales.");
}
if (IsFixed != null && parameters.Count != IsFixed.Count)
{
throw new ArgumentException("The parameters can't have different size from the IsFixed list.");
}
}
#region Numerical Derivatives
// Numerical derivatives by using the central or forward finite difference
private Matrix<double> NumericalJacobian(Vector<double> parameters, Vector<double> currentValues, int accuracyOrder = 2)
{
const double sqrtEpsilon = 1.4901161193847656250E-8; // sqrt(machineEpsilon)
Matrix<double> derivertives = Matrix<double>.Build.Dense(NumberOfObservations, NumberOfParameters);
var d = 0.000003 * parameters.PointwiseAbs().PointwiseMaximum(sqrtEpsilon);
var h = Vector<double>.Build.Dense(NumberOfParameters);
for (int i = 0; i < NumberOfObservations; i++)
{
var x = ObservedX[i];
for (int j = 0; j < NumberOfParameters; j++)
{
h[j] = d[j];
if (accuracyOrder >= 6)
{
// f'(x) = {- f(x - 3h) + 9f(x - 2h) - 45f(x - h) + 45f(x + h) - 9f(x + 2h) + f(x + 3h)} / 60h + O(h^6)
var f1 = userFunction(parameters - 3 * h, x);
var f2 = userFunction(parameters - 2 * h, x);
var f3 = userFunction(parameters - h, x);
var f4 = userFunction(parameters + h, x);
var f5 = userFunction(parameters + 2 * h, x);
var f6 = userFunction(parameters + 3 * h, x);
var prime = (-f1 + 9 * f2 - 45 * f3 + 45 * f4 - 9 * f5 + f6) / (60 * h[j]);
derivertives[i, j] = prime;
}
else if (accuracyOrder == 5)
{
// f'(x) = {-137f(x) + 300f(x + h) - 300f(x + 2h) + 200f(x + 3h) - 75f(x + 4h) + 12f(x + 5h)} / 60h + O(h^5)
var f1 = currentValues[i];
var f2 = userFunction(parameters + h, x);
var f3 = userFunction(parameters + 2 * h, x);
var f4 = userFunction(parameters + 3 * h, x);
var f5 = userFunction(parameters + 4 * h, x);
var f6 = userFunction(parameters + 5 * h, x);
var prime = (-137 * f1 + 300 * f2 - 300 * f3 + 200 * f4 - 75 * f5 + 12 * f6) / (60 * h[j]);
derivertives[i, j] = prime;
}
else if (accuracyOrder == 4)
{
// f'(x) = {f(x - 2h) - 8f(x - h) + 8f(x + h) - f(x + 2h)} / 12h + O(h^4)
var f1 = userFunction(parameters - 2 * h, x);
var f2 = userFunction(parameters - h, x);
var f3 = userFunction(parameters + h, x);
var f4 = userFunction(parameters + 2 * h, x);
var prime = (f1 - 8 * f2 + 8 * f3 - f4) / (12 * h[j]);
derivertives[i, j] = prime;
}
else if (accuracyOrder == 3)
{
// f'(x) = {-11f(x) + 18f(x + h) - 9f(x + 2h) + 2f(x + 3h)} / 6h + O(h^3)
var f1 = currentValues[i];
var f2 = userFunction(parameters + h, x);
var f3 = userFunction(parameters + 2 * h, x);
var f4 = userFunction(parameters + 3 * h, x);
var prime = (-11 * f1 + 18 * f2 - 9 * f3 + 2 * f4) / (6 * h[j]);
derivertives[i, j] = prime;
}
else if (accuracyOrder == 2)
{
// f'(x) = {f(x + h) - f(x - h)} / 2h + O(h^2)
var f1 = userFunction(parameters + h, x);
var f2 = userFunction(parameters - h, x);
var prime = (f1 - f2) / (2 * h[j]);
derivertives[i, j] = prime;
}
else
{
// f'(x) = {- f(x) + f(x + h)} / h + O(h)
var f1 = currentValues[i];
var f2 = userFunction(parameters + h, x);
var prime = (-f1 + f2) / h[j];
derivertives[i, j] = prime;
}
h[j] = 0;
}
}
return derivertives;
}
#endregion Numerical Derivatives
#region Projection
private Vector<double> ProjectParametersToInternal(Vector<double> Pext)
{
var Pint = Pext.Clone();
if (LowerBound != null && UpperBound != null)
{
for (int i = 0; i < Pext.Count; i++)
{
Pint[i] = Math.Asin((2.0 * (Pext[i] - LowerBound[i]) / (UpperBound[i] - LowerBound[i])) - 1.0);
}
return Pint;
}
else if (LowerBound != null && UpperBound == null)
{
for (int i = 0; i < Pext.Count; i++)
{
Pint[i] = Math.Sqrt(Math.Pow(Pext[i] - LowerBound[i] + 1.0, 2) - 1.0);
}
return Pint;
}
else if (LowerBound == null && UpperBound != null)
{
for (int i = 0; i < Pext.Count; i++)
{
Pint[i] = Math.Sqrt(Math.Pow(UpperBound[i] - Pext[i] + 1.0, 2) - 1.0);
}
return Pint;
}
else if (Scales != null)
{
for (int i = 0; i < Pext.Count; i++)
{
Pint[i] = Pext[i] / Scales[i];
}
return Pint;
}
return Pint;
}
private Vector<double> ProjectParametersToExternal(Vector<double> Pint)
{
var Pext = Pint.Clone();
if (LowerBound != null && UpperBound != null)
{
for (int i = 0; i < Pint.Count; i++)
{
Pext[i] = LowerBound[i] + (UpperBound[i] / 2.0 - LowerBound[i] / 2.0) * (Math.Sin(Pint[i]) + 1.0);
}
return Pext;
}
else if (LowerBound != null && UpperBound == null)
{
for (int i = 0; i < Pint.Count; i++)
{
Pext[i] = LowerBound[i] + Math.Sqrt(Pint[i] * Pint[i] + 1.0) - 1.0;
}
return Pext;
}
else if (LowerBound == null && UpperBound != null)
{
for (int i = 0; i < Pint.Count; i++)
{
Pext[i] = UpperBound[i] - Math.Sqrt(Pint[i] * Pint[i] + 1.0) + 1.0;
}
return Pext;
}
else if (Scales != null)
{
for (int i = 0; i < Pint.Count; i++)
{
Pext[i] = Pint[i] * Scales[i];
}
return Pext;
}
return Pext;
}
private Vector<double> ScaleFactorsOfJacobian(Vector<double> Pint)
{
var scale = Vector<double>.Build.Dense(Pint.Count, 1.0);
if (LowerBound != null && UpperBound != null)
{
for (int i = 0; i < Pint.Count; i++)
{
scale[i] = (UpperBound[i] - LowerBound[i]) / 2.0 * Math.Cos(Pint[i]);
}
return scale;
}
else if (LowerBound != null && UpperBound == null)
{
for (int i = 0; i < Pint.Count; i++)
{
scale[i] = Pint[i] / Math.Sqrt(Pint[i] * Pint[i] + 1.0);
}
return scale;
}
else if (LowerBound == null && UpperBound != null)
{
for (int i = 0; i < Pint.Count; i++)
{
scale[i] = -Pint[i] / Math.Sqrt(Pint[i] * Pint[i] + 1.0);
}
return scale;
}
else if (Scales != null)
{
return Scales;
}
return scale;
}
#endregion Projection
}
}

322
src/Numerics/Optimization/TrustRegionDogLegMinimizer.cs
View File

@ -0,0 +1,322 @@
using MathNet.Numerics.LinearAlgebra;
using System;
namespace MathNet.Numerics.Optimization
{
public sealed class TrustRegionDogLegMinimizer
{
#region Tolerances and options
/// <summary>
/// The stopping threshold for infinity norm of the gradient.
/// </summary>
public static double GradientTolerance { get; set; }
/// <summary>
/// The stopping threshold for L2 norm of the change of the parameters.
/// </summary>
public static double StepTolerance { get; set; }
/// <summary>
/// The stopping threshold for the function value or L2 norm of the residuals.
/// </summary>
public static double FunctionTolerance { get; set; }
/// <summary>
/// The stopping threshold for the trust region radius.
/// </summary>
public static double RadiusTolerance { get; set; }
/// <summary>
/// The maximum number of iterations.
/// </summary>
public int MaximumIterations { get; set; }
#endregion Tolerances and options
public TrustRegionDogLegMinimizer(double gradientTolerance = 1E-8, double stepTolerance = 1E-8, double functionTolerance = 1E-8, double radiusTolerance = 1E-8, int maximumIterations = -1)
{
FunctionTolerance = functionTolerance;
GradientTolerance = gradientTolerance;
StepTolerance = stepTolerance;
RadiusTolerance = radiusTolerance;
MaximumIterations = maximumIterations;
}
public ModelMinimizationResult FindMinimum(IObjectiveModel objective, Vector<double> initialGuess)
{
if (objective == null)
throw new ArgumentNullException("objective");
if (initialGuess == null)
throw new ArgumentNullException("initialGuess");
return Minimum(objective, initialGuess, GradientTolerance, StepTolerance, FunctionTolerance, RadiusTolerance, MaximumIterations);
}
public ModelMinimizationResult FindMinimum(IObjectiveModel objective, double[] initialGuess)
{
if (objective == null)
throw new ArgumentNullException("objective");
if (initialGuess == null)
throw new ArgumentNullException("initialGuess");
return Minimum(objective, CreateVector.DenseOfArray<double>(initialGuess), GradientTolerance, StepTolerance, FunctionTolerance, RadiusTolerance, MaximumIterations);
}
/// <summary>
/// Non-linear least square fitting by the trust-region dogleg algorithm.
/// </summary>
/// <param name="objective">The objective model, including function, jacobian, observations, and parameter bounds.</param>
/// <param name="initialGuess">The initial guess values.</param>
/// <param name="functionTolerance">The stopping threshold for L2 norm of the residuals.</param>
/// <param name="gradientTolerance">The stopping threshold for infinity norm of the gradient vector.</param>
/// <param name="stepTolerance">The stopping threshold for L2 norm of the change of parameters.</param>
/// <param name="radiusTolerance">The stopping threshold for trust region radius</param>
/// <param name="maximumIterations">The max iterations.</param>
/// <returns></returns>
public static ModelMinimizationResult Minimum(IObjectiveModel objective, Vector<double> initialGuess, double gradientTolerance = 1E-8, double stepTolerance = 1E-8, double functionTolerance = 1E-8, double radiusTolerance = 1E-18, int maximumIterations = -1)
{
// Non-linear least square fitting by the trust-region dogleg algorithm.
//
// The Powell's dogleg method is finding the minimum of a function F(p) that is a sum of squares of nonlinear functions.
//
// For given datum pair (x, y), uncertainties σ (or weighting W = 1 / σ^2) and model function f = f(x; p),
// let's find the parameters of the model so that the sum of the quares of the deviations is minimized.
//
// F(p) = 1/2 * ∑{ Wi * (yi - f(xi; p))^2 }
// pbest = argmin F(p)
//
// Here, we will use the following terms:
// Weighting W is the diagonal matrix and can be decomposed as LL', so L = 1/σ
// Residuals, R = L(y - f(x; p))
// Residual sum of squares, RSS = ||R||^2 = R.DotProduct(R)
// Jacobian J = df(x; p)/dp
// Gradient g = J'W(y − f(x; p)) = J'LR
// Approximated Hessian H = J'WJ
//
// The Powell's dogleg algorithm is summarized as follows:
// initially set trust-region radius, Δ
// repeat
// solve quadratic subproblem
// update Δ:
// let ρ = (RSS - RSSnew) / predRed
// if ρ > 0.75, Δ = 2Δ
// if ρ < 0.25, Δ = Δ/4
// if ρ > eta, P = P + ΔP
//
// References:
// [1]. Madsen, K., H. B. Nielsen, and O. Tingleff.
// "Methods for Non-Linear Least Squares Problems. Technical University of Denmark, 2004. Lecture notes." (2004).
// Available Online from: http://orbit.dtu.dk/files/2721358/imm3215.pdf
// [2]. SciPy
// Available Online from: https://github.com/scipy/scipy/blob/master/scipy/optimize/_trustregion_dogleg.py
double maxDelta = 1000;
double eta = 0;
if (objective == null)
throw new ArgumentNullException("objective");
if (initialGuess == null)
throw new ArgumentNullException("initialGuess");
ExitCondition exitCondition = ExitCondition.None;
// First, calculate function values and setup variables
objective.EvaluateFunction(initialGuess);
var P = objective.Parameters; // current parameters
var RSS = objective.Residue; // Residual Sum of Squares = R'R
var RSSinit = RSS; // RSS at initial gussing parameters
if (maximumIterations < 0)
{
maximumIterations = 200 * (initialGuess.Count + 1);
}
// if R == NaN, stop
if (double.IsNaN(RSS))
{
exitCondition = ExitCondition.InvalidValues;
return new ModelMinimizationResult(objective, -1, exitCondition);
}
// When only function evaluation is needed, set maximumIterations to zero,
if (maximumIterations == 0)
{
exitCondition = ExitCondition.ManuallyStopped;
}
// if ||R||^2 <= fTol, stop
if (RSS <= functionTolerance)
{
exitCondition = ExitCondition.Converged; // SmallRSS
}
// Evaluate projected Hessian, and gradient
objective.EvaluateJacobian(P);
var Hessian = objective.Hessian;
var Gradient = objective.Gradient;
// if ||g||_oo <= gtol, found and stop
if (Gradient.InfinityNorm() <= gradientTolerance)
{
exitCondition = ExitCondition.RelativeGradient; // SmallGradient
}
if (exitCondition != ExitCondition.None)
{
// finalize
objective.EvaluateCovariance(P);
return new ModelMinimizationResult(objective, -1, exitCondition);
}
// initialize trust-region radius, Δ
double delta = Gradient.DotProduct(Gradient) / (Hessian * Gradient).DotProduct(Gradient);
delta = Math.Max(1.0, Math.Min(delta, maxDelta));
int iterations = 0;
while (iterations < maximumIterations && exitCondition == ExitCondition.None)
{
iterations++;
// solve the subproblem
var subprogram = SolveQuadraticSubproblem(objective, delta);
var Pstep = subprogram.Item1;
var predictedReduction = subprogram.Item2;
var hitBoundary = subprogram.Item3;
if (Pstep.L2Norm() <= stepTolerance * (stepTolerance + P.L2Norm()))
{
exitCondition = ExitCondition.RelativePoints; // SmallRelativeParameters
break;
}
var Pnew = P + Pstep; // parameters to test
objective.EvaluateFunction(Pnew);
var RSSnew = objective.Residue;
// calculate the ratio of the actual to the predicted reduction.
double rho = (predictedReduction != 0)
? (RSS - RSSnew) / predictedReduction
: 0;
if (rho > 0.75 && hitBoundary)
{
delta = Math.Min(2.0 * delta, maxDelta);
}
else if (rho < 0.25)
{
delta = delta * 0.25;
if (delta <= radiusTolerance * (radiusTolerance + P.DotProduct(P)))
{
exitCondition = ExitCondition.RelativePoints; // SmallRelativeParameters
break;
}
}
if (rho > eta)
{
// accepted
Pnew.CopyTo(P);
RSS = RSSnew;
// update Jacobian, Hessian, and gradient
objective.EvaluateJacobian(P);
Gradient = objective.Gradient;
Hessian = objective.Hessian;
// if ||g||_oo <= gtol, found and stop
if (Gradient.InfinityNorm() <= gradientTolerance)
{
exitCondition = ExitCondition.RelativeGradient;
}
// if ||R||^2 < fTol, found and stop
if (RSS <= functionTolerance)
{
exitCondition = ExitCondition.Converged; // SmallRSS
}
}
}
if (iterations >= maximumIterations)
{
exitCondition = ExitCondition.ExceedIterations;
}
// finalize
objective.EvaluateCovariance(P);
return new ModelMinimizationResult(objective, iterations, exitCondition);
}
private static Tuple<Vector<double>, double, bool> SolveQuadraticSubproblem(IObjectiveModel objective, double delta)
{
Vector<double> Pstep;
double predictedReduction;
bool hitBoundary = false;
var Jacobian = objective.Jacobian;
var Gradient = objective.Gradient;
var Hessian = objective.Hessian;
var RSS = objective.Residue;
// newton point
// the Gauss–Newton step by solving the normal equations
var Pgn = Hessian.PseudoInverse() * Gradient; // Hessian.Solve(Gradient) fails so many times...
// cauchy point
// steepest descent direction is given by
var alpha = Gradient.DotProduct(Gradient) / (Hessian * Gradient).DotProduct(Gradient);
var Psd = alpha * Gradient;
// update step and prectted reduction
if (Pgn.L2Norm() <= delta)
{
// Pgn is inside trust region radius
hitBoundary = false;
Pstep = Pgn;
predictedReduction = RSS;
}
else if (alpha * Psd.L2Norm() >= delta)
{
// Psd is outside trust region radius
hitBoundary = true;
Pstep = delta / Psd.L2Norm() * Psd;
predictedReduction = delta * (2.0 * (alpha * Gradient).L2Norm() - delta) / 2.0 / alpha;
}
else
{
// Pstep is intersection of the trust region boundary
hitBoundary = true;
var beta = FindBeta(alpha, Psd, Pgn, delta);
Pstep = alpha * Psd + beta * (Pgn - alpha * Psd);
predictedReduction = 0.5 * alpha * (1 - beta) * (1 - beta) * Gradient.DotProduct(Gradient) + beta * (2 - beta) * RSS;
}
return new Tuple<Vector<double>, double, bool>(Pstep, predictedReduction, hitBoundary);
}
private static double FindBeta(double alpha, Vector<double> sd, Vector<double> gn, double delta)
{
// Pstep is intersection of the trust region boundary
// Pstep = α*Psd + β*(Pgn - α*Psd)
// find r so that ||Pstep|| = Δ
// z = α*Psd, d = (Pgn - z)
// (d^2)β^2 + (2*z*d)β + (z^2 - Δ^2) = 0
// get positive β by using the quadratic formula
var z = alpha * sd;
var d = gn - z;
var a = d.DotProduct(d);
var b = 2.0 * z.DotProduct(d);
var c = z.DotProduct(z) - delta * delta;
var aux = b + ((b >= 0) ? 1.0 : -1.0) * Math.Sqrt(b * b - 4.0 * a * c);
var beta = Math.Max(-aux / 2.0 / a, -2.0 * c / aux);
return beta;
}
}
}
Write Preview
Loading…
Cancel
Save