[R] optim bug/help?

[Ben Bolker]

 <rkevinburton <at> charter.net> writes:

> 
> In the documentation for 'optim' it gives the following function:
> 
> fr <- function(x) {   ## Rosenbrock Banana function
>     x1 <- x[1]
>     x2 <- x[2]
>     100 * (x2 - x1 * x1)^2 + (1 - x1)^2
> }
> optim(c(-1.2,1), fr)
> 
> When I run this code I get:
> 
> $par
> [1] 1.000260 1.000506
> 
> I am sure I am missing something but why isn't 1,1 a better answer? If I plug
1,1 in the function it seems that
> the function is zero. Whereas the answer given gives a function result of
8.82e-8. This was after 195 calls
> to the function (as reported by optim). The documentation indicates that the
'reltol' is about 1e-8. Is
> this the limit that I am bumping up against?
> 
> Kevin
> 

  Yes, this is basically just numeric fuzz, the
bulk of which probably comes from finite-difference
evaluation of the derivative.
As demonstrated below, you can get a lot closer
by defining an analytic gradient function.
May I ask why this level of accuracy is important?
(Insert Tukey quotation here about approximate answers
to the right question here ...)

  Ben Bolker
---------

fr <- function(x) {   ## Rosenbrock Banana function
    x1 <- x[1]
    x2 <- x[2]
    100 * (x2 - x1 * x1)^2 + (1 - x1)^2
}

## gradient function
frg <- function(x) {
  x1 <- x[1]
  x2 <- x[2]
  c(100*(4*x1^3-2*x2*2*x1)-2+2*x1,
    100*(2*x2-2*x1^2))
}

## use numericDeriv to double-check my calculus
x1 <- 1.5
x2 <- 1.7
numericDeriv(quote(fr(c(x1,x2))),c("x1","x2"))
frg(c(x1,x2))

##
optim(c(-1.2,1), fr)  ## Nelder-Mead
optim(c(-1.2,1), fr,method="BFGS")
optim(c(-1.2,1), fr, gr=frg,method="BFGS")  ## use gradient