From ripley at stats.ox.ac.uk  Thu Jan 19 11:49:07 2006
From: ripley at stats.ox.ac.uk (ripley@stats.ox.ac.uk)
Date: Thu, 19 Jan 2006 11:49:07 +0100 (CET)
Subject: [Rd] function 'eigen' (PR#8503)
Message-ID: <20060119104907.2B96826BC5@slim.kubism.ku.dk>

> I am reporting the mis-behaviour of the function 'eigen' in 'base', for the
> following input matrix:

All you can report is that the results do not accord with your
expectations, and I am afraid the problem is the latter.  However, I am 
not going be as impolite as to call that your `mis-behaviour'.

I am going to assume (because you did not tell us) that you are using a 
pre-built version of R for MacOS X linked against veclib.  We do ask you 
not to report issues with obselete versions of R, and 2.2.1 is current, 
not 2.1.1.  I believe that you will see different printed results on 2.2.1 
because of the NEWS item

     o	The printing of complex numbers has changed, handling numbers
 	as a whole rather than in two parts.  So both real and
 	imaginary parts are shown to the same accuracy, with the
 	'digits' parameter referring to the accuracy of the larger
 	component, and both components are shown in fixed or
 	scientific notation (unless one is entirely zero when it is
 	always shown in fixed notation).

J. F. Wilkinson wrote a large tome on the numerical analysis of the 
eigenproblem.  An over-simple precis is that all one can expect on a 
real-world computer is a solution to a problem somewhat close to the 
original one.  Since asymmetric matrices do not in general have a solution 
to the eigenproblem in the real domain, the analysis has to be conducted 
in the complex domain (where there is a solution, possibly degenerate). 
That is what you are seeing here: in particular for a theoretical solution 
with repeated eigenvalues the computed solution will almost certainly have 
close but not identical ones.  In general software makes internal 
decisions about whether to report repeated eigenvalues or not depending 
how close the computed values are.

The different results on MacOS 10.4.3/4 would reflect different versions 
of the OS services, in this case veclib. It looks as if Apple have already 
corrected one of their minor bugs.

eigen() for an asymmetric real matrix reports complex eigenvalues and 
eigenvectors whenever the computed eigenvalues are actually complex (that 
is, not all imaginary parts are identically zero).  So what you are seeing 
as '3+0i' probably does not have a zero imaginary part, but an imaginary 
part small compared to the modulus given that you asked for the result to 
7 significant figures.  Please distinguish between internal and printed 
representations of numbers.  (One could argue that the test should have a 
tolerance, but it is hard to know what it should be and I suspect it would 
be lower than your machine is giving.)

The problems you are seeing (on an obselete version of R) seem to indicate 
that Apple's OS services are less accurate than the LAPACK used by R on 
Windows, and indeed other LAPACKs I tried.  (My guess is that internally 
veclib is using some single-precision calculations for speed.)  If so, you 
need to adjust your expectations further to allow for the limited accuracy 
of your particular OS.


On Wed, 18 Jan 2006 Pierre.Legendre at umontreal.ca wrote:

> Full_Name: Pierre Legendre
> Version: 2.1.1
> OS: Mac OSX 10.4.3
> Submission from: (NULL) (132.204.120.81)
>
>
> I am reporting the mis-behaviour of the function 'eigen' in 'base', for the
> following input matrix:
>
> A <- matrix(c(2,3,4,-1,3,1,1,-2,0),3,3)
> eigen(A)
>
> I obtain the following results, which are incorrect for eigenvalues and
> eigenvectors 2 and 3 (incorrect imaginary portions):
> $values
> [1] 3+0.000000e+00i 1+1.869518e-08i 1-1.869518e-08i
>
> $vectors
>              [,1]                       [,2]                       [,3]
> [1,] -0.5345225+0i 4.720910e-17+2.643897e-09i 4.720910e-17-2.643897e-09i
> [2,] -0.2672612+0i 7.071068e-01+0.000000e+00i 7.071068e-01+0.000000e+00i
> [3,] -0.8017837+0i 7.071068e-01-2.643897e-09i 7.071068e-01+2.643897e-09i
>
> The eigenvalues of matrix A are 3, 1, and 1. The eigenvectors are real (not
> complex). Eigenvectors 2 and 3 of matrix A are identical.

More accurately, the geometric multiplicity of the repeated eigenvalue is 
1, so the eigenvectors span a one-dimensional space.  (This is one of 
those degenerate cases I mentioned earlier.)

> With R version 2.2.1 on a PC, 'eigen' produces the following (correct) results:
> Results obtained on a PC (WIndows) using R version 2.2.1:
>
> eigen(A)
> $values
> [1] 3 1 1
>
> $vectors
>           [,1]          [,2]          [,3]
> [1,] -0.5345225 -3.563430e-09 -3.563431e-09
> [2,] -0.2672612 -7.071068e-01  7.071068e-01
> [3,] -0.8017837 -7.071068e-01  7.071068e-01
>
> I obtained the  following results using R version 2.2.1 on a PowerMac G5 under
> Mac OSX 10.4.4:
>
> eigen(A)
> $values
> [1] 3+0i 1+0i 1-0i   # OK results, although there is no need to show the
> imaginary portion which is 0
>
> $vectors
>              [,1]                       [,2]                       [,3]
> [1,] -0.5345225+0i 4.720910e-17+2.643897e-09i 4.720910e-17-2.643897e-09i
> [2,] -0.2672612+0i 7.071068e-01+0.000000e+00i 7.071068e-01+0.000000e+00i
> [3,] -0.8017837+0i 7.071068e-01-2.643897e-09i 7.071068e-01+2.643897e-09i
>
> Incorrect imaginary portions, again, in eigenvectors 2 and 3.
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>
>

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595



[Rd] function 'eigen' (PR#8503)
ripley@stats.ox.ac.uk