Browse Source

Changed the sign of the gradient of the FittingObjectiveModel to match the scheme of the existing Minimizers.

arrays
diluculo 8 years ago
parent
commit
975cd21252
  1. 14
      src/Numerics/Optimization/LevenbergMarquardtMinimizer.cs
  2. 8
      src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs
  3. 4
      src/Numerics/Optimization/Subproblems/DogLegSubproblem.cs
  4. 5
      src/Numerics/Optimization/Subproblems/NewtonCGSubproblem.cs
  5. 6
      src/Numerics/Optimization/TrustRegionMinimizerBase.cs

14
src/Numerics/Optimization/LevenbergMarquardtMinimizer.cs

@ -92,14 +92,14 @@ namespace MathNet.Numerics.Optimization
// Residuals, R = L(y - f(x; p)) // Residuals, R = L(y - f(x; p))
// Residual sum of squares, RSS = ||R||^2 = R.DotProduct(R) // Residual sum of squares, RSS = ||R||^2 = R.DotProduct(R)
// Jacobian J = df(x; p)/dp // Jacobian J = df(x; p)/dp
// Gradient g = J'W(y − f(x; p)) = J'LR // Gradient g = -J'W(y − f(x; p)) = -J'LR
// Approximated Hessian H = J'WJ // Approximated Hessian H = J'WJ
// //
// The Levenberg-Marquardt algorithm is summarized as follows: // The Levenberg-Marquardt algorithm is summarized as follows:
// initially let μ = τ * max(diag(J'WJ)). // initially let μ = τ * max(diag(H)).
// repeat // repeat
// solve linear equations: (J'WJ + μI)ΔP = J'R // solve linear equations: (H + μI)ΔP = -g
// let ρ = (||R||^2 - ||Rnew||^2) / (Δp'(μΔp + J'R)). // let ρ = (||R||^2 - ||Rnew||^2) / (Δp'(μΔp - g)).
// if ρ > ε, P = P + ΔP; μ = μ * max(1/3, 1 - (2ρ - 1)^3); ν = 2; // if ρ > ε, P = P + ΔP; μ = μ * max(1/3, 1 - (2ρ - 1)^3); ν = 2;
// otherwise μ = μ*ν; ν = 2*ν; // otherwise μ = μ*ν; ν = 2*ν;
// //
@ -186,7 +186,7 @@ namespace MathNet.Numerics.Optimization
Hessian.SetDiagonal(Hessian.Diagonal() + mu); // hessian[i, i] = hessian[i, i] + mu; Hessian.SetDiagonal(Hessian.Diagonal() + mu); // hessian[i, i] = hessian[i, i] + mu;
// solve normal equations // solve normal equations
Pstep = Hessian.Solve(Gradient); Pstep = Hessian.Solve(-Gradient);
// if ||ΔP|| <= xTol * (||P|| + xTol), found and stop // if ||ΔP|| <= xTol * (||P|| + xTol), found and stop
if (Pstep.L2Norm() <= stepTolerance * (stepTolerance + P.DotProduct(P))) if (Pstep.L2Norm() <= stepTolerance * (stepTolerance + P.DotProduct(P)))
@ -207,8 +207,8 @@ namespace MathNet.Numerics.Optimization
} }
// calculate the ratio of the actual to the predicted reduction. // calculate the ratio of the actual to the predicted reduction.
// ρ = (RSS - RSSnew) / (Δp'(μΔp + g)) // ρ = (RSS - RSSnew) / (Δp'(μΔp - g))
var predictedReduction = Pstep.DotProduct(mu * Pstep + Gradient); var predictedReduction = Pstep.DotProduct(mu * Pstep - Gradient);
var rho = (predictedReduction != 0) var rho = (predictedReduction != 0)
? (RSS - RSSnew) / predictedReduction ? (RSS - RSSnew) / predictedReduction
: 0; : 0;

8
src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs

