Search code examples
roptimization

Problem minimising a function using L-BFGS-B method in R?

A similar question was asked here but it didn't receive any answers.

I have a function that I need to minimize. My function has 3 parameters, x, y, and w. Given w, I need to find the optimal pair of x and y that minimizes my function. However, x and y are restricted in the interval (-0.5 ,0.5]

I tried to use the optim function and it works when I don't have restrictions. However, when I apply the restrictions and use the L-BFGS-B method, I get an error. I cant figure out why I am getting this error.

Below is an example of my function and some of the methods I have tried. This is the function that Im trying to minimize:

fn <- function(par, w) {
  x <- par[1]
  y <- par[2]

  p1 <- w^x * (w + 1) * beta(x + 1, w + 1 - x)
  p2 <- w^y * (w + 1) * beta(y + 1, w + 1 - y)
  p3 <- w^(2 * x) * (w + 1) * (beta(2 * x + 1, w + 1 - 2 * x) - (w + 1) * (beta(x + 1, w + 1 - x)^2))
  p4 <- w^(2 * y) * (w + 1) * (beta(2 * y + 1, w + 1 - 2 * y) - (w + 1) * (beta(y + 1, w + 1 - y)^2))
  p5 <- w^(y + x) * (w + 1) * (beta(y + x + 1, w + 1 - y - x) - (w + 1) * beta(y + 1, w + 1 - y) * beta(x + 1, w + 1 - x))
  p6 <- (y^2) * p3 * p1^(2 * y - 2) * p2^(-2 * x)
  p7 <- y * x * p1^(2 * y - 1) * p2^(-2 * x - 1) * p5
  p8 <- (x^2) * p4 * p1^(2 * y) * p2^(-2 * x - 2)
  mu <- (p1^y) / (p2^x)
  sigma <- p6 - 2 * p7 + p8
  v <- (gamma(x + 1)^y / gamma(y + 1)^x)

  quant <- ceiling((qnorm(0.95)*sqrt(sigma)/(v-mu))^2)
  
  return(quant)
}

Here is my issue:

# Here is an example of it working when there are no restrictions
oo = optim(par = c(0.1, 0.2), fn, w = 0.1)
oo$par
# In this example, I apply restrictions and it returns an error
oo1 = optim(par = c(0.1, 0.2), fn, w = 0.1, lower = c(-0.49, -0.49), upper = c(0.5, 0.5), method="L-BFGS-B")
oo1$par

In the second example, I apply restrictions and get the following error:

Error in optim(par = c(0.1, 0.2), fn, w = 0.1, lower = c(-0.49, -0.49), : L-BFGS-B needs finite values of 'fn'

Additionally, I also tried the optimx package and it doesn't seem to work. It just returns NAs. For example:

library(optimx)
opt <- optimx(par = c(0.1, 0.2), fn = fn, lower = c(-0.4, -0.4), upper = c(0.5, 0.5), w = 0.1)
opt
> opt
         p1 p2         value fevals gevals niter convcode kkt1 kkt2 xtime
L-BFGS-B NA NA 8.988466e+307     NA     NA    NA     9999   NA   NA 0.001

Does anyone have any idea as to why I cant minimize my function when I apply the restrictions? (and why the optimx package returns NAs?)

Solution

  • Here's a start. I added a line

    cat(x,y, w, quant, "\n")

    right before return(quant). Your first (no-bounds) run gave:

    0.1 0.2 0.1 40 
0.12 0.2 0.1 40 
0.1 0.22 0.1 41 
0.12 0.18 0.1 39 
0.13 0.16 0.1 39 
0.14 0.18 0.1 40 
0.11 0.195 0.1 40 
0.11 0.19 0.1 39 
0.12 0.19 0.1 40 
0.11 0.18 0.1 39 
0.1125 0.1825 0.1 39

    Your second (with bounds) gave:

    0.1 0.2 0.1 40 
0.101 0.2 0.1 40 
0.099 0.2 0.1 39 
0.1 0.201 0.1 40 
0.1 0.199 0.1 39 
-0.49 -0.49 0.1 NaN

    in other words, R's attempt to evaluate the function at the lower bounds gave NaN.

    Poking around some more:

    library(emdbook)
cc <- curve3d(fn(c(x,y), w = 0.1), xlim = c(-0.5,0.5), ylim = c(-0.5, 0.5),
              sys3d = "image")

    image of objective function surface showing white (non-finite) values on x=y diagonal and at x=0 or y=0

    it looks like your function is giving non-finite values if x=y or x=0 or y=0 ? You might want to catch those cases and return a finite value ...

    In fact, that leads to a workaround, which is to not make the lower or upper bounds exactly equal, to stop R trying to evaluate the function in the bad places:

    oo1 = optim(par = c(0.1, 0.2), fn, w = 0.1, 
   lower = c(-0.49, -0.48), 
   upper = c(0.5, 0.49), method="L-BFGS-B")

    It probably would have been possible to diagnose this by looking at the objective function and thinking hard about where it would have non-finite values.