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

Merge pull request #614 from diluculo/LMTR 

Nonlinear least squares fitting
pull/627/merge
Christoph Ruegg 8 years ago
committed by GitHub
parent
7e45066368 19bd5f576d
commit
8aa6a56910
No known key found for this signature in database GPG Key ID: 4AEE18F83AFDEB23
16 changed files with 2300 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. 621
      src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs
  2. 6
      src/Numerics/Optimization/ExitCondition.cs
  3. 74
      src/Numerics/Optimization/IObjectiveModel.cs
  4. 12
      src/Numerics/Optimization/ITrustRegionSubProblem.cs
  5. 234
      src/Numerics/Optimization/LevenbergMarquardtMinimizer.cs
  6. 78
      src/Numerics/Optimization/NonlinearMinimizationResult.cs
  7. 302
      src/Numerics/Optimization/NonlinearMinimizerBase.cs
  8. 106
      src/Numerics/Optimization/ObjectiveFunction.cs
  9. 436
      src/Numerics/Optimization/ObjectiveFunctions/NonlinearObjectiveFunction.cs
  10. 45
      src/Numerics/Optimization/Subproblems/DogLegSubproblem.cs
  11. 65
      src/Numerics/Optimization/Subproblems/NewtonCGSubproblem.cs
  12. 35
      src/Numerics/Optimization/Subproblems/Util.cs
  13. 12
      src/Numerics/Optimization/TrustRegionDogLegMinimizer.cs
  14. 245
      src/Numerics/Optimization/TrustRegionMinimizerBase.cs
  15. 13
      src/Numerics/Optimization/TrustRegionNewtonCGMinimizer.cs
  16. 17
      src/Numerics/Optimization/TrustRegionSubProblem.cs

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

