Calling optim() within an optim() call leads to R crash

[valentingar]

Hello everyone!

R crashes when I try to call optim() within another optim() call. I am hoping that one of you may be able explain the reason for the cash or even hint toward a solution. I have attached a reprex below for illustration of the general problem. I am well aware that optim() can estimate multiple parameters at the same time, this is just for illustration.

The reason I am using optim() within optim() is the following: I have written code, lets call it fun_1(x_i, dat_i, y), that can automatically optimize some input parameters x_i to fit a very specific dataset dat_i. I have multiple datasets that all have different x_i but share a constant second variable y. Using a different dataset dat_effect I can now fit the output of fun_1() for all the first datasets dat_i to find the best value of y:

# pseudo code: ----------

fun_0 <- function(x, y){
# some output
}

fun_1 <- function(x_i, dat_i, y){

# get best x value:
optim(par = x0, fn = fun_0, y = y, ...) # calls f_0 and compares output to dat_i

# output fun_0  for optimized parameter set
}

fun_2 <- function(y_k, dat_i, dat_effect){
#using y_k, call fun_1 on all dat_i to get the x_i
# compare to dat_effect
}

# find best y value:
optim(par = y0, fn = fun_2, ...) # this will crash, presumably 
                                 # because you cannot call optim() within optim()?

In my case, fun_1() optimizes parameters at distinct time points, so that I get a time series for the output fun_1(x_i, dat_i, y) and the optimized parameters x_i . I would like to use fun_2() to then optimize yas well, that is assumed to be constant. I may have several hundred "separate" datasets dat_i each with distinct parameters x_i which makes combined optimization unfeasible.

Does anyone have an idea on how to get "recursive" optimization working? Do you have to run the different optimization in different environments or something like that? Or am I missing something obvious? I am also interested in knowing why it does not work as shown below (though I have absolutely no experience of the inner workings of optim()).

Thank you so much in advance!

Best,
Valentin

REPREX
Just for showing that calling optim() within itself crashes R.
!!this will crash R!!

# define datasets -----------------
x <- seq(0,1,0.1)
a <- 23
b <- 2

y1 <- a^x
y2 <- b*y1

# define objective functions -------
# inner function
obj_fun_1 <- function(X, 
                      x,
                      y1){
  yi <- X^x
  RMSE <- sqrt(mean((yi-y1)^2))
}

# outer function (optimises inner function first)
obj_fun_2 <- function(Y, 
                      x, 
                      y1,
                      y2){

# get an estimate for the first parameter
a <- 
optim(1,
      obj_fun_1, 
      method = "L-BFGS-B",
      upper = 100,
      lower = 1,
      x = x, 
      y1 = y1)$par

# use estimate to complete function
yi <- Y*a^x
RMSE <- sqrt(mean((yi-y2)^2))
}


optim(1,
      obj_fun_2, 
      method = "L-BFGS-B",
      upper = 100,
      lower = 1,
      x = x, 
      y1 = y1,
      y2 = y2)

^{Created on 2022-04-05 by the reprex package (v2.0.1)}

[arthur.t]

I don't know the answer to your question, but I do think you can change your approach to optimizing all parameters simultaneously. optim can work on a vector of values.

[nirgrahamuk]

That is unexpected behaviour. The easiest 'fix' is to use method 'Brent', at least in the outer loop, or remove the method setting on the outerloop altogether (alongside omitting setting of upper and lower on the outer loop)

[valentingar]

Hi,

Sorry for the late reply!

@arthur.t

I don't know the answer to your question, but I do think you can change your approach to optimizing all parameters simultaneously. optim can work on a vector of values.

Thanks for the reply! Yes, that is correct and definitely the way to go for the REPREX in my post. For my special case that I tried to explain above this would not be possible.

@nirgrahamuk

That is unexpected behaviour. The easiest 'fix' is to use method 'Brent', at least in the outer loop, or remove the method setting on the outerloop altogether (alongside omitting setting of upper and lower on the outer loop)

Thanks for the reply! Changing the method in the outer function helped. Only removing the upper and lower bounds but keeping the method "L-BFGS-B", it works once but leads to a crash when evaluated multiple times. So I guess the two instances of "L-BFGS-B" are somehow interfering each other.

This gave me an idea: what if I use callR to run the internal optimisation in another R session? This seems to work, however is really slow, as a new R session is created each time the outer optim() call evaluates its function (code below). Does anyone have another idea on how to untangle the "L-BFGS-B" calls?

x <- seq(0,1,0.1)
a <- 23
b <- 2

y1 <- a^x
y2 <- b*y1

# define objective functions -------
# inner function
obj_fun_1 <- function(X, 
                      x,
                      y1){
  yi <- X^x
  RMSE <- sqrt(mean((yi-y1)^2))
}

# optim call as a function:
opt_a <- 
  function(obj_fun, 
         x, 
         y1){
  optim(1,
        obj_fun, 
        method = "L-BFGS-B",
        upper = 100,
        lower = 1,
        x = x, 
        y1 = y1)
    }

# outer function (optimises inner function first)
obj_fun_2 <- function(Y, 
                      x, 
                      y1,
                      y2){
  
  # get an estimate for the first parameter
  opt_res <- 
  callr::r(opt_a,
    args = list(obj_fun = obj_fun_1,
                x = x, 
                y1 = y1))
  
  a <- opt_res$par
  
  # use estimate to complete function
  yi <- Y*a^x
  RMSE <- sqrt(mean((yi-y2)^2))
}


optim(1,
      obj_fun_2, 
      method = "L-BFGS-B",
      upper = 100,
      lower = 1,
      x = x, 
      y1 = y1,
      y2 = y2)
#> $par
#> [1] 2
#> 
#> $value
#> [1] 1.954476e-10
#> 
#> $counts
#> function gradient 
#>       38       38 
#> 
#> $convergence
#> [1] 52
#> 
#> $message
#> [1] "ERROR: ABNORMAL_TERMINATION_IN_LNSRCH"

^{Created on 2022-04-11 by the reprex package (v2.0.1)}