From 22d9041d59d48c667d0fb605ecbec9ad34e7e392 Mon Sep 17 00:00:00 2001 From: diluculo Date: Thu, 3 Jan 2019 21:52:52 +0900 Subject: [PATCH] Reworked ObjectiveModel. --- .../NonLinearCurveFittingTests.cs | 64 +-- src/Numerics/Optimization/IObjectiveModel.cs | 36 +- .../LevenbergMarquardtMinimizer.cs | 45 +- .../NonlinearMinimizationResult.cs | 54 +- src/Numerics/Optimization/ObjectiveModel.cs | 12 +- .../ObjectiveModels/FittingObjectiveModel.cs | 500 ++++++++---------- .../Subproblems/DogLegSubproblem.cs | 6 +- .../Subproblems/NewtonCGSubproblem.cs | 2 +- .../Optimization/TrustRegionMinimizerBase.cs | 64 ++- 9 files changed, 376 insertions(+), 407 deletions(-) diff --git a/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs b/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs index e1fc339f..faae726c 100644 --- a/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs +++ b/src/Numerics.Tests/OptimizationTests/NonLinearCurveFittingTests.cs @@ -55,10 +55,9 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests } // box constrained - obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockPrime, RosenbrockX, RosenbrockY, - lowerBound: RosebbrockLowerBound, upperBound: RosenbrockUpperBound); + obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockPrime, RosenbrockX, RosenbrockY); solver = new LevenbergMarquardtMinimizer(maximumIterations: 10000); - result = solver.FindMinimum(obj, RosenbrockStart1); + result = solver.FindMinimum(obj, RosenbrockStart1, lowerBound: RosebbrockLowerBound, upperBound: RosenbrockUpperBound); for (int i = 0; i < result.BestFitParameters.Count; i++) { @@ -80,11 +79,9 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests } // box constrained - obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockX, RosenbrockY, - lowerBound: RosebbrockLowerBound, upperBound: RosenbrockUpperBound, - accuracyOrder: 6); + obj = ObjectiveModel.FittingModel(RosenbrockModel, RosenbrockX, RosenbrockY, accuracyOrder: 6); solver = new LevenbergMarquardtMinimizer(maximumIterations: 10000); - result = solver.FindMinimum(obj, RosenbrockStart1); + result = solver.FindMinimum(obj, RosenbrockStart1, lowerBound: RosebbrockLowerBound, upperBound: RosenbrockUpperBound); for (int i = 0; i < result.BestFitParameters.Count; i++) { @@ -272,10 +269,9 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests // lower < parameters < upper // Note that in this case, scales have no effect. - obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY, - lowerBound: BoxBodLowerBound, upperBound: BoxBodUpperBound); + obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY); solver = new LevenbergMarquardtMinimizer(); - result = solver.FindMinimum(obj, BoxBodStart1); + result = solver.FindMinimum(obj, BoxBodStart1, lowerBound: BoxBodLowerBound, upperBound: BoxBodUpperBound); for (int i = 0; i < result.BestFitParameters.Count; i++) { @@ -285,10 +281,9 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests // lower < parameters, no scales - obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY, - lowerBound: BoxBodLowerBound); + obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY); solver = new LevenbergMarquardtMinimizer(); - result = solver.FindMinimum(obj, BoxBodStart1); + result = solver.FindMinimum(obj, BoxBodStart1, lowerBound: BoxBodLowerBound); for (int i = 0; i < result.BestFitParameters.Count; i++) { @@ -298,10 +293,9 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests // lower < parameters, scales - obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY, - lowerBound: BoxBodLowerBound, scales: BoxBodScales); + obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY); solver = new LevenbergMarquardtMinimizer(); - result = solver.FindMinimum(obj, BoxBodStart1); + result = solver.FindMinimum(obj, BoxBodStart1, lowerBound: BoxBodLowerBound, scales: BoxBodScales); for (int i = 0; i < result.BestFitParameters.Count; i++) { @@ -311,10 +305,9 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests // parameters < upper, no scales - obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY, - upperBound: BoxBodUpperBound); + obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY); solver = new LevenbergMarquardtMinimizer(); - result = solver.FindMinimum(obj, BoxBodStart1); + result = solver.FindMinimum(obj, BoxBodStart1, upperBound: BoxBodUpperBound); for (int i = 0; i < result.BestFitParameters.Count; i++) { @@ -324,10 +317,9 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests // parameters < upper, scales - obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY, - upperBound: BoxBodUpperBound, scales: BoxBodScales); + obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY); solver = new LevenbergMarquardtMinimizer(); - result = solver.FindMinimum(obj, BoxBodStart1); + result = solver.FindMinimum(obj, BoxBodStart1, upperBound: BoxBodUpperBound, scales: BoxBodScales); for (int i = 0; i < result.BestFitParameters.Count; i++) { @@ -337,10 +329,9 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests // only scales - obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY, - scales: BoxBodScales); + obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY); solver = new LevenbergMarquardtMinimizer(); - result = solver.FindMinimum(obj, BoxBodStart1); + result = solver.FindMinimum(obj, BoxBodStart1, scales: BoxBodScales); for (int i = 0; i < result.BestFitParameters.Count; i++) { @@ -364,11 +355,9 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests } // box constrained - obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodX, BoxBodY, - lowerBound: BoxBodLowerBound, upperBound: BoxBodUpperBound, - accuracyOrder: 6); + obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodX, BoxBodY, accuracyOrder: 6); solver = new LevenbergMarquardtMinimizer(); - result = solver.FindMinimum(obj, BoxBodStart1); + result = solver.FindMinimum(obj, BoxBodStart1, lowerBound: BoxBodLowerBound, upperBound: BoxBodUpperBound); for (int i = 0; i < result.BestFitParameters.Count; i++) { @@ -394,9 +383,10 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests [Test] public void BoxBod_TRNCG_Dif() { - var obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodX, BoxBodY, accuracyOrder: 6); + var obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodPrime, BoxBodX, BoxBodY); + //var obj = ObjectiveModel.FittingModel(BoxBodModel, BoxBodX, BoxBodY, accuracyOrder: 6); var solver = new TrustRegionNewtonCGMinimizer(); - var result = solver.FindMinimum(obj, BoxBodStart1); + var result = solver.FindMinimum(obj, BoxBodStart2); for (int i = 0; i < result.BestFitParameters.Count; i++) { @@ -528,11 +518,9 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests [Test] public void Thurber_TRDL_Dif() { - var obj = ObjectiveModel.FittingModel(ThurberModel, ThurberX, ThurberY, - scales: ThurberScales, - accuracyOrder: 6); + var obj = ObjectiveModel.FittingModel(ThurberModel, ThurberX, ThurberY, accuracyOrder: 6); var solver = new TrustRegionDogLegMinimizer(); - var result = solver.FindMinimum(obj, ThurberStart); + var result = solver.FindMinimum(obj, ThurberStart, scales: ThurberScales); for (int i = 0; i < result.BestFitParameters.Count; i++) { @@ -544,11 +532,9 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests [Test] public void Thurber_TRNCG_Dif() { - var obj = ObjectiveModel.FittingModel(ThurberModel, ThurberX, ThurberY, - scales: ThurberScales, - accuracyOrder: 6); + var obj = ObjectiveModel.FittingModel(ThurberModel, ThurberX, ThurberY, accuracyOrder: 6); var solver = new TrustRegionNewtonCGMinimizer(); - var result = solver.FindMinimum(obj, ThurberStart); + var result = solver.FindMinimum(obj, ThurberStart, scales: ThurberScales); for (int i = 0; i < result.BestFitParameters.Count; i++) { diff --git a/src/Numerics/Optimization/IObjectiveModel.cs b/src/Numerics/Optimization/IObjectiveModel.cs index ff117091..b106a347 100644 --- a/src/Numerics/Optimization/IObjectiveModel.cs +++ b/src/Numerics/Optimization/IObjectiveModel.cs @@ -1,4 +1,5 @@ using MathNet.Numerics.LinearAlgebra; +using System.Collections.Generic; namespace MathNet.Numerics.Optimization { @@ -19,44 +20,33 @@ namespace MathNet.Numerics.Optimization /// /// Get the y-values of the fitted model that correspond to the independent values. /// - Vector Values { get; } + Vector ModelValues { get; } /// /// Get the values of the parameters. /// - Vector Parameters { get; } + Vector Point { get; } /// /// Get the residual sum of squares. /// - double Residue { get; } + double Value { get; } - /// - /// Get the Jacobian matrix, J(x; p) = df(x; p)/dp. - /// - Matrix Jacobian { get; } /// /// Get the Gradient vector. G = J'(y - f(x; p)) /// Vector Gradient { get; } + /// /// Get the approximated Hessian matrix. H = J'J /// Matrix Hessian { get; } - /// - /// Get the covariance matrix. - /// - Matrix Covariance { get; } - - /// - /// Get the correlation matrix. - /// - Matrix Correlation { get; } /// /// Get the number of calls to function. /// int FunctionEvaluations { get; set; } + /// /// Get the number of calls to jacobian. /// @@ -67,17 +57,17 @@ namespace MathNet.Numerics.Optimization /// int DegreeOfFreedom { get; } - /// - /// Get whether or not the analytical jacobian is supported. - /// - bool IsJacobianSupported { get; } + bool IsGradientSupported { get; } + bool IsHessianSupported { get; } + + bool IsFinished { get; set; } } public interface IObjectiveModel : IObjectiveModelEvaluation { - void EvaluateFunction(Vector parameters); - void EvaluateJacobian(Vector parameters); - void EvaluateCovariance(Vector parameters); + void SetParameters(Vector initialGuess, Vector lowerBound = null, Vector upperBound = null, Vector scales = null, List isFixed = null); + + void EvaluateAt(Vector parameters); /// Create a new independent copy of this objective function, evaluated at the same point. IObjectiveModel Fork(); diff --git a/src/Numerics/Optimization/LevenbergMarquardtMinimizer.cs b/src/Numerics/Optimization/LevenbergMarquardtMinimizer.cs index ab6bb421..c39b0119 100644 --- a/src/Numerics/Optimization/LevenbergMarquardtMinimizer.cs +++ b/src/Numerics/Optimization/LevenbergMarquardtMinimizer.cs @@ -1,5 +1,6 @@ using MathNet.Numerics.LinearAlgebra; using System; +using System.Collections.Generic; using System.Linq; namespace MathNet.Numerics.Optimization @@ -44,24 +45,31 @@ namespace MathNet.Numerics.Optimization MaximumIterations = maximumIterations; } - public NonlinearMinimizationResult FindMinimum(IObjectiveModel objective, Vector initialGuess) + 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, InitialMu, FunctionTolerance, GradientTolerance, StepTolerance, MaximumIterations); + return Minimum(objective, initialGuess, lowerBound, upperBound, scales, isFixed, InitialMu, FunctionTolerance, GradientTolerance, StepTolerance, MaximumIterations); } - public NonlinearMinimizationResult FindMinimum(IObjectiveModel objective, double[] initialGuess) + 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"); - return Minimum(objective, CreateVector.DenseOfArray(initialGuess), InitialMu, GradientTolerance, StepTolerance, FunctionTolerance, MaximumIterations); + var lb = (lowerBound == null) ? null : CreateVector.Dense(lowerBound); + var ub = (upperBound == null) ? null : CreateVector.Dense(upperBound); + var sc = (scales == null) ? null : CreateVector.Dense(scales); + var fx = (isFixed == null) ? null : isFixed.ToList(); + + return Minimum(objective, CreateVector.DenseOfArray(initialGuess), lb, ub, sc, fx, InitialMu, GradientTolerance, StepTolerance, FunctionTolerance, MaximumIterations); } /// @@ -75,7 +83,9 @@ namespace MathNet.Numerics.Optimization /// The stopping threshold for L2 norm of the residuals. /// The max iterations. /// The result of the Levenberg-Marquardt minimization - public static NonlinearMinimizationResult Minimum(IObjectiveModel objective, Vector initialGuess, double initialMu = 1E-3, double gradientTolerance = 1E-18, double stepTolerance = 1E-18, double functionTolerance = 1E-18, int maximumIterations = -1) + 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) { // Non-linear least square fitting by the Levenberg-Marduardt algorithm. // @@ -118,23 +128,24 @@ namespace MathNet.Numerics.Optimization if (initialGuess == null) throw new ArgumentNullException("initialGuess"); + objective.SetParameters(initialGuess, lowerBound, upperBound, scales, isFixed); + ExitCondition exitCondition = ExitCondition.None; // Initialize objective objective.FunctionEvaluations = 0; objective.JacobianEvaluations = 0; + objective.IsFinished = false; // First, calculate function values and setup variables - objective.EvaluateFunction(initialGuess); - var P = objective.Parameters; // current parameters + objective.EvaluateAt(initialGuess); + var P = objective.Point; // current parameters var Pstep = Vector.Build.Dense(P.Count); // the change of parameters - var RSS = objective.Residue; // Residual Sum of Squares = R'R + var RSS = objective.Value; // Residual Sum of Squares = R'R if (maximumIterations < 0) { - maximumIterations = (objective.IsJacobianSupported) - ? 100 * (initialGuess.Count + 1) - : 200 * (initialGuess.Count + 1); + maximumIterations = 200 * (initialGuess.Count + 1); } // if RSS == NaN, stop @@ -157,9 +168,8 @@ namespace MathNet.Numerics.Optimization } // Evaluate gradient and Hessian - objective.EvaluateJacobian(P); var Gradient = objective.Gradient; - var Hessian = objective.Hessian; + var Hessian = objective.Hessian; var diagonalOfHessian = Hessian.Diagonal(); // diag(H) // if ||g||oo <= gtol, found and stop @@ -170,7 +180,6 @@ namespace MathNet.Numerics.Optimization if (exitCondition != ExitCondition.None) { - objective.EvaluateCovariance(P); return new NonlinearMinimizationResult(objective, -1, exitCondition); } @@ -197,8 +206,8 @@ namespace MathNet.Numerics.Optimization var Pnew = P + Pstep; // new parameters to test - objective.EvaluateFunction(Pnew); - var RSSnew = objective.Residue; + objective.EvaluateAt(Pnew); + var RSSnew = objective.Value; if (double.IsNaN(RSSnew)) { @@ -220,7 +229,6 @@ namespace MathNet.Numerics.Optimization RSS = RSSnew; // update gradient and Hessian - objective.EvaluateJacobian(P); Gradient = objective.Gradient; Hessian = objective.Hessian; diagonalOfHessian = Hessian.Diagonal(); @@ -258,9 +266,6 @@ namespace MathNet.Numerics.Optimization exitCondition = ExitCondition.ExceedIterations; } - // finalize - objective.EvaluateCovariance(P); - return new NonlinearMinimizationResult(objective, iterations, exitCondition); } } diff --git a/src/Numerics/Optimization/NonlinearMinimizationResult.cs b/src/Numerics/Optimization/NonlinearMinimizationResult.cs index a88569ba..12381bbf 100644 --- a/src/Numerics/Optimization/NonlinearMinimizationResult.cs +++ b/src/Numerics/Optimization/NonlinearMinimizationResult.cs @@ -13,31 +13,25 @@ namespace MathNet.Numerics.Optimization /// /// Returns the best fit parameters. /// - public Vector BestFitParameters { get { return ModelInfoAtMinimum.Parameters; } } + public Vector BestFitParameters { get { return ModelInfoAtMinimum.Point; } } /// /// Returns the standard errors of the corresponding parameters /// - public Vector StandardErrors - { - get - { - if (ModelInfoAtMinimum.Covariance == null) - return null; - return ModelInfoAtMinimum.Covariance.Diagonal().PointwiseSqrt(); - } - } + public Vector StandardErrors { get; private set; } /// /// Returns the y-values of the fitted model that correspond to the independent values. /// - public Vector BestFitValues { get { return ModelInfoAtMinimum.Values; } } + public Vector BestFitValues { get { return ModelInfoAtMinimum.ModelValues; } } /// /// Returns the residual sum of squares. /// - public double Residue { get { return ModelInfoAtMinimum.Residue; } } - public double DegreeOfFreedom { get { return ModelInfoAtMinimum.DegreeOfFreedom; } } + public double Residue { get { return ModelInfoAtMinimum.Value; } } + public double DegreeOfFreedom { get { return ModelInfoAtMinimum.DegreeOfFreedom; } } + public Matrix Covariance { get; private set; } + public Matrix Correlation { get; private set; } public int Iterations { get; private set; } public ExitCondition ReasonForExit { get; private set; } @@ -47,6 +41,40 @@ namespace MathNet.Numerics.Optimization ModelInfoAtMinimum = modelInfo; Iterations = iterations; ReasonForExit = reasonForExit; + + AnalyzeResult(modelInfo); + } + + private void AnalyzeResult(IObjectiveModel objective) + { + objective.IsFinished = true; + objective.EvaluateAt(objective.Point); + + var Hessian = objective.Hessian; + if (Hessian == null || DegreeOfFreedom < 1) + { + Covariance = null; + Correlation = null; + StandardErrors = null; + return; + } + + Covariance = Hessian.PseudoInverse() * objective.Value / DegreeOfFreedom; + + if (Covariance != null) + { + StandardErrors = Covariance.Diagonal().PointwiseSqrt(); + + var correlation = Covariance.Clone(); + var d = correlation.Diagonal().PointwiseSqrt(); + var dd = d.OuterProduct(d); + Correlation = correlation.PointwiseDivide(dd); + } + else + { + StandardErrors = null; + Correlation = null; + } } } } diff --git a/src/Numerics/Optimization/ObjectiveModel.cs b/src/Numerics/Optimization/ObjectiveModel.cs index 3d786c18..77cae337 100644 --- a/src/Numerics/Optimization/ObjectiveModel.cs +++ b/src/Numerics/Optimization/ObjectiveModel.cs @@ -1,9 +1,6 @@ using MathNet.Numerics.LinearAlgebra; using MathNet.Numerics.Optimization.ObjectiveModels; using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; namespace MathNet.Numerics.Optimization { @@ -13,27 +10,22 @@ namespace MathNet.Numerics.Optimization /// Fitting model with a user supplied jacobian for non-linear least squares regression. /// public static IObjectiveModel FittingModel(Func, double, double> function, Func, double, Vector> derivatives, - Vector observedX, Vector observedY, Vector weight = null, - Vector lowerBound = null, Vector upperBound = null, Vector scales = null, List isFixed = null) + Vector observedX, Vector observedY, Vector weight = null) { var objective = new FittingObjectiveModel(function, derivatives); objective.SetObserved(observedX, observedY, weight); - objective.SetParameters(lowerBound, upperBound, scales, isFixed); return objective; } /// /// Fitting model for non-linear least squares regression. - /// The numerical jacobian with accuracy order is used. /// public static IObjectiveModel FittingModel(Func, double, double> function, Vector observedX, Vector observedY, Vector weight = null, - Vector lowerBound = null, Vector upperBound = null, Vector scales = null, List isFixed = null, int accuracyOrder = 2) { - var objective = new FittingObjectiveModel(function, null, accuracyOrder: accuracyOrder); + var objective = new FittingObjectiveModel(function, accuracyOrder: accuracyOrder); objective.SetObserved(observedX, observedY, weight); - objective.SetParameters(lowerBound, upperBound, scales, isFixed); return objective; } diff --git a/src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs b/src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs index 22401f63..eb30f357 100644 --- a/src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs +++ b/src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs @@ -8,10 +8,30 @@ 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 + #region Public Variables - Observed Data /// /// Set or get the values of the independent variable. @@ -25,178 +45,183 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels /// /// Set or get the values of the weights for the observations. - /// inverse of the standard measurement errors - /// If null, unity weighting is used. /// public Matrix Weights { get; private set; } - // W = LL' - private Vector L; + private Vector L; // Weights = LL' /// - /// Set or get the values of the parameters. + /// Get the number of observations. /// - public Vector Parameters { get; private set; } + public int NumberOfObservations { get { return (ObservedY == null) ? 0 : ObservedY.Count; } } - /// - /// Set or get the values of the parameters. - /// - public List IsFixed { get; set; } + #endregion Public Variables - Observed Data - /// - /// Set or get the values of the parameters. - /// - public Vector LowerBound { get; set; } + #region Public Variables - Bounds of Parameter /// - /// Set or get the values of the parameters. + /// Get the values of the parameters. /// - public Vector UpperBound { get; set; } + public List IsFixed { get; private set; } /// - /// Set or get the scale factor of the parameters. + /// Get the values of the parameters. /// - public Vector Scales { get; set; } + public Vector LowerBound { get; private set; } /// - /// Set of get whether or not the parameters are bounded. + /// Get the values of the parameters. /// - public bool IsBounded { get; set; } + public Vector UpperBound { get; private set; } /// - /// Get the y-values of the fitted model that correspond to the independent values. + /// Get the scale factor of the parameters. /// - public Vector Values { get; private set; } + public Vector Scales { get; private set; } /// - /// Get the error values, R(x; p) = L * (y - f(x; p)) where L = sqrt(W) + /// Get the number of unknown parameters. /// - private Vector Residuals; + public int NumberOfParameters { get { return (Point == null) ? 0 : Point.Count; } } - /// - /// Get the residual sum of squares, R.DotProduct(R) - /// - public double Residue { get; private set; } + #endregion Public Variables - Bounds of Parameter + + #region Public Variables - Others /// - /// Get the Jacobian matrix of x and p, J(x; p). + /// Get the number of calls to function. /// - public Matrix Jacobian { get; private set; } - + public int FunctionEvaluations { get; set; } /// - /// Get the Gradient vector of x and p, J'WR + /// Get the number of calls to jacobian. /// - public Vector Gradient { get; private set; } + public int JacobianEvaluations { get; set; } - /// - /// Get the Hessian matrix of x and p, J'WJ - /// - public Matrix Hessian { get; private 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); + } /// - /// Get the number of observations. + /// Set or get the values of the parameters. /// - public int NumberOfObservations { get { return (ObservedY == null) ? 0 : ObservedY.Count; } } + public Vector Point { get { return coefficients; } } /// - /// Get the number of unknown parameters. + /// Get the y-values of the fitted model that correspond to the independent values. /// - public int NumberOfParameters { get { return (Parameters == null) ? 0 : Parameters.Count; } } + public Vector ModelValues { get; private set; } /// - /// Get the degree of freedom + /// Get the residual sum of squares. /// - public int DegreeOfFreedom + public double Value { get { - var dof = NumberOfObservations - NumberOfParameters; - if (IsFixed != null) + if (!hasFunctionValue) { - dof = dof + IsFixed.Count(p => p == true); + EvaluateFunction(); + hasFunctionValue = true; } - return dof; + return functionValue; } } /// - /// Get the covariance matrix. - /// - public Matrix Covariance { get; private set; } - - /// - /// Get the correlation matrix. - /// - public Matrix Correlation { get; private set; } - - /// - /// Get the number of calls to function. - /// - public int FunctionEvaluations { get; set; } - /// - /// Get the number of calls to jacobian. - /// - public int JacobianEvaluations { get; set; } - - /// - /// Set or get the desired accuracy order of the numerical jacobian. + /// Get the Gradient vector of x and p. /// - public int AccuracyOrder { get; set; } + public Vector Gradient + { + get + { + if (!hasJacobianValue) + { + EvaluateJacobian(); + hasJacobianValue = true; + } + return gradientValue; + } + } /// - /// Get whether or not the analytical jacobian is supported. + /// Get the Hessian matrix of x and p, J'WJ /// - public bool IsJacobianSupported { get { return userDerivatives != null; } } - - #endregion Public Variables - - public FittingObjectiveModel(Func, double, double>function, Func, double, Vector> derivatives, int accuracyOrder = 2) + public Matrix Hessian { - userFunction = function; - userDerivatives = derivatives; - AccuracyOrder = Math.Min(6, Math.Max(1, accuracyOrder)); + get + { + if (!hasJacobianValue) + { + EvaluateJacobian(); + hasJacobianValue = true; + } + return hessianValue; + } } - public IObjectiveModel Fork() + /// + /// Get the degree of freedom + /// + public int DegreeOfFreedom { - return new FittingObjectiveModel(userFunction, userDerivatives, AccuracyOrder) + get { - ObservedX = ObservedX, - ObservedY = ObservedY, - Weights = Weights, - - Parameters = Parameters, - LowerBound = LowerBound, - UpperBound = UpperBound, - IsFixed = IsFixed, - Scales = Scales, - IsBounded = IsBounded, - - Residue = Residue, - Jacobian = Jacobian - }; + var df = NumberOfObservations - NumberOfParameters; + if (IsFixed != null) + { + df = df + IsFixed.Count(p => p == true); + } + return df; + } } - public IObjectiveModel CreateNew() - { - return new FittingObjectiveModel(userFunction, userDerivatives, AccuracyOrder); - } + 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) { - EvaluateFunction(point); - EvaluateJacobian(point); + EvaluateAt(point); - return new Tuple, Matrix>(Residue, Gradient, Hessian); + return new Tuple, Matrix>(Value, Gradient, Hessian); } - LowerBound = null; - UpperBound = null; - Scales = null; - IsFixed = null; - IsBounded = false; - var objective = new GradientHessianObjectiveFunction(function); return objective; } @@ -244,50 +269,37 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels } /// - /// Set observed data to fit. - /// - public void SetObserved(double[] observedX, double[] observedY, double[] weights = null) - { - if (observedX == null || observedY == null) - { - throw new ArgumentNullException("The data set can't be null."); - } - if (observedX.Length != observedY.Length) - { - throw new ArgumentException("The observed x data can't have different from observed y data."); - } - - var wVector = (weights == null) - ? null - : Vector.Build.DenseOfArray(weights); - SetObserved(Vector.Build.DenseOfArray(observedX), Vector.Build.DenseOfArray(observedY), wVector); - } - - /// - /// Set parameters. - /// - /// If bounded, the paramneters will be projected to unconstrained range by the mapping rule from the MINPACK. - /// If the projection is not needed, set IsBounded = false befre calling the Minimization method. + /// Set parameters and bounds. /// /// The lower bounds of parameters. /// The upper bounds of parameters. - /// /// The scaling constants of parameters + /// The scaling constants of parameters /// The list to the parameters fix or free. - public void SetParameters(Vector lowerBound = null, Vector upperBound = null, Vector scales = null, List isFixed = null) + 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 && lowerBound != null && upperBound.Count != lowerBound.Count) + if (upperBound != null && upperBound.Count != initialGuess.Count) { - throw new ArgumentException("The upper bounds can't have different elements from the lower bounds."); + throw new ArgumentException("The upper bounds can't have different elements from the initial guess."); } UpperBound = upperBound; @@ -295,13 +307,9 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels { throw new ArgumentException("The scales must be finite."); } - if (scales != null && lowerBound != null && scales.Count != lowerBound.Count) - { - throw new ArgumentException("The upper bounds can't have different elements from the lower bounds."); - } - if (scales != null && upperBound != null && scales.Count != upperBound.Count) + if (scales != null && scales.Count != initialGuess.Count) { - throw new ArgumentException("The upper bounds can't have different elements from the upper bounds."); + throw new ArgumentException("The upper bounds can't have different elements from the initial guess."); } if (scales != null && scales.Count(x => x < 0) > 0) { @@ -309,48 +317,20 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels } Scales = scales; - IsBounded = (LowerBound != null || UpperBound != null || Scales != null); - - if (isFixed != null && lowerBound != null && isFixed.Count != lowerBound.Count) - { - throw new ArgumentException("The initial guess can't have different elements from the lower bounds."); - } - if (isFixed != null && upperBound != null && isFixed.Count != upperBound.Count) + if (isFixed != null && isFixed.Count != initialGuess.Count) { - throw new ArgumentException("The initial guess can't have different elements from the upper bounds."); - } - if (isFixed != null && scales != null && isFixed.Count != scales.Count) - { - throw new ArgumentException("The initial guess can't have different elements from the scales."); + 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; - } - - /// - /// Set parameters. - /// - /// If bounded, the paramneters will be projected to unconstrained range by the mapping rule. - /// If the projection is not needed, set IsBounded = false befre calling the Minimization method. - /// - /// The lower bounds of parameters. - /// The upper bounds of parameters. - /// The scaling constants of parameters - /// The list to the parameters fix or free. - public void SetParameters(double[] lowerBound = null, double[] upperBound = null, double[] scales = null, bool[] isFixed = null) - { - var lb = (lowerBound == null) ? null : Vector.Build.DenseOfArray(lowerBound); - var ub = (upperBound == null) ? null : Vector.Build.DenseOfArray(upperBound); - var sc = (scales == null) ? null : Vector.Build.DenseOfArray(scales); - var fp = (isFixed == null) ? null : isFixed.ToList(); - SetParameters(lb, ub, sc, fp); + isBounded = LowerBound != null || UpperBound != null || Scales != null; } - public void EvaluateFunction(Vector parameters) + public void EvaluateAt(Vector parameters) { ValidateParameters(parameters); @@ -386,69 +366,84 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels // // Except when it is initial guess, the parameters argument is always internal parameter. // So, first map the parameters argument to the external parameters in order to calculate function values. - var Pext = (FunctionEvaluations > 0 && this.IsBounded) - ? ProjectParametersToExternal(parameters) - : parameters.Clone(); - // Project parameters, now this.Parameters are the internal parameters. - Parameters = (this.IsBounded) + 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 (Values == null) + if (ModelValues == null) { - Values = Vector.Build.Dense(NumberOfObservations); + ModelValues = Vector.Build.Dense(NumberOfObservations); } for (int i = 0; i < NumberOfObservations; i++) { - Values[i] = userFunction(Pext, ObservedX[i]); + ModelValues[i] = userFunction(Pext, ObservedX[i]); } FunctionEvaluations++; // calculate the weighted residuals - Residuals = (Weights == null) - ? ObservedY - Values - : (ObservedY - Values).PointwiseMultiply(L); + residuals = (Weights == null) + ? ObservedY - ModelValues + : (ObservedY - ModelValues).PointwiseMultiply(L); // Calculate the residual sum of squares - Residue = Residuals.DotProduct(Residuals); + functionValue = residuals.DotProduct(residuals); return; } - public void EvaluateJacobian(Vector parameters) + private void EvaluateJacobian() { - var Pext = (IsBounded) - ? ProjectParametersToExternal(parameters) - : parameters.Clone(); - // Calculates the jacobian of x and p. if (userDerivatives != null) { // analytical jacobian - if (Jacobian == null) + if (jacobianValue == null) { - Jacobian = Matrix.Build.Dense(NumberOfObservations, NumberOfParameters); + jacobianValue = Matrix.Build.Dense(NumberOfObservations, NumberOfParameters); } for (int i = 0; i < NumberOfObservations; i++) { - Jacobian.SetRow(i, userDerivatives(Pext, ObservedX[i])); + jacobianValue.SetRow(i, userDerivatives(Pext, ObservedX[i])); } JacobianEvaluations++; } else { // numerical jacobian - Jacobian = NumericalJacobian(Pext, Values, AccuracyOrder); - FunctionEvaluations += AccuracyOrder; + jacobianValue = NumericalJacobian(Pext, ModelValues, accuracyOrder); + FunctionEvaluations += accuracyOrder; } - var scaleFactors = (this.IsBounded) - ? ScaleFactorsOfJacobian(Parameters) - : Vector.Build.Dense(Parameters.Count, 1.0); + var scaleFactors = (isBounded && !IsFinished) + ? ScaleFactorsOfJacobian(Pint) + : Vector.Build.Dense(Pint.Count, 1.0); - // Jint(x; Pint) = Jext(x; Pext) * scale where scale = dPext/dPint + // 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++) @@ -456,58 +451,24 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels if (IsFixed != null && IsFixed[j]) { // if j-th parameter is fixed, set J[i, j] = 0 - Jacobian[i, j] = 0.0; + jacobianValue[i, j] = 0.0; } else { - Jacobian[i, j] = Jacobian[i, j] * scaleFactors[j]; + jacobianValue[i, j] = jacobianValue[i, j] * scaleFactors[j]; } } } // Gradient, g = -J'W(y − f(x; p)) = -J'L(L'E) = -J'LR - Gradient = (Weights == null) - ? -Jacobian.Transpose() * (ObservedY - Values) - : -Jacobian.Transpose() * Weights * (ObservedY - Values); + gradientValue = (Weights == null) + ? -jacobianValue.Transpose() * (ObservedY - ModelValues) + : -jacobianValue.Transpose() * Weights * (ObservedY - ModelValues); // approximated Hessian, H = J'WJ + ∑LRiHi ~ J'WJ near the minimum - Hessian = (Weights == null) - ? Jacobian.Transpose() * Jacobian - : Jacobian.Transpose() * Weights * Jacobian; - } - - public void EvaluateCovariance(Vector parameters) - { - // convert to bounded(external) parameters - var Pext = (IsBounded) - ? ProjectParametersToExternal(parameters) - : parameters.Clone(); - - // set IsBounded = false to get external Parameters and covariance matrix - this.IsBounded = false; - - EvaluateFunction(Pext); - EvaluateJacobian(Pext); - - // restore isBounded - this.IsBounded = (LowerBound != null || UpperBound != null); - - if (Hessian == null || Residuals == null || DegreeOfFreedom < 1) - { - Covariance = null; - Correlation = null; - return; - } - - var covariance = Hessian.PseudoInverse() * Residuals.DotProduct(Residuals) / DegreeOfFreedom; - Covariance = covariance; - - var correlation = covariance.Clone(); - var d = correlation.Diagonal().PointwiseSqrt(); - var dd = d.OuterProduct(d); - Correlation = correlation.PointwiseDivide(dd); - - return; + hessianValue = (Weights == null) + ? jacobianValue.Transpose() * jacobianValue + : jacobianValue.Transpose() * Weights * jacobianValue; } private void ValidateParameters(Vector parameters) @@ -538,16 +499,13 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels } } - #region Numerical Derivatives - - // Numerical derivatives by using the central or forward finite difference - private Matrix NumericalJacobian(Vector parameters, Vector currentValues, int accuracyOrder = 2) + 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 * parameters.PointwiseAbs().PointwiseMaximum(sqrtEpsilon); + var d = 0.000003 * Pext.PointwiseAbs().PointwiseMaximum(sqrtEpsilon); var h = Vector.Build.Dense(NumberOfParameters); for (int i = 0; i < NumberOfObservations; i++) @@ -560,12 +518,12 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels if (accuracyOrder >= 6) { // f'(x) = {- f(x - 3h) + 9f(x - 2h) - 45f(x - h) + 45f(x + h) - 9f(x + 2h) + f(x + 3h)} / 60h + O(h^6) - var f1 = userFunction(parameters - 3 * h, x); - var f2 = userFunction(parameters - 2 * h, x); - var f3 = userFunction(parameters - h, x); - var f4 = userFunction(parameters + h, x); - var f5 = userFunction(parameters + 2 * h, x); - var f6 = userFunction(parameters + 3 * h, x); + var 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; @@ -574,11 +532,11 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels { // f'(x) = {-137f(x) + 300f(x + h) - 300f(x + 2h) + 200f(x + 3h) - 75f(x + 4h) + 12f(x + 5h)} / 60h + O(h^5) var f1 = currentValues[i]; - var f2 = userFunction(parameters + h, x); - var f3 = userFunction(parameters + 2 * h, x); - var f4 = userFunction(parameters + 3 * h, x); - var f5 = userFunction(parameters + 4 * h, x); - var f6 = userFunction(parameters + 5 * h, x); + var 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; @@ -586,10 +544,10 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels else if (accuracyOrder == 4) { // f'(x) = {f(x - 2h) - 8f(x - h) + 8f(x + h) - f(x + 2h)} / 12h + O(h^4) - var f1 = userFunction(parameters - 2 * h, x); - var f2 = userFunction(parameters - h, x); - var f3 = userFunction(parameters + h, x); - var f4 = userFunction(parameters + 2 * h, x); + var 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; @@ -598,9 +556,9 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels { // f'(x) = {-11f(x) + 18f(x + h) - 9f(x + 2h) + 2f(x + 3h)} / 6h + O(h^3) var f1 = currentValues[i]; - var f2 = userFunction(parameters + h, x); - var f3 = userFunction(parameters + 2 * h, x); - var f4 = userFunction(parameters + 3 * h, x); + var 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; @@ -608,8 +566,8 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels else if (accuracyOrder == 2) { // f'(x) = {f(x + h) - f(x - h)} / 2h + O(h^2) - var f1 = userFunction(parameters + h, x); - var f2 = userFunction(parameters - h, x); + var f1 = userFunction(Pext + h, x); + var f2 = userFunction(Pext - h, x); var prime = (f1 - f2) / (2 * h[j]); derivertives[i, j] = prime; @@ -618,7 +576,7 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels { // f'(x) = {- f(x) + f(x + h)} / h + O(h) var f1 = currentValues[i]; - var f2 = userFunction(parameters + h, x); + var f2 = userFunction(Pext + h, x); var prime = (-f1 + f2) / h[j]; derivertives[i, j] = prime; @@ -631,10 +589,6 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels return derivertives; } - #endregion Numerical Derivatives - - #region Projection - private Vector ProjectParametersToInternal(Vector Pext) { var Pint = Pext.Clone(); @@ -771,6 +725,6 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels return scale; } - #endregion Projection + #endregion Private Methods } } diff --git a/src/Numerics/Optimization/Subproblems/DogLegSubproblem.cs b/src/Numerics/Optimization/Subproblems/DogLegSubproblem.cs index eb62acef..ce1ac838 100644 --- a/src/Numerics/Optimization/Subproblems/DogLegSubproblem.cs +++ b/src/Numerics/Optimization/Subproblems/DogLegSubproblem.cs @@ -13,12 +13,10 @@ namespace MathNet.Numerics.Optimization.Subproblems var Gradient = objective.Gradient; var Hessian = objective.Hessian; - // newton point - // the Gauss–Newton step by solving the normal equations + // newton point, the Gauss–Newton step by solving the normal equations var Pgn = -Hessian.PseudoInverse() * Gradient; // Hessian.Solve(Gradient) fails so many times... - // cauchy point - // steepest descent direction is given by + // cauchy point, steepest descent direction is given by var alpha = Gradient.DotProduct(Gradient) / (Hessian * Gradient).DotProduct(Gradient); var Psd = -alpha * Gradient; diff --git a/src/Numerics/Optimization/Subproblems/NewtonCGSubproblem.cs b/src/Numerics/Optimization/Subproblems/NewtonCGSubproblem.cs index ec26b34a..73149818 100644 --- a/src/Numerics/Optimization/Subproblems/NewtonCGSubproblem.cs +++ b/src/Numerics/Optimization/Subproblems/NewtonCGSubproblem.cs @@ -58,7 +58,7 @@ namespace MathNet.Numerics.Optimization.Subproblems z = znext; r = rnext; - d = -rnext + rnext_sq / r_sq * d; ; + d = -rnext + rnext_sq / r_sq * d; } } } diff --git a/src/Numerics/Optimization/TrustRegionMinimizerBase.cs b/src/Numerics/Optimization/TrustRegionMinimizerBase.cs index 22de97d8..d32c68f8 100644 --- a/src/Numerics/Optimization/TrustRegionMinimizerBase.cs +++ b/src/Numerics/Optimization/TrustRegionMinimizerBase.cs @@ -1,5 +1,7 @@ using MathNet.Numerics.LinearAlgebra; using System; +using System.Collections.Generic; +using System.Linq; namespace MathNet.Numerics.Optimization { @@ -46,39 +48,50 @@ namespace MathNet.Numerics.Optimization MaximumIterations = maximumIterations; } - public NonlinearMinimizationResult FindMinimum(IObjectiveModel objective, Vector initialGuess) + 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, Subproblem, GradientTolerance, StepTolerance, FunctionTolerance, RadiusTolerance, MaximumIterations); + return Minimum(Subproblem, objective, initialGuess, lowerBound, upperBound, scales, isFixed, + GradientTolerance, StepTolerance, FunctionTolerance, RadiusTolerance, MaximumIterations); } - public NonlinearMinimizationResult FindMinimum(IObjectiveModel objective, double[] initialGuess) + 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"); - return Minimum(objective, CreateVector.DenseOfArray(initialGuess), Subproblem, GradientTolerance, StepTolerance, FunctionTolerance, RadiusTolerance, MaximumIterations); + var lb = (lowerBound == null) ? null : CreateVector.Dense(lowerBound); + var ub = (upperBound == null) ? null : CreateVector.Dense(upperBound); + var sc = (scales == null) ? null : CreateVector.Dense(scales); + var fx = (isFixed == null) ? null : isFixed.ToList(); + + return Minimum(Subproblem, objective, CreateVector.DenseOfArray(initialGuess), lb, ub, sc, fx, + GradientTolerance, StepTolerance, FunctionTolerance, RadiusTolerance, MaximumIterations); } /// /// Non-linear least square fitting by the trust-region algorithm. /// /// The objective model, including function, jacobian, observations, and parameter bounds. + /// The subproblem /// The initial guess values. - /// The subproblem /// The stopping threshold for L2 norm of the residuals. /// The stopping threshold for infinity norm of the gradient vector. /// The stopping threshold for L2 norm of the change of parameters. /// The stopping threshold for trust region radius /// The max iterations. /// - public static NonlinearMinimizationResult Minimum(IObjectiveModel objective, Vector initialGuess, ITrustRegionSubproblem subproblem, double gradientTolerance = 1E-8, double stepTolerance = 1E-8, double functionTolerance = 1E-8, double radiusTolerance = 1E-18, int maximumIterations = -1) + public static NonlinearMinimizationResult Minimum(ITrustRegionSubproblem subproblem, IObjectiveModel objective, Vector initialGuess, + Vector lowerBound = null, Vector upperBound = null, Vector scales = null, List isFixed = null, + double gradientTolerance = 1E-8, double stepTolerance = 1E-8, double functionTolerance = 1E-8, double radiusTolerance = 1E-18, int maximumIterations = -1) { // Non-linear least square fitting by the trust-region algorithm. // @@ -123,16 +136,19 @@ namespace MathNet.Numerics.Optimization if (initialGuess == null) throw new ArgumentNullException("initialGuess"); + objective.SetParameters(initialGuess, lowerBound, upperBound, scales, isFixed); + ExitCondition exitCondition = ExitCondition.None; // Initialize objective objective.FunctionEvaluations = 0; objective.JacobianEvaluations = 0; + objective.IsFinished = false; // First, calculate function values and setup variables - objective.EvaluateFunction(initialGuess); - var P = objective.Parameters; // current parameters - var RSS = objective.Residue; // Residual Sum of Squares = R'R + 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 if (maximumIterations < 0) @@ -140,7 +156,7 @@ namespace MathNet.Numerics.Optimization maximumIterations = 200 * (initialGuess.Count + 1); } - // if R == NaN, stop + // if RSS == NaN, stop if (double.IsNaN(RSS)) { exitCondition = ExitCondition.InvalidValues; @@ -159,10 +175,9 @@ namespace MathNet.Numerics.Optimization exitCondition = ExitCondition.Converged; // SmallRSS } - // Evaluate projected Hessian, and gradient - objective.EvaluateJacobian(P); - var Hessian = objective.Hessian; + // evaluate projected gradient and Hessian var Gradient = objective.Gradient; + var Hessian = objective.Hessian; // if ||g||_oo <= gtol, found and stop if (Gradient.InfinityNorm() <= gradientTolerance) @@ -172,8 +187,6 @@ namespace MathNet.Numerics.Optimization if (exitCondition != ExitCondition.None) { - // finalize - objective.EvaluateCovariance(P); return new NonlinearMinimizationResult(objective, -1, exitCondition); } @@ -191,7 +204,7 @@ namespace MathNet.Numerics.Optimization var Pstep = subproblem.Pstep; var hitBoundary = subproblem.HitBoundary; // predicted reduction = L(0) - L(Δp) = -Δp'g - 1/2 * Δp'HΔp - var predictedReduction = -objective.Gradient.DotProduct(Pstep) - 0.5 * Pstep.DotProduct(objective.Hessian * Pstep); + var predictedReduction = -Gradient.DotProduct(Pstep) - 0.5 * Pstep.DotProduct(Hessian * Pstep); if (Pstep.L2Norm() <= stepTolerance * (stepTolerance + P.L2Norm())) { @@ -201,13 +214,20 @@ namespace MathNet.Numerics.Optimization var Pnew = P + Pstep; // parameters to test - objective.EvaluateFunction(Pnew); - var RSSnew = objective.Residue; + objective.EvaluateAt(Pnew); + var RSSnew = objective.Value; + + // if RSS == NaN, stop + if (double.IsNaN(RSSnew)) + { + exitCondition = ExitCondition.InvalidValues; + break; + } // calculate the ratio of the actual to the predicted reduction. double rho = (predictedReduction != 0) ? (RSS - RSSnew) / predictedReduction - : 0; + : 0.0; if (rho > 0.75 && hitBoundary) { @@ -229,8 +249,7 @@ namespace MathNet.Numerics.Optimization Pnew.CopyTo(P); RSS = RSSnew; - // update Jacobian, Hessian, and gradient - objective.EvaluateJacobian(P); + // evaluate projected gradient and Hessian Gradient = objective.Gradient; Hessian = objective.Hessian; @@ -253,9 +272,6 @@ namespace MathNet.Numerics.Optimization exitCondition = ExitCondition.ExceedIterations; } - // finalize - objective.EvaluateCovariance(P); - return new NonlinearMinimizationResult(objective, iterations, exitCondition); } }