Re: [Rd] eigen(symmetric=TRUE) for complex matrices
On Jun 18, 2013, at 03:30 , robin hankin wrote: > R-3.0.1 rev 62743, binary downloaded from CRAN just now; macosx 10.8.3 > > Hello, > > eigen(symmetric=TRUE) behaves strangely when given complex matrices. > > > The following two lines define 'A', a 100x100 (real) symmetric matrix > which theoretical considerations [Bochner's theorem] show to be positive > definite: > > jj <- matrix(0,100,100) > A <- exp(-0.1*(row(jj)-col(jj))^2) > > > A's being positive-definite is important to me: > > >> min(eigen(A,T,T)$values) > [1] 2.521153e-10 >> > > Coercing A to a complex matrix should make no difference, but makes > eigen() return the wrong answer: > >> min(eigen(A+0i,T,T)$values) > [1] -0.359347 >> > > This is very, very wrong. Yep. I see this also on 10.6/7 (Snow Leopard, Lion) and 3.0.x, but NOT with a MacPorts build of 2.15.3 that I had lying around. So this sits somewhere between Mac builds, R versions, and possibly LAPACK issues. Can anyone reproduce on non-Mac? It's not the first time complex matrices have acted up. We may need to put in a regression test to catch it earlier. > > I would expect these two commands to return identical values, up to > numerical precision. Compare svd(): > > >> dput(min(eigen(A,T,T)$values)) > 2.52115250343783e-10 >> dput(min(eigen(A+0i,T,T)$values)) > -0.359346984206908 >> dput(min(svd(A)$d)) > 2.52115166468044e-10 >> dput(min(svd(A+0i)$d)) > 2.52115166468044e-10 >> > > So svd() doesn't care about the coercion to complex. The 'A' matrix > isn't particularly badly conditioned: > > >> eigen(A,T)$vectors -> e >> crossprod(e)[1:4,1:4] > > also: > >> crossprod(A,solve(A)) > > > [and the associated commands with A+0i in place of A], give errors of > order 1e-14 or less. > > > I think the eigenvectors are misbehaving too: > >> eigen(A,T)$vectors -> ev1 >> eigen(A+0i,T)$vectors -> ev2 >> range(Re((A %*% ev1[,100])/ev1[,100])) > [1] 2.497662e-10 2.566555e-10 # min=max mathematically; > differences due to numerics >> range(Re((A %*% ev2[,100])/ev2[,100])) > [1] -19.407290 4.412938 # off the scale errors > [note the difference in sign] >> > > > FWIW, these problems do not appear to occur if symmetric=FALSE: > >> min(Re(eigen(A+0i,F,T)$values)) > [1] 2.521153e-10 >> min(Re(eigen(A,F,T)$values)) > [1] 2.521153e-10 >> > > and the eigenvectors appear to behave themselves too. > > > Also, can I raise a doco? The documentation for eigen() is not > entirely transparent with regard to the 'symmetric' argument. For > complex matrices, 'symmetric' should read 'Hermitian': > > >> B <- matrix(c(2,1i,-1i,2),2,2) # 'B' is Hermitian >> eigen(B,F,T)$values > [1] 3+0i 1+0i >> eigen(B,T,T)$values# answers agree as expected if 'symmetric' means > 'Hermitian' > [1] 3 1 > > > >> C <- matrix(c(2,1i,1i,2),2,2)# 'C' is symmetric >> eigen(C,F,T)$values > [1] 2-1i 2+1i >> eigen(C,T,T)$values # answers disagree because 'C' is not Hermitian > [1] 3 1 >> > This issue, however, I do not understand; ?eigen already has: symmetric: if ‘TRUE’, the matrix is assumed to be symmetric (or Hermitian if complex) and only its lower triangle (diagonal included) is used. If ‘symmetric’ is not specified, the matrix is inspected for symmetry. Which part of "Hermitian if complex" is unclear? -- Peter Dalgaard, Professor, Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd@cbs.dk Priv: pda...@gmail.com __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] eigen(symmetric=TRUE) for complex matrices
On some further digging, on the Lion machine, the problem seems absent from builds against the veclib framework. I strongly suspect that the issue is the same as PR#14964, which also affects Hermitian matrices of dimension 33x33 and above. At smaller dimensions, a curious effect is visible: > jj <- matrix(0,33,33) > A <- exp(-0.1*(row(jj)-col(jj))^2) > eigen(A,,T) $values [1] 5.497434e+00 5.187095e+00 4.708628e+00 4.112589e+00 3.456602e+00 [6] 2.796253e+00 2.177674e+00 1.633085e+00 1.179630e+00 8.209787e-01 [11] 5.506744e-01 3.560858e-01 2.220308e-01 1.335197e-01 7.744384e-02 [16] 4.332359e-02 2.337139e-02 1.215408e-02 6.089878e-03 2.937761e-03 [21] 1.363015e-03 6.073985e-04 2.595197e-04 1.060717e-04 4.134966e-05 [26] 1.531412e-05 5.360284e-06 1.760457e-06 5.369131e-07 1.496798e-07 [31] 3.715172e-08 7.807420e-09 1.238432e-09 $vectors NULL > eigen(A+0i,,T) $values [1] 5.497434e+00 5.187095e+00 4.708628e+00 4.112589e+00 3.456602e+00 [6] 2.796253e+00 2.177674e+00 1.633085e+00 1.179630e+00 1.00e+00 [11] 8.209787e-01 5.506744e-01 3.560858e-01 2.220308e-01 1.335197e-01 [16] 7.744384e-02 4.332359e-02 2.337139e-02 1.215408e-02 6.089878e-03 [21] 2.937761e-03 1.363015e-03 6.073985e-04 2.595197e-04 1.060717e-04 [26] 4.134966e-05 1.531412e-05 5.360284e-06 1.760457e-06 5.369131e-07 [31] 1.496793e-07 3.707403e-08 6.664029e-09 $vectors NULL Notice that a bogus eigenvalue of 1.000 has sneaked into the 10th position in the complex case. This seems to be happening with increasing frequency as the dimension is increased, e.g. > jj <- matrix(0,39,39) > A <- exp(-0.1*(row(jj)-col(jj))^2) > eigen(A+0i,,T) $values [1] 5.525963e+00 5.295552e+00 4.932840e+00 4.466692e+00 4.313221e+00 [6] 3.931940e+00 3.365090e+00 2.800272e+00 2.266053e+00 1.992364e+00 [11] 1.783467e+00 1.365369e+00 1.016934e+00 7.369943e-01 5.730933e-01 [16] 5.197947e-01 3.568288e-01 2.384533e-01 1.551331e-01 1.070839e-01 [21] 9.826254e-02 6.059788e-02 3.638230e-02 2.126343e-02 1.318810e-02 [26] 1.209480e-02 6.693571e-03 3.602773e-03 1.884985e-03 1.011671e-03 [31] 9.580352e-04 4.725998e-04 2.260430e-04 1.046915e-04 4.687560e-05 [36] 3.794374e-05 2.020133e-05 7.751127e-06 1.472817e-06 $vectors NULL > eigen(A,,T) $values [1] 5.525963e+00 5.295552e+00 4.932840e+00 4.466692e+00 3.931940e+00 [6] 3.365090e+00 2.800272e+00 2.266053e+00 1.783467e+00 1.365369e+00 [11] 1.016934e+00 7.369943e-01 5.197947e-01 3.568288e-01 2.384533e-01 [16] 1.551331e-01 9.826254e-02 6.059788e-02 3.638230e-02 2.126343e-02 [21] 1.209480e-02 6.693571e-03 3.602773e-03 1.884985e-03 9.580352e-04 [26] 4.725998e-04 2.260430e-04 1.046915e-04 4.687623e-05 2.025048e-05 [31] 8.418654e-06 3.356873e-06 1.278239e-06 4.620504e-07 1.572170e-07 [36] 4.971980e-08 1.431462e-08 3.613385e-09 7.510902e-10 $vectors NULL in which most of the rightmost column of the complex case appear to be insertions. -pd On Jun 18, 2013, at 09:57 , peter dalgaard wrote: > > On Jun 18, 2013, at 03:30 , robin hankin wrote: > >> R-3.0.1 rev 62743, binary downloaded from CRAN just now; macosx 10.8.3 >> >> Hello, >> >> eigen(symmetric=TRUE) behaves strangely when given complex matrices. >> >> >> The following two lines define 'A', a 100x100 (real) symmetric matrix >> which theoretical considerations [Bochner's theorem] show to be positive >> definite: >> >> jj <- matrix(0,100,100) >> A <- exp(-0.1*(row(jj)-col(jj))^2) >> >> >> A's being positive-definite is important to me: >> >> >>> min(eigen(A,T,T)$values) >> [1] 2.521153e-10 >>> >> >> Coercing A to a complex matrix should make no difference, but makes >> eigen() return the wrong answer: >> >>> min(eigen(A+0i,T,T)$values) >> [1] -0.359347 >>> >> >> This is very, very wrong. > > Yep. I see this also on 10.6/7 (Snow Leopard, Lion) and 3.0.x, but NOT with > a MacPorts build of 2.15.3 that I had lying around. > > So this sits somewhere between Mac builds, R versions, and possibly LAPACK > issues. Can anyone reproduce on non-Mac? > > It's not the first time complex matrices have acted up. We may need to put in > a regression test to catch it earlier. > >> >> I would expect these two commands to return identical values, up to >> numerical precision. Compare svd(): >> >> >>> dput(min(eigen(A,T,T)$values)) >> 2.52115250343783e-10 >>> dput(min(eigen(A+0i,T,T)$values)) >> -0.359346984206908 >>> dput(min(svd(A)$d)) >> 2.52115166468044e-10 >>> dput(min(svd(A+0i)$d)) >> 2.52115166468044e-10 >>> >> >> So svd() doesn't care about the coercion to complex. The 'A' matrix >> isn't particularly badly conditioned: >> >> >>> eigen(A,T)$vectors -> e >>> crossprod(e)[1:4,1:4] >> >> also: >> >>> crossprod(A,solve(A)) >> >> >> [and the associated commands with A+0i in place of A], give errors of >> order 1e-14 or less. >> >> >> I think the eigenvectors are misbehaving too: >> >>> eigen(A,T)$vectors -> ev1 >>> eigen(A+0i,T)$vectors -> ev2 >>> range(Re((A %*% ev1[,100])/ev1[,100])) >> [1] 2.497662e-10 2.566555e-10 # min
Re: [Rd] eigen(symmetric=TRUE) for complex matrices
On 13-06-18 12:57 AM, peter dalgaard wrote: Coercing A to a complex matrix should make no difference, but makes >eigen() return the wrong answer: > >>min(eigen(A+0i,T,T)$values) >[1] -0.359347 >> > >This is very, very wrong. Yep. I see this also on 10.6/7 (Snow Leopard, Lion) and 3.0.x, but NOT with a MacPorts build of 2.15.3 that I had lying around. So this sits somewhere between Mac builds, R versions, and possibly LAPACK issues. Can anyone reproduce on non-Mac? It doesn't happen with the CRAN binary of 2.15.3 built for Ubuntu Precise. Davor __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] eigen(symmetric=TRUE) for complex matrices
On 18-06-2013, at 09:57, peter dalgaard wrote: > > On Jun 18, 2013, at 03:30 , robin hankin wrote: > >> R-3.0.1 rev 62743, binary downloaded from CRAN just now; macosx 10.8.3 >> >> Hello, >> >> eigen(symmetric=TRUE) behaves strangely when given complex matrices. >> >> >> The following two lines define 'A', a 100x100 (real) symmetric matrix >> which theoretical considerations [Bochner's theorem] show to be positive >> definite: >> >> jj <- matrix(0,100,100) >> A <- exp(-0.1*(row(jj)-col(jj))^2) >> >> >> A's being positive-definite is important to me: >> >> >>> min(eigen(A,T,T)$values) >> [1] 2.521153e-10 >>> >> >> Coercing A to a complex matrix should make no difference, but makes >> eigen() return the wrong answer: >> >>> min(eigen(A+0i,T,T)$values) >> [1] -0.359347 >>> >> >> This is very, very wrong. > > Yep. I see this also on 10.6/7 (Snow Leopard, Lion) and 3.0.x, but NOT with > a MacPorts build of 2.15.3 that I had lying around. > > So this sits somewhere between Mac builds, R versions, and possibly LAPACK > issues. Can anyone reproduce on non-Mac? > The problem does not occur with the Cran binary of R-3.0.1 Kubuntu 12.04 64-bit. That R uses the system provided Blas (libblas 1.2.30110419) and Lapack 3.3.1; I don't know if these have been patched. I have been able to reproduce the problem on a self compiled version of R-3.0.1 using Rlapack and Rblas on Ubuntu 10.04 64-bit. Berend __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] eigen(symmetric=TRUE) for complex matrices
1. OpenSuse 12.1 R compiled against ACML 5.3.1 > sessionInfo() R version 3.0.1 RC (2013-05-08 r62723) Platform: x86_64-unknown-linux-gnu (64-bit) > A <- exp(-0.1*(row(jj)-col(jj))^2) > min(eigen(A,T,T)$values) [1] 2.521151e-10 > min(eigen(A+0i,T,T)$values) [1] 2.521154e-10 2. OpenSuse 12.3 R compiled against Intel MKL 10.2 R version 3.0.0 Patched (2013-04-23 r62650) Platform: x86_64-unknown-linux-gnu (64-bit) > jj <- matrix(0,100,100) > A <- exp(-0.1*(row(jj)-col(jj))^2) > min(eigen(A,T,T)$values) [1] 2.521153e-10 > min(eigen(A+0i,T,T)$values) [1] 2.521154e-10 3. Windows 7 64, R 3.0.0p (binary, built-in libraries) R version 3.0.0 Patched (2013-04-03 r62488) Platform: x86_64-w64-mingw32/x64 (64-bit) > jj <- matrix(0,100,100) > A <- exp(-0.1*(row(jj)-col(jj))^2) > min(eigen(A,T,T)$values) [1] 2.521153e-10 > min(eigen(A+0i,T,T)$values) [1] -0.359347 same behaviour with the 32 bit binaries. On Tue, Jun 18, 2013 at 9:57 AM, peter dalgaard wrote: > > > On Jun 18, 2013, at 03:30 , robin hankin wrote: > > > R-3.0.1 rev 62743, binary downloaded from CRAN just now; macosx 10.8.3 > > > > Hello, > > > > eigen(symmetric=TRUE) behaves strangely when given complex matrices. > > > > > > The following two lines define 'A', a 100x100 (real) symmetric matrix > > which theoretical considerations [Bochner's theorem] show to be positive > > definite: > > > > jj <- matrix(0,100,100) > > A <- exp(-0.1*(row(jj)-col(jj))^2) > > > > > > A's being positive-definite is important to me: > > > > > >> min(eigen(A,T,T)$values) > > [1] 2.521153e-10 > >> > > > > Coercing A to a complex matrix should make no difference, but makes > > eigen() return the wrong answer: > > > >> min(eigen(A+0i,T,T)$values) > > [1] -0.359347 > >> > > > > This is very, very wrong. > > Yep. I see this also on 10.6/7 (Snow Leopard, Lion) and 3.0.x, but NOT with > a MacPorts build of 2.15.3 that I had lying around. > > So this sits somewhere between Mac builds, R versions, and possibly LAPACK > issues. Can anyone reproduce on non-Mac? > > It's not the first time complex matrices have acted up. We may need to put in > a regression test to catch it earlier. > > > > > I would expect these two commands to return identical values, up to > > numerical precision. Compare svd(): > > > > > >> dput(min(eigen(A,T,T)$values)) > > 2.52115250343783e-10 > >> dput(min(eigen(A+0i,T,T)$values)) > > -0.359346984206908 > >> dput(min(svd(A)$d)) > > 2.52115166468044e-10 > >> dput(min(svd(A+0i)$d)) > > 2.52115166468044e-10 > >> > > > > So svd() doesn't care about the coercion to complex. The 'A' matrix > > isn't particularly badly conditioned: > > > > > >> eigen(A,T)$vectors -> e > >> crossprod(e)[1:4,1:4] > > > > also: > > > >> crossprod(A,solve(A)) > > > > > > [and the associated commands with A+0i in place of A], give errors of > > order 1e-14 or less. > > > > > > I think the eigenvectors are misbehaving too: > > > >> eigen(A,T)$vectors -> ev1 > >> eigen(A+0i,T)$vectors -> ev2 > >> range(Re((A %*% ev1[,100])/ev1[,100])) > > [1] 2.497662e-10 2.566555e-10 # min=max mathematically; > > differences due to numerics > >> range(Re((A %*% ev2[,100])/ev2[,100])) > > [1] -19.407290 4.412938 # off the scale errors > > [note the difference in sign] > >> > > > > > > FWIW, these problems do not appear to occur if symmetric=FALSE: > > > >> min(Re(eigen(A+0i,F,T)$values)) > > [1] 2.521153e-10 > >> min(Re(eigen(A,F,T)$values)) > > [1] 2.521153e-10 > >> > > > > and the eigenvectors appear to behave themselves too. > > > > > > Also, can I raise a doco? The documentation for eigen() is not > > entirely transparent with regard to the 'symmetric' argument. For > > complex matrices, 'symmetric' should read 'Hermitian': > > > > > >> B <- matrix(c(2,1i,-1i,2),2,2) # 'B' is Hermitian > >> eigen(B,F,T)$values > > [1] 3+0i 1+0i > >> eigen(B,T,T)$values# answers agree as expected if 'symmetric' means > > 'Hermitian' > > [1] 3 1 > > > > > > > >> C <- matrix(c(2,1i,1i,2),2,2)# 'C' is symmetric > >> eigen(C,F,T)$values > > [1] 2-1i 2+1i > >> eigen(C,T,T)$values # answers disagree because 'C' is not Hermitian > > [1] 3 1 > >> > > > > This issue, however, I do not understand; ?eigen already has: > > > symmetric: if ‘TRUE’, the matrix is assumed to be symmetric (or > Hermitian if complex) and only its lower triangle (diagonal > included) is used. If ‘symmetric’ is not specified, the > matrix is inspected for symmetry. > > Which part of "Hermitian if complex" is unclear? > > -- > Peter Dalgaard, Professor, > Center for Statistics, Copenhagen Business School > Solbjerg Plads 3, 2000 Frederiksberg, Denmark > Phone: (+45)38153501 > Email: pd@cbs.dk Priv: pda...@gmail.com > > __ > R-devel@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel --
Re: [Rd] eigen(symmetric=TRUE) for complex matrices
On 18-06-2013, at 13:23, peter dalgaard wrote:
> On some further digging, on the Lion machine, the problem seems absent from
> builds against the veclib framework. I strongly suspect that the issue is the
> same as PR#14964, which also affects Hermitian matrices of dimension 33x33
> and above.
>
With the CRAN distribution of R-3.0.1 on OS X 10.8.4 I encountered the
problem too.
I've done a small experiment to see if I could find a possible cause.
This version of R uses the libRblas and libRlapack.
For complex matrices I'm assuming that eigen uses the Lapack routine zheev and
that R first calls zheev to determine the optimal workspace.
I've made an alternative eigen that also uses Lapack's zheev but calls it with
minimal workspace to force zheev to use the unblocked algorithm. To determine
if the blocked algorithm that Lapack uses could be the cause of the problem.
This will work when eigen is called first so that the Lapack library is
dyn.loaded.
Code:
# quick and dirty
zeigen <- function(A,symmetric,only.values=FALSE) {
# SUBROUTINE ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
# * INFO )
#.. Scalar Arguments ..
# * CHARACTER JOBZ, UPLO
# * INTEGERINFO, LDA, LWORK, N
# * ..
# * .. Array Arguments ..
# * DOUBLE PRECISION RWORK( * ), W( * )
# * COMPLEX*16 A( LDA, * ), WORK( * )
n <- dim(A)[1]
# use minimal workspace
lwork <- as.integer(2*n-1)
# warning: "N" and "L" work on OSX but you may have to write an interface
routine since some systems
# give Fortran runtime errors when passing character arguments to Fortran
routines.
z <- .Fortran("zheev", "N", "L", n, A=A, n, w=numeric(n),
work=complex(lwork), lwork,
numeric(3*n-2), info=integer(1L))
list(values=z$w, vectors=z$A)
}
N <- 100
jj <- matrix(0,N,N)
A <- exp(-0.1*(row(jj)-col(jj))^2)
min(eigen(A,only.values=TRUE,symmetric=TRUE)$values)
## [1] 2.521153e-10
# Coercing A to a complex matrix should make no difference, but makes
# eigen() return the wrong answer:
min(eigen(A+0i,only.values=TRUE,symmetric=TRUE)$values)
## [1] -0.359347
min(zeigen(A+0i,only.values=TRUE,symmetric=TRUE)$values)
##[1] 2.521154e-10
So it seems that the blocked algorithm is the cause of the error and that using
the (possibly slow) unblocked algorithm gives the correct result.
Berend
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] eigen(symmetric=TRUE) for complex matrices
On Jun 18, 2013, at 21:23 , Berend Hasselman wrote: > > So it seems that the blocked algorithm is the cause of the error and that > using the (possibly slow) unblocked algorithm gives the correct result. Thanks Berend, The finger seems to be pointing to the internal libRlapack/libRblas. The apparent pattern is that things work when R is linked against external libraries and not when the internal ones are used. So it could be time to start looking for differences between R's lapack module and the original LAPACK code. -pd -- Peter Dalgaard, Professor, Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd@cbs.dk Priv: pda...@gmail.com __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
