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.

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
>





-- 
Robin Hankin
Uncertainty Analyst
hankin.robin@gmail.com

	[[alternative HTML version deleted]]

