diff --git a/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs b/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs index faae726c..2b7b553d 100644 --- a/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs +++ b/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs @@ -19,16 +19,23 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests // best fitted parameters: // a = 1 // b = 1 - private double RosenbrockModel(Vector p, double x) + private Vector RosenbrockModel(Vector p, Vector x) { - var y = Math.Pow(1.0 - p[0], 2) + 100.0 * Math.Pow(p[1] - p[0] * p[0], 2); + 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 Vector RosenbrockPrime(Vector p, double x) + private Matrix RosenbrockPrime(Vector p, Vector 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]); + 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); @@ -39,29 +46,27 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests private Vector RosebbrockLowerBound = new DenseVector(new double[] { -5.0, -5.0 }); private Vector RosenbrockUpperBound = new DenseVector(new double[] { 5.0, 5.0 }); - #endregion Rosenbrock - [Test] public void Rosenbrock_LM_Der() { // unconstrained - var obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockPrime, RosenbrockX, RosenbrockY); + 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.BestFitParameters.Count; i++) + for (int i = 0; i < result.MinimizingPoint.Count; i++) { - AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.BestFitParameters[i], 3); + AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.MinimizingPoint[i], 2); } // box constrained - obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockPrime, RosenbrockX, RosenbrockY); + 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.BestFitParameters.Count; i++) + for (int i = 0; i < result.MinimizingPoint.Count; i++) { - AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.BestFitParameters[i], 3); + AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.MinimizingPoint[i], 2); } } @@ -69,52 +74,54 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests public void Rosenbrock_LM_Dif() { // unconstrained - var obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockX, RosenbrockY, accuracyOrder:2); + 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.BestFitParameters.Count; i++) + for (int i = 0; i < result.MinimizingPoint.Count; i++) { - AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.BestFitParameters[i], 3); + AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.MinimizingPoint[i], 2); } // box constrained - obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockX, RosenbrockY, accuracyOrder: 6); + 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.BestFitParameters.Count; i++) + for (int i = 0; i < result.MinimizingPoint.Count; i++) { - AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.BestFitParameters[i], 3); + AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.MinimizingPoint[i], 2); } } [Test] public void Rosenbrock_Bfgs_Dif() { - var obj = ObjectiveModel.FittingFunction(RosenbrockModel, RosenbrockX, RosenbrockY, accuracyOrder: 6); - var solver = new BfgsMinimizer(1e-10, 1e-10, 1e-10, 1000); + 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], 3); + AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.MinimizingPoint[i], 2); } } [Test] public void Rosenbrock_LBfgs_Dif() { - var obj = ObjectiveModel.FittingFunction(RosenbrockModel, RosenbrockX, RosenbrockY, accuracyOrder: 6); - var solver = new LimitedMemoryBfgsMinimizer(1e-10, 1e-10, 1e-10, 1000); + 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], 3); + 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) @@ -124,9 +131,13 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests // 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) + private Vector Rat43Model(Vector p, Vector x) { - var y = p[0] / Math.Pow(1.0 + Math.Exp(p[1] - p[2] * x), 1.0 / p[3]); + 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[] { @@ -147,18 +158,16 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests private Vector Rat43Start1 = new DenseVector(new double[] { 100, 10, 1, 1 }); private Vector Rat43Start2 = new DenseVector(new double[] { 700, 5, 0.75, 1.3 }); - #endregion Rat43 - [Test] public void Rat43_LM_Dif() { - var obj = ObjectiveModel.FittingModel(Rat43Model, Rat43X, Rat43Y, accuracyOrder: 6); + 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.BestFitParameters.Count; i++) + for (int i = 0; i < result.MinimizingPoint.Count; i++) { - AssertHelpers.AlmostEqualRelative(Rat43Pbest[i], result.BestFitParameters[i], 6); + AssertHelpers.AlmostEqualRelative(Rat43Pbest[i], result.MinimizingPoint[i], 6); AssertHelpers.AlmostEqualRelative(Rat43Pstd[i], result.StandardErrors[i], 6); } } @@ -166,13 +175,13 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests [Test] public void Rat43_TRDL_Dif() { - var obj = ObjectiveModel.FittingModel(Rat43Model, Rat43X, Rat43Y, accuracyOrder: 6); + 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.BestFitParameters.Count; i++) + for (int i = 0; i < result.MinimizingPoint.Count; i++) { - AssertHelpers.AlmostEqualRelative(Rat43Pbest[i], result.BestFitParameters[i], 2); + AssertHelpers.AlmostEqualRelative(Rat43Pbest[i], result.MinimizingPoint[i], 2); AssertHelpers.AlmostEqualRelative(Rat43Pstd[i], result.StandardErrors[i], 2); } } @@ -180,13 +189,13 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests [Test] public void Rat43_TRNCG_Dif() { - var obj = ObjectiveModel.FittingModel(Rat43Model, Rat43X, Rat43Y, accuracyOrder: 6); + 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.BestFitParameters.Count; i++) + for (int i = 0; i < result.MinimizingPoint.Count; i++) { - AssertHelpers.AlmostEqualRelative(Rat43Pbest[i], result.BestFitParameters[i], 2); + AssertHelpers.AlmostEqualRelative(Rat43Pbest[i], result.MinimizingPoint[i], 2); AssertHelpers.AlmostEqualRelative(Rat43Pstd[i], result.StandardErrors[i], 2); } } @@ -194,7 +203,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests [Test] public void Rat43_Bfgs_Dif() { - var obj = ObjectiveModel.FittingFunction(Rat43Model, Rat43X, Rat43Y, accuracyOrder: 6); + 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); @@ -207,7 +216,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests [Test] public void Rat43_LBfgs_Dif() { - var obj = ObjectiveModel.FittingFunction(Rat43Model, Rat43X, Rat43Y, accuracyOrder: 6); + 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); @@ -217,6 +226,8 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests } } + #endregion Rat43 + #region BoxBod // model: BoxBod (https://www.itl.nist.gov/div898/strd/nls/data/boxbod.shtml) @@ -227,16 +238,23 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests // best fitted parameters: // a = 2.1380940889E+02 +/- 1.2354515176E+01 // b = 5.4723748542E-01 +/- 1.0455993237E-01 - private double BoxBodModel(Vector p, double x) + private Vector BoxBodModel(Vector p, Vector x) { - var y = p[0] * (1.0 - Math.Exp(-p[1] * 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 Vector BoxBodPrime(Vector p, double x) + private Matrix BoxBodPrime(Vector p, Vector 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); + 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 }); @@ -250,92 +268,90 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests private Vector BoxBodUpperBound = new DenseVector(new double[] { 1000.0, 100 }); private Vector BoxBodScales = new DenseVector(new double[] { 100.0, 0.1 }); - #endregion BoxBod - [Test] public void BoxBod_LM_Der() { // unconstrained - var obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY); + var obj = ObjectiveFunction.NonlinearModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY); var solver = new LevenbergMarquardtMinimizer(); var result = solver.FindMinimum(obj, BoxBodStart1); - for (int i = 0; i < result.BestFitParameters.Count; i++) + for (int i = 0; i < result.MinimizingPoint.Count; i++) { - AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6); + 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 = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY); + 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.BestFitParameters.Count; i++) + for (int i = 0; i < result.MinimizingPoint.Count; i++) { - AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6); + AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.MinimizingPoint[i], 6); AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6); } // lower < parameters, no scales - obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY); + obj = ObjectiveFunction.NonlinearModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY); solver = new LevenbergMarquardtMinimizer(); result = solver.FindMinimum(obj, BoxBodStart1, lowerBound: BoxBodLowerBound); - for (int i = 0; i < result.BestFitParameters.Count; i++) + for (int i = 0; i < result.MinimizingPoint.Count; i++) { - AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6); + AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.MinimizingPoint[i], 6); AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6); } // lower < parameters, scales - obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY); + 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.BestFitParameters.Count; i++) + for (int i = 0; i < result.MinimizingPoint.Count; i++) { - AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6); + AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.MinimizingPoint[i], 6); AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6); } // parameters < upper, no scales - obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY); + obj = ObjectiveFunction.NonlinearModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY); solver = new LevenbergMarquardtMinimizer(); result = solver.FindMinimum(obj, BoxBodStart1, upperBound: BoxBodUpperBound); - for (int i = 0; i < result.BestFitParameters.Count; i++) + for (int i = 0; i < result.MinimizingPoint.Count; i++) { - AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6); + AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.MinimizingPoint[i], 6); AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6); } // parameters < upper, scales - obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY); + 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.BestFitParameters.Count; i++) + for (int i = 0; i < result.MinimizingPoint.Count; i++) { - AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6); + AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.MinimizingPoint[i], 6); AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6); } // only scales - obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY); + obj = ObjectiveFunction.NonlinearModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY); solver = new LevenbergMarquardtMinimizer(); result = solver.FindMinimum(obj, BoxBodStart1, scales: BoxBodScales); - for (int i = 0; i < result.BestFitParameters.Count; i++) + for (int i = 0; i < result.MinimizingPoint.Count; i++) { - AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6); + AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.MinimizingPoint[i], 6); AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6); } } @@ -344,24 +360,24 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests public void BoxBod_LM_Dif() { // unconstrained - var obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodX, BoxBodY, accuracyOrder:6); + 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.BestFitParameters.Count; i++) + for (int i = 0; i < result.MinimizingPoint.Count; i++) { - AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6); + AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.MinimizingPoint[i], 6); AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6); } // box constrained - obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodX, BoxBodY, accuracyOrder: 6); + 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.BestFitParameters.Count; i++) + for (int i = 0; i < result.MinimizingPoint.Count; i++) { - AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6); + AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.MinimizingPoint[i], 6); AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6); } } @@ -369,13 +385,13 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests [Test] public void BoxBod_TRDL_Dif() { - var obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodX, BoxBodY, accuracyOrder: 6); + 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.BestFitParameters.Count; i++) + for (int i = 0; i < result.MinimizingPoint.Count; i++) { - AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 3); + AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.MinimizingPoint[i], 3); AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 3); } } @@ -383,14 +399,13 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests [Test] public void BoxBod_TRNCG_Dif() { - var obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY); - //var obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodX, BoxBodY, accuracyOrder: 6); + 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.BestFitParameters.Count; i++) + for (int i = 0; i < result.MinimizingPoint.Count; i++) { - AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 3); + AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.MinimizingPoint[i], 3); AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 3); } } @@ -398,7 +413,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests [Test] public void BoxBod_Bfgs_Der() { - var obj = ObjectiveModel.FittingFunction(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY); + 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); @@ -408,6 +423,21 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests } } + [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) @@ -428,32 +458,38 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests // 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) + private Vector ThurberModel(Vector p, Vector x) { - var xSq = x * x; - var xCb = xSq * 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]; - var y = (p[0] + p[1] * x + p[2] * xSq + p[3] * xCb) - / (1 + p[4] * x + p[5] * xSq + p[6] * xCb); + 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 Vector ThurberPrime(Vector p, double x) + private Matrix ThurberPrime(Vector p, Vector 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; + 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[] { @@ -485,18 +521,16 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests 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 }); - #endregion Thurber - [Test] public void Thurber_LM_Der() { - var obj = ObjectiveModel.FittingModel(ThurberModel, ThurberPrime, ThurberX, ThurberY); + var obj = ObjectiveFunction.NonlinearModel(ThurberModel, ThurberPrime, ThurberX, ThurberY); var solver = new LevenbergMarquardtMinimizer(); var result = solver.FindMinimum(obj, ThurberStart); - for (int i = 0; i < result.BestFitParameters.Count; i++) + for (int i = 0; i < result.MinimizingPoint.Count; i++) { - AssertHelpers.AlmostEqualRelative(ThurberPbest[i], result.BestFitParameters[i], 6); + AssertHelpers.AlmostEqualRelative(ThurberPbest[i], result.MinimizingPoint[i], 6); AssertHelpers.AlmostEqualRelative(ThurberPstd[i], result.StandardErrors[i], 6); } } @@ -504,13 +538,13 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests [Test] public void Thurber_LM_Dif() { - var obj = ObjectiveModel.FittingModel(ThurberModel, ThurberX, ThurberY, accuracyOrder: 6); + 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.BestFitParameters.Count; i++) + for (int i = 0; i < result.MinimizingPoint.Count; i++) { - AssertHelpers.AlmostEqualRelative(ThurberPbest[i], result.BestFitParameters[i], 6); + AssertHelpers.AlmostEqualRelative(ThurberPbest[i], result.MinimizingPoint[i], 6); AssertHelpers.AlmostEqualRelative(ThurberPstd[i], result.StandardErrors[i], 6); } } @@ -518,13 +552,13 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests [Test] public void Thurber_TRDL_Dif() { - var obj = ObjectiveModel.FittingModel(ThurberModel, ThurberX, ThurberY, accuracyOrder: 6); + 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.BestFitParameters.Count; i++) + for (int i = 0; i < result.MinimizingPoint.Count; i++) { - AssertHelpers.AlmostEqualRelative(ThurberPbest[i], result.BestFitParameters[i], 3); + AssertHelpers.AlmostEqualRelative(ThurberPbest[i], result.MinimizingPoint[i], 3); AssertHelpers.AlmostEqualRelative(ThurberPstd[i], result.StandardErrors[i], 3); } } @@ -532,13 +566,13 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests [Test] public void Thurber_TRNCG_Dif() { - var obj = ObjectiveModel.FittingModel(ThurberModel, ThurberX, ThurberY, accuracyOrder: 6); + 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.BestFitParameters.Count; i++) + for (int i = 0; i < result.MinimizingPoint.Count; i++) { - AssertHelpers.AlmostEqualRelative(ThurberPbest[i], result.BestFitParameters[i], 3); + AssertHelpers.AlmostEqualRelative(ThurberPbest[i], result.MinimizingPoint[i], 3); AssertHelpers.AlmostEqualRelative(ThurberPstd[i], result.StandardErrors[i], 3); } } @@ -546,7 +580,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests [Test] public void Thurber_Bfgs_Dif() { - var obj = ObjectiveModel.FittingFunction(ThurberModel, ThurberX, ThurberY, accuracyOrder: 6); + 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); @@ -559,7 +593,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests [Test] public void Thurber_BfgsB_Dif() { - var obj = ObjectiveModel.FittingFunction(ThurberModel, ThurberX, ThurberY, accuracyOrder: 6); + 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); @@ -572,7 +606,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests [Test] public void Thurber_LBfgs_Dif() { - var obj = ObjectiveModel.FittingFunction(ThurberModel, ThurberX, ThurberY, accuracyOrder: 6); + 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); @@ -581,5 +615,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests AssertHelpers.AlmostEqualRelative(ThurberPbest[i], result.MinimizingPoint[i], 6); } } + + #endregion Thurber } } diff --git a/src/Numerics/Optimization/IObjectiveModel.cs b/src/Numerics/Optimization/IObjectiveModel.cs index b106a347..c9a7c66f 100644 --- a/src/Numerics/Optimization/IObjectiveModel.cs +++ b/src/Numerics/Optimization/IObjectiveModel.cs @@ -59,17 +59,14 @@ namespace MathNet.Numerics.Optimization bool IsGradientSupported { get; } bool IsHessianSupported { get; } - - bool IsFinished { get; set; } } public interface IObjectiveModel : IObjectiveModelEvaluation { - void SetParameters(Vector initialGuess, Vector lowerBound = null, Vector upperBound = null, Vector scales = null, List isFixed = null); + void SetParameters(Vector initialGuess, List isFixed = null); void EvaluateAt(Vector parameters); - /// Create a new independent copy of this objective function, evaluated at the same point. IObjectiveModel Fork(); IObjectiveFunction ToObjectiveFunction(); diff --git a/src/Numerics/Optimization/LevenbergMarquardtMinimizer.cs b/src/Numerics/Optimization/LevenbergMarquardtMinimizer.cs index c39b0119..f340515d 100644 --- a/src/Numerics/Optimization/LevenbergMarquardtMinimizer.cs +++ b/src/Numerics/Optimization/LevenbergMarquardtMinimizer.cs @@ -5,55 +5,23 @@ using System.Linq; namespace MathNet.Numerics.Optimization { - public sealed class LevenbergMarquardtMinimizer + public class LevenbergMarquardtMinimizer : NonlinearMinimizerBase { - #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) + 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; - GradientTolerance = gradientTolerance; - StepTolerance = stepTolerance; - FunctionTolerance = functionTolerance; - MaximumIterations = maximumIterations; } public NonlinearMinimizationResult FindMinimum(IObjectiveModel objective, Vector initialGuess, Vector lowerBound = null, Vector upperBound = null, Vector scales = null, List isFixed = null) { - if (objective == null) - throw new ArgumentNullException("objective"); - if (initialGuess == null) - throw new ArgumentNullException("initialGuess"); - - return Minimum(objective, initialGuess, lowerBound, upperBound, scales, isFixed, InitialMu, FunctionTolerance, GradientTolerance, StepTolerance, MaximumIterations); + return Minimum(objective, initialGuess, lowerBound, upperBound, scales, isFixed, InitialMu, GradientTolerance, StepTolerance, FunctionTolerance, MaximumIterations); } public NonlinearMinimizationResult FindMinimum(IObjectiveModel objective, double[] initialGuess, @@ -85,7 +53,7 @@ namespace MathNet.Numerics.Optimization /// 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-18, double stepTolerance = 1E-18, double functionTolerance = 1E-18, int maximumIterations = -1) + 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. // @@ -125,23 +93,16 @@ namespace MathNet.Numerics.Optimization if (objective == null) throw new ArgumentNullException("objective"); - if (initialGuess == null) - throw new ArgumentNullException("initialGuess"); - - objective.SetParameters(initialGuess, lowerBound, upperBound, scales, isFixed); + ValidateBounds(initialGuess, lowerBound, upperBound, scales); + + objective.SetParameters(initialGuess, isFixed); ExitCondition exitCondition = ExitCondition.None; - // Initialize objective - objective.FunctionEvaluations = 0; - objective.JacobianEvaluations = 0; - objective.IsFinished = false; - // First, calculate function values and setup variables - objective.EvaluateAt(initialGuess); - var P = objective.Point; // current parameters + var P = ProjectToInternalParameters(initialGuess); // current internal parameters var Pstep = Vector.Build.Dense(P.Count); // the change of parameters - var RSS = objective.Value; // Residual Sum of Squares = R'R + var RSS = EvaluateFunction(objective, P); // Residual Sum of Squares = R'R if (maximumIterations < 0) { @@ -168,8 +129,9 @@ namespace MathNet.Numerics.Optimization } // Evaluate gradient and Hessian - var Gradient = objective.Gradient; - var Hessian = objective.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 @@ -200,14 +162,13 @@ namespace MathNet.Numerics.Optimization // if ||ΔP|| <= xTol * (||P|| + xTol), found and stop if (Pstep.L2Norm() <= stepTolerance * (stepTolerance + P.DotProduct(P))) { - exitCondition = ExitCondition.RelativePoints; // SmallRelativeParameters + exitCondition = ExitCondition.RelativePoints; break; } var Pnew = P + Pstep; // new parameters to test - - objective.EvaluateAt(Pnew); - var RSSnew = objective.Value; + // evaluate function at Pnew + var RSSnew = EvaluateFunction(objective, Pnew); if (double.IsNaN(RSSnew)) { @@ -229,8 +190,9 @@ namespace MathNet.Numerics.Optimization RSS = RSSnew; // update gradient and Hessian - Gradient = objective.Gradient; - Hessian = objective.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 diff --git a/src/Numerics/Optimization/NonlinearMinimizationResult.cs b/src/Numerics/Optimization/NonlinearMinimizationResult.cs index 12381bbf..a67d3c85 100644 --- a/src/Numerics/Optimization/NonlinearMinimizationResult.cs +++ b/src/Numerics/Optimization/NonlinearMinimizationResult.cs @@ -1,8 +1,4 @@ using MathNet.Numerics.LinearAlgebra; -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; namespace MathNet.Numerics.Optimization { @@ -13,7 +9,7 @@ namespace MathNet.Numerics.Optimization /// /// Returns the best fit parameters. /// - public Vector BestFitParameters { get { return ModelInfoAtMinimum.Point; } } + public Vector MinimizingPoint { get { return ModelInfoAtMinimum.Point; } } /// /// Returns the standard errors of the corresponding parameters @@ -23,17 +19,20 @@ namespace MathNet.Numerics.Optimization /// /// Returns the y-values of the fitted model that correspond to the independent values. /// - public Vector BestFitValues { get { return ModelInfoAtMinimum.ModelValues; } } + public Vector MinimizedValues { get { return ModelInfoAtMinimum.ModelValues; } } /// - /// Returns the residual sum of squares. + /// Returns the covariance matrix at minimizing point. /// - public double Residue { get { return ModelInfoAtMinimum.Value; } } - public double DegreeOfFreedom { get { return ModelInfoAtMinimum.DegreeOfFreedom; } } 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) @@ -42,16 +41,15 @@ namespace MathNet.Numerics.Optimization Iterations = iterations; ReasonForExit = reasonForExit; - AnalyzeResult(modelInfo); + EvaluateCovariance(modelInfo); } - private void AnalyzeResult(IObjectiveModel objective) + private void EvaluateCovariance(IObjectiveModel objective) { - objective.IsFinished = true; - objective.EvaluateAt(objective.Point); + objective.EvaluateAt(objective.Point); // Hessian may be not yet updated. var Hessian = objective.Hessian; - if (Hessian == null || DegreeOfFreedom < 1) + if (Hessian == null || objective.DegreeOfFreedom < 1) { Covariance = null; Correlation = null; @@ -59,7 +57,7 @@ namespace MathNet.Numerics.Optimization return; } - Covariance = Hessian.PseudoInverse() * objective.Value / DegreeOfFreedom; + Covariance = Hessian.PseudoInverse() * objective.Value / objective.DegreeOfFreedom; if (Covariance != 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/ObjectiveModel.cs b/src/Numerics/Optimization/ObjectiveModel.cs deleted file mode 100644 index 77cae337..00000000 --- a/src/Numerics/Optimization/ObjectiveModel.cs +++ /dev/null @@ -1,56 +0,0 @@ -using MathNet.Numerics.LinearAlgebra; -using MathNet.Numerics.Optimization.ObjectiveModels; -using System; - -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) - { - var objective = new FittingObjectiveModel(function, derivatives); - objective.SetObserved(observedX, observedY, weight); - return objective; - } - - /// - /// Fitting model for non-linear least squares regression. - /// - public static IObjectiveModel FittingModel(Func, double, double> function, - Vector observedX, Vector observedY, Vector weight = null, - int accuracyOrder = 2) - { - var objective = new FittingObjectiveModel(function, accuracyOrder: accuracyOrder); - objective.SetObserved(observedX, observedY, weight); - return objective; - } - - /// - /// Fitting function with a user supplied jacobian for nonlinear least squares regression by the line search algorithm. - /// - public static IObjectiveFunction FittingFunction(Func, double, double> function, Func, double, Vector> derivatives, - Vector observedX, Vector observedY, Vector weight = null) - { - var objective = new FittingObjectiveModel(function, derivatives); - objective.SetObserved(observedX, observedY, weight); - return objective.ToObjectiveFunction(); - } - - /// - /// Fitting function for nonlinear least squares regression by the line search algorithm. - /// The numerical jacobian with accuracy order is used. - /// - public static IObjectiveFunction FittingFunction(Func, double, double> function, - Vector observedX, Vector observedY, Vector weight = null, - int accuracyOrder = 2) - { - var objective = new FittingObjectiveModel(function, null, accuracyOrder: accuracyOrder); - objective.SetObserved(observedX, observedY, weight); - return objective.ToObjectiveFunction(); - } - } -} diff --git a/src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs b/src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs deleted file mode 100644 index eb30f357..00000000 --- a/src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs +++ /dev/null @@ -1,730 +0,0 @@ -using MathNet.Numerics.LinearAlgebra; -using MathNet.Numerics.Optimization.ObjectiveFunctions; -using System; -using System.Collections.Generic; -using System.Linq; - -namespace MathNet.Numerics.Optimization.ObjectiveModels -{ - internal class FittingObjectiveModel : IObjectiveModel - { - #region Private Variables - - readonly Func, double, double> userFunction; // (p, x) => f(x; p) - readonly Func, double, Vector> userDerivatives; // (p, x) => df(x; p)/dp - readonly int accuracyOrder; // the desired accuracy order to evaluate the jacobian by numerical approximaiton. - - Vector coefficients; - Vector Pint; // internal(unbounded) coefficients - public Vector Pext; // external(bounded) coefficients - - bool hasFunctionValue; - double functionValue; // the residual sum of squares. Residuals * Residuals - Vector residuals; // the error values - - bool hasJacobianValue; - Matrix jacobianValue; // the Jacobian matrix. - Vector gradientValue; // the Gradient vector. - Matrix hessianValue; // the Hessian matrix. - - bool isBounded; - - #endregion Private Variables - - #region Public Variables - Observed Data - - /// - /// 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 the number of observations. - /// - public int NumberOfObservations { get { return (ObservedY == null) ? 0 : ObservedY.Count; } } - - #endregion Public Variables - Observed Data - - #region Public Variables - Bounds of Parameter - - /// - /// Get the values of the parameters. - /// - public List IsFixed { get; private set; } - - /// - /// Get the values of the parameters. - /// - public Vector LowerBound { get; private set; } - - /// - /// Get the values of the parameters. - /// - public Vector UpperBound { get; private set; } - - /// - /// Get the scale factor of the parameters. - /// - public Vector Scales { get; private set; } - - /// - /// Get the number of unknown parameters. - /// - public int NumberOfParameters { get { return (Point == null) ? 0 : Point.Count; } } - - #endregion Public Variables - Bounds of Parameter - - #region Public Variables - Others - - /// - /// 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 - Others - - public FittingObjectiveModel(Func, double, double>function, Func, double, Vector> derivatives = null, int accuracyOrder = 2) - { - this.userFunction = function; - this.userDerivatives = derivatives; - this.accuracyOrder = Math.Min(6, Math.Max(1, accuracyOrder)); - - IsFinished = false; - } - - public IObjectiveModel Fork() - { - return new FittingObjectiveModel(userFunction, userDerivatives, accuracyOrder) - { - ObservedX = ObservedX, - ObservedY = ObservedY, - Weights = Weights, - - coefficients = coefficients, - Pint = Pint, - Pext = Pext, - - hasFunctionValue = hasFunctionValue, - functionValue = functionValue, - - hasJacobianValue = hasJacobianValue, - jacobianValue = jacobianValue, - gradientValue = gradientValue, - hessianValue = hessianValue - }; - } - - public IObjectiveModel CreateNew() - { - return new FittingObjectiveModel(userFunction, userDerivatives, 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; - } - } - - /// - /// 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; - } - } - - public bool IsGradientSupported { get { return true; } } - public bool IsHessianSupported { get { return true; } } - - public bool IsFinished { get; set; } - - public IObjectiveFunction ToObjectiveFunction() - { - Tuple, Matrix> function(Vector point) - { - EvaluateAt(point); - - return new Tuple, Matrix>(Value, Gradient, Hessian); - } - - var objective = new GradientHessianObjectiveFunction(function); - return objective; - } - - /// - /// 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 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 initialGuess, Vector lowerBound = null, Vector upperBound = null, Vector scales = null, List isFixed = null) - { - if (initialGuess == null) - { - throw new ArgumentNullException("initialGuess"); - } - coefficients = initialGuess; - - 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 != initialGuess.Count) - { - throw new ArgumentException("The upper bounds can't have different elements from the initial guess."); - } - 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 != initialGuess.Count) - { - throw new ArgumentException("The upper bounds can't have different elements from the initial guess."); - } - 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 != initialGuess.Count) - { - throw new ArgumentException("The upper bounds can't have different elements from the initial guess."); - } - if (scales != null && scales.Count(x => x < 0) > 0) - { - scales.PointwiseAbs(); - } - Scales = scales; - - if (isFixed != null && isFixed.Count != initialGuess.Count) - { - throw new ArgumentException("The isFixed can't have different elements 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; - - isBounded = LowerBound != null || UpperBound != null || Scales != null; - } - - public void EvaluateAt(Vector parameters) - { - ValidateParameters(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. - - Pext = (FunctionEvaluations > 0 && isBounded) - ? ProjectParametersToExternal(parameters) - : parameters.Clone(); - Pint = (isBounded) - ? ProjectParametersToInternal(Pext) - : Pext; - - this.coefficients = Pint; - - if (IsFinished) - { - this.coefficients = Pext; - } - - hasFunctionValue = false; - hasJacobianValue = false; - - // don't keep references unnecessarily - jacobianValue = null; - gradientValue = null; - hessianValue = null; - } - - #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); - } - for (int i = 0; i < NumberOfObservations; i++) - { - ModelValues[i] = userFunction(Pext, ObservedX[i]); - } - 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 (userDerivatives != null) - { - // analytical jacobian - if (jacobianValue == null) - { - jacobianValue = Matrix.Build.Dense(NumberOfObservations, NumberOfParameters); - } - for (int i = 0; i < NumberOfObservations; i++) - { - jacobianValue.SetRow(i, userDerivatives(Pext, ObservedX[i])); - } - JacobianEvaluations++; - } - else - { - // numerical jacobian - jacobianValue = NumericalJacobian(Pext, ModelValues, accuracyOrder); - FunctionEvaluations += accuracyOrder; - } - - var scaleFactors = (isBounded && !IsFinished) - ? ScaleFactorsOfJacobian(Pint) - : Vector.Build.Dense(Pint.Count, 1.0); - - // project jacobian: 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 - jacobianValue[i, j] = 0.0; - } - else - { - jacobianValue[i, j] = jacobianValue[i, j] * scaleFactors[j]; - } - } - } - - // Gradient, g = -J'W(y − f(x; p)) = -J'L(L'E) = -J'LR - gradientValue = (Weights == null) - ? -jacobianValue.Transpose() * (ObservedY - ModelValues) - : -jacobianValue.Transpose() * Weights * (ObservedY - ModelValues); - - // approximated Hessian, H = J'WJ + ∑LRiHi ~ J'WJ near the minimum - hessianValue = (Weights == null) - ? jacobianValue.Transpose() * jacobianValue - : jacobianValue.Transpose() * Weights * jacobianValue; - } - - 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."); - } - } - - private Matrix NumericalJacobian(Vector Pext, Vector currentValues, int accuracyOrder = 2) - { - const double sqrtEpsilon = 1.4901161193847656250E-8; // sqrt(machineEpsilon) - - Matrix derivertives = Matrix.Build.Dense(NumberOfObservations, NumberOfParameters); - - var d = 0.000003 * Pext.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(Pext - 3 * h, x); - var f2 = userFunction(Pext - 2 * h, x); - var f3 = userFunction(Pext - h, x); - var f4 = userFunction(Pext + h, x); - var f5 = userFunction(Pext + 2 * h, x); - var f6 = userFunction(Pext + 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(Pext + h, x); - var f3 = userFunction(Pext + 2 * h, x); - var f4 = userFunction(Pext + 3 * h, x); - var f5 = userFunction(Pext + 4 * h, x); - var f6 = userFunction(Pext + 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(Pext - 2 * h, x); - var f2 = userFunction(Pext - h, x); - var f3 = userFunction(Pext + h, x); - var f4 = userFunction(Pext + 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(Pext + h, x); - var f3 = userFunction(Pext + 2 * h, x); - var f4 = userFunction(Pext + 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(Pext + h, x); - var f2 = userFunction(Pext - 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(Pext + h, x); - - var prime = (-f1 + f2) / h[j]; - derivertives[i, j] = prime; - } - - h[j] = 0; - } - } - - return derivertives; - } - - 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] = (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; - } - - 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] = (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; - } - - 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] = (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 Private Methods - } -} diff --git a/src/Numerics/Optimization/TrustRegionMinimizerBase.cs b/src/Numerics/Optimization/TrustRegionMinimizerBase.cs index d32c68f8..1ab8eb36 100644 --- a/src/Numerics/Optimization/TrustRegionMinimizerBase.cs +++ b/src/Numerics/Optimization/TrustRegionMinimizerBase.cs @@ -5,57 +5,32 @@ using System.Linq; namespace MathNet.Numerics.Optimization { - public abstract class TrustRegionMinimizerBase + public abstract class TrustRegionMinimizerBase : NonlinearMinimizerBase { - public static ITrustRegionSubproblem Subproblem; - - /// - /// 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. + /// The trust region subproblem. /// - 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; } + public static ITrustRegionSubproblem Subproblem; /// /// The stopping threshold for the trust region radius. /// public static double RadiusTolerance { get; set; } - /// - /// The maximum number of iterations. - /// - public int MaximumIterations { 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; - FunctionTolerance = functionTolerance; - GradientTolerance = gradientTolerance; - StepTolerance = stepTolerance; RadiusTolerance = radiusTolerance; - MaximumIterations = maximumIterations; } public NonlinearMinimizationResult FindMinimum(IObjectiveModel objective, Vector initialGuess, Vector lowerBound = null, Vector upperBound = null, Vector scales = null, List isFixed = null) { - if (objective == null) - throw new ArgumentNullException("objective"); - if (initialGuess == null) - throw new ArgumentNullException("initialGuess"); - return Minimum(Subproblem, objective, initialGuess, lowerBound, upperBound, scales, isFixed, GradientTolerance, StepTolerance, FunctionTolerance, RadiusTolerance, MaximumIterations); } @@ -63,11 +38,6 @@ namespace MathNet.Numerics.Optimization 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); @@ -133,23 +103,17 @@ namespace MathNet.Numerics.Optimization if (objective == null) throw new ArgumentNullException("objective"); - if (initialGuess == null) - throw new ArgumentNullException("initialGuess"); - objective.SetParameters(initialGuess, lowerBound, upperBound, scales, isFixed); + ValidateBounds(initialGuess, lowerBound, upperBound, scales); - ExitCondition exitCondition = ExitCondition.None; + objective.SetParameters(initialGuess, isFixed); - // Initialize objective - objective.FunctionEvaluations = 0; - objective.JacobianEvaluations = 0; - objective.IsFinished = false; + ExitCondition exitCondition = ExitCondition.None; // First, calculate function values and setup variables - objective.EvaluateAt(initialGuess); - var P = objective.Point; // current parameters - var RSS = objective.Value; // Residual Sum of Squares = R'R - var RSSinit = RSS; // RSS at initial gussing parameters + 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) { @@ -176,8 +140,9 @@ namespace MathNet.Numerics.Optimization } // evaluate projected gradient and Hessian - var Gradient = objective.Gradient; - var Hessian = objective.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) @@ -195,14 +160,16 @@ namespace MathNet.Numerics.Optimization 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); - var Pstep = subproblem.Pstep; - var hitBoundary = subproblem.HitBoundary; + 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); @@ -213,10 +180,9 @@ namespace MathNet.Numerics.Optimization } var Pnew = P + Pstep; // parameters to test - - objective.EvaluateAt(Pnew); - var RSSnew = objective.Value; - + // evaluate function at Pnew + var RSSnew = EvaluateFunction(objective, Pnew); + // if RSS == NaN, stop if (double.IsNaN(RSSnew)) { @@ -238,7 +204,7 @@ namespace MathNet.Numerics.Optimization delta = delta * 0.25; if (delta <= radiusTolerance * (radiusTolerance + P.DotProduct(P))) { - exitCondition = ExitCondition.RelativePoints; // SmallRelativeParameters + exitCondition = ExitCondition.LackOfProgress; break; } } @@ -250,8 +216,9 @@ namespace MathNet.Numerics.Optimization RSS = RSSnew; // evaluate projected gradient and Hessian - Gradient = objective.Gradient; - Hessian = objective.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)