@ -0,0 +1,621 @@
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
{
#region Rosenbrock
// 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 Vector<double> RosenbrockModel(Vector<double> p, Vector<double> x)
{
var y = CreateVector.Dense<double>(x.Count);
for (int i = 0; i < x.Count; i++)
{
y[i] = Math.Pow(1.0 - p[0], 2) + 100.0 * Math.Pow(p[1] - p[0] * p[0], 2);
}
return y;
}
private Matrix<double> RosenbrockPrime(Vector<double> p, Vector<double> x)
{
var prime = Matrix<double>.Build.Dense(x.Count, p.Count);
for (int i = 0; i < x.Count; i++)
{
prime[i, 0] = 400.0 * p[0] * p[0] * p[0] - 400.0 * p[0] * p[1] + 2.0 * p[0] - 2.0;
prime[i, 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 Rosenbrock_LM_Der()
{
// unconstrained
var obj = ObjectiveFunction.NonlinearModel(RosenbrockModel, RosenbrockPrime, RosenbrockX, RosenbrockY);
var solver = new LevenbergMarquardtMinimizer(maximumIterations: 10000);
var result = solver.FindMinimum(obj, RosenbrockStart1);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.MinimizingPoint[i], 2);
}
// box constrained
obj = ObjectiveFunction.NonlinearModel(RosenbrockModel, RosenbrockPrime, RosenbrockX, RosenbrockY);
solver = new LevenbergMarquardtMinimizer(maximumIterations: 10000);
result = solver.FindMinimum(obj, RosenbrockStart1, lowerBound: RosebbrockLowerBound, upperBound: RosenbrockUpperBound);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.MinimizingPoint[i], 2);
}
}
[Test]
public void Rosenbrock_LM_Dif()
{
// unconstrained
var obj = ObjectiveFunction.NonlinearModel(RosenbrockModel, RosenbrockX, RosenbrockY, accuracyOrder:2);
var solver = new LevenbergMarquardtMinimizer(maximumIterations: 10000);
var result = solver.FindMinimum(obj, RosenbrockStart1);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.MinimizingPoint[i], 2);
}
// box constrained
obj = ObjectiveFunction.NonlinearModel(RosenbrockModel, RosenbrockX, RosenbrockY, accuracyOrder: 6);
solver = new LevenbergMarquardtMinimizer(maximumIterations: 10000);
result = solver.FindMinimum(obj, RosenbrockStart1, lowerBound: RosebbrockLowerBound, upperBound: RosenbrockUpperBound);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.MinimizingPoint[i], 2);
}
}
[Test]
public void Rosenbrock_Bfgs_Dif()
{
var obj = ObjectiveFunction.NonlinearFunction(RosenbrockModel, RosenbrockX, RosenbrockY, accuracyOrder: 6);
var solver = new BfgsMinimizer(1e-8, 1e-8, 1e-8, 1000);
var result = solver.FindMinimum(obj, RosenbrockStart1);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.MinimizingPoint[i], 2);
}
}
[Test]
public void Rosenbrock_LBfgs_Dif()
{
var obj = ObjectiveFunction.NonlinearFunction(RosenbrockModel, RosenbrockX, RosenbrockY, accuracyOrder: 6);
var solver = new LimitedMemoryBfgsMinimizer(1e-8, 1e-8, 1e-8, 1000);
var result = solver.FindMinimum(obj, RosenbrockStart1);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.MinimizingPoint[i], 2);
}
}
#endregion Rosenbrock
#region Rat43
// 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 Vector<double> Rat43Model(Vector<double> p, Vector<double> x)
{
var y = CreateVector.Dense<double>(x.Count);
for (int i = 0; i < x.Count; i++)
{
y[i] = p[0] / Math.Pow(1.0 + Math.Exp(p[1] - p[2] * x[i]), 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 Rat43_LM_Dif()
{
var obj = ObjectiveFunction.NonlinearModel(Rat43Model, Rat43X, Rat43Y, accuracyOrder: 6);
var solver = new LevenbergMarquardtMinimizer();
var result = solver.FindMinimum(obj, Rat43Start1);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(Rat43Pbest[i], result.MinimizingPoint[i], 6);
AssertHelpers.AlmostEqualRelative(Rat43Pstd[i], result.StandardErrors[i], 6);
}
}
[Test]
public void Rat43_TRDL_Dif()
{
var obj = ObjectiveFunction.NonlinearModel(Rat43Model, Rat43X, Rat43Y, accuracyOrder: 6);
var solver = new TrustRegionDogLegMinimizer();
var result = solver.FindMinimum(obj, Rat43Start2);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(Rat43Pbest[i], result.MinimizingPoint[i], 2);
AssertHelpers.AlmostEqualRelative(Rat43Pstd[i], result.StandardErrors[i], 2);
}
}
[Test]
public void Rat43_TRNCG_Dif()
{
var obj = ObjectiveFunction.NonlinearModel(Rat43Model, Rat43X, Rat43Y, accuracyOrder: 6);
var solver = new TrustRegionNewtonCGMinimizer();
var result = solver.FindMinimum(obj, Rat43Start2);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(Rat43Pbest[i], result.MinimizingPoint[i], 2);
AssertHelpers.AlmostEqualRelative(Rat43Pstd[i], result.StandardErrors[i], 2);
}
}
[Test]
public void Rat43_Bfgs_Dif()
{
var obj = ObjectiveFunction.NonlinearFunction(Rat43Model, Rat43X, Rat43Y, accuracyOrder: 6);
var solver = new BfgsMinimizer(1e-10, 1e-10, 1e-10, 1000);
var result = solver.FindMinimum(obj, Rat43Start2);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(Rat43Pbest[i], result.MinimizingPoint[i], 2);
}
}
[Test]
public void Rat43_LBfgs_Dif()
{
var obj = ObjectiveFunction.NonlinearFunction(Rat43Model, Rat43X, Rat43Y, accuracyOrder: 6);
var solver = new LimitedMemoryBfgsMinimizer(1e-10, 1e-10, 1e-10, 1000);
var result = solver.FindMinimum(obj, Rat43Start2);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(Rat43Pbest[i], result.MinimizingPoint[i], 2);
}
}
#endregion Rat43
#region BoxBod
// 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 Vector<double> BoxBodModel(Vector<double> p, Vector<double> x)
{
var y = CreateVector.Dense<double>(x.Count);
for (int i = 0; i < x.Count; i++)
{
y[i] = p[0] * (1.0 - Math.Exp(-p[1] * x[i]));
}
return y;
}
private Matrix<double> BoxBodPrime(Vector<double> p, Vector<double> x)
{
var prime = Matrix<double>.Build.Dense(x.Count, p.Count);
for (int i = 0; i < x.Count; i++)
{
prime[i, 0] = 1.0 - Math.Exp(-p[1] * x[i]);
prime[i, 1] = p[0] * x[i] * Math.Exp(-p[1] * x[i]);
}
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[] { -1000, -100 });
private Vector<double> BoxBodUpperBound = new DenseVector(new double[] { 1000.0, 100 });
private Vector<double> BoxBodScales = new DenseVector(new double[] { 100.0, 0.1 });
[Test]
public void BoxBod_LM_Der()
{
// unconstrained
var obj = ObjectiveFunction.NonlinearModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY);
var solver = new LevenbergMarquardtMinimizer();
var result = solver.FindMinimum(obj, BoxBodStart1);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.MinimizingPoint[i], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6);
}
// lower < parameters < upper
// Note that in this case, scales have no effect.
obj = ObjectiveFunction.NonlinearModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY);
solver = new LevenbergMarquardtMinimizer();
result = solver.FindMinimum(obj, BoxBodStart1, lowerBound: BoxBodLowerBound, upperBound: BoxBodUpperBound);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.MinimizingPoint[i], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6);
}
// lower < parameters, no scales
obj = ObjectiveFunction.NonlinearModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY);
solver = new LevenbergMarquardtMinimizer();
result = solver.FindMinimum(obj, BoxBodStart1, lowerBound: BoxBodLowerBound);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.MinimizingPoint[i], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6);
}
// lower < parameters, scales
obj = ObjectiveFunction.NonlinearModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY);
solver = new LevenbergMarquardtMinimizer();
result = solver.FindMinimum(obj, BoxBodStart1, lowerBound: BoxBodLowerBound, scales: BoxBodScales);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.MinimizingPoint[i], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6);
}
// parameters < upper, no scales
obj = ObjectiveFunction.NonlinearModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY);
solver = new LevenbergMarquardtMinimizer();
result = solver.FindMinimum(obj, BoxBodStart1, upperBound: BoxBodUpperBound);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.MinimizingPoint[i], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6);
}
// parameters < upper, scales
obj = ObjectiveFunction.NonlinearModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY);
solver = new LevenbergMarquardtMinimizer();
result = solver.FindMinimum(obj, BoxBodStart1, upperBound: BoxBodUpperBound, scales: BoxBodScales);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.MinimizingPoint[i], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6);
}
// only scales
obj = ObjectiveFunction.NonlinearModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY);
solver = new LevenbergMarquardtMinimizer();
result = solver.FindMinimum(obj, BoxBodStart1, scales: BoxBodScales);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.MinimizingPoint[i], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6);
}
}
[Test]
public void BoxBod_LM_Dif()
{
// unconstrained
var obj = ObjectiveFunction.NonlinearModel(BoxBodModel, BoxBodX, BoxBodY, accuracyOrder:6);
var solver = new LevenbergMarquardtMinimizer();
var result = solver.FindMinimum(obj, BoxBodStart1);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.MinimizingPoint[i], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6);
}
// box constrained
obj = ObjectiveFunction.NonlinearModel(BoxBodModel, BoxBodX, BoxBodY, accuracyOrder: 6);
solver = new LevenbergMarquardtMinimizer();
result = solver.FindMinimum(obj, BoxBodStart1, lowerBound: BoxBodLowerBound, upperBound: BoxBodUpperBound);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.MinimizingPoint[i], 6);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6);
}
}
[Test]
public void BoxBod_TRDL_Dif()
{
var obj = ObjectiveFunction.NonlinearModel(BoxBodModel, BoxBodX, BoxBodY, accuracyOrder: 6);
var solver = new TrustRegionDogLegMinimizer();
var result = solver.FindMinimum(obj, BoxBodStart1);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.MinimizingPoint[i], 3);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 3);
}
}
[Test]
public void BoxBod_TRNCG_Dif()
{
var obj = ObjectiveFunction.NonlinearModel(BoxBodModel, BoxBodX, BoxBodY, accuracyOrder: 6);
var solver = new TrustRegionNewtonCGMinimizer();
var result = solver.FindMinimum(obj, BoxBodStart2);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.MinimizingPoint[i], 3);
AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 3);
}
}
[Test]
public void BoxBod_Bfgs_Der()
{
var obj = ObjectiveFunction.NonlinearFunction(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY);
var solver = new BfgsMinimizer(1e-10, 1e-10, 1e-10, 100);
var result = solver.FindMinimum(obj, BoxBodStart2);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.MinimizingPoint[i], 6);
}
}
[Test]
public void BoxBod_Newton_Der()
{
var obj = ObjectiveFunction.NonlinearFunction(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY);
var solver = new NewtonMinimizer(1e-10, 100);
var result = solver.FindMinimum(obj, BoxBodStart2);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.MinimizingPoint[i], 6);
}
}
#endregion BoxBod
#region Thurber
// 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 Vector<double> ThurberModel(Vector<double> p, Vector<double> x)
{
var y = CreateVector.Dense<double>(x.Count);
for (int i = 0; i < x.Count; i++)
{
var xSq = x[i] * x[i];
var xCb = xSq * x[i];
y[i] = (p[0] + p[1] * x[i] + p[2] * xSq + p[3] * xCb)
/ (1 + p[4] * x[i] + p[5] * xSq + p[6] * xCb);
}
return y;
}
private Matrix<double> ThurberPrime(Vector<double> p, Vector<double> x)
{
var prime = Matrix<double>.Build.Dense(x.Count, p.Count);
for (int i = 0; i < x.Count; i++)
{
var xSq = x[i] * x[i];
var xCb = xSq * x[i];
var num = p[0] + x[i] * (p[1] + x[i] * (p[2] + p[3] * x[i]));
var den = p[4] * x[i] + p[5] * xSq + p[6] * xCb + 1.0;
var denSq = den * den;
prime[i, 0] = 1.0 / den;
prime[i, 1] = x[i] / den;
prime[i, 2] = xSq / den;
prime[i, 3] = xCb / den;
prime[i, 4] = -(x[i] * num) / denSq;
prime[i, 5] = -(xSq * num) / denSq;
prime[i, 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> ThurberStart = new DenseVector(new double[] { 1000.0, 1000.0, 400.0, 40.0, 0.7, 0.3, 0.03 });
private Vector<double> ThurberLowerBound = new DenseVector(new double[] { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 });
private Vector<double> ThurberUpperBound = new DenseVector(new double[] { 1E6, 1E6, 1E6, 1E6, 1E6, 1E6, 1E6 });
private Vector<double> ThurberScales = new DenseVector(new double[7] { 1000, 1000, 400, 40, 0.7, 0.3, 0.03 });
[Test]
public void Thurber_LM_Der()
{
var obj = ObjectiveFunction.NonlinearModel(ThurberModel, ThurberPrime, ThurberX, ThurberY);
var solver = new LevenbergMarquardtMinimizer();
var result = solver.FindMinimum(obj, ThurberStart);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(ThurberPbest[i], result.MinimizingPoint[i], 6);
AssertHelpers.AlmostEqualRelative(ThurberPstd[i], result.StandardErrors[i], 6);
}
}
[Test]
public void Thurber_LM_Dif()
{
var obj = ObjectiveFunction.NonlinearModel(ThurberModel, ThurberX, ThurberY, accuracyOrder: 6);
var solver = new LevenbergMarquardtMinimizer();
var result = solver.FindMinimum(obj, ThurberStart);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(ThurberPbest[i], result.MinimizingPoint[i], 6);
AssertHelpers.AlmostEqualRelative(ThurberPstd[i], result.StandardErrors[i], 6);
}
}
[Test]
public void Thurber_TRDL_Dif()
{
var obj = ObjectiveFunction.NonlinearModel(ThurberModel, ThurberX, ThurberY, accuracyOrder: 6);
var solver = new TrustRegionDogLegMinimizer();
var result = solver.FindMinimum(obj, ThurberStart, scales: ThurberScales);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(ThurberPbest[i], result.MinimizingPoint[i], 3);
AssertHelpers.AlmostEqualRelative(ThurberPstd[i], result.StandardErrors[i], 3);
}
}
[Test]
public void Thurber_TRNCG_Dif()
{
var obj = ObjectiveFunction.NonlinearModel(ThurberModel, ThurberX, ThurberY, accuracyOrder: 6);
var solver = new TrustRegionNewtonCGMinimizer();
var result = solver.FindMinimum(obj, ThurberStart, scales: ThurberScales);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(ThurberPbest[i], result.MinimizingPoint[i], 3);
AssertHelpers.AlmostEqualRelative(ThurberPstd[i], result.StandardErrors[i], 3);
}
}
[Test]
public void Thurber_Bfgs_Dif()
{
var obj = ObjectiveFunction.NonlinearFunction(ThurberModel, ThurberX, ThurberY, accuracyOrder: 6);
var solver = new BfgsMinimizer(1e-10, 1e-10, 1e-10, 1000);
var result = solver.FindMinimum(obj, ThurberStart);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(ThurberPbest[i], result.MinimizingPoint[i], 6);
}
}
[Test]
public void Thurber_BfgsB_Dif()
{
var obj = ObjectiveFunction.NonlinearFunction(ThurberModel, ThurberX, ThurberY, accuracyOrder: 6);
var solver = new BfgsBMinimizer(1e-10, 1e-10, 1e-10, 1000);
var result = solver.FindMinimum(obj, ThurberLowerBound, ThurberUpperBound, ThurberStart);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(ThurberPbest[i], result.MinimizingPoint[i], 6);
}
}
[Test]
public void Thurber_LBfgs_Dif()
{
var obj = ObjectiveFunction.NonlinearFunction(ThurberModel, ThurberX, ThurberY, accuracyOrder: 6);
var solver = new LimitedMemoryBfgsMinimizer(1e-10, 1e-10, 1e-10, 1000);
var result = solver.FindMinimum(obj, ThurberStart);
for (int i = 0; i < result.MinimizingPoint.Count; i++)
{
AssertHelpers.AlmostEqualRelative(ThurberPbest[i], result.MinimizingPoint[i], 6);
}
}
#endregion Thurber
}
}

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
}
}

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

