Re: [R] Calculating large determinants

From: Ravi Varadhan <>
Date: Wed, 5 Dec 2007 09:30:56 -0500

Hi Murray,

A likely reason for the observed information matrix not to be positive definite is the inaccuracies in the numerical estimation of the scores. You might want to try more accurate methods (e.g. Richardson extrapolation) for approximating numerical derivatives, such as available in the package "numDeriv". This would likely alleviate your problem.

Please disregard my advice, if you are using analytical, closed-form expressions for the scores. In this case, you might try the approaches suggested by Martin.


Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625



-----Original Message-----
From: [] On Behalf Of Martin Maechler
Sent: Wednesday, December 05, 2007 4:42 AM To:
Subject: Re: [R] Calculating large determinants

>>>>> "MJ" == maj <>
>>>>> on Wed, 5 Dec 2007 14:18:23 +1300 (NZDT) writes:

    MJ> I apologise for not including a reproducible example
    MJ> with this query but I hope that I can make things clear
    MJ> without one.

    MJ> I am fitting some finite mixture models to data. Each
    MJ> mixture component has p parameters (p=29 in my
    MJ> application) and there are q components to the     MJ> mixture. The number of data points is n ~ 1500.

    MJ> I need to select a good q and I have been considering model selection

    MJ> methods suggested in Chapter 6 of
    MJ> @BOOK{mp01,
    MJ> author    = {McLachlan, G. J. and Peel, D.},
    MJ> title     = {Finite Mixture Models},
    MJ> publisher = {Wiley},
    MJ> address   = {New York},
    MJ> year      = {2001}
    MJ> }

    MJ> One of these methods involves an "empirical information
    MJ> matrix" which is the matrix of products of parameter
    MJ> scores at the observation level evaluated at the MLE and
    MJ> summed over all observations. For example a six-component
    MJ> mixture will have 6 - 1 + 29*6 = 179 parameters. So for
    MJ> observation i I form the 179 by 179 matrix of products of     MJ> scores and sum these up over all 1500-odd observations.

    MJ> Actually it is the log of the determinant of the resultant matrix that I

    MJ> really need rather than the matrix itself. I am seeking advice on what may

    MJ> be the best way to evaluate this log(det()). I have been encountering

    MJ> problems using
    MJ> determinant(SS,logarithm=TRUE)

    MJ> and eigen(SS,only.values = TRUE)$values

    MJ> shows some negative eigenvalues.

which is a problem?
In that case I guess your problem is that you want to estimate a positive definite matrix S but your estimate S^ is not quite positive definite.

Function posdefify() in CRAN package "sfsmisc" provides an old cheap solution to this problem,
where nearPD() in package 'Matrix' (based on a donation from Jens Oehlschlaegel) provides a more sophisticated algorithm for this problem.

If you really only need the eigenvalues of the "corrected" matrix, you might want to abbreviate the nearPD() function by just returning the final 'd' vector of eigenvalues.

    MJ> Advice will be gratefully received!

I'll be glad to hear if and how you'd use one of these two functions.

Martin Maechler, ETH Zurich

    MJ> Murray Jorgensen mailing list PLEASE do read the posting guide and provide commented, minimal, self-contained, reproducible code. mailing list PLEASE do read the posting guide and provide commented, minimal, self-contained, reproducible code. Received on Wed 05 Dec 2007 - 14:38:14 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Wed 05 Dec 2007 - 15:30:17 GMT.

Mailing list information is available at Please read the posting guide before posting to the list.