R help - May 2005 - eigen() may fail for some symmetric matrices, affects mvrnorm()

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

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

Hi all,

Professor Ripley is right. This problem does not affect newer versions 
of mvrnorm. We had an older version of mvrnorm which did not set 
EISPACK=TRUE as it should. Sorry for the confusion. Thank you for the 
correction.

Jouni Kerman

On May 1, 2005, at 6:08 AM, Prof Brian Ripley <ripley at stats.ox.ac.uk> 
wrote:
> 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.