@ -0,0 +1,74 @@
using MathNet.Numerics.LinearAlgebra;
using System.Collections.Generic;
namespace MathNet.Numerics.Optimization
{
public interface IObjectiveModelEvaluation
{
IObjectiveModel CreateNew();
/// <summary>
/// Get the y-values of the observations.
/// </summary>
Vector<double> ObservedY { get; }
/// <summary>
/// Get the values of the weights for the observations.
/// </summary>
Matrix<double> Weights { get; }
/// <summary>
/// Get the y-values of the fitted model that correspond to the independent values.
/// </summary>
Vector<double> ModelValues { get; }
/// <summary>
/// Get the values of the parameters.
/// </summary>
Vector<double> Point { get; }
/// <summary>
/// Get the residual sum of squares.
/// </summary>
double Value { 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 number of calls to function.
/// </summary>
int FunctionEvaluations { get; set; }
/// <summary>
/// Get the number of calls to jacobian.
/// </summary>
int JacobianEvaluations { get; set; }
/// <summary>
/// Get the degree of freedom.
/// </summary>
int DegreeOfFreedom { get; }
bool IsGradientSupported { get; }
bool IsHessianSupported { get; }
}
public interface IObjectiveModel : IObjectiveModelEvaluation
{
void SetParameters(Vector<double> initialGuess, List<bool> isFixed = null);
void EvaluateAt(Vector<double> parameters);
IObjectiveModel Fork();
IObjectiveFunction ToObjectiveFunction();
}
}

12
src/Numerics/Optimization/ITrustRegionSubProblem.cs
View File

@ -0,0 +1,12 @@
using MathNet.Numerics.LinearAlgebra;
namespace MathNet.Numerics.Optimization
{
public interface ITrustRegionSubproblem
{
Vector<double> Pstep { get; }
bool HitBoundary { get; }
void Solve(IObjectiveModel objective, double radius);
}
}

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

@ -0,0 +1,234 @@
using MathNet.Numerics.LinearAlgebra;
using System;
using System.Collections.Generic;
using System.Linq;
namespace MathNet.Numerics.Optimization
{
public class LevenbergMarquardtMinimizer : NonlinearMinimizerBase
{
/// <summary>
/// The scale factor for initial mu
/// </summary>
public static double InitialMu { get; set; }
public LevenbergMarquardtMinimizer(double initialMu = 1E-3, double gradientTolerance = 1E-15, double stepTolerance = 1E-15, double functionTolerance = 1E-15, int maximumIterations = -1)
: base(gradientTolerance, stepTolerance, functionTolerance, maximumIterations)
{
InitialMu = initialMu;
}
public NonlinearMinimizationResult FindMinimum(IObjectiveModel objective, Vector<double> initialGuess,
Vector<double> lowerBound = null, Vector<double> upperBound = null, Vector<double> scales = null, List<bool> isFixed = null)
{
return Minimum(objective, initialGuess, lowerBound, upperBound, scales, isFixed, InitialMu, GradientTolerance, StepTolerance, FunctionTolerance, MaximumIterations);
}
public NonlinearMinimizationResult FindMinimum(IObjectiveModel objective, double[] initialGuess,
double[] lowerBound = null, double[] upperBound = null, double[] scales = null, bool[] isFixed = null)
{
if (objective == null)
throw new ArgumentNullException("objective");
if (initialGuess == null)
throw new ArgumentNullException("initialGuess");
var lb = (lowerBound == null) ? null : CreateVector.Dense<double>(lowerBound);
var ub = (upperBound == null) ? null : CreateVector.Dense<double>(upperBound);
var sc = (scales == null) ? null : CreateVector.Dense<double>(scales);
var fx = (isFixed == null) ? null : isFixed.ToList();
return Minimum(objective, CreateVector.DenseOfArray<double>(initialGuess), lb, ub, sc, fx, 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 NonlinearMinimizationResult Minimum(IObjectiveModel objective, Vector<double> initialGuess,
Vector<double> lowerBound = null, Vector<double> upperBound = null, Vector<double> scales = null, List<bool> isFixed = null,
double initialMu = 1E-3, double gradientTolerance = 1E-15, double stepTolerance = 1E-15, double functionTolerance = 1E-15, 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(H)).
// repeat
// solve linear equations: (H + μI)ΔP = -g
// let ρ = (||R||^2 - ||Rnew||^2) / (Δp'(μΔp - g)).
// 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");
ValidateBounds(initialGuess, lowerBound, upperBound, scales);
objective.SetParameters(initialGuess, isFixed);
ExitCondition exitCondition = ExitCondition.None;
// First, calculate function values and setup variables
var P = ProjectToInternalParameters(initialGuess); // current internal parameters
var Pstep = Vector<double>.Build.Dense(P.Count); // the change of parameters
var RSS = EvaluateFunction(objective, P); // Residual Sum of Squares = R'R
if (maximumIterations < 0)
{
maximumIterations = 200 * (initialGuess.Count + 1);
}
// if RSS == NaN, stop
if (double.IsNaN(RSS))
{
exitCondition = ExitCondition.InvalidValues;
return new NonlinearMinimizationResult(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
var jac = EvaluateJacobian(objective, P);
var Gradient = jac.Item1; // objective.Gradient;
var Hessian = jac.Item2; // 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)
{
return new NonlinearMinimizationResult(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;
break;
}
var Pnew = P + Pstep; // new parameters to test
// evaluate function at Pnew
var RSSnew = EvaluateFunction(objective, Pnew);
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
jac = EvaluateJacobian(objective, P);
Gradient = jac.Item1; // objective.Gradient;
Hessian = jac.Item2; // 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;
}
return new NonlinearMinimizationResult(objective, iterations, exitCondition);
}
}
}

78
src/Numerics/Optimization/NonlinearMinimizationResult.cs
View File

@ -0,0 +1,78 @@
using MathNet.Numerics.LinearAlgebra;
namespace MathNet.Numerics.Optimization
{
public class NonlinearMinimizationResult
{
public IObjectiveModel ModelInfoAtMinimum { get; private set; }
/// <summary>
/// Returns the best fit parameters.
/// </summary>
public Vector<double> MinimizingPoint { get { return ModelInfoAtMinimum.Point; } }
/// <summary>
/// Returns the standard errors of the corresponding parameters
/// </summary>
public Vector<double> StandardErrors { get; private set; }
/// <summary>
/// Returns the y-values of the fitted model that correspond to the independent values.
/// </summary>
public Vector<double> MinimizedValues { get { return ModelInfoAtMinimum.ModelValues; } }
/// <summary>
/// Returns the covariance matrix at minimizing point.
/// </summary>
public Matrix<double> Covariance { get; private set; }
/// <summary>
/// Returns the correlation matrix at minimizing point.
/// </summary>
public Matrix<double> Correlation { get; private set; }
public int Iterations { get; private set; }
public ExitCondition ReasonForExit { get; private set; }
public NonlinearMinimizationResult(IObjectiveModel modelInfo, int iterations, ExitCondition reasonForExit)
{
ModelInfoAtMinimum = modelInfo;
Iterations = iterations;
ReasonForExit = reasonForExit;
EvaluateCovariance(modelInfo);
}
private void EvaluateCovariance(IObjectiveModel objective)
{
objective.EvaluateAt(objective.Point); // Hessian may be not yet updated.
var Hessian = objective.Hessian;
if (Hessian == null || objective.DegreeOfFreedom < 1)
{
Covariance = null;
Correlation = null;
StandardErrors = null;
return;
}
Covariance = Hessian.PseudoInverse() * objective.Value / objective.DegreeOfFreedom;
if (Covariance != null)
{
StandardErrors = Covariance.Diagonal().PointwiseSqrt();
var correlation = Covariance.Clone();
var d = correlation.Diagonal().PointwiseSqrt();
var dd = d.OuterProduct(d);
Correlation = correlation.PointwiseDivide(dd);
}
else
{
StandardErrors = null;
Correlation = null;
}
}
}
}

302
src/Numerics/Optimization/NonlinearMinimizerBase.cs
View File

@ -0,0 +1,302 @@
using MathNet.Numerics.LinearAlgebra;
using System;
using System.Linq;
namespace MathNet.Numerics.Optimization
{
public abstract class NonlinearMinimizerBase
{
/// <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 L2 norm of the change of the parameters.
/// </summary>
public static double StepTolerance { get; set; }
/// <summary>
/// The stopping threshold for infinity norm of the gradient.
/// </summary>
public static double GradientTolerance { get; set; }
/// <summary>
/// The maximum number of iterations.
/// </summary>
public static int MaximumIterations { get; set; }
/// <summary>
/// The lower bound of the parameters.
/// </summary>
public static Vector<double> LowerBound { get; private set; }
/// <summary>
/// The upper bound of the parameters.
/// </summary>
public static Vector<double> UpperBound { get; private set; }
/// <summary>
/// The scale factors for the parameters.
/// </summary>
public static Vector<double> Scales { get; private set; }
private static bool IsBounded { get { return LowerBound != null || UpperBound != null || Scales != null; } }
protected NonlinearMinimizerBase(double gradientTolerance = 1E-18, double stepTolerance = 1E-18, double functionTolerance = 1E-18, int maximumIterations = -1)
{
GradientTolerance = gradientTolerance;
StepTolerance = stepTolerance;
FunctionTolerance = functionTolerance;
MaximumIterations = maximumIterations;
}
protected static void ValidateBounds(Vector<double> parameters, Vector<double> lowerBound = null, Vector<double> upperBound = null, Vector<double> scales = null)
{
if (parameters == null)
{
throw new ArgumentNullException("parameters");
}
if (lowerBound != null && lowerBound.Count(x => double.IsInfinity(x) || double.IsNaN(x)) > 0)
{
throw new ArgumentException("The lower bounds must be finite.");
}
if (lowerBound != null && lowerBound.Count != parameters.Count)
{
throw new ArgumentException("The lower bounds can't have different size from the parameters.");
}
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 && upperBound.Count != parameters.Count)
{
throw new ArgumentException("The upper bounds can't have different size from the parameetrs.");
}
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 && scales.Count != parameters.Count)
{
throw new ArgumentException("The scales can't have different size from the parameters.");
}
if (scales != null && scales.Count(x => x < 0) > 0)
{
scales.PointwiseAbs();
}
Scales = scales;
}
protected static double EvaluateFunction(IObjectiveModel objective, Vector<double> Pint)
{
var Pext = ProjectToExternalParameters(Pint);
objective.EvaluateAt(Pext);
return objective.Value;
}
protected static Tuple<Vector<double>, Matrix<double>> EvaluateJacobian(IObjectiveModel objective, Vector<double> Pint)
{
var gradient = objective.Gradient;
var hessian = objective.Hessian;
if (IsBounded)
{
var scaleFactors = ScaleFactorsOfJacobian(Pint); // the parameters argument is always internal.
for (int i = 0; i < gradient.Count; i++)
{
gradient[i] = gradient[i] * scaleFactors[i];
}
for (int i = 0; i < hessian.RowCount; i++)
{
for (int j = 0; j < hessian.ColumnCount; j++)
{
hessian[i, j] = hessian[i, j] * scaleFactors[i] * scaleFactors[j];
}
}
}
return new Tuple<Vector<double>, Matrix<double>>(gradient, hessian);
}
#region Projection of Parameters
// To handle the box constrained minimization as the unconstrained minimization,
// the parameters are mapping by the following rules,
// which are modified the rules shown in the ref[1] in order to introduce scales.
//
// 1. lower < Pext < 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 < Pext
// Pint = sqrt((Pext/scale - lower/scale + 1)^2 - 1)
// Pext = lower + scale * (sqrt(Pint^2 + 1) - 1)
// dPext/dPint = scale * Pint / sqrt(Pint^2 + 1)
//
// 3. Pext < upper
// Pint = sqrt((upper / scale - Pext / scale + 1)^2 - 1)
// Pext = upper + scale - scale * sqrt(Pint^2 + 1)
// dPext/dPint = - scale * Pint / sqrt(Pint^2 + 1)
//
// 4. no bounds, but scales
// Pint = Pext / scale
// Pext = Pint * scale
// dPext/dPint = scale
//
// The rules are applied in ProjectParametersToInternal, ProjectParametersToExternal, and ScaleFactorsOfJacobian methods.
//
// References:
// [1] https://lmfit.github.io/lmfit-py/bounds.html
// [2] MINUIT User's Guide, https://root.cern.ch/download/minuit.pdf
//
// 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.
protected static Vector<double> ProjectToInternalParameters(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] = (Scales == null)
? Math.Sqrt(Math.Pow(Pext[i] - LowerBound[i] + 1.0, 2) - 1.0)
: Math.Sqrt(Math.Pow((Pext[i] - LowerBound[i]) / Scales[i] + 1.0, 2) - 1.0);
}
return Pint;
}
else if (LowerBound == null && UpperBound != null)
{
for (int i = 0; i < Pext.Count; i++)
{
Pint[i] = (Scales == null)
? Math.Sqrt(Math.Pow(UpperBound[i] - Pext[i] + 1.0, 2) - 1.0)
: Math.Sqrt(Math.Pow((UpperBound[i] - Pext[i]) / Scales[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;
}
protected static Vector<double> ProjectToExternalParameters(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] = (Scales == null)
? LowerBound[i] + Math.Sqrt(Pint[i] * Pint[i] + 1.0) - 1.0
: LowerBound[i] + Scales[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] = (Scales == null)
? UpperBound[i] - Math.Sqrt(Pint[i] * Pint[i] + 1.0) + 1.0
: UpperBound[i] - Scales[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;
}
protected static 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] = (Scales == null)
? Pint[i] / Math.Sqrt(Pint[i] * Pint[i] + 1.0)
: Scales[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] = (Scales == null)
? -Pint[i] / Math.Sqrt(Pint[i] * Pint[i] + 1.0)
: -Scales[i] * Pint[i] / Math.Sqrt(Pint[i] * Pint[i] + 1.0);
}
return scale;
}
else if (Scales != null)
{
return Scales;
}
return scale;
}
#endregion Projection of Parameters
}
}

