diff --git a/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs b/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs index 4918a263..f0787922 100644 --- a/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs +++ b/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs @@ -93,8 +93,9 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests } [Test] - public void TRLMDER_FindMinimum_Rosenbrock_Unconstrained() + public void TRDLDER_FindMinimum_Rosenbrock_Unconstrained() { + // DogLeg Minimizer var obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockPrime, RosenbrockX, RosenbrockY); var solver = new TrustRegionDogLegMinimizer(maximumIterations: 10000); var result = solver.FindMinimum(obj, RosenbrockStart1); @@ -103,10 +104,20 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests { AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.BestFitParameters[i], 1); } + + // NewtonCG Minimizer + obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockPrime, RosenbrockX, RosenbrockY); + var solverNCG = new TrustRegionNewtonCGMinimizer(maximumIterations: 10000); + result = solverNCG.FindMinimum(obj, RosenbrockStart1); + + for (int i = 0; i < result.BestFitParameters.Count; i++) + { + AssertHelpers.AlmostEqualRelative(RosenbrockPbest[i], result.BestFitParameters[i], 1); + } } [Test] - public void TRLMDIF_FindMinimum_Rosenbrock_Unconstrained() + public void TRDLDIF_FindMinimum_Rosenbrock_Unconstrained() { var obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockPrime, RosenbrockX, RosenbrockY); var solver = new TrustRegionDogLegMinimizer(maximumIterations: 10000); @@ -118,6 +129,19 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests } } + [Test] + public void TRNCGDER_FindMinimum_Rosenbrock_Unconstrained() + { + var obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockPrime, RosenbrockX, RosenbrockY); + var solver = new TrustRegionNewtonCGMinimizer(maximumIterations: 10000); + var result = solver.FindMinimum(obj, RosenbrockStart1); + + 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) // f(x; a, b, c, d) = a / ((1 + exp(b - c * x))^(1 / d)) // best fitted parameters: @@ -335,7 +359,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests } [Test] - public void TRLMDIF_FindMinimum_BoxBod_Unconstrained() + public void TRDLDIF_FindMinimum_BoxBod_Unconstrained() { var obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodX, BoxBodY, accuracyOrder: 6); var solver = new TrustRegionDogLegMinimizer(); @@ -348,6 +372,20 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests } } + [Test] + public void TRNCGDIF_FindMinimum_BoxBod_Unconstrained() + { + var obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodX, BoxBodY, accuracyOrder: 6); + var solver = new TrustRegionNewtonCGMinimizer(); + var result = solver.FindMinimum(obj, BoxBodStart1); + + 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) // f(x; b1 ... b7) = (b1 + b2*x + b3*x^2 + b4*x^3) / (1 + b5*x + b6*x^2 + b7*x^3) // derivatives: @@ -450,7 +488,7 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests } [Test] - public void TRLMDIF_FindMinimum_Thurber_Scaled() + public void TRDLDIF_FindMinimum_Thurber_Scaled() { var obj = ObjectiveModel.FittingModel(ThurberModel, ThurberX, ThurberY, scales: ThurberScales, @@ -464,5 +502,21 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests AssertHelpers.AlmostEqualRelative(ThurberPstd[i], result.StandardErrors[i], 3); } } + + [Test] + public void TRNCGDIF_FindMinimum_Thurber_Scaled() + { + var obj = ObjectiveModel.FittingModel(ThurberModel, ThurberX, ThurberY, + scales: ThurberScales, + accuracyOrder: 6); + var solver = new TrustRegionNewtonCGMinimizer(); + var result = solver.FindMinimum(obj, ThurberInitialGuess); + + for (int i = 0; i < result.BestFitParameters.Count; i++) + { + AssertHelpers.AlmostEqualRelative(ThurberPbest[i], result.BestFitParameters[i], 3); + AssertHelpers.AlmostEqualRelative(ThurberPstd[i], result.StandardErrors[i], 3); + } + } } } diff --git a/src/Numerics/Optimization/Subproblems/DogLegSubproblem.cs b/src/Numerics/Optimization/Subproblems/DogLegSubproblem.cs index e3f06d15..8d41288d 100644 --- a/src/Numerics/Optimization/Subproblems/DogLegSubproblem.cs +++ b/src/Numerics/Optimization/Subproblems/DogLegSubproblem.cs @@ -13,7 +13,6 @@ namespace MathNet.Numerics.Optimization.Subproblems var Jacobian = objective.Jacobian; var Gradient = objective.Gradient; var Hessian = objective.Hessian; - var RSS = objective.Residue; // newton point // the Gauss–Newton step by solving the normal equations diff --git a/src/Numerics/Optimization/Subproblems/NewtonCGSubproblem.cs b/src/Numerics/Optimization/Subproblems/NewtonCGSubproblem.cs new file mode 100644 index 00000000..ba6f9a6c --- /dev/null +++ b/src/Numerics/Optimization/Subproblems/NewtonCGSubproblem.cs @@ -0,0 +1,65 @@ +using MathNet.Numerics.LinearAlgebra; +using System; + +namespace MathNet.Numerics.Optimization.Subproblems +{ + internal class NewtonCGSubproblem : ITrustRegionSubproblem + { + public Vector Pstep { get; private set; } + + public bool HitBoundary { get; private set; } + + public void Solve(IObjectiveModel objective, double delta) + { + var Jacobian = objective.Jacobian; + var Gradient = objective.Gradient; + var Hessian = objective.Hessian; + + // define tolerance + var tolerance = Math.Min(0.5, Math.Sqrt(Gradient.L2Norm())) * Gradient.L2Norm(); + + // initialize internal variables + var z = Vector.Build.Dense(Hessian.RowCount); + var r = -Gradient; + var d = -r; + + while (true) + { + var Bd = Hessian * d; + var dBd = d.DotProduct(Bd); + + if (dBd <= 0) + { + var t = Util.FindBeta(1, z, d, delta); + Pstep = z + t.Item1 * d; + HitBoundary = true; + return; + } + + var r_sq = r.DotProduct(r); + var alpha = r_sq / dBd; + var znext = z + alpha * d; + if(znext.L2Norm() >= delta) + { + var t = Util.FindBeta(1, z, d, delta); + Pstep = z + t.Item2 * d; + HitBoundary = true; + return; + } + + var rnext = r + alpha * Bd; + var rnext_sq = rnext.DotProduct(rnext); + if (Math.Sqrt(rnext_sq) < tolerance) + { + Pstep = znext; + HitBoundary = false; + return; + } + + z = znext; + r = rnext; + d = -rnext + rnext_sq / r_sq * d; ; + } + } + } +} diff --git a/src/Numerics/Optimization/TrustRegionDogLegMinimizer.cs b/src/Numerics/Optimization/TrustRegionDogLegMinimizer.cs index 596f7637..77a220f0 100644 --- a/src/Numerics/Optimization/TrustRegionDogLegMinimizer.cs +++ b/src/Numerics/Optimization/TrustRegionDogLegMinimizer.cs @@ -1,10 +1,10 @@ -using MathNet.Numerics.LinearAlgebra; -using System; - -namespace MathNet.Numerics.Optimization +namespace MathNet.Numerics.Optimization { public sealed class TrustRegionDogLegMinimizer : TrustRegionMinimizerBase { + /// + /// Non-linear least square fitting by the trust region dogleg algorithm. + /// public TrustRegionDogLegMinimizer(double gradientTolerance = 1E-8, double stepTolerance = 1E-8, double functionTolerance = 1E-8, double radiusTolerance = 1E-8, int maximumIterations = -1) : base(TrustRegionSubproblem.DogLeg(), gradientTolerance, stepTolerance, functionTolerance, radiusTolerance, maximumIterations) { } diff --git a/src/Numerics/Optimization/TrustRegionNewtonCGMinimizer.cs b/src/Numerics/Optimization/TrustRegionNewtonCGMinimizer.cs new file mode 100644 index 00000000..7e6d5f5a --- /dev/null +++ b/src/Numerics/Optimization/TrustRegionNewtonCGMinimizer.cs @@ -0,0 +1,13 @@ +namespace MathNet.Numerics.Optimization +{ + public sealed class TrustRegionNewtonCGMinimizer : TrustRegionMinimizerBase + { + /// + /// Non-linear least square fitting by the trust region Newton-Conjugate-Gradient algorithm. + /// + public TrustRegionNewtonCGMinimizer(double gradientTolerance = 1E-8, double stepTolerance = 1E-8, double functionTolerance = 1E-8, double radiusTolerance = 1E-8, int maximumIterations = -1) + : base(TrustRegionSubproblem.NewtonCG(), gradientTolerance, stepTolerance, functionTolerance, radiusTolerance, maximumIterations) + { } + } +} + diff --git a/src/Numerics/Optimization/TrustRegionSubProblem.cs b/src/Numerics/Optimization/TrustRegionSubProblem.cs index cde8ecf8..3018f8b6 100644 --- a/src/Numerics/Optimization/TrustRegionSubProblem.cs +++ b/src/Numerics/Optimization/TrustRegionSubProblem.cs @@ -8,5 +8,10 @@ namespace MathNet.Numerics.Optimization { return new DogLegSubproblem(); } + + public static ITrustRegionSubproblem NewtonCG() + { + return new NewtonCGSubproblem(); + } } }