difficult optimization

An example showing that a seemingly simple optimization problem can be tricky.

f <- function (x) {
  y = -cos(3*pi*x)/x
  return(y)
}

This is our fake log-probability function.

plotf <- function() {
  xs <- (0:2000)/1000
  plot(xs, f(xs), t="l", lwd=2, xlab="x", ylab="log-probability", ylim=c(-6,6))
}

plotf()

brent1.5 <- optim(par=1.5, fn=f, method="Brent", lower=0, upper=2, control=list(fnscale=-1))

brent0.5 <- optim(par=0.5, fn=f, method="Brent", lower=0, upper=2, control=list(fnscale=-1))

NelderMead1.5 <- optim(par=1.5, fn=f, method="Nelder-Mead", control=list(fnscale=-1))

## Warning in optim(par = 1.5, fn = f, method = "Nelder-Mead", control = list(fnscale = -1)): l'optimisation de Nelder-Mead unidimensionnelle est instable :
## use "Brent" or optimize() directement

NelderMead0.5 <- optim(par=0.5, fn=f, method="Nelder-Mead", control=list(fnscale=-1))

## Warning in optim(par = 0.5, fn = f, method = "Nelder-Mead", control = list(fnscale = -1)): l'optimisation de Nelder-Mead unidimensionnelle est instable :
## use "Brent" or optimize() directement

plotf <- function() {
  jit_factor = 100
  xs <- (0:2000)/1000
  plot(xs, f(xs), t="l", lwd=2, xlab="x", ylab="log-probability", ylim=c(-6,6))
  segments(x0=1.5, x1=1.5,  y0=-6, y1=f(1.5), col="grey", lwd=2, lty=2)
  segments(x0=0.5, x1=0.5,  y0=-6, y1=f(0.5), col="darkgrey", lwd=2, lty=2)
  }

plotf()

plotf <- function() {
  jit_factor = 100
  xs <- (0:2000)/1000
  plot(xs, f(xs), t="l", lwd=2, xlab="x", ylab="log-probability", ylim=c(-6,6))
  segments(x0=1.5, x1=1.5,  y0=-6, y1=f(1.5), col="grey", lwd=2, lty=2)
  segments(x0=0.5, x1=0.5,  y0=-6, y1=f(0.5), col="darkgrey", lwd=2, lty=2)
  points((brent1.5$par), jitter(brent1.5$value, factor=jit_factor), col="red", pch=20, cex=2)
  points((brent0.5$par), jitter(brent1.5$value, factor=jit_factor), col="darkred", pch=20, cex=2)
  points((NelderMead1.5$par), jitter(NelderMead1.5$value, factor=jit_factor), col="blue", pch=20, cex=2)
  points((NelderMead0.5$par), NelderMead0.5$value, col="darkblue", pch=20, cex=2)
  legend(legend=c("Brent starting at 1.5", "Brent starting at 0.5", "Nelder-Mead starting at 1.5", "Nelder-Mead starting at 0.5"), col=c("red", "darkred", "blue", "darkblue"), "topright", pch=20)
}

plotf()