106
src/Numerics/Optimization/ObjectiveFunction.cs
View File

@ -114,5 +114,111 @@ namespace MathNet.Numerics.Optimization
{
return new ScalarObjectiveFunction(function, derivative, secondDerivative);
}
/// <summary>
/// objective model with a user supplied jacobian for non-linear least squares regression.
/// </summary>
public static IObjectiveModel NonlinearModel(Func<Vector<double>, Vector<double>, Vector<double>> function,
Func<Vector<double>, Vector<double>, Matrix<double>> derivatives,
Vector<double> observedX, Vector<double> observedY, Vector<double> weight = null)
{
var objective = new NonlinearObjectiveFunction(function, derivatives);
objective.SetObserved(observedX, observedY, weight);
return objective;
}
/// <summary>
/// Objective model for non-linear least squares regression.
/// </summary>
public static IObjectiveModel NonlinearModel(Func<Vector<double>, Vector<double>, Vector<double>> function,
Vector<double> observedX, Vector<double> observedY, Vector<double> weight = null,
int accuracyOrder = 2)
{
var objective = new NonlinearObjectiveFunction(function, accuracyOrder: accuracyOrder);
objective.SetObserved(observedX, observedY, weight);
return objective;
}
/// <summary>
/// Objective model with a user supplied jacobian for non-linear least squares regression.
/// </summary>
public static IObjectiveModel NonlinearModel(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> func(Vector<double> point, Vector<double> x)
{
var functionValues = CreateVector.Dense<double>(x.Count);
for (int i = 0; i < x.Count; i++)
{
functionValues[i] = function(point, x[i]);
}
return functionValues;
}
Matrix<double> prime(Vector<double> point, Vector<double> x)
{
var derivativeValues = CreateMatrix.Dense<double>(x.Count, point.Count);
for (int i = 0; i < x.Count; i++)
{
derivativeValues.SetRow(i, derivatives(point, x[i]));
}
return derivativeValues;
}
var objective = new NonlinearObjectiveFunction(func, prime);
objective.SetObserved(observedX, observedY, weight);
return objective;
}
/// <summary>
/// Objective model for non-linear least squares regression.
/// </summary>
public static IObjectiveModel NonlinearModel(Func<Vector<double>, double, double> function,
Vector<double> observedX, Vector<double> observedY, Vector<double> weight = null,
int accuracyOrder = 2)
{
Vector<double> func(Vector<double> point, Vector<double> x)
{
var functionValues = CreateVector.Dense<double>(x.Count);
for (int i = 0; i < x.Count; i++)
{
functionValues[i] = function(point, x[i]);
}
return functionValues;
}
var objective = new NonlinearObjectiveFunction(func, accuracyOrder: accuracyOrder);
objective.SetObserved(observedX, observedY, weight);
return objective;
}
/// <summary>
/// Objective function with a user supplied jacobian for nonlinear least squares regression.
/// </summary>
public static IObjectiveFunction NonlinearFunction(Func<Vector<double>, Vector<double>, Vector<double>> function,
Func<Vector<double>, Vector<double>, Matrix<double>> derivatives,
Vector<double> observedX, Vector<double> observedY, Vector<double> weight = null)
{
var objective = new NonlinearObjectiveFunction(function, derivatives);
objective.SetObserved(observedX, observedY, weight);
return objective.ToObjectiveFunction();
}
/// <summary>
/// Objective function for nonlinear least squares regression.
/// The numerical jacobian with accuracy order is used.
/// </summary>
public static IObjectiveFunction NonlinearFunction(Func<Vector<double>, Vector<double>, Vector<double>> function,
Vector<double> observedX, Vector<double> observedY, Vector<double> weight = null,
int accuracyOrder = 2)
{
var objective = new NonlinearObjectiveFunction(function, null, accuracyOrder: accuracyOrder);
objective.SetObserved(observedX, observedY, weight);
return objective.ToObjectiveFunction();
}
}
}

