Skip to content
Prev 4411 / 63421 Next

bounds violations, infinite loops in optim/L-BFGS-B (PR#671)

Different day, more or less the same story -- turning off optimization
(and leaving out -ffloat-store) leads to hanging in different places, but
not any fewer hanging cases out of 1000 bootstrap samples.  I could try
the no optimization/-ffloat-store combination, but I strongly suspect the
results would be the same.  (As I mentioned before, I also get hanging in
different places with and without adding parscale to the optim() control
list.)
  These aren't pathological samples, as far as I can tell.
  I'm afraid that fussing with compiler options isn't going to do it ...
given enough time I might be able to dig into the lbfgsb.c code and figure
out what ABNORMAL_TERMINATION_IN_LNSRCH means and (at least) how to modify
the iteration counter so that optim() bails out when it gets into this
kind of loop rather than looping indefinitely.
  On the other hand, I might choose to work on adapting a different set of
code -- the "nlminb" function in Splus5 appears to be based on a different
set of Netlib FORTRAN code (dmnfb,dmngb,dmnhb) -- don't know the
particular reasons for using a different algorithm here.

  I managed to strip down the full code that I use for doing the
bootstraps to about 58 lines (and I removed the embarrassing "<<-"
constructs), so that it if anyone wants to work on this they don't have to
come back to me for particular examples that hang for a particular
platform/set of compiler options/etc..  (It's still not pretty.)

  Ben Bolker

# Stripped-down file to find hanging optim() 

total <- rep(c(100,150,200,300,500,900),rep(2,6))
pred <- c(83,73,22,47,53,31,56,83,30,73,44,91)

nllfun2g.boot <- function(q,debug=FALSE,data.pred,data.total) {
  N1 <- q[1]
  N2 <- q[2]
  G <- q[3]
  maxval <- 1e8
  L <- cfunc2(data.total,N1,N2,G)
  r <- -sum(dbinom(data.pred, data.total, L, log=TRUE))
  r[!is.finite(r)] <- maxval*sign(r[!is.finite(r)])
  r
}

cfunc2 <- function(input,N1,N2,G=1.0,debug=FALSE){
  if (debug)cat(c(N1,N2,G),"\n")
  pmin(1,1/((1+((input - N1)/N2))^G))
}

fuzz <- 0.001

do.boot <- function(tot=1000,seed=101) {
  bootmat <- matrix(nrow = tot, ncol = 5)
  set.seed(seed)
  i <- 1
  while (i<=tot) {
    print(i)
    boot.ind <- sample(1:length(pred), rep = TRUE)
    boot.pred <- pred[boot.ind]  ## kluge!
    boot.total <- total[boot.ind]
    n <- try(optim(c(min(boot.total)-1,100,1),nllfun2g.boot,
                   lower=rep(fuzz,3),upper=c(min(boot.total)-fuzz,Inf,Inf),
                   method="L-BFGS-B",control=list(trace=FALSE,
                                       parscale=c(100,100,1)),
                   data.pred=boot.pred,
                   data.total=boot.total))
     if (is.null(class(n))) {
      bootmat[i,] <- c(n$val,n$par, n$convergence)
      i <- i+1
    }
  }
  bootmat
}

bootmat <- do.boot(5)

## find looping bootstrap sample

find.boot <- function(n,seed=101,cat=FALSE) {
  set.seed(seed)
  for (i in 1:n)
    boot.ind <- sample(1:length(pred), rep = TRUE)
  boot.pred <- pred[boot.ind]  ## kluge!
  boot.total <- total[boot.ind]
  if (cat) {
    cat(paste("boot.pred <- c(",paste(boot.pred,collapse=","),")",sep=""),"\n")
    cat(paste("boot.total <- c(",paste(boot.total,collapse=","),")",sep=""),"\n")
  }
  else list(boot.pred,boot.total)
}
On 3 Oct 2000, Peter Dalgaard BSA wrote: