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