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

[Simone Giannerini]

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 <pdalgd at 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.mes at cbs.dk  Priv: PDalgd at gmail.com
>
> ______________________________________________
> R-devel at 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/