436
src/Numerics/Optimization/ObjectiveFunctions/NonlinearObjectiveFunction.cs
View File

@ -0,0 +1,436 @@
using MathNet.Numerics.LinearAlgebra;
using System;
using System.Collections.Generic;
using System.Linq;
namespace MathNet.Numerics.Optimization.ObjectiveFunctions
{
internal class NonlinearObjectiveFunction : IObjectiveModel
{
#region Private Variables
readonly Func<Vector<double>, Vector<double>, Vector<double>> userFunction; // (p, x) => f(x; p)
readonly Func<Vector<double>, Vector<double>, Matrix<double>> userDerivative; // (p, x) => df(x; p)/dp
readonly int accuracyOrder; // the desired accuracy order to evaluate the jacobian by numerical approximaiton.
Vector<double> coefficients;
bool hasFunctionValue;
double functionValue; // the residual sum of squares, residuals * residuals.
Vector<double> residuals; // the weighted error values
bool hasJacobianValue;
Matrix<double> jacobianValue; // the Jacobian matrix.
Vector<double> gradientValue; // the Gradient vector.
Matrix<double> hessianValue; // the Hessian matrix.
#endregion Private Variables
#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.
/// </summary>
public Matrix<double> Weights { get; private set; }
private Vector<double> L; // Weights = LL'
/// <summary>
/// Get whether parameters are fixed or free.
/// </summary>
public List<bool> IsFixed { 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 (Point == null) ? 0 : Point.Count; } }
/// <summary>
/// Get the degree of freedom
/// </summary>
public int DegreeOfFreedom
{
get
{
var df = NumberOfObservations - NumberOfParameters;
if (IsFixed != null)
{
df = df + IsFixed.Count(p => p == true);
}
return df;
}
}
/// <summary>
/// Get the number of calls to function.
/// </summary>
public int FunctionEvaluations { get; set; }
/// <summary>
/// Get the number of calls to jacobian.
/// </summary>
public int JacobianEvaluations { get; set; }
#endregion Public Variables
public NonlinearObjectiveFunction(Func<Vector<double>, Vector<double>, Vector<double>> function,
Func<Vector<double>, Vector<double>, Matrix<double>> derivative = null, int accuracyOrder = 2)
{
this.userFunction = function;
this.userDerivative = derivative;
this.accuracyOrder = Math.Min(6, Math.Max(1, accuracyOrder));
}
public IObjectiveModel Fork()
{
return new NonlinearObjectiveFunction(userFunction, userDerivative, accuracyOrder)
{
ObservedX = ObservedX,
ObservedY = ObservedY,
Weights = Weights,
coefficients = coefficients,
hasFunctionValue = hasFunctionValue,
functionValue = functionValue,
hasJacobianValue = hasJacobianValue,
jacobianValue = jacobianValue,
gradientValue = gradientValue,
hessianValue = hessianValue
};
}
public IObjectiveModel CreateNew()
{
return new NonlinearObjectiveFunction(userFunction, userDerivative, accuracyOrder);
}
/// <summary>
/// Set or get the values of the parameters.
/// </summary>
public Vector<double> Point { get { return coefficients; } }
/// <summary>
/// Get the y-values of the fitted model that correspond to the independent values.
/// </summary>
public Vector<double> ModelValues { get; private set; }
/// <summary>
/// Get the residual sum of squares.
/// </summary>
public double Value
{
get
{
if (!hasFunctionValue)
{
EvaluateFunction();
hasFunctionValue = true;
}
return functionValue;
}
}
/// <summary>
/// Get the Gradient vector of x and p.
/// </summary>
public Vector<double> Gradient
{
get
{
if (!hasJacobianValue)
{
EvaluateJacobian();
hasJacobianValue = true;
}
return gradientValue;
}
}
/// <summary>
/// Get the Hessian matrix of x and p, J'WJ
/// </summary>
public Matrix<double> Hessian
{
get
{
if (!hasJacobianValue)
{
EvaluateJacobian();
hasJacobianValue = true;
}
return hessianValue;
}
}
public bool IsGradientSupported { get { return true; } }
public bool IsHessianSupported { get { return true; } }
/// <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 parameters and bounds.
/// </summary>
/// <param name="initialGuess">The initial values of parameters.</param>
/// <param name="isFixed">The list to the parameters fix or free.</param>
public void SetParameters(Vector<double> initialGuess, List<bool> isFixed = null)
{
if (initialGuess == null)
{
throw new ArgumentNullException("initialGuess");
}
coefficients = initialGuess;
if (isFixed != null && isFixed.Count != initialGuess.Count)
{
throw new ArgumentException("The isFixed can't have different size from the initial guess.");
}
if (isFixed != null && isFixed.Count(p => p == true) == isFixed.Count)
{
throw new ArgumentException("All the parameters can't be fixed.");
}
IsFixed = isFixed;
}
public void EvaluateAt(Vector<double> parameters)
{
if (parameters == null)
{
throw new ArgumentNullException("parameters");
}
if (parameters.Count(p => double.IsNaN(p) || double.IsInfinity(p)) > 0)
{
throw new ArgumentException("The parameters must be finite.");
}
coefficients = parameters;
hasFunctionValue = false;
hasJacobianValue = false;
jacobianValue = null;
gradientValue = null;
hessianValue = null;
}
public IObjectiveFunction ToObjectiveFunction()
{
Tuple<double, Vector<double>, Matrix<double>> function(Vector<double> point)
{
EvaluateAt(point);
return new Tuple<double, Vector<double>, Matrix<double>>(Value, Gradient, Hessian);
}
var objective = new GradientHessianObjectiveFunction(function);
return objective;
}
#region Private Methods
private void EvaluateFunction()
{
// Calculates the residuals, (y[i] - f(x[i]; p)) * L[i]
if (ModelValues == null)
{
ModelValues = Vector<double>.Build.Dense(NumberOfObservations);
}
ModelValues = userFunction(Point, ObservedX);
FunctionEvaluations++;
// calculate the weighted residuals
residuals = (Weights == null)
? ObservedY - ModelValues
: (ObservedY - ModelValues).PointwiseMultiply(L);
// Calculate the residual sum of squares
functionValue = residuals.DotProduct(residuals);
return;
}
private void EvaluateJacobian()
{
// Calculates the jacobian of x and p.
if (userDerivative != null)
{
// analytical jacobian
jacobianValue = userDerivative(Point, ObservedX);
JacobianEvaluations++;
}
else
{
// numerical jacobian
jacobianValue = NumericalJacobian(Point, ModelValues, accuracyOrder);
FunctionEvaluations += accuracyOrder;
}
// weighted jacobian
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
jacobianValue[i, j] = 0.0;
}
else
{
jacobianValue[i, j] = (Weights == null)
? jacobianValue[i, j]
: jacobianValue[i, j] * L[j];
}
}
}
// Gradient, g = -J'W(y − f(x; p)) = -J'L(L'E) = -J'LR
gradientValue = -jacobianValue.Transpose() * residuals;
// approximated Hessian, H = J'WJ + ∑LRiHi ~ J'WJ near the minimum
hessianValue = jacobianValue.Transpose() * jacobianValue;
}
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 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, ObservedX);
var f2 = userFunction(parameters - 2 * h, ObservedX);
var f3 = userFunction(parameters - h, ObservedX);
var f4 = userFunction(parameters + h, ObservedX);
var f5 = userFunction(parameters + 2 * h, ObservedX);
var f6 = userFunction(parameters + 3 * h, ObservedX);
var prime = (-f1 + 9 * f2 - 45 * f3 + 45 * f4 - 9 * f5 + f6) / (60 * h[j]);
derivertives.SetColumn(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;
var f2 = userFunction(parameters + h, ObservedX);
var f3 = userFunction(parameters + 2 * h, ObservedX);
var f4 = userFunction(parameters + 3 * h, ObservedX);
var f5 = userFunction(parameters + 4 * h, ObservedX);
var f6 = userFunction(parameters + 5 * h, ObservedX);
var prime = (-137 * f1 + 300 * f2 - 300 * f3 + 200 * f4 - 75 * f5 + 12 * f6) / (60 * h[j]);
derivertives.SetColumn(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, ObservedX);
var f2 = userFunction(parameters - h, ObservedX);
var f3 = userFunction(parameters + h, ObservedX);
var f4 = userFunction(parameters + 2 * h, ObservedX);
var prime = (f1 - 8 * f2 + 8 * f3 - f4) / (12 * h[j]);
derivertives.SetColumn(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;
var f2 = userFunction(parameters + h, ObservedX);
var f3 = userFunction(parameters + 2 * h, ObservedX);
var f4 = userFunction(parameters + 3 * h, ObservedX);
var prime = (-11 * f1 + 18 * f2 - 9 * f3 + 2 * f4) / (6 * h[j]);
derivertives.SetColumn(j, prime);
}
else if (accuracyOrder == 2)
{
// f'(x) = {f(x + h) - f(x - h)} / 2h + O(h^2)
var f1 = userFunction(parameters + h, ObservedX);
var f2 = userFunction(parameters - h, ObservedX);
var prime = (f1 - f2) / (2 * h[j]);
derivertives.SetColumn(j, prime);
}
else
{
// f'(x) = {- f(x) + f(x + h)} / h + O(h)
var f1 = currentValues;
var f2 = userFunction(parameters + h, ObservedX);
var prime = (-f1 + f2) / h[j];
derivertives.SetColumn(j, prime);
}
h[j] = 0;
}
return derivertives;
}
#endregion Private Methods
}
}

