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

[Prof Brian Ripley]

I do not begin to understand this.

mvrnorm() uses EISPACK=TRUE and always has done, explicitly since there
was such an argument to eigen() (since 15 Jan 2003: the argument was
introduced in R 1.7.0 in April 2003).  Just take a look at the code.

The claimed `patch' (which is not a patch at all but a replacement) is not 
based on any remotely current version of MASS.

The puported example against eigen(EISPACK=FALSE) (the default) works 
correctly (no NaNs) on all my systems, including the distributed version 
of rw2010.  It seems that students (and `Jouni Kerman' is listed as a 
student) at the Department of Statistics, Columbia University are provided 
with a deficient computing environment (unstated despite the posting 
guide): my surmise is that a broken external BLAS or LAPACK is in use.


On Sat, 30 Apr 2005, Jouni Kerman wrote:

>
> 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
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>

-- 
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