[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 20:30:49 CEST 2004
At 11:25 28/07/2004, Roger D. Peng wrote:
>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 : 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
>>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 : stephane.dray at umontreal.ca
--------------------------------------------------------------------------------------------------
Web http://www.steph280.freesurf.fr/
More information about the R-help
mailing list