[R] eigen() may fail for some symmetric matrices, affects mvrnorm()
Jouni Kerman
jpk2005 at columbia.edu
Sun May 1 04:19:30 CEST 2005
Hi all,
Recently our statistics students noticed that their Gibbs samplers were
crashing due to some NaNs in some parameters. The NaNs came from
mvrnorm (Ripley & Venables' MASS package multivariate normal sampling
function) and with some more investigation it turned out that they were
generated by function eigen, the eigenvalue computing function. The
problem did not seem to happen when the option "symmetric" was set to
FALSE; at first I thought that only matrices that were not completely
symmetric (such as ones with small errors of the magnitude 1e-17 or so)
were affected but then we encountered matrices that were perfectly
symmetric but eigen still generated NaNs when symmetric=TRUE.
I wrote a trivial patch to mvrnorm to detect this problem and to
recompute the eigenvalues with symmetric=TRUE and EISPACK=TRUE. This
seems to work. symmetric=FALSE also seems to work, but is somewhat
slower.
If you use mvrnorm then you may want to try this. Details and the patch
are on my web page, along with one example of a problematic symmetric
matrix (to be loaded with the 'load' function). The NaNs should appear
in the eigenvectors. This particular matrix may not produce an error on
all platforms.
http://www.stat.columbia.edu/~kerman/Software/rprog.html
Thanks to Radford Neal who suggested on our research web log (link
below) that setting EISPACK=TRUE may help.
http://www.stat.columbia.edu/~cook/movabletype/mlm/
Jouni Kerman
Department of Statistics
Columbia University
New York
More information about the R-help
mailing list