[Rd] eigen of a real pd symmetric matrix gives NaNs in $vector (PR#7987)

[Patrick Burns]

I would presume this is another manifestation of what I reported
(reproduced below) on 2003-12-01.

cajo.terbraak at wur.nl wrote:

>Full_Name: cajo ter Braak
>Version: 2.1.1
>OS: Windows
>Submission from: (NULL) (137.224.10.105)
>
>
># I would like to attach the matrix C in the Rdata file; it is 50x50 and comes
>from a geostatistical problem (spherical covariogram)
>
>  
>
>>rm(list=ls(all=TRUE))
>>load(file= "test.eigen.Rdata")
>>ls()
>>    
>>
>[1] "C"  "eW"
>  
>
>>sym.check = max(abs(C - t(C)))    # should be 0 for symmetry
>>sym.check
>>    
>>
>[1] 0
>  
>
>>eW <-eigen(C, symmetric = TRUE)
>>l_eW <- eW$values
>>print(eW$values)
>>    
>>
> [1] 4.5730646 4.5730646 3.3066738 3.3066738 3.3066738 3.3066738 2.3935268
> [8] 2.3935268 1.9367508 1.9367508 1.9347787 1.9347787 1.4276845 1.4276845
>[15] 1.4276845 1.4276845 0.9858318 0.9858318 0.9858318 0.9858318 0.9123115
>[22] 0.9123115 0.7945283 0.7945283 0.7880493 0.7880493 0.6047920 0.6047920
>[29] 0.6047920 0.6047920 0.5689609 0.5689609 0.5681210 0.5681210 0.5440676
>[36] 0.5440676 0.5440676 0.5440676 0.5224040 0.5224040 0.5139844 0.5139844
>[43] 0.5077485 0.5077485 0.5008249 0.5008249 0.5008249 0.5008249 0.4960220
>[50] 0.4960220
>  
>
>>#print(eW$vector)
>>
>>
>>#library(MASS)
>>#n = nrow(C)
>>#y = runif(n)
>>
>>#lm1 = lm.gls(y~1, W = C)
>>#summary(lm1)
>>
>>eW <-eigen(C, symmetric = TRUE)
>>l_eW <- eW$values
>># the thirdteens eigenvector contrains NaN
>>print(eW$values)
>>    
>>
> [1] 4.5730646 4.5730646 3.3066738 3.3066738 3.3066738 3.3066738 2.3935268
> [8] 2.3935268 1.9367508 1.9367508 1.9347787 1.9347787 1.4276845 1.4276845
>[15] 1.4276845 1.4276845 0.9858318 0.9858318 0.9858318 0.9858318 0.9123115
>[22] 0.9123115 0.7945283 0.7945283 0.7880493 0.7880493 0.6047920 0.6047920
>[29] 0.6047920 0.6047920 0.5689609 0.5689609 0.5681210 0.5681210 0.5440676
>[36] 0.5440676 0.5440676 0.5440676 0.5224040 0.5224040 0.5139844 0.5139844
>[43] 0.5077485 0.5077485 0.5008249 0.5008249 0.5008249 0.5008249 0.4960220
>[50] 0.4960220
>  
>
>>print(eW$vector[,13])
>>    
>>
> [1]   0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
>NaN
>[20] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
>NaN
>[39] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
>  
>

This is undoubtedly a bug, but I doubt that it is really down to R.

Synopsis: a certain matrix causes eigen(symmetric=TRUE)
to produce NaN's in some of the returned eigenvectors.
This happens using SuSe 8.2 Professional and the precompiled
R rpm (happens in both 1.8.1 and 1.7.1).  I don't see it under
Windows.

To reproduce the bug:
The matrix (75 by 75, about 45K) is in
http://www.burns-stat.com/pages/Flotsam/eigenbugmatrix.RData

Load that into R.

sum(is.na(eigen(eigenbugmatrix, symmetric=TRUE)$vectors))

is non-zero (600 in my experience) when the bug is exhibited and
is zero when the bug is not.

It is quite sensitive to the numbers.  The bug is still there with some
scaling of the matrix (from about divide by 2 to multiply by 8).  The
bug disappears if the matrix is dumped and sourced back in again.

The only clue that I can offer is that it is vectors 43:50 (I think) that
are NaN's and the matrix is logically of rank 50.  That is, it is a
covariance matrix on 75 variables using 50 observations.

Hopefully, someone has the experience and tenacity to figure out
what is going on here.


Patrick Burns

Burns Statistics
patrick at burns-stat.com
+44 (0)20 8525 0696
http://www.burns-stat.com
(home of S Poetry and "A Guide for the Unwilling S User")

>______________________________________________
>R-devel at r-project.org mailing list
>https://stat.ethz.ch/mailman/listinfo/r-devel
>
>
>
>  
>