From 975cd212520d1540607275d2d9adb2366ae10876 Mon Sep 17 00:00:00 2001 From: diluculo Date: Thu, 3 Jan 2019 10:13:21 +0900 Subject: [PATCH] Changed the sign of the gradient of the FittingObjectiveModel to match the scheme of the existing Minimizers. --- .../Optimization/LevenbergMarquardtMinimizer.cs | 14 +++++++------- .../ObjectiveModels/FittingObjectiveModel.cs | 8 ++++---- .../Optimization/Subproblems/DogLegSubproblem.cs | 4 ++-- .../Optimization/Subproblems/NewtonCGSubproblem.cs | 5 +++-- .../Optimization/TrustRegionMinimizerBase.cs | 6 +++--- 5 files changed, 19 insertions(+), 18 deletions(-) diff --git a/src/Numerics/Optimization/LevenbergMarquardtMinimizer.cs b/src/Numerics/Optimization/LevenbergMarquardtMinimizer.cs index ddae4d72..ab6bb421 100644 --- a/src/Numerics/Optimization/LevenbergMarquardtMinimizer.cs +++ b/src/Numerics/Optimization/LevenbergMarquardtMinimizer.cs @@ -92,14 +92,14 @@ namespace MathNet.Numerics.Optimization // Residuals, R = L(y - f(x; p)) // Residual sum of squares, RSS = ||R||^2 = R.DotProduct(R) // 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 // // The Levenberg-Marquardt algorithm is summarized as follows: - // initially let μ = τ * max(diag(J'WJ)). + // initially let μ = τ * max(diag(H)). // repeat - // solve linear equations: (J'WJ + μI)ΔP = J'R - // let ρ = (||R||^2 - ||Rnew||^2) / (Δp'(μΔp + J'R)). + // solve linear equations: (H + μI)ΔP = -g + // let ρ = (||R||^2 - ||Rnew||^2) / (Δp'(μΔp - g)). // if ρ > ε, P = P + ΔP; μ = μ * max(1/3, 1 - (2ρ - 1)^3); ν = 2; // otherwise μ = μ*ν; ν = 2*ν; // @@ -186,7 +186,7 @@ namespace MathNet.Numerics.Optimization Hessian.SetDiagonal(Hessian.Diagonal() + mu); // hessian[i, i] = hessian[i, i] + mu; // solve normal equations - Pstep = Hessian.Solve(Gradient); + Pstep = Hessian.Solve(-Gradient); // if ||ΔP|| <= xTol * (||P|| + xTol), found and stop 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. - // ρ = (RSS - RSSnew) / (Δp'(μΔp + g)) - var predictedReduction = Pstep.DotProduct(mu * Pstep + Gradient); + // ρ = (RSS - RSSnew) / (Δp'(μΔp - g)) + var predictedReduction = Pstep.DotProduct(mu * Pstep - Gradient); var rho = (predictedReduction != 0) ? (RSS - RSSnew) / predictedReduction : 0; diff --git a/src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs b/src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs index f1204189..22401f63 100644 --- a/src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs +++ b/src/Numerics/Optimization/ObjectiveModels/FittingObjectiveModel.cs @@ -188,7 +188,7 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels EvaluateFunction(point); EvaluateJacobian(point); - return new Tuple, Matrix>(Residue, -Gradient, Hessian); + return new Tuple, Matrix>(Residue, Gradient, Hessian); } 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) - ? Jacobian.Transpose() * (ObservedY - Values) - : Jacobian.Transpose() * Weights * (ObservedY - Values); + ? -Jacobian.Transpose() * (ObservedY - Values) + : -Jacobian.Transpose() * Weights * (ObservedY - Values); // approximated Hessian, H = J'WJ + ∑LRiHi ~ J'WJ near the minimum Hessian = (Weights == null) diff --git a/src/Numerics/Optimization/Subproblems/DogLegSubproblem.cs b/src/Numerics/Optimization/Subproblems/DogLegSubproblem.cs index 83671ce7..eb62acef 100644 --- a/src/Numerics/Optimization/Subproblems/DogLegSubproblem.cs +++ b/src/Numerics/Optimization/Subproblems/DogLegSubproblem.cs @@ -15,12 +15,12 @@ namespace MathNet.Numerics.Optimization.Subproblems // newton point // 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 // steepest descent direction is given by var alpha = Gradient.DotProduct(Gradient) / (Hessian * Gradient).DotProduct(Gradient); - var Psd = alpha * Gradient; + var Psd = -alpha * Gradient; // update step and prectted reduction if (Pgn.L2Norm() <= delta) diff --git a/src/Numerics/Optimization/Subproblems/NewtonCGSubproblem.cs b/src/Numerics/Optimization/Subproblems/NewtonCGSubproblem.cs index 04c98214..ec26b34a 100644 --- a/src/Numerics/Optimization/Subproblems/NewtonCGSubproblem.cs +++ b/src/Numerics/Optimization/Subproblems/NewtonCGSubproblem.cs @@ -15,11 +15,12 @@ namespace MathNet.Numerics.Optimization.Subproblems var Hessian = objective.Hessian; // 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 var z = Vector.Build.Dense(Hessian.RowCount); - var r = -Gradient; + var r = Gradient; var d = -r; while (true) diff --git a/src/Numerics/Optimization/TrustRegionMinimizerBase.cs b/src/Numerics/Optimization/TrustRegionMinimizerBase.cs index 50711a69..22de97d8 100644 --- a/src/Numerics/Optimization/TrustRegionMinimizerBase.cs +++ b/src/Numerics/Optimization/TrustRegionMinimizerBase.cs @@ -93,7 +93,7 @@ namespace MathNet.Numerics.Optimization // Residuals, R = L(y - f(x; p)) // Residual sum of squares, RSS = ||R||^2 = R.DotProduct(R) // 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 // // The trust region algorithm is summarized as follows: @@ -190,8 +190,8 @@ namespace MathNet.Numerics.Optimization subproblem.Solve(objective, delta); 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); + // 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); if (Pstep.L2Norm() <= stepTolerance * (stepTolerance + P.L2Norm())) {