From 67e1b2a277b38f5991e8d12c23b669b3e63e7339 Mon Sep 17 00:00:00 2001
From: Jesper Andersson <jesper.andersson@ndcn.ox.ac.uk>
Date: Mon, 24 Feb 2020 12:00:31 +0000
Subject: [PATCH] Made changes to Levenberg-Marquardt (levmar) so that there is
now a pure Gauss-Newton option
---
nonlin.cpp | 10 +++++++---
nonlin.h | 6 ++++--
2 files changed, 11 insertions(+), 5 deletions(-)
diff --git a/nonlin.cpp b/nonlin.cpp
index 16e5641..6f7c3c5 100644
--- a/nonlin.cpp
+++ b/nonlin.cpp
@@ -359,7 +359,6 @@ NonlinOut levmar(const NonlinParam& p, const NonlinCF& cfo)
H = cfo.hess(p.Par(),H); // Hessian evaluated at current parameters
prof.EndEntry(hess_key);
}
- double nudge_key = prof.StartEntry("Nudge hessian");
for (int i=1; i<=p.NPar(); i++) { // Nudge it
if (p.GaussNewtonType() == LM_LM) { // If Levenberg-Marquardt
// H->AddTo(i,i,(p.Lambda()-olambda)*H->Peek(i,i));
@@ -368,8 +367,8 @@ NonlinOut levmar(const NonlinParam& p, const NonlinCF& cfo)
else if (p.GaussNewtonType() == LM_L) { // If Levenberg
H->AddTo(i,i,p.Lambda()-olambda);
}
+ else if (p.GaussNewtonType() == LM_GN) {} // If it is pure Gauss-Newton we don't nudge at all
}
- prof.EndEntry(nudge_key);
ColumnVector step;
double ncf = 0.0;
bool inv_fail = false; // Signals failure of equation solving
@@ -384,7 +383,7 @@ NonlinOut levmar(const NonlinParam& p, const NonlinCF& cfo)
prof.EndEntry(fail_key);
}
prof.EndEntry(solve_key);
- if (!inv_fail && (success = (ncf < p.CF()))) { // If last step successful
+ if (!inv_fail && (success = (ncf < p.CF()))) { // If last step successful
olambda = 0.0; // Pristine Hessian, so no need to undo old lambda
p.SetPar(p.Par()+step); // Set attempt as new parameters
p.SetLambda(p.Lambda()/10.0); // Decrease nudge factor
@@ -395,6 +394,11 @@ NonlinOut levmar(const NonlinParam& p, const NonlinCF& cfo)
p.SetCF(ncf); // Store value of cost-function
}
else { // If last step was unsuccesful
+ if (p.GaussNewtonType() == LM_GN) { // For pure Gauss-Newton we stop if unsuccessful in reducing cf
+ if (inv_fail) p.SetStatus(NL_INV_FAIL); // Inversion of Hessian failed
+ else if (!(ncf < p.CF())) p.SetStatus(NL_CF_INC_CONV);// Cost-function failed to decrease
+ return(p.Status());
+ }
success = false;
olambda = p.Lambda(); // Returning to same H, so must undo old lambda
p.SetLambda(10.0*p.Lambda()); // Increase nudge factor
diff --git a/nonlin.h b/nonlin.h
index 183e63c..63c7f02 100644
--- a/nonlin.h
+++ b/nonlin.h
@@ -34,7 +34,7 @@ enum NLMethod {NL_VM, // Variable-Metric (see NRin
NL_GD, // Gradient Descent
NL_NM}; // Nelder-Mead simplex method (see NRinC)
-enum LMType {LM_L, LM_LM}; // Levenberg or Levenberg-Marquardt
+enum LMType {LM_GN, LM_L, LM_LM}; // Gauss-Newton, Levenberg or Levenberg-Marquardt
enum VMUpdateType {VM_DFP, VM_BFGS}; // See NRinC chapter 10.
@@ -50,10 +50,12 @@ enum LinOut {LM_MAXITER, // Too many iterations in li
enum NonlinOut {NL_UNDEFINED, // Initial value before minimisation
NL_MAXITER, // Too many iterations
- NL_LM_MAXITER, // To many iterations during a line-minimisation
+ NL_LM_MAXITER, // Too many iterations during a line-minimisation
+ NL_INV_FAIL, // Failure to invert Hessian
NL_PARCONV, // Convergence. Step in parameter space small
NL_GRADCONV, // Convergence. Gradient small
NL_CFCONV, // Convergence. Change in cost-function small
+ NL_CF_INC_CONV, // Sort of convergence. Gauss-Newton step caused increase in cost
NL_LCONV}; // Convergence, lambda very large
const double EPS = 2.0e-16; // Losely based on NRinC 20.1
--
GitLab