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] 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.
Thanks in advance. ++++++++++++++++++++++++++++++++++++
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
Note that I have do the same with random number and never find this kind of problems
> 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"
St�phane DRAY
--------------------------------------------------------------------------------------------------
D�partement des Sciences Biologiques Universit� de Montr�al, C.P. 6128, succursale centre-ville Montr�al, Qu�bec H3C 3J7, Canada
Tel : (514) 343-6111 poste 1233 Fax : (514) 343-2293
E-mail : [EMAIL PROTECTED]
--------------------------------------------------------------------------------------------------
Web http://www.steph280.freesurf.fr/
______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
