Implementing a fixed-point solver using ADjoint using optimizers

[Phillip Vetter]

Hi,

Im trying to understand the capabilities of ADjoint, here by trying to implement a fixed-point solver with stats::nlminb - towards more robustness (by using the hessian) than what MakeTape$newton seems to provide(?), but I'm running into some problems. 

My understanding of the capabilities of ADjoint is clearly limited - what am I (not) doing wrong here? I provide a "minimal" example below with 2 ADjoint attempts that fail.

library(RTMB)

# We consider test function f(x) = theta * (x - mu) with root x = mu
mu <- 10 #root


# 1. Use built-in newton() but less robust!:
# ----------------------------------------------------
.f <- function(x) x[1] * (x[2] - x[3])
f <- MakeTape(function(x) .f(x)^2, numeric(3))
g <- MakeTape(function(x) f$newton(2)(x), numeric(2))
# For small theta (<= 1e-5), $newton() fails
g(c(scale=1e-5, root=mu)) # fail
g(c(scale=1e-4, root=mu)) # success
# ----------------------------------------------------


ff <- function(scale, root) function(x) f(c(scale,x,root))
tape <- MakeTape(function(x) ff(scale=1e-5, root=mu)(x), 1)
grad <- tape$jacfun()
hess <- grad$jacfun()
# 2. Use nlminb - more robust, but adjoint approach fails
# ----------------------------------------------------
stats::nlminb(0, tape, grad, hess)$par # works
# but using nlminb will require specifying the adjoint:
minimizer <- ADjoint(
  function(x) stats::nlminb(x, tape, grad, hess)$par,
  function(x,y,dy) x # wrong adjoint, shape correct, for test
)
# but the tape fails
try(MakeADFun(function(p) minimizer(p), 1))
# Error in stats::nlminb(x, tape, grad, hess) : REAL() can only be applied to a 'numeric', not a 'complex'
# In addition:
# Warning messages:
# 1: In stats::nlminb(x, tape, grad, hess) :
#   imaginary parts discarded in coercion
# 2: In stats::nlminb(x, tape, grad, hess) :
#   imaginary parts discarded in coercion
# ----------------------------------------------------


# 3. use ad-friendly (no if's) minimizer
# ----------------------------------------------------
quasi.newton.minimizer <- function(x, grad, hess, iter = 10) {
  for (i in seq_len(iter)) {
    g <- grad(x)
    dim(g) <- rev(dim(g)) # row-vec to col-vec
    H <- hess(x)
    x <- x - solve(H,g)
    }
  x
}
minimizer <- ADjoint(
  function(x) quasi.newton.minimizer(x, grad, hess),
  function(x,y,dy) x # wrong adjoint, shape correct, for test
)
# but the tape fails
try(MakeADFun(function(p) minimizer(p), 1))
# Error: EvalOp: Function must return 'real' or 'integer'
# In addition:
# Warning message:
# In TapedEval(F, x) : imaginary parts discarded in coercion
# ----------------------------------------------------


# 4. Differentiating the self-written minimizer directly works
# ----------------------------------------------------
obj <- MakeADFun(function(p) quasi.newton.minimizer(p, grad, hess), 1)
obj$fn(0) # finds root mu regardless of input
# ----------------------------------------------------

[Kasper Kristensen]

1.

The tape$newton has a number of control options:

https://kaskr.github.io/adcomp/structnewton_1_1newton__config.html

I used the 'trace' argument to realize that a lower gradient tolerance 'grad_tol' was needed by your example:

g <- MakeTape(function(x) f$newton(2, grad_tol=1e-12)(x), numeric(2))

2.

The failing ADjoint was a bug in RTMB (now fixed!).

[Phillip Vetter]

Thanks a lot for the quick fix, and the insight into newton behaviour / options!