Retain Newtown update only if objective function improved #455

Author: seananderson

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: sdmTMB
 Title: Spatial and Spatiotemporal SPDE-Based GLMMs with 'TMB'
-Version: 0.7.2.9003
+Version: 0.7.2.9004
 Authors@R: c(
     person(c("Sean", "C."), "Anderson", , "sean@seananderson.ca",
            role = c("aut", "cre"),

NEWS.md

@@ -2,15 +2,15 @@
 
 ## Minor improvements and fixes
 
+* Only retain Newton update parameters if they improve the objective function.
+  #455
+  
 * Suppress `nlminb()` warnings by default, which can usually be ignored by the
   user and may be confusing. This can be controlled via
   `sdmTMB(..., control = sdmTMBcontrol(suppress_nlminb_warnings = FALSE))`.
+  This option now mirrors tinyVAST.
 
-* Only run Newton optimization steps if maximum absolute gradient is >= -1e6
-  to save time.
-
-* Fix reported gradient if Newton loops used in optimization. Previously,
-  the gradients would have looked overly pessimistic for some models. #455
+* Only run Newton updates if maximum absolute gradient is >= 1e-9 to save time.
 
 * Round time-varying AR(1) rho to 2 decimals in model printing/summary.
 

R/extra-optimization.R

@@ -50,26 +50,37 @@ run_extra_optimization <- function(object,
     }
   }
 
-  for (i in seq_len(newton_loops)) {
-    g <- as.numeric(new_obj$tmb_obj$gr(tmb_opt$par))
-    h <- stats::optimHess(
-      tmb_opt$par,
-      fn = new_obj$tmb_obj$fn,
-      gr = new_obj$tmb_obj$gr,
-      lower = object$lower,
-      upper = object$upper
-    )
-    tmb_opt$par <- tmb_opt$par - solve(h, g)
-    tmb_opt$objective <- new_obj$tmb_obj$fn(tmb_opt$par)
-    if (i == newton_loops) {
-      new_obj$tmb_obj$fn(tmb_opt$par) # call once to update environment
-    }
-  }
+  tmb_opt <- run_newton_loops(newton_loops = newton_loops, tmb_opt, new_obj, silent = FALSE)
   new_obj$model <- tmb_opt
-  new_obj$sd_report <- TMB::sdreport(new_obj$tmb_obj, par.fixed = tmb_opt$par,
+  new_obj$sd_report <- TMB::sdreport(new_obj$tmb_obj,
     getJointPrecision = "jointPrecision" %in% names(object$sd_report))
   conv <- get_convergence_diagnostics(new_obj$sd_report)
   new_obj$gradients <- conv$final_grads
   new_obj$bad_eig <- conv$bad_eig
   new_obj
 }
+
+run_newton_loops <- function(newton_loops, opt, obj, silent = TRUE) {
+  if (newton_loops > 0) {
+    if (!silent) cli_inform("attempting to improve convergence with a Newton update\n")
+    for (i in seq_len(newton_loops)) {
+      g <- as.numeric(obj$gr(opt$par))
+      if (max(abs(g)) < 1e-9) {
+        if (!silent) cli_inform(c("maximum absolute gradient is already < 1e-9;",
+          "skipping any remaining Newton updates for speed\n"))
+        break
+      }
+      h <- stats::optimHess(opt$par, fn = obj$fn, gr = obj$gr)
+      new_par <- opt$par - solve(h, g)
+      new_objective <- obj$fn(new_par) # also updates obj$env$last.par and obj$env$last.par.best!
+      if (new_objective < opt$objective) {
+        if (!silent) cli_inform("accepting parameters from Newton update\n")
+        opt$par <- new_par
+        opt$objective <- new_objective
+      } else {
+        if (!silent) cli_inform("retaining parameters from before Newton update\n")
+      }
+    }
+  }
+  opt
+}

R/fit.R

@@ -1609,21 +1609,12 @@ sdmTMB <- function(
       cli_inform("Upper or lower limits were set. `stats::optimHess()` will ignore these limits. Set `control = sdmTMBcontrol(newton_loops = 0)` to avoid the `stats::optimHess()` optimization if desired.")
     }
   }
-  if (newton_loops > 0) {
-    if (!silent) cli_inform("attempting to improve convergence with optimHess\n")
-    for (i in seq_len(newton_loops)) {
-      g <- as.numeric(tmb_obj$gr(tmb_opt$par))
-      if (max(abs(g)) < 1e-6) break
-      h <- stats::optimHess(tmb_opt$par, fn = tmb_obj$fn, gr = tmb_obj$gr)
-      tmb_opt$par <- tmb_opt$par - solve(h, g)
-      tmb_opt$objective <- tmb_obj$fn(tmb_opt$par)
-    }
-  }
+
+  tmb_opt <- run_newton_loops(newton_loops = newton_loops, tmb_opt, tmb_obj, silent)
   check_bounds(tmb_opt$par, lim$lower, lim$upper)
 
   if (!silent) cli_inform("running TMB sdreport\n")
-  sd_report <- TMB::sdreport(tmb_obj, par.fixed = tmb_opt$par,
-    getJointPrecision = get_joint_precision)
+  sd_report <- TMB::sdreport(tmb_obj, getJointPrecision = get_joint_precision)
   conv <- get_convergence_diagnostics(sd_report)
 
   ## save params that families need to grab from environments: