[R] convergence=0 in optim and nlminb is real?

[Adelchi Azzalini]

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