# [R] a bug with LAPACK ? non orthogonal vectors obtained with eigen of a symmetric matrix

Stephane DRAY stephane.dray at umontreal.ca
Tue Jul 20 18:57:18 CEST 2004

```I have continue my experiments in changing the size of my matrix :
(3^2 by 3^2, 4^2 by 4^2... 20^2 by 20^2)

EISPACK is always correct but LINPACK provide very strange results:

> for(i in 3:20){
+ x=dbcenter(nb2mat(cell2nb(i,i),style="B"))
+ res=eigen(x,EIS=T)
+ eq0 <- apply(as.matrix(res\$values),1,function(x) identical(all.equal(x,
0), TRUE))
+ res0=res\$vec[,-which(eq0)]
+ print(sum((diag(1,ncol(res0))-crossprod(res0))^2))
+ }
[1] 7.939371e-30
[1] 2.268788e-29
[1] 9.237286e-29
[1] 1.806393e-28
[1] 3.24619e-28
[1] 5.239195e-28
[1] 9.78079e-28
[1] 1.315542e-27
[1] 1.838600e-27
[1] 3.114150e-27
[1] 5.499297e-27
[1] 5.471782e-27
[1] 1.075098e-26
[1] 1.534822e-26
[1] 1.771326e-26
[1] 2.342404e-26
[1] 3.462522e-26
[1] 4.310143e-26
> for(i in 3:20){
+ x=dbcenter(nb2mat(cell2nb(i,i),style="B"))
+ res=eigen(x)
+ eq0 <- apply(as.matrix(res\$values),1,function(x) identical(all.equal(x,
0), TRUE))
+ res0=res\$vec[,-which(eq0)]
+ print(sum((diag(1,ncol(res0))-crossprod(res0))^2))
+ }
[1] 1.515139e-30
[1] 1.054286e-27
[1] 9.553017e-29
[1] 2.263455e-28
[1] 5.641993e-27
[1] 4.442088e-26
[1] 3.996714
[1] 3.986387
[1] 3.996545
[1] 7.396718
[1] NaN
[1] 7.980621
[1] 7.996769
[1] 3.984399
[1] NaN
[1] NaN
[1] NaN
[1] NaN

> R.Version()
\$platform
[1] "i386-pc-mingw32"

\$arch
[1] "i386"

\$os
[1] "mingw32"

\$system
[1] "i386, mingw32"

\$status
[1] ""

\$major
[1] "1"

\$minor
[1] "9.1"

\$year
[1] "2004"

\$month
[1] "06"

\$day
[1] "21"

\$language
[1] "R"

At 10:44 20/07/2004, Stephane DRAY wrote:
>Hello,
>I have obtained strange results using eigen on a symmetric matrix:
>
># this function perform a double centering of a matrix
>(xij-rowmean(i)-colmean(j)+meantot)
>dbcenter=function(mat){
>rmean=apply(mat,1,mean)
>cmean=apply(mat,2,mean)
>newmat=sweep(mat,1,rmean,"-")
>newmat=sweep(newmat,2,cmean,"-")
>newmat=newmat+mean(mat)
>newmat}
>
># i use spdep package to create a spatial contiguity matrix
>library(spdep)
>x=dbcenter(nb2mat(cell2nb(3,3),style="B"))
>
>#compute eigenvalues of a 9 by 9 matrix
>res=eigen(x)
>
># some eigenvalues are equal to 0
>eq0 <- apply(as.matrix(res\$values),1,function(x) identical(all.equal(x,
>0), TRUE))
>
># I remove the corresponding eigenvectors
>res0=res\$vec[,-which(eq0)]
>
># then I compute the Froebenius norm with the identity matrix
>sum((diag(1,ncol(res0))-crossprod(res0))^2)
>
>[1] 1.515139e-30
>
># The results are correct,
># then I do it again with a biggest matrix(100 by 100)
>
>x=dbcenter(nb2mat(cell2nb(10,10),style="B"))
>res=eigen(x)
>eq0 <- apply(as.matrix(res\$values),1,function(x) identical(all.equal(x,
>0), TRUE))
>res0=res\$vec[,-which(eq0)]
>sum((diag(1,ncol(res0))-crossprod(res0))^2)
>
>[1] 3.986387
>
>
>I have try the same with res=eigen(x,EISPACK=T):
>
>  x=dbcenter(nb2mat(cell2nb(10,10),style="B"))
>res=eigen(x,EISPACK=T)
>eq0 <- apply(as.matrix(res\$values),1,function(x) identical(all.equal(x,
>0), TRUE))
>res0=res\$vec[,-which(eq0)]
>sum((diag(1,ncol(res0))-crossprod(res0))^2)
>[1] 1.315542e-27
>
>
>So I wonder I there is a bug in the LAPACK algorithm or if there are some
>things that I have not understood. Note that my matrix has some pairs of
>equal eigenvalues.
>
>
>Stéphane DRAY
>--------------------------------------------------------------------------------------------------
>
>Département des Sciences Biologiques
>Université de Montréal, C.P. 6128, succursale centre-ville
>
>Tel : (514) 343-6111 poste 1233         Fax : (514) 343-2293
>E-mail : stephane.dray at umontreal.ca
>--------------------------------------------------------------------------------------------------
>
>Web                                          http://www.steph280.freesurf.fr/
>
>______________________________________________
>R-help at stat.math.ethz.ch mailing list
>https://www.stat.math.ethz.ch/mailman/listinfo/r-help

Stéphane DRAY
--------------------------------------------------------------------------------------------------

Département des Sciences Biologiques
Université de Montréal, C.P. 6128, succursale centre-ville