[Rd] eigen(symmetric=TRUE) for complex matrices

[peter dalgaard]

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