From kkr at dfu.min.dk  Thu Oct 27 18:50:49 2005
From: kkr at dfu.min.dk (Kasper Kristensen)
Date: Thu, 27 Oct 2005 12:50:49 -0400
Subject: [Rd] do_optimhess: ndeps parscale
Message-ID: <1130431849.1233.47.camel@box.dfu.local>

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



[Rd] do_optimhess: ndeps parscale
Kasper Kristensen