This is interesting. I can reproduce your results but cannot come up with an explanation. However, using svd(LINPACK = FALSE) seems to work every time. Might you consider trying that instead?
-roger
I have try to do it with svd but here the problem is that the matrices under study have some equal singular values and in this case, I do not have u=v for a symmetric matrix (I supose that a linear combination of u corresponding to equal singular values allow to obtain the v for the same singular values).
i.e. if d2=d3=d4
u2=av2+bv3+cv4
... and not necessarily u2=v2
But here, I found another problem:
#This is for svd with LINPACK. It works, u vectors are orthogonal (same for v)
or(i in 3:20){
x=dbcenter(nb2mat(cell2nb(i,i),style="B"))
res=svd(x,LIN=T)
eq0 <- apply(as.matrix(res$d),1,function(x) identical(all.equal(x, 0), TRUE))
resu=res$u[,-which(eq0)]
resv=res$v[,-which(eq0)]
print(paste("Grid size",i))
print(paste("u ortho ?",all.equal(diag(1,ncol(resu)),t(resu)%*%resu)))
print(paste("v ortho ?",all.equal(diag(1,ncol(resv)),t(resv)%*%resv)))
print(paste("u == v ?",all.equal(abs(resv),abs(resu))))}
#But without LINPACK:
for(i in 3:20){
x=dbcenter(nb2mat(cell2nb(i,i),style="B"))
res=svd(x,LIN=F)
eq0 <- apply(as.matrix(res$d),1,function(x) identical(all.equal(x, 0), TRUE))
resu=res$u[,-which(eq0)]
resv=res$v[,-which(eq0)]
print(paste("Grid size",i))
print(paste("u ortho ?",all.equal(diag(1,ncol(resu)),t(resu)%*%resu)))
print(paste("v ortho ?",all.equal(diag(1,ncol(resv)),t(resv)%*%resv)))
print(paste("u == v ?",all.equal(abs(resv),abs(resu))))
}
All is fine, except for i=14 where vectors u (and v) are not orthogonal [1] "Grid size 14" [1] "u ortho ? Mean relative difference: 269.2564" [1] "v ortho ? Mean relative difference: 2097.224" [1] "u == v ? Mean relative difference: 0.7612647"
I wonder if it is a general problem, a problem due to the structure of my matrix (I don't think so) , a windoze problem or a bug in my head !!!
Sincerely.
Stephane DRAY wrote: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] 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
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
