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);
}
}