thr3ads.net - R help - [R] a bug with LAPACK ? non orthogonal vectors obtained with eigen of a symmetric matrix [Jul 2004]

If this information is useful, please help other people find it:
Share via:

Stephane DRAY

2004-Jul-20 14:44 UTC

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

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.

Thanks in advance.

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/

Stephane DRAY

2004-Jul-20 16:57 UTC

head link

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

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.
>
>Thanks in advance.
>
>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/

Stephane DRAY

2004-Jul-28 14:57 UTC

head link

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

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/

Roger D. Peng

2004-Jul-28 15:25 UTC

head link

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

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

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
>

R help - Jul 2004 - a bug with LAPACK ? non orthogonal vectors obtained with eigen of a symmetric matrix

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

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

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

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