[R] optim() and resultant hessian

[Prof Brian Ripley]

On Fri, 16 Feb 2007, Brian Healey wrote:

> R users;
>
> A question about optimization within R.
>
> I've been using both optim() and nlminb() to estimate parameters and all
> seems to be working fine. For context (but without getting into specifics -
> sorry), I'm working with a problem that is known to have correlated
> parameters, and parameter estimation can be difficult. I have a question on
> optim() - I'm using method="L-BFGS-B" to accommodate box constraints.
>
> For my dataset, I obtain parameter estimates using a few iterations of
> optim() - iterations in that I'm simply taking the results from a previous
> optim() call and using these as starting values in the next function call.
>
> The final call to optim() returns the following:
> $par
> [1] 0.2272361 0.8037642 26.8591998 3.0631280 0.2224566
> $value
> [1] -46.13906
> $counts
> function gradient
> 4 4
> $convergence
> [1] 0
> $message
> [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
> $hessian
> [,1] [,2] [,3] [,4] [,5]
> [1,] 1.267070e+17 1.012691e+17 1.348054e+15 625551.58724 9.359559e+07
> [2,] 1.012691e+17 8.189877e+16 1.144248e+15 569562.44945 8.699072e+07
> [3,] 1.348054e+15 1.144248e+15 2.457323e+05 3426.60293 -2.297009e+03
> [4,] 6.255516e+05 5.695624e+05 3.426603e+03 99.06880 -6.750806e+01
> [5,] 9.359559e+07 8.699072e+07 -2.297009e+03 -67.50806 1.905247e+03
>
> i.e. convergence and message report that things look "ok".
>
> However; if I take the hessian and compute eigenvalues and eigenvectors, the
> result is:
> $values
> [1] 2.080357e+17 5.889416e+14 1.907746e+03 9.648828e+01 -1.886641e+13
> $vectors
> [,1] [,2] [,3] [,4] [,5]
> [1,] 7.797175e-01 6.246383e-01 -5.024932e-09 -2.214316e-10 -4.321721e-02
> [2,] 6.260739e-01 -7.768497e-01 5.917821e-09 2.637892e-10 6.734993e-02
> [3,] 8.496067e-03 -7.957109e-02 9.710924e-08 4.049691e-09 -9.967930e-01
> [4,] 4.058779e-12 -8.828316e-11 3.729724e-02 -9.993042e-01 -4.192523e-10
> [5,] 6.125907e-10 -1.547717e-08 -9.993042e-01 -3.729724e-02 -9.626469e-08
>
> Optim() indicates convergence when one of the eigenvalues is negative !?!?
> Any reason to not be concerned?

I don't think you have scaled your problem as recommended on the help 
page, or the Hessian would not have elements differing by 15 orders of 
magnitude.  So, yes, there is reason to be worried about the 
finite-difference approximations used, as well as the numerical accuracy 
of such extreme differences.  (If you did scale, you don't want to do this 
on the rescaled Hessian.)

> Possibly the Hessian is simply computed at
> the solution and not used to arrive at the estimate?)

Please do read what it says about the method on the help page and the 
references.

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595