45
src/Numerics/Optimization/Subproblems/DogLegSubproblem.cs
View File

@ -0,0 +1,45 @@
using MathNet.Numerics.LinearAlgebra;
namespace MathNet.Numerics.Optimization.Subproblems
{
internal class DogLegSubproblem : ITrustRegionSubproblem
{
public Vector<double> Pstep { get; private set; }
public bool HitBoundary { get; private set; }
public void Solve(IObjectiveModel objective, double delta)
{
var Gradient = objective.Gradient;
var Hessian = objective.Hessian;
// 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;
}
else if (alpha * Psd.L2Norm() >= delta)
{
// Psd is outside trust region radius
HitBoundary = true;
Pstep = delta / Psd.L2Norm() * Psd;
}
else
{
// Pstep is intersection of the trust region boundary
HitBoundary = true;
var beta = Util.FindBeta(alpha, Psd, Pgn, delta).Item2;
Pstep = alpha * Psd + beta * (Pgn - alpha * Psd);
}
}
}
}

65
src/Numerics/Optimization/Subproblems/NewtonCGSubproblem.cs
View File

@ -0,0 +1,65 @@
using MathNet.Numerics.LinearAlgebra;
using System;
namespace MathNet.Numerics.Optimization.Subproblems
{
internal class NewtonCGSubproblem : ITrustRegionSubproblem
{
public Vector<double> Pstep { get; private set; }
public bool HitBoundary { get; private set; }
public void Solve(IObjectiveModel objective, double delta)
{
var Gradient = objective.Gradient;
var Hessian = objective.Hessian;
// define tolerance
var gnorm = Gradient.L2Norm();
var tolerance = Math.Min(0.5, Math.Sqrt(gnorm)) * gnorm;
// initialize internal variables
var z = Vector<double>.Build.Dense(Hessian.RowCount);
var r = Gradient;
var d = -r;
while (true)
{
var Bd = Hessian * d;
var dBd = d.DotProduct(Bd);
if (dBd <= 0)
{
var t = Util.FindBeta(1, z, d, delta);
Pstep = z + t.Item1 * d;
HitBoundary = true;
return;
}
var r_sq = r.DotProduct(r);
var alpha = r_sq / dBd;
var znext = z + alpha * d;
if(znext.L2Norm() >= delta)
{
var t = Util.FindBeta(1, z, d, delta);
Pstep = z + t.Item2 * d;
HitBoundary = true;
return;
}
var rnext = r + alpha * Bd;
var rnext_sq = rnext.DotProduct(rnext);
if (Math.Sqrt(rnext_sq) < tolerance)
{
Pstep = znext;
HitBoundary = false;
return;
}
z = znext;
r = rnext;
d = -rnext + rnext_sq / r_sq * d;
}
}
}
}

