It must be the case that this issue has already been rised before,
but I did not manage to find it in past posting.
In some cases, optim() and nlminb() declare a successful convergence,
but the corresponding Hessian is not positive-definite. A simplified
version of the original problem is given in the code which for
readability is placed below this text. The example is built making use
of package 'sn', but this is only required to set-up the example: the
question is about the outcome of the optimizers. At the end of the run,
a certain point is declared to correspont to a minimum since
'convergence=0' is reported, but the eigenvalues of the (numerically
evaluated) Hessian matrix at that point are not all positive.
Any views on the cause of the problem? (i) the point does not
correspong to a real minimum, (ii) it does dive a minimum but the
Hessian matrix is wrong, (iii) the eigenvalues are not right.
...and, in case, how to get the real solution.
Adelchi Azzalini
#------------------
library(sn) # version 0.4-18
data(ais, package="sn")
attach(ais)
X <- cbind(1,Ht,Wt)
y <- cbind(bmi, lbm)
dettach(ais)
negLogLik <- function(vdp, x, y)
{
d <- ncol(y)
p <- ncol(x)
beta <- matrix(vdp[1:(d*p)], p, d)
Omega <- matrix(0, d, d)
Omega[lower.tri(Omega, TRUE)] <- vdp[(p*d+1):(p*d+d*(d+1)/2)]
Omega <- Omega + t(Omega) - diag(diag(Omega), d)
if(any(eigen(Omega)$values <= 0)) return(Inf)
omega <- sqrt(diag(Omega))
alpha <- vdp[(p*d+d*(d+1)/2+1):(p*d+d*(d+1)/2+d)]
nu <- vdp[p*d+d*(d+1)/2+d+1]
if(nu <= 0) return(Inf)
logL <- sum(dmst(y, x %*% beta, Omega, alpha, nu, log=TRUE))
return(-logL)
}
# run 1
vdp <- c(44, 0, 0, -4, 0, 0, 0.05, -0.5, 35, 0.5, -20, 3.5)
opt <- optim(vdp, negLogLik, method="BFGS", hessian=TRUE, x=X, y=y)
opt$value
# [1] 625.3
opt$convergence
# [1] 0
eigen(opt$hessian)$values
# [1] 7.539e+07 1.523e+06 .... 5.684e-02 -4.516e-01
#---
# run 2
vdp <- c(44.17, -0.2441, 0.303, -3.620, 0.04044, 0.8906,
0.0487, -0.5072, 36.33, 0.4445, -20.87, 3.5)
opt <- optim(vdp, negLogLik, method="BFGS", hessian=TRUE, x=X, y=y)
opt$value
# [1] 599.7
opt$convergence
# [1] 0
eigen(opt$hessian)$values
# [1] 1.235e+08 .... 3.845e-02 -1.311e-03 -6.701e+02
#---
# run 3
vdp <- c(44.17, -0.2441, 0.303, -3.620, 0.04044, 0.8906,
0.0487, -0.5072, 36.33, 0.4445, -20.87, 3.5)
opt <- optim(vdp, negLogLik, method="SANN", hessian=TRUE, x=X, y=y)
opt$value
# [1] 599.8
opt$convergence
# [1] 0
eigen(opt$hessian)$values
# [1] 1.232e+08 .... 3.225e-02 -6.681e-02 -7.513e+02
#--
# run 4
vdp <- c(44.17, -0.2441, 0.303, -3.620, 0.04044, 0.8906,
0.0487, -0.5072, 36.33, 0.4445, -20.87, 3.5)
opt <- optim(vdp, negLogLik, method="Nelder", hessian=TRUE, x=X, y=y)
opt$value
# [1] 599.7
opt$convergence
# [1] 1
#--
# run 5
vdp <- c(44.17, -0.2441, 0.303, -3.620, 0.04044, 0.8906,
0.0487, -0.5072, 36.33, 0.4445, -20.87, 3.5)
opt <- optim(vdp, negLogLik, method="CG", hessian=TRUE, x=X, y=y)
opt$value
# [1] 599.7
opt$convergence
# [1] 0
eigen(opt$hessian)$values
# [1] 1.236e+08 3.026e+06 .... 3.801e-02 -2.348e-04 -7.344e+02
#--
# run 6
vdp <- c(44.17, -0.2441, 0.303, -3.620, 0.04044, 0.8906,
0.0487, -0.5072, 36.33, 0.4445, -20.87, 3.5)
opt <- nlminb(vdp, negLogLik, x=X, y=y)
opt$obj
# [1] 599.7
H <- optimHess(opt$par, negLogLik, x=X, y=y)
eigen(H)$values
# [1] 1.236e+08 3.041e+06 ... 4.090e-05 -7.176e+02
=============
_
platform x86_64-unknown-linux-gnu
arch x86_64
os linux-gnu
system x86_64, linux-gnu
status
major 3
minor 0.2
year 2013
month 09
day 25
svn rev 63987
language R
version.string R version 3.0.2 (2013-09-25)
nickname Frisbee Sailing
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