@ -188,7 +188,7 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels
EvaluateFunction(point); EvaluateFunction(point);
EvaluateJacobian(point); EvaluateJacobian(point);
return new Tuple<double, Vector<double>, Matrix<double>>(Residue, -Gradient, Hessian); return new Tuple<double, Vector<double>, Matrix<double>>(Residue, Gradient, Hessian);
} }
LowerBound = null; LowerBound = null;
@ -465,10 +465,10 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels
} }
} }
// Gradient, g = J'W(y − f(x; p)) = J'L(L'E) = J'LR // Gradient, g = -J'W(y − f(x; p)) = -J'L(L'E) = -J'LR
Gradient = (Weights == null) Gradient = (Weights == null)
? Jacobian.Transpose() * (ObservedY - Values) ? -Jacobian.Transpose() * (ObservedY - Values)
: Jacobian.Transpose() * Weights * (ObservedY - Values); : -Jacobian.Transpose() * Weights * (ObservedY - Values);
// approximated Hessian, H = J'WJ + ∑LRiHi ~ J'WJ near the minimum // approximated Hessian, H = J'WJ + ∑LRiHi ~ J'WJ near the minimum
Hessian = (Weights == null) Hessian = (Weights == null)

4
src/Numerics/Optimization/Subproblems/DogLegSubproblem.cs

@ -15,12 +15,12 @@ namespace MathNet.Numerics.Optimization.Subproblems
// newton point // newton point
// the Gauss–Newton step by solving the normal equations // the Gauss–Newton step by solving the normal equations
var Pgn = Hessian.PseudoInverse() * Gradient; // Hessian.Solve(Gradient) fails so many times... var Pgn = -Hessian.PseudoInverse() * Gradient; // Hessian.Solve(Gradient) fails so many times...
// cauchy point // cauchy point
// steepest descent direction is given by // steepest descent direction is given by
var alpha = Gradient.DotProduct(Gradient) / (Hessian * Gradient).DotProduct(Gradient); var alpha = Gradient.DotProduct(Gradient) / (Hessian * Gradient).DotProduct(Gradient);
var Psd = alpha * Gradient; var Psd = -alpha * Gradient;
// update step and prectted reduction // update step and prectted reduction
if (Pgn.L2Norm() <= delta) if (Pgn.L2Norm() <= delta)

5
src/Numerics/Optimization/Subproblems/NewtonCGSubproblem.cs

@ -15,11 +15,12 @@ namespace MathNet.Numerics.Optimization.Subproblems
var Hessian = objective.Hessian; var Hessian = objective.Hessian;
// define tolerance // define tolerance
var tolerance = Math.Min(0.5, Math.Sqrt(Gradient.L2Norm())) * Gradient.L2Norm(); var gnorm = Gradient.L2Norm();
var tolerance = Math.Min(0.5, Math.Sqrt(gnorm)) * gnorm;
// initialize internal variables // initialize internal variables
var z = Vector<double>.Build.Dense(Hessian.RowCount); var z = Vector<double>.Build.Dense(Hessian.RowCount);
var r = -Gradient; var r = Gradient;
var d = -r; var d = -r;
while (true) while (true)

6
src/Numerics/Optimization/TrustRegionMinimizerBase.cs

@ -93,7 +93,7 @@ namespace MathNet.Numerics.Optimization
// Residuals, R = L(y - f(x; p)) // Residuals, R = L(y - f(x; p))
// Residual sum of squares, RSS = ||R||^2 = R.DotProduct(R) // Residual sum of squares, RSS = ||R||^2 = R.DotProduct(R)
// Jacobian J = df(x; p)/dp // Jacobian J = df(x; p)/dp
// Gradient g = J'W(y − f(x; p)) = J'LR // Gradient g = -J'W(y − f(x; p)) = -J'LR
// Approximated Hessian H = J'WJ // Approximated Hessian H = J'WJ
// //
// The trust region algorithm is summarized as follows: // The trust region algorithm is summarized as follows:
@ -190,8 +190,8 @@ namespace MathNet.Numerics.Optimization
subproblem.Solve(objective, delta); subproblem.Solve(objective, delta);
var Pstep = subproblem.Pstep; var Pstep = subproblem.Pstep;
var hitBoundary = subproblem.HitBoundary; var hitBoundary = subproblem.HitBoundary;
// predicted reduction = L(0) - L(Δp) = Δp'g - 1/2 * Δp'HΔp // 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 = -objective.Gradient.DotProduct(Pstep) - 0.5 * Pstep.DotProduct(objective.Hessian * Pstep);
if (Pstep.L2Norm() <= stepTolerance * (stepTolerance + P.L2Norm())) if (Pstep.L2Norm() <= stepTolerance * (stepTolerance + P.L2Norm()))
{ {

Loading…
Cancel
Save