35
src/Numerics/Optimization/Subproblems/Util.cs
View File

@ -0,0 +1,35 @@
using MathNet.Numerics.LinearAlgebra;
using System;
namespace MathNet.Numerics.Optimization.Subproblems
{
internal static class Util
{
public static Tuple<double, 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
//
// positive β is used for 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 beta1 = -aux / 2.0 / a;
var beta2 = -2.0 * c / aux;
// return sorted beta
return (beta1 < beta2)
? new Tuple<double, double>(beta1, beta2)
: new Tuple<double, double>(beta2, beta1);
}
}
}

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

@ -0,0 +1,12 @@
namespace MathNet.Numerics.Optimization
{
public sealed class TrustRegionDogLegMinimizer : TrustRegionMinimizerBase
{
/// <summary>
/// Non-linear least square fitting by the trust region dogleg algorithm.
/// </summary>
public TrustRegionDogLegMinimizer(double gradientTolerance = 1E-8, double stepTolerance = 1E-8, double functionTolerance = 1E-8, double radiusTolerance = 1E-8, int maximumIterations = -1)
: base(TrustRegionSubproblem.DogLeg(), gradientTolerance, stepTolerance, functionTolerance, radiusTolerance, maximumIterations)
{ }
}
}

245
src/Numerics/Optimization/TrustRegionMinimizerBase.cs
View File

