From ravi.varadhan at jhu.edu  Tue Jun 18 16:16:11 2013
From: ravi.varadhan at jhu.edu (Ravi Varadhan)
Date: Tue, 18 Jun 2013 14:16:11 +0000
Subject: [Rd] eigen(symmetric=TRUE) for complex matrices
In-Reply-To: <267A30A5-1865-4D23-8B9A-F04F9BBED1EE@gmail.com>
References: <CAHHjBM4z8d58QNb1HuqT4-k7BKaTo_=zAeuFeDr4pZdjTP2-Tg@mail.gmail.com>
	<5CEBFA58-9D40-42F0-87B8-F006EC3A848B@gmail.com>
	<267A30A5-1865-4D23-8B9A-F04F9BBED1EE@gmail.com>
Message-ID: <2F9EA67EF9AE1C48A147CB41BE2E15C34621A1@DOM-MTW-MAIL2.win.ad.jhu.edu>

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



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