From 220ab5aadb1550cedb1b50aadef70fe4b3cccd0f Mon Sep 17 00:00:00 2001 From: diluculo Date: Mon, 31 Dec 2018 17:04:41 +0900 Subject: [PATCH] Updated the box-constrained mapping rule in the FittingObjectiveModel. --- .../NonLinearCurveFittingTests.cs | 178 ++++++++++++------ .../ObjectiveModels/FittingObjectiveModel.cs | 52 +++-- 2 files changed, 155 insertions(+), 75 deletions(-) diff --git a/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs b/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs index fbbd3448..4918a263 100644 --- a/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs +++ b/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs @@ -42,11 +42,12 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests { var obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockPrime, RosenbrockX, RosenbrockY); var solver = new LevenbergMarquardtMinimizer(maximumIterations: 10000); - var result = solver.FindMinimum(obj, RosenbrockStart1); - AssertHelpers.AlmostEqualRelative(RosenbrockPbest[0], result.BestFitParameters[0], 3); - AssertHelpers.AlmostEqualRelative(RosenbrockPbest[1], result.BestFitParameters[1], 3); + for (int i = 0; i < result.BestFitParameters.Count; i++) + { + AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.BestFitParameters[i], 3); + } } [Test] @@ -54,11 +55,12 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests { var obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockX, RosenbrockY, accuracyOrder:2); var solver = new LevenbergMarquardtMinimizer(maximumIterations: 10000); - var result = solver.FindMinimum(obj, RosenbrockStart1); - - AssertHelpers.AlmostEqualRelative(RosenbrockPbest[0], result.BestFitParameters[0], 3); - AssertHelpers.AlmostEqualRelative(RosenbrockPbest[1], result.BestFitParameters[1], 3); + + for (int i = 0; i < result.BestFitParameters.Count; i++) + { + AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.BestFitParameters[i], 3); + } } [Test] @@ -67,11 +69,12 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests var obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockPrime, RosenbrockX, RosenbrockY, lowerBound : RosebbrockLowerBound, upperBound : RosenbrockUpperBound); var solver = new LevenbergMarquardtMinimizer(maximumIterations: 10000); - var result = solver.FindMinimum(obj, RosenbrockStart1); - AssertHelpers.AlmostEqualRelative(RosenbrockPbest[0], result.BestFitParameters[0], 3); - AssertHelpers.AlmostEqualRelative(RosenbrockPbest[1], result.BestFitParameters[1], 3); + for (int i = 0; i < result.BestFitParameters.Count; i++) + { + AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.BestFitParameters[i], 3); + } } [Test] @@ -79,13 +82,14 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests { var obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockX, RosenbrockY, lowerBound: RosebbrockLowerBound, upperBound: RosenbrockUpperBound, - accuracyOrder: 2); + accuracyOrder: 6); var solver = new LevenbergMarquardtMinimizer(maximumIterations: 10000); - var result = solver.FindMinimum(obj, RosenbrockStart1); - AssertHelpers.AlmostEqualRelative(RosenbrockPbest[0], result.BestFitParameters[0], 3); - AssertHelpers.AlmostEqualRelative(RosenbrockPbest[1], result.BestFitParameters[1], 3); + for (int i = 0; i < result.BestFitParameters.Count; i++) + { + AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.BestFitParameters[i], 3); + } } [Test] @@ -93,11 +97,12 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests { var obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockPrime, RosenbrockX, RosenbrockY); var solver = new TrustRegionDogLegMinimizer(maximumIterations: 10000); - var result = solver.FindMinimum(obj, RosenbrockStart1); - AssertHelpers.AlmostEqualRelative(RosenbrockPbest[0], result.BestFitParameters[0], 1); - AssertHelpers.AlmostEqualRelative(RosenbrockPbest[1], result.BestFitParameters[1], 1); + for (int i = 0; i < result.BestFitParameters.Count; i++) + { + AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.BestFitParameters[i], 1); + } } [Test] @@ -105,11 +110,12 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests { var obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockPrime, RosenbrockX, RosenbrockY); var solver = new TrustRegionDogLegMinimizer(maximumIterations: 10000); - var result = solver.FindMinimum(obj, RosenbrockStart1); - AssertHelpers.AlmostEqualRelative(RosenbrockPbest[0], result.BestFitParameters[0], 1); - AssertHelpers.AlmostEqualRelative(RosenbrockPbest[1], result.BestFitParameters[1], 1); + for (int i = 0; i < result.BestFitParameters.Count; i++) + { + AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.BestFitParameters[i], 1); + } } // model: Rat43 (https://www.itl.nist.gov/div898/strd/nls/data/ratkowsky3.shtml) @@ -147,7 +153,6 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests { var obj = ObjectiveModel.FittingModel(Rat43Model, Rat43X, Rat43Y, accuracyOrder: 6); var solver = new LevenbergMarquardtMinimizer(); - var result = solver.FindMinimum(obj, Rat43Start1); for (int i = 0; i < result.BestFitParameters.Count; i++) @@ -162,7 +167,6 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests { var obj = ObjectiveModel.FittingModel(Rat43Model, Rat43X, Rat43Y, accuracyOrder: 6); var solver = new TrustRegionDogLegMinimizer(); - var result = solver.FindMinimum(obj, Rat43Start2); for (int i = 0; i < result.BestFitParameters.Count; i++) @@ -199,23 +203,22 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests private Vector BoxBodStart1 = new DenseVector(new double[] { 1.0, 1.0 }); private Vector BoxBodStart2 = new DenseVector(new double[] { 100.0, 0.75 }); - private Vector BoxBodLowerBound = new DenseVector(new double[] { 0, 0 }); - private Vector BoxBodUpperBound = new DenseVector(new double[] { 500.0, 10 }); - private Vector BoxBodScales = new DenseVector(new double[] { 100.0, 1 }); + 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 LMDER_FindMinimum_BoxBod_Unconstrained() { var obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY); var solver = new LevenbergMarquardtMinimizer(); - var result = solver.FindMinimum(obj, BoxBodStart1); - AssertHelpers.AlmostEqualRelative(BoxBodPbest[0], result.BestFitParameters[0], 6); - AssertHelpers.AlmostEqualRelative(BoxBodPbest[1], result.BestFitParameters[1], 6); - - AssertHelpers.AlmostEqualRelative(BoxBodPstd[0], result.StandardErrors[0], 6); - AssertHelpers.AlmostEqualRelative(BoxBodPstd[1], result.StandardErrors[1], 6); + for (int i = 0; i < result.BestFitParameters.Count; i++) + { + AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6); + AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6); + } } [Test] @@ -223,30 +226,96 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests { var obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodX, BoxBodY, accuracyOrder:6); var solver = new LevenbergMarquardtMinimizer(); - var result = solver.FindMinimum(obj, BoxBodStart1); - AssertHelpers.AlmostEqualRelative(BoxBodPbest[0], result.BestFitParameters[0], 6); - AssertHelpers.AlmostEqualRelative(BoxBodPbest[1], result.BestFitParameters[1], 6); - - AssertHelpers.AlmostEqualRelative(BoxBodPstd[0], result.StandardErrors[0], 6); - AssertHelpers.AlmostEqualRelative(BoxBodPstd[1], result.StandardErrors[1], 6); + for (int i = 0; i < result.BestFitParameters.Count; i++) + { + AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6); + AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6); + } } [Test] public void LMDER_FindMinimum_BoxBod_BoxConstrained() { + // lower < parameters < upper + // Note that in this case, scales have no effect. + var obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY, lowerBound: BoxBodLowerBound, upperBound: BoxBodUpperBound); var solver = new LevenbergMarquardtMinimizer(); - var result = solver.FindMinimum(obj, BoxBodStart1); - AssertHelpers.AlmostEqualRelative(BoxBodPbest[0], result.BestFitParameters[0], 6); - AssertHelpers.AlmostEqualRelative(BoxBodPbest[1], result.BestFitParameters[1], 6); + for (int i = 0; i < result.BestFitParameters.Count; i++) + { + AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6); + AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6); + } - AssertHelpers.AlmostEqualRelative(BoxBodPstd[0], result.StandardErrors[0], 6); - AssertHelpers.AlmostEqualRelative(BoxBodPstd[1], result.StandardErrors[1], 6); + // lower < parameters, no scales + + obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY, + lowerBound: BoxBodLowerBound); + solver = new LevenbergMarquardtMinimizer(); + result = solver.FindMinimum(obj, BoxBodStart1); + + for (int i = 0; i < result.BestFitParameters.Count; i++) + { + AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6); + AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6); + } + + // lower < parameters, scales + + obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY, + lowerBound: BoxBodLowerBound, scales: BoxBodScales); + solver = new LevenbergMarquardtMinimizer(); + result = solver.FindMinimum(obj, BoxBodStart1); + + for (int i = 0; i < result.BestFitParameters.Count; i++) + { + AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6); + AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6); + } + + // parameters < upper, no scales + + obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY, + upperBound: BoxBodUpperBound); + solver = new LevenbergMarquardtMinimizer(); + result = solver.FindMinimum(obj, BoxBodStart1); + + for (int i = 0; i < result.BestFitParameters.Count; i++) + { + AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6); + AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6); + } + + // parameters < upper, scales + + obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY, + upperBound: BoxBodUpperBound, scales: BoxBodScales); + solver = new LevenbergMarquardtMinimizer(); + result = solver.FindMinimum(obj, BoxBodStart1); + + for (int i = 0; i < result.BestFitParameters.Count; i++) + { + AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6); + AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6); + } + + // only scales + + obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY, + scales: BoxBodScales); + solver = new LevenbergMarquardtMinimizer(); + result = solver.FindMinimum(obj, BoxBodStart1); + + for (int i = 0; i < result.BestFitParameters.Count; i++) + { + AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6); + AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6); + } } [Test] @@ -256,14 +325,13 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests lowerBound: BoxBodLowerBound, upperBound: BoxBodUpperBound, accuracyOrder: 6); var solver = new LevenbergMarquardtMinimizer(); - var result = solver.FindMinimum(obj, BoxBodStart1); - AssertHelpers.AlmostEqualRelative(BoxBodPbest[0], result.BestFitParameters[0], 6); - AssertHelpers.AlmostEqualRelative(BoxBodPbest[1], result.BestFitParameters[1], 6); - - AssertHelpers.AlmostEqualRelative(BoxBodPstd[0], result.StandardErrors[0], 6); - AssertHelpers.AlmostEqualRelative(BoxBodPstd[1], result.StandardErrors[1], 6); + for (int i = 0; i < result.BestFitParameters.Count; i++) + { + AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 6); + AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 6); + } } [Test] @@ -271,14 +339,13 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests { var obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodX, BoxBodY, accuracyOrder: 6); var solver = new TrustRegionDogLegMinimizer(); - var result = solver.FindMinimum(obj, BoxBodStart1); - AssertHelpers.AlmostEqualRelative(BoxBodPbest[0], result.BestFitParameters[0], 3); - AssertHelpers.AlmostEqualRelative(BoxBodPbest[1], result.BestFitParameters[1], 3); - - AssertHelpers.AlmostEqualRelative(BoxBodPstd[0], result.StandardErrors[0], 3); - AssertHelpers.AlmostEqualRelative(BoxBodPstd[1], result.StandardErrors[1], 3); + for (int i = 0; i < result.BestFitParameters.Count; i++) + { + AssertHelpers.AlmostEqualRelative(BoxBodPbest[i], result.BestFitParameters[i], 3); + AssertHelpers.AlmostEqualRelative(BoxBodPstd[i], result.StandardErrors[i], 3); + } } // model : Thurber (https://www.itl.nist.gov/div898/strd/nls/data/thurber.shtml) @@ -359,7 +426,6 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests { var obj = ObjectiveModel.FittingModel(ThurberModel, ThurberPrime, ThurberX, ThurberY); var solver = new LevenbergMarquardtMinimizer(); - var result = solver.FindMinimum(obj, ThurberInitialGuess); for (int i = 0; i < result.BestFitParameters.Count; i++) @@ -374,7 +440,6 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests { var obj = ObjectiveModel.FittingModel(ThurberModel, ThurberX, ThurberY, accuracyOrder: 6); var solver = new LevenbergMarquardtMinimizer(); - var result = solver.FindMinimum(obj, ThurberInitialGuess); for (int i = 0; i < result.BestFitParameters.Count; i++) @@ -391,7 +456,6 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests scales: ThurberScales, accuracyOrder: 6); var solver = new TrustRegionDogLegMinimizer(); - var result = solver.FindMinimum(obj, ThurberInitialGuess); for (int i = 0; i < result.BestFitParameters.Count; i++) diff --git a/src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs b/src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs index a1ea2beb..ce81fb54 100644 --- a/src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs +++ b/src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs @@ -329,30 +329,34 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels ValidateParameters(parameters); // To handle the box constrained minimization as the unconstrained minimization, - // parameters are mapping by the following rule. + // 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 < P < upper + // 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 < P - // Pint = sqrt((Pext - lower + 1)^2 - 1) - // Pext = lower - 1 + sqrt(Pint^2 + 1) - // dPext/dPint = Pint / sqrt(Pint^2 + 1) - // 3. P < upper - // Pint = sqrt((upper - Pext + 1)^2 - 1) - // Pext = upper + 1 - sqrt(Pint^2 + 1) - // dPext/dPint = - Pint / sqrt(Pint^2 + 1) + // + // 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 // - // see ProjectParametersToInternal(Pext), ProjectParametersToExternal(Pint), ScaleFactorsOfJacobian(Pint) methods. + // 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. @@ -617,7 +621,9 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels { for (int i = 0; i < Pext.Count; i++) { - Pint[i] = Math.Sqrt(Math.Pow(Pext[i] - LowerBound[i] + 1.0, 2) - 1.0); + 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; @@ -626,7 +632,9 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels { for (int i = 0; i < Pext.Count; i++) { - Pint[i] = Math.Sqrt(Math.Pow(UpperBound[i] - Pext[i] + 1.0, 2) - 1.0); + 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; @@ -661,7 +669,9 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels { for (int i = 0; i < Pint.Count; i++) { - Pext[i] = LowerBound[i] + Math.Sqrt(Pint[i] * Pint[i] + 1.0) - 1.0; + 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; @@ -670,7 +680,9 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels { for (int i = 0; i < Pint.Count; i++) { - Pext[i] = UpperBound[i] - Math.Sqrt(Pint[i] * Pint[i] + 1.0) + 1.0; + 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; @@ -704,7 +716,9 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels { for (int i = 0; i < Pint.Count; i++) { - scale[i] = Pint[i] / Math.Sqrt(Pint[i] * Pint[i] + 1.0); + 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; } @@ -712,7 +726,9 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels { for (int i = 0; i < Pint.Count; i++) { - scale[i] = -Pint[i] / Math.Sqrt(Pint[i] * Pint[i] + 1.0); + 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; }