MakeTape$newton's (sometimes?) has NaN effect on obj$he()

[phillipbvetter]

Hi.

As always, thanks for a great package!

I believe I've found some inconsistent behaviour with obj$he() when the objective function includes calls to MakeTape$newton(). For certain parameters columns of the hessian are NaNs. The behaviour seems to be affected by previous calls to obj$he.

I provide a code example below taken from https://github.com/kaskr/RTMB/blob/master/tmb_examples/linreg.R, with the addition of each parameter wrapped in a MakeTape$newton() function - which is the identity function.

library(RTMB)
set.seed(14)

g <- RTMB::MakeTape(function(x) ((x[1])-x[2])^2, rep(1,2))
f_newton <- g$newton(2)

data <- list(Y = rnorm(10) + 1:10, x=1:10)
parameters <- list(a=0, b=0, logSigma=0)

f <- function(parms) {
    Y <- data$Y
    x <- data$x
    a <- f_newton(parms$a)
    b <- f_newton(parms$b)
    logSigma <- f_newton(parms$logSigma)
    nll = -sum(dnorm(Y, a+b*x, exp(logSigma), TRUE))
    nll
}
obj <- MakeADFun(f, parameters, silent=TRUE)
obj$he(obj$par) #works fine
obj$he(obj$par + 4.0) #works fine
obj$he(obj$par + 2.0) # gives NaN
obj$he(obj$par + 4.0) # now this gives NaN
obj$he(obj$par) #works fine

# Optimize
opt <- do.call("optim", obj)
# nlminb(obj$par,obj$fn,obj$gr,obj$he) #doesn't work because some hessian calls will give NaN eventually...

# Check hessian
obj$he(opt$par) #gives NaN's

rep <- sdreport(obj)
summary(rep) # but here we get std. errors!

# A small perturbation gives answers
obj$he(opt$par+1e-15) #does work
obj$he(opt$par-1e-15) # does work

# Check inversion of neg. log like hessian against sdreport result
std.error <- sqrt(diag(solve(obj$he(opt$par+1e-15))))
cbind(
  sdreport = sqrt(diag(rep$cov.fixed)),
  inverted_hessian = std.error,
  difference  = sqrt(diag(rep$cov.fixed)) - std.error
)

[kaskr]

For future reference, I'm adding here some info on how to debug this issue:
Once we have an unexpected NaN such as

obj$he(obj$par + 2.0) # gives NaN

we can inspect the internal tapes by

GetTape(obj, "ADFun")$print() ## Function value tape - No NaNs here
GetTape(obj, "ADGrad")$print() ## Gradient tape - NaNs in the derivative column!

We then search for the first NaN running backwards. It appears right after the PowOp operator. A closer inspection tells us that the problem is the reverse derivatives of pow(x,y) at y=0. This operation is generated by 'mistake' because the TMBad::pow operator lacks the constant operand optimization, i.e. pow(x,2) should only update the first operand in the reverse pass.

This should obviously be fixed in TMB, but an obvious workaround for this case is to replace x^2 by x*x:

g <- RTMB::MakeTape(function(x) {dx<-x[1]-x[2]; dx*dx}, rep(1,2))

[phillipbvetter]

Thanks for the fast reply @kaskr, and the quick fix.

I can confirm that this has solved my own issue with vectors/matrices where MakeTape$newton() is used in an Implicit Euler ODE solver.

    # function for MakeTape
    maketape_ode_integrator <- function(x){
      # extract elements of x
      ...

      # implicit euler step
      # vector
      X1 = stateVec_rootfind - (stateVec + f__(stateVec_rootfind, parVec, inputVec) * dt)
      # matrix
      P1 = covMat_rootfind - (covMat + cov_ode_1step(covMat_rootfind, stateVec_rootfind, parVec, inputVec) * dt)

      # squared sum for minimization (newton)
      # cost <- sum(X1^2) + sum(c(P1)^2) # this did not work
      cost <- sum(X1*X1) + sum(c(P1*P1)) #this does work!
      return(cost) 
    }

ode_integrator0 <- RTMB::MakeTape(maketape_ode_integrator, ...)$newton(...)