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 <pda...@gmail.com> 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
--
___________________________________________________
Simone Giannerini
Dipartimento di Scienze Statistiche "Paolo Fortunati"
Universita' di Bologna
Via delle belle arti 41 - 40126 Bologna, ITALY
Tel: +39 051 2098262 Fax: +39 051 232153
http://www2.stat.unibo.it/giannerini/
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
