[Rd] eigen(symmetric=TRUE) for complex matrices
Ravi Varadhan
ravi.varadhan at jhu.edu
Tue Jun 18 16:16:11 CEST 2013
I am also able to reproduce this problem on Windows:
library(eiginv)
n <- 33 # problem does not arise for n <= 32
evals <- sort(rnorm(n))
system.time(A <- eiginv(evals, symmetric=TRUE))
all.equal(evals, sort(eigen(A+0i,,TRUE)$val))
cbind(evals, sort(eigen(A+0i,,TRUE)$val))
Best,
Ravi
> version
_
platform i386-w64-mingw32
arch i386
os mingw32
system i386, mingw32
status
major 3
minor 0.1
year 2013
month 05
day 16
svn rev 62743
language R
version.string R version 3.0.1 (2013-05-16)
nickname Good Sport
>
-----Original Message-----
From: r-devel-bounces at r-project.org [mailto:r-devel-bounces at r-project.org] On Behalf Of peter dalgaard
Sent: Tuesday, June 18, 2013 7:23 AM
To: robin hankin
Cc: r-devel at r-project.org
Subject: Re: [Rd] eigen(symmetric=TRUE) for complex matrices
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
______________________________________________
R-devel at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
More information about the R-devel
mailing list