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....@cbs.dk Priv: pda...@gmail.com > > > > > > > > -- Peter Dalgaard, Professor Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd....@cbs.dk Priv: pda...@gmail.com ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
