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

Stephane DRAY dray at biomserv.univ-lyon1.fr
Wed Jul 28 16:57:20 CEST 2004

```Hello,
I have send send this message one week ago but I have receive no answer.
Perhaps, some of RListers were in holidays and do not read my message. I
try again..
My problem is that I obtained non orthonormal eigenvectors with some
matrices with LAPACK while EISPACK seems to provide "good" results. Is
there some restrictions with the use of LAPACK ? Is it a bug ? I did not
find the answer. Here is my experiment:

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

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

++++++++++++++++++++++++++++++++++++

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))
+ }
 7.939371e-30
 2.268788e-29
 9.237286e-29
 1.806393e-28
 3.24619e-28
 5.239195e-28
 9.78079e-28
 1.315542e-27
 1.838600e-27
 3.114150e-27
 5.499297e-27
 5.471782e-27
 1.075098e-26
 1.534822e-26
 1.771326e-26
 2.342404e-26
 3.462522e-26
 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.515139e-30
 1.054286e-27
 9.553017e-29
 2.263455e-28
 5.641993e-27
 4.442088e-26
 3.996714
 3.986387
 3.996545
 7.396718
 NaN
 7.980621
 7.996769
 3.984399
 NaN
 NaN
 NaN
 NaN

Note that I have do the same with random number and never find this kind of
problems

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

\$arch
 "i386"

\$os
 "mingw32"

\$system
 "i386, mingw32"

\$status
 ""

\$major
 "1"

\$minor
 "9.1"

\$year
 "2004"

\$month
 "06"

\$day
 "21"

\$language
 "R"

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

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