@ -0,0 +1,245 @@
using MathNet.Numerics.LinearAlgebra;
using System;
using System.Collections.Generic;
using System.Linq;
namespace MathNet.Numerics.Optimization
{
public abstract class TrustRegionMinimizerBase : NonlinearMinimizerBase
{
/// <summary>
/// The trust region subproblem.
/// </summary>
public static ITrustRegionSubproblem Subproblem;
/// <summary>
/// The stopping threshold for the trust region radius.
/// </summary>
public static double RadiusTolerance { get; set; }
public TrustRegionMinimizerBase(ITrustRegionSubproblem subproblem,
double gradientTolerance = 1E-8, double stepTolerance = 1E-8, double functionTolerance = 1E-8, double radiusTolerance = 1E-8, int maximumIterations = -1)
: base(gradientTolerance, stepTolerance, functionTolerance, maximumIterations)
{
if (subproblem == null)
throw new ArgumentNullException("subproblem");
Subproblem = subproblem;
RadiusTolerance = radiusTolerance;
}
public NonlinearMinimizationResult FindMinimum(IObjectiveModel objective, Vector<double> initialGuess,
Vector<double> lowerBound = null, Vector<double> upperBound = null, Vector<double> scales = null, List<bool> isFixed = null)
{
return Minimum(Subproblem, objective, initialGuess, lowerBound, upperBound, scales, isFixed,
GradientTolerance, StepTolerance, FunctionTolerance, RadiusTolerance, MaximumIterations);
}
public NonlinearMinimizationResult FindMinimum(IObjectiveModel objective, double[] initialGuess,
double[] lowerBound = null, double[] upperBound = null, double[] scales = null, bool[] isFixed = null)
{
var lb = (lowerBound == null) ? null : CreateVector.Dense<double>(lowerBound);
var ub = (upperBound == null) ? null : CreateVector.Dense<double>(upperBound);
var sc = (scales == null) ? null : CreateVector.Dense<double>(scales);
var fx = (isFixed == null) ? null : isFixed.ToList();
return Minimum(Subproblem, objective, CreateVector.DenseOfArray<double>(initialGuess), lb, ub, sc, fx,
GradientTolerance, StepTolerance, FunctionTolerance, RadiusTolerance, MaximumIterations);
}
/// <summary>
/// Non-linear least square fitting by the trust-region algorithm.
/// </summary>
/// <param name="objective">The objective model, including function, jacobian, observations, and parameter bounds.</param>
/// <param name="subproblem">The subproblem</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 NonlinearMinimizationResult Minimum(ITrustRegionSubproblem subproblem, IObjectiveModel objective, Vector<double> initialGuess,
Vector<double> lowerBound = null, Vector<double> upperBound = null, Vector<double> scales = null, List<bool> isFixed = null,
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 algorithm.
//
// 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 trust region algorithm is summarized as follows:
// initially set trust-region radius, Δ
// repeat
// solve 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]. Nocedal, Jorge, and Stephen J. Wright.
// Numerical optimization (2006): 101-134.
// [3]. SciPy
// Available Online from: https://github.com/scipy/scipy/blob/master/scipy/optimize/_trustregion.py
double maxDelta = 1000;
double eta = 0;
if (objective == null)
throw new ArgumentNullException("objective");
ValidateBounds(initialGuess, lowerBound, upperBound, scales);
objective.SetParameters(initialGuess, isFixed);
ExitCondition exitCondition = ExitCondition.None;
// First, calculate function values and setup variables
var P = ProjectToInternalParameters(initialGuess); // current internal parameters
var Pstep = Vector<double>.Build.Dense(P.Count); // the change of parameters
var RSS = EvaluateFunction(objective, initialGuess); // Residual Sum of Squares
if (maximumIterations < 0)
{
maximumIterations = 200 * (initialGuess.Count + 1);
}
// if RSS == NaN, stop
if (double.IsNaN(RSS))
{
exitCondition = ExitCondition.InvalidValues;
return new NonlinearMinimizationResult(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 gradient and Hessian
var jac = EvaluateJacobian(objective, P);
var Gradient = jac.Item1; // objective.Gradient;
var Hessian = jac.Item2; // objective.Hessian;
// if ||g||_oo <= gtol, found and stop
if (Gradient.InfinityNorm() <= gradientTolerance)
{
exitCondition = ExitCondition.RelativeGradient; // SmallGradient
}
if (exitCondition != ExitCondition.None)
{
return new NonlinearMinimizationResult(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;
bool hitBoundary = false;
while (iterations < maximumIterations && exitCondition == ExitCondition.None)
{
iterations++;
// solve the subproblem
subproblem.Solve(objective, delta);
Pstep = subproblem.Pstep;
hitBoundary = subproblem.HitBoundary;
// predicted reduction = L(0) - L(Δp) = -Δp'g - 1/2 * Δp'HΔp
var predictedReduction = -Gradient.DotProduct(Pstep) - 0.5 * Pstep.DotProduct(Hessian * Pstep);
if (Pstep.L2Norm() <= stepTolerance * (stepTolerance + P.L2Norm()))
{
exitCondition = ExitCondition.RelativePoints; // SmallRelativeParameters
break;
}
var Pnew = P + Pstep; // parameters to test
// evaluate function at Pnew
var RSSnew = EvaluateFunction(objective, Pnew);
// if RSS == NaN, stop
if (double.IsNaN(RSSnew))
{
exitCondition = ExitCondition.InvalidValues;
break;
}
// calculate the ratio of the actual to the predicted reduction.
double rho = (predictedReduction != 0)
? (RSS - RSSnew) / predictedReduction
: 0.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.LackOfProgress;
break;
}
}
if (rho > eta)
{
// accepted
Pnew.CopyTo(P);
RSS = RSSnew;
// evaluate projected gradient and Hessian
jac = EvaluateJacobian(objective, P);
Gradient = jac.Item1; // objective.Gradient;
Hessian = jac.Item2; // 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;
}
return new NonlinearMinimizationResult(objective, iterations, exitCondition);
}
}
}

13
src/Numerics/Optimization/TrustRegionNewtonCGMinimizer.cs
View File

@ -0,0 +1,13 @@
namespace MathNet.Numerics.Optimization
{
public sealed class TrustRegionNewtonCGMinimizer : TrustRegionMinimizerBase
{
/// <summary>
/// Non-linear least square fitting by the trust region Newton-Conjugate-Gradient algorithm.
/// </summary>
public TrustRegionNewtonCGMinimizer(double gradientTolerance = 1E-8, double stepTolerance = 1E-8, double functionTolerance = 1E-8, double radiusTolerance = 1E-8, int maximumIterations = -1)
: base(TrustRegionSubproblem.NewtonCG(), gradientTolerance, stepTolerance, functionTolerance, radiusTolerance, maximumIterations)
{ }
}
}

17
src/Numerics/Optimization/TrustRegionSubProblem.cs
View File

@ -0,0 +1,17 @@
using MathNet.Numerics.Optimization.Subproblems;
namespace MathNet.Numerics.Optimization
{
public static class TrustRegionSubproblem
{
public static ITrustRegionSubproblem DogLeg()
{
return new DogLegSubproblem();
}
public static ITrustRegionSubproblem NewtonCG()
{
return new NewtonCGSubproblem();
}
}
}
Write Preview
Loading…
Cancel
Save