pirls doesn't work with bobyqa in nloptwrap

[stevencarlislewalker]

Background on this problem is given in issue #2, which is now closed.

Things we know/suspect (see #2 for more info on each topic):

1. minqa bobyqa works
2. nloptwrap bobyqa gives a strange error message
3. nloptwrap bobyqa seem to be handling environments properly
4. we doubt the problem is with NLopt itself

[stevencarlislewalker]

Here is the test example (FYI @bbolker):

library(lme4)
library(nloptwrap) # or library(minqa) for an example that works
library(lme4pureR)

# model structure
form <- cbind(incidence, size - incidence) ~ period + (1 | herd)
glmod <- glFormula(form, cbpp, binomial)

# get initial values
gm1 <- glm(nobars(form), binomial, cbpp)
weights <- weights(gm1)
beta0 <- coef(gm1)
ratios <- gm1$y
eta <- gm1$linear.predictors

# create deviance function with `lme4pureR`
devf <- pirls(glmod, ratios, binomial(), weights=weights, eta=eta, nAGQ=1)

# run `bobyqa`
bobyqa(c(1,beta0), devf, lower=c(0,rep.int(-Inf,length(beta0))))$par

The main result is that bobyqa returns the initial beta (beta0) in the nloptwrap implementation but optimizes to the correct value in the minqa version:
With nloptwrap

> bobyqa(c(1,beta0), devf, lower=c(0,rep.int(-Inf,length(beta0))))$par
[1]  1.000000 -1.269023 -1.170763 -1.301405 -1.782279

With minqa:

> bobyqa(c(1,beta0), devf, lower=c(0,rep.int(-Inf,length(beta0))))$par
[1]  0.6422594 -1.3985243 -0.9923337 -1.1286733 -1.5803150

[bbolker]

I can narrow this down a lot, it looks like a bug (or at least a misfeature??) in the nlopt wrapper. In particular, it doesn't like -Inf lower bounds:

goodpars <- c(0.6422594,-1.3985243,-0.9923337,-1.1286733,-1.5803150)
testfun <- function(minval=-Inf,  ## lower bounds on beta
                    bigval=stop("failed")) {     ## value returned on try-error
    devf2 <- function(x) {
        v <- try(devf(x),silent=TRUE)
        r <- if (is(v,"try-error")) bigval else v
        r
    }
    b <- bobyqa(c(1,beta0), devf2, lower=c(0,rep.int(minval,length(beta0))))$par
    cat("value: ",b,"\n")
    all.equal(b,goodpars)
}

testfun()  ## default: bogus answer, stops at iteration 1
testfun(minval=-5)  ## works
testfun(minval=-10)  ## works
testfun(minval=-20)  ## works
testfun(minval=-1e20)  ## fails
## try to address failure by dropping in a 'big' answer 
testfun(minval=-1e20,bigval=1e20)  ## gets to a bogus answer
testfun(minval=-1e20,bigval=NA)  ## fails
testfun(minval=-1e20,bigval=100)  ## bogus
testfun(minval=-1000,bigval=100)  ## bogus

I wonder if we can make up a trivial test case that lets the nlopt[wrap] developers see the problem without having to install lme4 etc. ???

[stevencarlislewalker]

Thanks Ben. This makes sense given that we were getting errors about bounds. Here's a simple example that seems to isolate the problem:

library(nloptwrap)
f <- function(x) sum(x^2)
bobyqa(1, f) # works
bobyqa(1, f, lower = -Inf) # works
bobyqa(c(1, -1), f, lower = c(-Inf, -Inf)) # fails
optim(c(1, -1), f, lower = c(-Inf, -Inf), method = "L-BFGS-B") # works

And here's the output:

> library(nloptwrap)
> f <- function(x) sum(x^2)
> bobyqa(1, f) # works
$par
[1] 0

$value
[1] 0

$iter
[1] 13

$convergence
[1] 4

$message
[1] "NLOPT_XTOL_REACHED: Optimization stopped because xtol_rel or xtol_abs (above) was reached."

> bobyqa(1, f, lower = -Inf) # works
$par
[1] 0

$value
[1] 0

$iter
[1] 13

$convergence
[1] 4

$message
[1] "NLOPT_XTOL_REACHED: Optimization stopped because xtol_rel or xtol_abs (above) was reached."

> bobyqa(c(1, -1), f, lower = c(-Inf, -Inf)) # fails
$par
[1]  1 -1

$value
[1] Inf

$iter
[1] 0

$convergence
[1] -2

$message
[1] "NLOPT_INVALID_ARGS: Invalid arguments (e.g. lower bounds are bigger than upper bounds, an unknown algorithm was specified, etcetera)."

> optim(c(1, -1), f, lower = c(-Inf, -Inf), method = "L-BFGS-B") # works
$par
[1] -4.235165e-21  4.235165e-21

$value
[1] 3.587324e-41

$counts
function gradient 
       4        4 

$convergence
[1] 0

$message
[1] "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL"

[stevencarlislewalker]

I think the above conclusively proves that this is not an lme4pureR problem.

[bbolker]

But (1) how shall we resolve it (e.g. use minqa::bobyqa until such time as nloptwrap::bobyqa is fixed)? (2) Should we submit this as a bug report to the nloptwrap folks? (It's not 100% clear to me whether this is an NLoptr or an nloptwrap problem.) (3) Doug, does the equivalent test work with the Julia wrapper? (4) it's interesting/suprising to me that setting bad boundaries affects the results even in cases where the bounds shouldn't be active ...

[stevencarlislewalker]

In my view the goal of lme4pureR is not to be perfectly compatible with all possible optimizers, and so I consider the matter closed because I don't mind using minqa over nloptwrap. I've started an email to the nloptwrap maintainer about the issues we've encountered. After that I think there's nothing more to be done on this issue. But I could be missing reasons for why its important to have nloptwrap::bobyqa working with lme4pureR.

[bbolker]

IMO using minqa::bobyqa is fine. Just wanted to make sure we had done whatever follow-up we intended before dropping this ...

[dmbates]

I was unable to get the example to run in Julia and have created an issue in stevengj/NLopt.jl