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.ro...@gmail.com
[[alternative HTML version deleted]]
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
