[R] eigen() may fail for some symmetric matrices, affects mvrnorm()

[Jouni Kerman]

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