diff --git a/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs b/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs new file mode 100644 index 00000000..fbbd3448 --- /dev/null +++ b/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs @@ -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 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 RosenbrockPrime(Vector p, double x) + { + var prime = Vector.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 RosenbrockX = Vector.Build.Dense(2); + private Vector RosenbrockY = Vector.Build.Dense(2); + private Vector RosenbrockPbest = new DenseVector(new double[] { 1.0, 1.0 }); + + private Vector RosenbrockStart1 = new DenseVector(new double[] { -1.2, 1.0 }); + private Vector RosebbrockLowerBound = new DenseVector(new double[] { -5.0, -5.0 }); + private Vector 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 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 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 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 Rat43Pbest = new DenseVector(new double[] { + 6.9964151270E+02, 5.2771253025E+00, 7.5962938329E-01, 1.2792483859E+00 + }); + private Vector Rat43Pstd = new DenseVector(new double[]{ + 1.6302297817E+01, 2.0828735829E+00, 1.9566123451E-01, 6.8761936385E-01 + }); + + private Vector Rat43Start1 = new DenseVector(new double[] { 100, 10, 1, 1 }); + private Vector 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 p, double x) + { + var y = p[0] * (1.0 - Math.Exp(-p[1] * x)); + return y; + } + private Vector BoxBodPrime(Vector p, double x) + { + var prime = Vector.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 BoxBodX = new DenseVector(new double[] { 1, 2, 3, 5, 7, 10 }); + private Vector BoxBodY = new DenseVector(new double[] { 109, 149, 149, 191, 213, 224 }); + private Vector BoxBodPbest = new DenseVector(new double[] { 2.1380940889E+02, 5.4723748542E-01 }); + private Vector BoxBodPstd = new DenseVector(new double[] { 1.2354515176E+01, 1.0455993237E-01 }); + + private Vector BoxBodStart1 = new DenseVector(new double[] { 1.0, 1.0 }); + private Vector BoxBodStart2 = new DenseVector(new double[] { 100.0, 0.75 }); + private Vector BoxBodLowerBound = new DenseVector(new double[] { 0, 0 }); + private Vector BoxBodUpperBound = new DenseVector(new double[] { 500.0, 10 }); + private Vector 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 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 ThurberPrime(Vector p, double x) + { + var prime = Vector.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 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 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 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 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 ThurberInitialGuess = new DenseVector(new double[] { 1000.0, 1000.0, 400.0, 40.0, 0.7, 0.3, 0.03 }); + private Vector 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); + } + } + } +} diff --git a/src/Numerics/Optimization/ExitCondition.cs b/src/Numerics/Optimization/ExitCondition.cs index 90f30714..dc83feb4 100644 --- a/src/Numerics/Optimization/ExitCondition.cs +++ b/src/Numerics/Optimization/ExitCondition.cs @@ -32,12 +32,16 @@ namespace MathNet.Numerics.Optimization public enum ExitCondition { None, + InvalidValues, + ExceedIterations, + RelativePoints, RelativeGradient, LackOfProgress, AbsoluteGradient, WeakWolfeCriteria, BoundTolerance, StrongWolfeCriteria, - Converged + Converged, + ManuallyStopped } } diff --git a/src/Numerics/Optimization/IObjectiveModel.cs b/src/Numerics/Optimization/IObjectiveModel.cs new file mode 100644 index 00000000..45f94409 --- /dev/null +++ b/src/Numerics/Optimization/IObjectiveModel.cs @@ -0,0 +1,70 @@ +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization +{ + public interface IObjectiveModelEvaluation + { + IObjectiveModel CreateNew(); + + /// + /// Get the y-values of the fitted model that correspond to the independent values. + /// + Vector Values { get; } + + /// + /// Get the values of the parameters. + /// + Vector Parameters { get; } + + /// + /// Get the residual sum of squares. + /// + double Residue { get; } + + /// + /// Get the Jacobian matrix, J(x; p) = df(x; p)/dp. + /// + Matrix Jacobian { get; } + /// + /// Get the Gradient vector. G = J'(y - f(x; p)) + /// + Vector Gradient { get; } + /// + /// Get the approximated Hessian matrix. H = J'J + /// + Matrix Hessian { get; } + /// + /// Get the covariance matrix. + /// + Matrix Covariance { get; } + + /// + /// Get the number of calls to function. + /// + int FunctionEvaluations { get; } + /// + /// Get the number of calls to jacobian. + /// + int JacobianEvaluations { get; } + + /// + /// Get the degree of freedom. + /// + int DegreeOfFreedom { get; } + + /// + /// Get whether or not the analytical jacobian is supported. + /// + bool IsJacobianSupported { get; } + } + + public interface IObjectiveModel : IObjectiveModelEvaluation + { + void EvaluateFunction(Vector parameters); + void EvaluateJacobian(Vector parameters); + void EvaluateCovariance(Vector parameters); + + /// Create a new independent copy of this objective function, evaluated at the same point. + IObjectiveModel Fork(); + } +} diff --git a/src/Numerics/Optimization/LevenbergMarquardtMinimizer.cs b/src/Numerics/Optimization/LevenbergMarquardtMinimizer.cs new file mode 100644 index 00000000..681eb3dd --- /dev/null +++ b/src/Numerics/Optimization/LevenbergMarquardtMinimizer.cs @@ -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 + + /// + /// The scale factor for initial mu + /// + public static double InitialMu { get; set; } + + /// + /// The stopping threshold for infinity norm of the gradient. + /// + public static double GradientTolerance { get; set; } + + /// + /// The stopping threshold for L2 norm of the change of the parameters. + /// + public static double StepTolerance { get; set; } + + /// + /// The stopping threshold for the function value or L2 norm of the residuals. + /// + public static double FunctionTolerance { get; set; } + + /// + /// The maximum number of iterations. + /// + 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 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(initialGuess), InitialMu, GradientTolerance, StepTolerance, FunctionTolerance, MaximumIterations); + } + + /// + /// Non-linear least square fitting by the Levenberg-Marduardt algorithm. + /// + /// The objective function, including model, observations, and parameter bounds. + /// The initial guess values. + /// The initial damping parameter of mu. + /// The stopping threshold for infinity norm of the gradient vector. + /// The stopping threshold for L2 norm of the change of parameters. + /// The stopping threshold for L2 norm of the residuals. + /// The max iterations. + /// The result of the Levenberg-Marquardt minimization + public static ModelMinimizationResult Minimum(IObjectiveModel objective, Vector 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.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); + } + } +} diff --git a/src/Numerics/Optimization/ModelMinimizationResult.cs b/src/Numerics/Optimization/ModelMinimizationResult.cs new file mode 100644 index 00000000..baaf3710 --- /dev/null +++ b/src/Numerics/Optimization/ModelMinimizationResult.cs @@ -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; } + + /// + /// Returns the best fit parameters. + /// + public Vector BestFitParameters { get { return ModelInfoAtMinimum.Parameters; } } + + /// + /// Returns the standard errors of the corresponding parameters + /// + public Vector StandardErrors + { + get + { + if (ModelInfoAtMinimum.Covariance == null) + return null; + return ModelInfoAtMinimum.Covariance.Diagonal().PointwiseSqrt(); + } + } + + /// + /// Returns the y-values of the fitted model that correspond to the independent values. + /// + public Vector BestFitValues { get { return ModelInfoAtMinimum.Values; } } + + /// + /// Returns the residual sum of squares. + /// + 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; + } + } +} diff --git a/src/Numerics/Optimization/ObjectiveModel.cs b/src/Numerics/Optimization/ObjectiveModel.cs new file mode 100644 index 00000000..631defae --- /dev/null +++ b/src/Numerics/Optimization/ObjectiveModel.cs @@ -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 + { + /// + /// Fitting model with a user supplied jacobian for non-linear least squares regression. + /// + public static IObjectiveModel FittingModel(Func, double, double> function, Func, double, Vector> derivatives, + Vector observedX, Vector observedY, Vector weight = null, + Vector lowerBound = null, Vector upperBound = null, Vector scales = null, List isFixed = null) + { + var objective = new FittingObjectiveModel(function, derivatives); + objective.SetObserved(observedX, observedY, weight); + objective.SetParameters(lowerBound, upperBound, scales, isFixed); + return objective; + } + + /// + /// Fitting model for non-linear least squares regression. + /// The numerical jacobian with accuracy order is used. + /// + public static IObjectiveModel FittingModel(Func, double, double> function, + Vector observedX, Vector observedY, Vector weight = null, + Vector lowerBound = null, Vector upperBound = null, Vector scales = null, List 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; + } + } +} diff --git a/src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs b/src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs new file mode 100644 index 00000000..a1ea2beb --- /dev/null +++ b/src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs @@ -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, double, double> userFunction; // (p, x) => f(x; p) + readonly Func, double, Vector> userDerivatives; // (p, x) => df(x; p)/dp + + #region Public Variables + + /// + /// Set or get the values of the independent variable. + /// + public Vector ObservedX { get; private set; } + + /// + /// Set or get the values of the observations. + /// + public Vector ObservedY { get; private set; } + + /// + /// Set or get the values of the weights for the observations. + /// inverse of the standard measurement errors + /// If null, unity weighting is used. + /// + public Matrix Weights { get; private set; } + // W = LL' + private Vector L; + + /// + /// Set or get the values of the parameters. + /// + public Vector Parameters { get; private set; } + + /// + /// Set or get the values of the parameters. + /// + public List IsFixed { get; set; } + + /// + /// Set or get the values of the parameters. + /// + public Vector LowerBound { get; set; } + + /// + /// Set or get the values of the parameters. + /// + public Vector UpperBound { get; set; } + + /// + /// Set or get the scale factor of the parameters. + /// + public Vector Scales { get; set; } + + /// + /// Set of get whether or not the parameters are bounded. + /// + public bool IsBounded { get; set; } + + /// + /// Get the y-values of the fitted model that correspond to the independent values. + /// + public Vector Values { get; private set; } + + /// + /// Get the error values, R(x; p) = L * (y - f(x; p)) where L = sqrt(W) + /// + private Vector Residuals; + + /// + /// Get the residual sum of squares, R.DotProduct(R) + /// + public double Residue { get; private set; } + + /// + /// Get the Jacobian matrix of x and p, J(x; p). + /// + public Matrix Jacobian { get; private set; } + + /// + /// Get the Gradient vector of x and p, J'WR + /// + public Vector Gradient { get; private set; } + + /// + /// Get the Hessian matrix of x and p, J'WJ + /// + public Matrix Hessian { get; private set; } + + /// + /// Get the number of observations. + /// + public int NumberOfObservations { get { return (ObservedY == null) ? 0 : ObservedY.Count; } } + + /// + /// Get the number of unknown parameters. + /// + public int NumberOfParameters { get { return (Parameters == null) ? 0 : Parameters.Count; } } + + /// + /// Get the degree of freedom + /// + public int DegreeOfFreedom + { + get + { + var dof = NumberOfObservations - NumberOfParameters; + if (IsFixed != null) + { + dof = dof + IsFixed.Count(p => p == true); + } + return dof; + } + } + + /// + /// Get the covariance matrix. + /// + public Matrix Covariance { get; private set; } + + /// + /// Get the number of calls to function. + /// + public int FunctionEvaluations { get; private set; } + /// + /// Get the number of calls to jacobian. + /// + public int JacobianEvaluations { get; private set; } + + /// + /// Set or get the desired accuracy order of the numerical jacobian. + /// + public int AccuracyOrder { get; set; } + + /// + /// Get whether or not the analytical jacobian is supported. + /// + public bool IsJacobianSupported { get { return userDerivatives != null; } } + + #endregion Public Variables + + public FittingObjectiveModel(Func, double, double>function, Func, double, Vector> 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); + } + + /// + /// Set observed data to fit. + /// + public void SetObserved(Vector observedX, Vector observedY, Vector 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.Build.DenseOfDiagonalVector(weights); + + L = (weights == null) + ? null + : Weights.Diagonal().PointwiseSqrt(); + } + + /// + /// Set observed data to fit. + /// + 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.Build.DenseOfArray(weights); + SetObserved(Vector.Build.DenseOfArray(observedX), Vector.Build.DenseOfArray(observedY), wVector); + } + + /// + /// Set parameters. + /// + /// 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. + /// + /// The lower bounds of parameters. + /// The upper bounds of parameters. + /// /// The scaling constants of parameters + /// The list to the parameters fix or free. + public void SetParameters(Vector lowerBound = null, Vector upperBound = null, Vector scales = null, List 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; + } + + /// + /// Set parameters. + /// + /// 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. + /// + /// The lower bounds of parameters. + /// The upper bounds of parameters. + /// The scaling constants of parameters + /// The list to the parameters fix or free. + public void SetParameters(double[] lowerBound = null, double[] upperBound = null, double[] scales = null, bool[] isFixed = null) + { + var lb = (lowerBound == null) ? null : Vector.Build.DenseOfArray(lowerBound); + var ub = (upperBound == null) ? null : Vector.Build.DenseOfArray(upperBound); + var sc = (scales == null) ? null : Vector.Build.DenseOfArray(scales); + var fp = (isFixed == null) ? null : isFixed.ToList(); + + SetParameters(lb, ub, sc, fp); + } + + public void EvaluateFunction(Vector 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.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 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.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.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 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 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 NumericalJacobian(Vector parameters, Vector currentValues, int accuracyOrder = 2) + { + const double sqrtEpsilon = 1.4901161193847656250E-8; // sqrt(machineEpsilon) + + Matrix derivertives = Matrix.Build.Dense(NumberOfObservations, NumberOfParameters); + + var d = 0.000003 * parameters.PointwiseAbs().PointwiseMaximum(sqrtEpsilon); + + var h = Vector.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 ProjectParametersToInternal(Vector 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 ProjectParametersToExternal(Vector 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 ScaleFactorsOfJacobian(Vector Pint) + { + var scale = Vector.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 + } +} diff --git a/src/Numerics/Optimization/TrustRegionDogLegMinimizer.cs b/src/Numerics/Optimization/TrustRegionDogLegMinimizer.cs new file mode 100644 index 00000000..a8ab570f --- /dev/null +++ b/src/Numerics/Optimization/TrustRegionDogLegMinimizer.cs @@ -0,0 +1,322 @@ +using MathNet.Numerics.LinearAlgebra; +using System; + +namespace MathNet.Numerics.Optimization +{ + public sealed class TrustRegionDogLegMinimizer + { + #region Tolerances and options + + /// + /// The stopping threshold for infinity norm of the gradient. + /// + public static double GradientTolerance { get; set; } + + /// + /// The stopping threshold for L2 norm of the change of the parameters. + /// + public static double StepTolerance { get; set; } + + /// + /// The stopping threshold for the function value or L2 norm of the residuals. + /// + public static double FunctionTolerance { get; set; } + + /// + /// The stopping threshold for the trust region radius. + /// + public static double RadiusTolerance { get; set; } + + /// + /// The maximum number of iterations. + /// + 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 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(initialGuess), GradientTolerance, StepTolerance, FunctionTolerance, RadiusTolerance, MaximumIterations); + } + + /// + /// Non-linear least square fitting by the trust-region dogleg algorithm. + /// + /// The objective model, including function, jacobian, observations, and parameter bounds. + /// The initial guess values. + /// The stopping threshold for L2 norm of the residuals. + /// The stopping threshold for infinity norm of the gradient vector. + /// The stopping threshold for L2 norm of the change of parameters. + /// The stopping threshold for trust region radius + /// The max iterations. + /// + public static ModelMinimizationResult Minimum(IObjectiveModel objective, Vector 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, double, bool> SolveQuadraticSubproblem(IObjectiveModel objective, double delta) + { + Vector 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, double, bool>(Pstep, predictedReduction, hitBoundary); + } + + private static double FindBeta(double alpha, Vector sd, Vector 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; + } + } +}