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

Tue Jun 18 13:23:02 CEST 2013

On some further digging, on the Lion machine, the problem seems absent from builds against the veclib framework. I strongly suspect that the issue is the same as PR#14964, which also affects Hermitian matrices of dimension 33x33 and above.

At smaller dimensions, a curious effect is visible:

> jj <- matrix(0,33,33)
> A <- exp(-0.1*(row(jj)-col(jj))^2)
> eigen(A,,T)
$values
 [1] 5.497434e+00 5.187095e+00 4.708628e+00 4.112589e+00 3.456602e+00
 [6] 2.796253e+00 2.177674e+00 1.633085e+00 1.179630e+00 8.209787e-01
[11] 5.506744e-01 3.560858e-01 2.220308e-01 1.335197e-01 7.744384e-02
[16] 4.332359e-02 2.337139e-02 1.215408e-02 6.089878e-03 2.937761e-03
[21] 1.363015e-03 6.073985e-04 2.595197e-04 1.060717e-04 4.134966e-05
[26] 1.531412e-05 5.360284e-06 1.760457e-06 5.369131e-07 1.496798e-07
[31] 3.715172e-08 7.807420e-09 1.238432e-09

$vectors
NULL

> eigen(A+0i,,T)
$values
 [1] 5.497434e+00 5.187095e+00 4.708628e+00 4.112589e+00 3.456602e+00
 [6] 2.796253e+00 2.177674e+00 1.633085e+00 1.179630e+00 1.000000e+00
[11] 8.209787e-01 5.506744e-01 3.560858e-01 2.220308e-01 1.335197e-01
[16] 7.744384e-02 4.332359e-02 2.337139e-02 1.215408e-02 6.089878e-03
[21] 2.937761e-03 1.363015e-03 6.073985e-04 2.595197e-04 1.060717e-04
[26] 4.134966e-05 1.531412e-05 5.360284e-06 1.760457e-06 5.369131e-07
[31] 1.496793e-07 3.707403e-08 6.664029e-09

$vectors
NULL

Notice that a bogus eigenvalue of 1.000 has sneaked into the 10th position in the complex case. This seems to be happening with increasing frequency as the dimension is increased, e.g.

> jj <- matrix(0,39,39)
> A <- exp(-0.1*(row(jj)-col(jj))^2)
> eigen(A+0i,,T)
$values
 [1] 5.525963e+00 5.295552e+00 4.932840e+00 4.466692e+00 4.313221e+00
 [6] 3.931940e+00 3.365090e+00 2.800272e+00 2.266053e+00 1.992364e+00
[11] 1.783467e+00 1.365369e+00 1.016934e+00 7.369943e-01 5.730933e-01
[16] 5.197947e-01 3.568288e-01 2.384533e-01 1.551331e-01 1.070839e-01
[21] 9.826254e-02 6.059788e-02 3.638230e-02 2.126343e-02 1.318810e-02
[26] 1.209480e-02 6.693571e-03 3.602773e-03 1.884985e-03 1.011671e-03
[31] 9.580352e-04 4.725998e-04 2.260430e-04 1.046915e-04 4.687560e-05
[36] 3.794374e-05 2.020133e-05 7.751127e-06 1.472817e-06

$vectors
NULL

> eigen(A,,T)
$values
 [1] 5.525963e+00 5.295552e+00 4.932840e+00 4.466692e+00 3.931940e+00
 [6] 3.365090e+00 2.800272e+00 2.266053e+00 1.783467e+00 1.365369e+00
[11] 1.016934e+00 7.369943e-01 5.197947e-01 3.568288e-01 2.384533e-01
[16] 1.551331e-01 9.826254e-02 6.059788e-02 3.638230e-02 2.126343e-02
[21] 1.209480e-02 6.693571e-03 3.602773e-03 1.884985e-03 9.580352e-04
[26] 4.725998e-04 2.260430e-04 1.046915e-04 4.687623e-05 2.025048e-05
[31] 8.418654e-06 3.356873e-06 1.278239e-06 4.620504e-07 1.572170e-07
[36] 4.971980e-08 1.431462e-08 3.613385e-09 7.510902e-10

$vectors
NULL

in which most of the rightmost column of the complex case appear to be insertions.

-pd

On Jun 18, 2013, at 09:57 , peter dalgaard 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
> 
> 
> 
> 
> 
> 
> 
> 

-- 
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