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.mes at cbs.dk  Priv: PDalgd at gmail.com