[Rd] do_optimhess: ndeps parscale

Kasper Kristensen kkr at dfu.min.dk
Thu Oct 27 18:50:49 CEST 2005


Dear list

I think there might be an inconsistency in the function "do_optimhess"
in src/main/optim.c. 
Consider changing the line
"eps = OS->ndeps[i]/(OS->parscale[i]);" 
to
"eps = OS->ndeps[i];"
Then "ndeps" is the finite-difference step on "par/parscale" as it
should be according to the documentation.

Motivation:

> ## EXAMPLE
> f <- function(x){sum(exp(.5*(a*x)^2))}
> a <- c(1e-3,1e3)
>
> ## On "par/parscale"-domain we want g(y)=sum(exp(.5*y^2))
> parscale <- 1/a
>
> ## Optimization works only when using a good parscale:
> xstart <- 1/a
> optim(xstart,f,hessian=FALSE,method="CG")$par
[1] 1.000000e+03 7.726989e-09
>
optim(xstart,f,hessian=FALSE,method="CG",control=list(parscale=parscale))$par
[1] -2.321611e-09 -2.277797e-15
>
> ## But even with a good parscale the hessian computation is wrong:
> xhat <- c(0,0)
> optim(xhat,f,hessian=TRUE,method="CG")$hessian
             [,1]    [,2]
[1,] 1.000089e-06       0
[2,] 0.000000e+00 3194528
>
optim(xhat,f,hessian=TRUE,method="CG",control=list(parscale=parscale))$hessian
             [,1]    [,2]
[1,] 1.000000e-06       0
[2,] 0.000000e+00 1648722

##### RECOMPILING WITH THE SUGGESTED CHANGE GIVES CORRECT HESSIAN:
>
optim(xhat,f,hessian=TRUE,method="CG",control=list(parscale=parscale))$hessian
             [,1]    [,2]
[1,] 1.000001e-06       0
[2,] 0.000000e+00 1000001


Best regards
Kasper Kristensen



More information about the R-devel mailing list