[Rd] Bug in optim for specific orders of magnitude

[J C Nash]

Extreme scaling quite often ruins optimization calculations. If you think available methods
are capable of doing this, there's a bridge I can sell you in NYC.

I've been trying for some years to develop a good check on scaling so I can tell users
who provide functions like this to send (lots of) money and I'll give them the best answer
there is (generally no answer at all). Or, more seriously, to inform them that they should
not expect results unless they scale. Richard Varga once said some decades ago that any
problem was trivially solvable in the right scale, and he was mostly right. Scaling is
important.

To see the range of answers from a number of methods, the script below is helpful. I had
to remove lbfgsb3c from the mix as it stopped mid-calculation in unrecoverable way. Note
that I use my development version of optimx, so some methods might not be included in
CRAN offering. Just remove the methods from the ameth and bmeth lists if necessary.

Cheers, John Nash

# CErickson221223.R
# optim(c(0,0), function(x) {x[1]*1e-317}, lower=c(-1,-1), upper=c(1,1),
#      method='L-BFGS-B')

tfun <- function(x, xpnt=317){
       if ((length(x)) != 2) {stop("Must have length 2")}
       scl <- 10^(-xpnt)
       val <- x[1]*scl # note that x[2] unused. May be an issue!
       val
}
gtfun <- function(x, xpnt=317){ # gradient
   scl <- 10^(-xpnt)
   gg <- c(scl, 0.0)
   gg
}


xx <- c(0,0)
lo <- c(-1,-1)
up <- c(1,1)
print(tfun(xx))
library(optimx)
ameth <- c("BFGS", "CG", "Nelder-Mead", "L-BFGS-B", "nlm", "nlminb",
              "Rcgmin", "Rtnmin", "Rvmmin", "spg", "ucminf", "newuoa", "bobyqa",
              "nmkb", "hjkb",  "hjn", "lbfgs", "subplex", "ncg", "nvm", "mla",
              "slsqp", "anms")

bmeth <- c("L-BFGS-B", "nlminb", "Rcgmin", "Rtnmin", "nvm",
             "bobyqa", "nmkb", "hjkb", "hjn", "ncg", "slsqp")

tstu <- opm(x<-c(0,0), fn=tfun, gr=gtfun, method=ameth, control=list(trace=0))
summary(tstu, order=value)

tstb <- opm(x<-c(0,0), fn=tfun, gr=gtfun, method=bmeth, lower=lo, upper=up,
             control=list(trace=0))
summary(tstb, order=value)


On 2022-12-23 13:41, Rui Barradas wrote:
> Às 17:30 de 23/12/2022, Collin Erickson escreveu:
>> Hello,
>>
>> I've come across what seems to be a bug in optim that has become a nuisance
>> for me.
>>
>> To recreate the bug, run:
>>
>> optim(c(0,0), function(x) {x[1]*1e-317}, lower=c(-1,-1), upper=c(1,1),
>> method='L-BFGS-B')
>>
>> The error message says:
>>
>> Error in optim(c(0, 0), function(x) { :
>>    non-finite value supplied by optim
>>
>> What makes this particularly treacherous is that this error only occurs for
>> specific powers. By running the following code you will find that the error
>> only occurs when the power is between -309 and -320; above and below that
>> work fine.
>>
>> p <- 1:1000
>> giveserror <- rep(NA, length(p))
>> for (i in seq_along(p)) {
>>    tryout <- try({
>>      optim(c(0,0), function(x) {x[1]*10^-p[i]}, lower=c(-1,-1),
>> upper=c(1,1), method='L-BFGS-B')
>>    })
>>    giveserror[i] <- inherits(tryout, "try-error")
>> }
>> p[giveserror]
>>
>> Obviously my function is much more complex than this and usually doesn't
>> fail, but this reprex demonstrates that this is a problem. To avoid the
>> error I may multiply by a factor or take the log, but it seems like a
>> legitimate bug that should be fixed.
>>
>> I tried to look inside of optim to track down the error, but the error lies
>> within the external C code:
>>
>> .External2(C_optim, par, fn1, gr1, method, con, lower,
>>          upper)
>>
>> For reference, I am running R 4.2.2, but was also able to recreate this bug
>> on another device running R 4.1.2 and another running 4.0.3.
>>
>> Thanks,
>> Collin Erickson
>>
>>     [[alternative HTML version deleted]]
>>
>> ______________________________________________
>> R-devel using r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
> 
> Hello,
> 
> See if this R-Help thread [1] earlier this year is relevant.
> In particular, the post by R Core team member Luke Tierney [2], that answers much better than what I could.
> 
> The very small numbers in your question seem to have hit a limit and this limit is not R related.
> 
> 
> [1] https://stat.ethz.ch/pipermail/r-help/2022-February/473840.html
> [2] https://stat.ethz.ch/pipermail/r-help/2022-February/473844.html
> 
> 
> Hope this helps,
> 
> Rui Barradas
> 
> ______________________________________________
> R-devel using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel