diff --git a/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs b/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs new file mode 100644 index 00000000..2b7b553d --- /dev/null +++ b/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs @@ -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 RosenbrockModel(Vector p, Vector x) + { + var y = CreateVector.Dense(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 RosenbrockPrime(Vector p, Vector x) + { + var prime = Matrix.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 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 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 Rat43Model(Vector p, Vector x) + { + var y = CreateVector.Dense(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 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 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 BoxBodModel(Vector p, Vector x) + { + var y = CreateVector.Dense(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 BoxBodPrime(Vector p, Vector x) + { + var prime = Matrix.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 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[] { -1000, -100 }); + private Vector BoxBodUpperBound = new DenseVector(new double[] { 1000.0, 100 }); + private Vector 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 ThurberModel(Vector p, Vector x) + { + var y = CreateVector.Dense(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 ThurberPrime(Vector p, Vector x) + { + var prime = Matrix.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 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 ThurberStart = new DenseVector(new double[] { 1000.0, 1000.0, 400.0, 40.0, 0.7, 0.3, 0.03 }); + private Vector ThurberLowerBound = new DenseVector(new double[] { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }); + private Vector ThurberUpperBound = new DenseVector(new double[] { 1E6, 1E6, 1E6, 1E6, 1E6, 1E6, 1E6 }); + private Vector 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 + } +} 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..c9a7c66f --- /dev/null +++ b/src/Numerics/Optimization/IObjectiveModel.cs @@ -0,0 +1,74 @@ +using MathNet.Numerics.LinearAlgebra; +using System.Collections.Generic; + +namespace MathNet.Numerics.Optimization +{ + public interface IObjectiveModelEvaluation + { + IObjectiveModel CreateNew(); + + /// + /// Get the y-values of the observations. + /// + Vector ObservedY { get; } + + /// + /// Get the values of the weights for the observations. + /// + Matrix Weights { get; } + + /// + /// Get the y-values of the fitted model that correspond to the independent values. + /// + Vector ModelValues { get; } + + /// + /// Get the values of the parameters. + /// + Vector Point { get; } + + /// + /// Get the residual sum of squares. + /// + double Value { 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 number of calls to function. + /// + int FunctionEvaluations { get; set; } + + /// + /// Get the number of calls to jacobian. + /// + int JacobianEvaluations { get; set; } + + /// + /// Get the degree of freedom. + /// + int DegreeOfFreedom { get; } + + bool IsGradientSupported { get; } + bool IsHessianSupported { get; } + } + + public interface IObjectiveModel : IObjectiveModelEvaluation + { + void SetParameters(Vector initialGuess, List isFixed = null); + + void EvaluateAt(Vector parameters); + + IObjectiveModel Fork(); + + IObjectiveFunction ToObjectiveFunction(); + } +} diff --git a/src/Numerics/Optimization/ITrustRegionSubProblem.cs b/src/Numerics/Optimization/ITrustRegionSubProblem.cs new file mode 100644 index 00000000..5305fc44 --- /dev/null +++ b/src/Numerics/Optimization/ITrustRegionSubProblem.cs @@ -0,0 +1,12 @@ +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization +{ + public interface ITrustRegionSubproblem + { + Vector Pstep { get; } + bool HitBoundary { get; } + + void Solve(IObjectiveModel objective, double radius); + } +} diff --git a/src/Numerics/Optimization/LevenbergMarquardtMinimizer.cs b/src/Numerics/Optimization/LevenbergMarquardtMinimizer.cs new file mode 100644 index 00000000..f340515d --- /dev/null +++ b/src/Numerics/Optimization/LevenbergMarquardtMinimizer.cs @@ -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 + { + /// + /// The scale factor for initial mu + /// + 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 initialGuess, + Vector lowerBound = null, Vector upperBound = null, Vector scales = null, List 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(lowerBound); + var ub = (upperBound == null) ? null : CreateVector.Dense(upperBound); + var sc = (scales == null) ? null : CreateVector.Dense(scales); + var fx = (isFixed == null) ? null : isFixed.ToList(); + + return Minimum(objective, CreateVector.DenseOfArray(initialGuess), lb, ub, sc, fx, 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 NonlinearMinimizationResult Minimum(IObjectiveModel objective, Vector initialGuess, + Vector lowerBound = null, Vector upperBound = null, Vector scales = null, List 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.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); + } + } +} diff --git a/src/Numerics/Optimization/NonlinearMinimizationResult.cs b/src/Numerics/Optimization/NonlinearMinimizationResult.cs new file mode 100644 index 00000000..a67d3c85 --- /dev/null +++ b/src/Numerics/Optimization/NonlinearMinimizationResult.cs @@ -0,0 +1,78 @@ +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization +{ + public class NonlinearMinimizationResult + { + public IObjectiveModel ModelInfoAtMinimum { get; private set; } + + /// + /// Returns the best fit parameters. + /// + public Vector MinimizingPoint { get { return ModelInfoAtMinimum.Point; } } + + /// + /// Returns the standard errors of the corresponding parameters + /// + public Vector StandardErrors { get; private set; } + + /// + /// Returns the y-values of the fitted model that correspond to the independent values. + /// + public Vector MinimizedValues { get { return ModelInfoAtMinimum.ModelValues; } } + + /// + /// Returns the covariance matrix at minimizing point. + /// + public Matrix Covariance { get; private set; } + + /// + /// Returns the correlation matrix at minimizing point. + /// + public Matrix 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; + } + } + } +} diff --git a/src/Numerics/Optimization/NonlinearMinimizerBase.cs b/src/Numerics/Optimization/NonlinearMinimizerBase.cs new file mode 100644 index 00000000..1ff68382 --- /dev/null +++ b/src/Numerics/Optimization/NonlinearMinimizerBase.cs @@ -0,0 +1,302 @@ +using MathNet.Numerics.LinearAlgebra; +using System; +using System.Linq; + +namespace MathNet.Numerics.Optimization +{ + public abstract class NonlinearMinimizerBase + { + /// + /// The stopping threshold for the function value or L2 norm of the residuals. + /// + public static double FunctionTolerance { get; set; } + + /// + /// The stopping threshold for L2 norm of the change of the parameters. + /// + public static double StepTolerance { get; set; } + + /// + /// The stopping threshold for infinity norm of the gradient. + /// + public static double GradientTolerance { get; set; } + + /// + /// The maximum number of iterations. + /// + public static int MaximumIterations { get; set; } + + /// + /// The lower bound of the parameters. + /// + public static Vector LowerBound { get; private set; } + + /// + /// The upper bound of the parameters. + /// + public static Vector UpperBound { get; private set; } + + /// + /// The scale factors for the parameters. + /// + public static Vector 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 parameters, Vector lowerBound = null, Vector upperBound = null, Vector 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 Pint) + { + var Pext = ProjectToExternalParameters(Pint); + objective.EvaluateAt(Pext); + return objective.Value; + } + + protected static Tuple, Matrix> EvaluateJacobian(IObjectiveModel objective, Vector 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, Matrix>(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 ProjectToInternalParameters(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] = (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 ProjectToExternalParameters(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] = (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 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] = (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 + } +} diff --git a/src/Numerics/Optimization/ObjectiveFunction.cs b/src/Numerics/Optimization/ObjectiveFunction.cs index a0c981d9..eaf1bfed 100644 --- a/src/Numerics/Optimization/ObjectiveFunction.cs +++ b/src/Numerics/Optimization/ObjectiveFunction.cs @@ -114,5 +114,111 @@ namespace MathNet.Numerics.Optimization { return new ScalarObjectiveFunction(function, derivative, secondDerivative); } + + /// + /// objective model with a user supplied jacobian for non-linear least squares regression. + /// + public static IObjectiveModel NonlinearModel(Func, Vector, Vector> function, + Func, Vector, Matrix> derivatives, + Vector observedX, Vector observedY, Vector weight = null) + { + var objective = new NonlinearObjectiveFunction(function, derivatives); + objective.SetObserved(observedX, observedY, weight); + return objective; + } + + /// + /// Objective model for non-linear least squares regression. + /// + public static IObjectiveModel NonlinearModel(Func, Vector, Vector> function, + Vector observedX, Vector observedY, Vector weight = null, + int accuracyOrder = 2) + { + var objective = new NonlinearObjectiveFunction(function, accuracyOrder: accuracyOrder); + objective.SetObserved(observedX, observedY, weight); + return objective; + } + + /// + /// Objective model with a user supplied jacobian for non-linear least squares regression. + /// + public static IObjectiveModel NonlinearModel(Func, double, double> function, + Func, double, Vector> derivatives, + Vector observedX, Vector observedY, Vector weight = null) + { + Vector func(Vector point, Vector x) + { + var functionValues = CreateVector.Dense(x.Count); + for (int i = 0; i < x.Count; i++) + { + functionValues[i] = function(point, x[i]); + } + + return functionValues; + } + + Matrix prime(Vector point, Vector x) + { + var derivativeValues = CreateMatrix.Dense(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; + } + + /// + /// Objective model for non-linear least squares regression. + /// + public static IObjectiveModel NonlinearModel(Func, double, double> function, + Vector observedX, Vector observedY, Vector weight = null, + int accuracyOrder = 2) + { + Vector func(Vector point, Vector x) + { + var functionValues = CreateVector.Dense(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; + } + + /// + /// Objective function with a user supplied jacobian for nonlinear least squares regression. + /// + public static IObjectiveFunction NonlinearFunction(Func, Vector, Vector> function, + Func, Vector, Matrix> derivatives, + Vector observedX, Vector observedY, Vector weight = null) + { + var objective = new NonlinearObjectiveFunction(function, derivatives); + objective.SetObserved(observedX, observedY, weight); + return objective.ToObjectiveFunction(); + } + + /// + /// Objective function for nonlinear least squares regression. + /// The numerical jacobian with accuracy order is used. + /// + public static IObjectiveFunction NonlinearFunction(Func, Vector, Vector> function, + Vector observedX, Vector observedY, Vector weight = null, + int accuracyOrder = 2) + { + var objective = new NonlinearObjectiveFunction(function, null, accuracyOrder: accuracyOrder); + objective.SetObserved(observedX, observedY, weight); + return objective.ToObjectiveFunction(); + } } } diff --git a/src/Numerics/Optimization/ObjectiveFunctions/NonlinearObjectiveFunction.cs b/src/Numerics/Optimization/ObjectiveFunctions/NonlinearObjectiveFunction.cs new file mode 100644 index 00000000..7c5e36b0 --- /dev/null +++ b/src/Numerics/Optimization/ObjectiveFunctions/NonlinearObjectiveFunction.cs @@ -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, Vector> userFunction; // (p, x) => f(x; p) + readonly Func, Vector, Matrix> userDerivative; // (p, x) => df(x; p)/dp + readonly int accuracyOrder; // the desired accuracy order to evaluate the jacobian by numerical approximaiton. + + Vector coefficients; + + bool hasFunctionValue; + double functionValue; // the residual sum of squares, residuals * residuals. + Vector residuals; // the weighted error values + + bool hasJacobianValue; + Matrix jacobianValue; // the Jacobian matrix. + Vector gradientValue; // the Gradient vector. + Matrix hessianValue; // the Hessian matrix. + + #endregion Private Variables + + #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. + /// + public Matrix Weights { get; private set; } + private Vector L; // Weights = LL' + + /// + /// Get whether parameters are fixed or free. + /// + public List IsFixed { 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 (Point == null) ? 0 : Point.Count; } } + + /// + /// Get the degree of freedom + /// + public int DegreeOfFreedom + { + get + { + var df = NumberOfObservations - NumberOfParameters; + if (IsFixed != null) + { + df = df + IsFixed.Count(p => p == true); + } + return df; + } + } + + /// + /// Get the number of calls to function. + /// + public int FunctionEvaluations { get; set; } + + /// + /// Get the number of calls to jacobian. + /// + public int JacobianEvaluations { get; set; } + + #endregion Public Variables + + public NonlinearObjectiveFunction(Func, Vector, Vector> function, + Func, Vector, Matrix> 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); + } + + /// + /// Set or get the values of the parameters. + /// + public Vector Point { get { return coefficients; } } + + /// + /// Get the y-values of the fitted model that correspond to the independent values. + /// + public Vector ModelValues { get; private set; } + + /// + /// Get the residual sum of squares. + /// + public double Value + { + get + { + if (!hasFunctionValue) + { + EvaluateFunction(); + hasFunctionValue = true; + } + return functionValue; + } + } + + /// + /// Get the Gradient vector of x and p. + /// + public Vector Gradient + { + get + { + if (!hasJacobianValue) + { + EvaluateJacobian(); + hasJacobianValue = true; + } + return gradientValue; + } + } + + /// + /// Get the Hessian matrix of x and p, J'WJ + /// + public Matrix Hessian + { + get + { + if (!hasJacobianValue) + { + EvaluateJacobian(); + hasJacobianValue = true; + } + return hessianValue; + } + } + + public bool IsGradientSupported { get { return true; } } + public bool IsHessianSupported { get { return true; } } + + /// + /// 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 parameters and bounds. + /// + /// The initial values of parameters. + /// The list to the parameters fix or free. + public void SetParameters(Vector initialGuess, List 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 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, Matrix> function(Vector point) + { + EvaluateAt(point); + + return new Tuple, Matrix>(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.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 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 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 + } +} diff --git a/src/Numerics/Optimization/Subproblems/DogLegSubproblem.cs b/src/Numerics/Optimization/Subproblems/DogLegSubproblem.cs new file mode 100644 index 00000000..ce1ac838 --- /dev/null +++ b/src/Numerics/Optimization/Subproblems/DogLegSubproblem.cs @@ -0,0 +1,45 @@ +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Optimization.Subproblems +{ + internal class DogLegSubproblem : ITrustRegionSubproblem + { + public Vector 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); + } + } + } +} diff --git a/src/Numerics/Optimization/Subproblems/NewtonCGSubproblem.cs b/src/Numerics/Optimization/Subproblems/NewtonCGSubproblem.cs new file mode 100644 index 00000000..73149818 --- /dev/null +++ b/src/Numerics/Optimization/Subproblems/NewtonCGSubproblem.cs @@ -0,0 +1,65 @@ +using MathNet.Numerics.LinearAlgebra; +using System; + +namespace MathNet.Numerics.Optimization.Subproblems +{ + internal class NewtonCGSubproblem : ITrustRegionSubproblem + { + public Vector 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.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; + } + } + } +} diff --git a/src/Numerics/Optimization/Subproblems/Util.cs b/src/Numerics/Optimization/Subproblems/Util.cs new file mode 100644 index 00000000..4ea986f8 --- /dev/null +++ b/src/Numerics/Optimization/Subproblems/Util.cs @@ -0,0 +1,35 @@ +using MathNet.Numerics.LinearAlgebra; +using System; + +namespace MathNet.Numerics.Optimization.Subproblems +{ + internal static class Util + { + public static Tuple 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 + // + // 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(beta1, beta2) + : new Tuple(beta2, beta1); + } + } +} diff --git a/src/Numerics/Optimization/TrustRegionDogLegMinimizer.cs b/src/Numerics/Optimization/TrustRegionDogLegMinimizer.cs new file mode 100644 index 00000000..77a220f0 --- /dev/null +++ b/src/Numerics/Optimization/TrustRegionDogLegMinimizer.cs @@ -0,0 +1,12 @@ +namespace MathNet.Numerics.Optimization +{ + public sealed class TrustRegionDogLegMinimizer : TrustRegionMinimizerBase + { + /// + /// Non-linear least square fitting by the trust region dogleg algorithm. + /// + 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) + { } + } +} diff --git a/src/Numerics/Optimization/TrustRegionMinimizerBase.cs b/src/Numerics/Optimization/TrustRegionMinimizerBase.cs new file mode 100644 index 00000000..1ab8eb36 --- /dev/null +++ b/src/Numerics/Optimization/TrustRegionMinimizerBase.cs @@ -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 + { + /// + /// The trust region subproblem. + /// + public static ITrustRegionSubproblem Subproblem; + + /// + /// The stopping threshold for the trust region radius. + /// + 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 initialGuess, + Vector lowerBound = null, Vector upperBound = null, Vector scales = null, List 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(lowerBound); + var ub = (upperBound == null) ? null : CreateVector.Dense(upperBound); + var sc = (scales == null) ? null : CreateVector.Dense(scales); + var fx = (isFixed == null) ? null : isFixed.ToList(); + + return Minimum(Subproblem, objective, CreateVector.DenseOfArray(initialGuess), lb, ub, sc, fx, + GradientTolerance, StepTolerance, FunctionTolerance, RadiusTolerance, MaximumIterations); + } + + /// + /// Non-linear least square fitting by the trust-region algorithm. + /// + /// The objective model, including function, jacobian, observations, and parameter bounds. + /// The subproblem + /// 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 NonlinearMinimizationResult Minimum(ITrustRegionSubproblem subproblem, IObjectiveModel objective, Vector initialGuess, + Vector lowerBound = null, Vector upperBound = null, Vector scales = null, List 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.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); + } + } +} diff --git a/src/Numerics/Optimization/TrustRegionNewtonCGMinimizer.cs b/src/Numerics/Optimization/TrustRegionNewtonCGMinimizer.cs new file mode 100644 index 00000000..7e6d5f5a --- /dev/null +++ b/src/Numerics/Optimization/TrustRegionNewtonCGMinimizer.cs @@ -0,0 +1,13 @@ +namespace MathNet.Numerics.Optimization +{ + public sealed class TrustRegionNewtonCGMinimizer : TrustRegionMinimizerBase + { + /// + /// Non-linear least square fitting by the trust region Newton-Conjugate-Gradient algorithm. + /// + 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) + { } + } +} + diff --git a/src/Numerics/Optimization/TrustRegionSubProblem.cs b/src/Numerics/Optimization/TrustRegionSubProblem.cs new file mode 100644 index 00000000..3018f8b6 --- /dev/null +++ b/src/Numerics/Optimization/TrustRegionSubProblem.cs @@ -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(); + } + } +}