From p.dalgaard at biostat.ku.dk  Tue Jan 10 16:18:27 2006
From: p.dalgaard at biostat.ku.dk (Peter Dalgaard)
Date: 10 Jan 2006 16:18:27 +0100
Subject: [Rd] eigen()
In-Reply-To: <Pine.LNX.4.61.0601101451540.4195@gannet.stats>
References: <3B7754B8-4982-49BE-B424-A0D9CFC6F056@soc.soton.ac.uk>
	<x2zmm49qb3.fsf@viggo.kubism.ku.dk> <43C3B9DF.1090803@cimr.cam.ac.uk>
	<x2mzi49nl7.fsf@viggo.kubism.ku.dk>
	<70356E1A-618F-4610-8E17-ABCB10485282@soc.soton.ac.uk>
	<Pine.LNX.4.61.0601101451540.4195@gannet.stats>
Message-ID: <x2irss9kn0.fsf@viggo.kubism.ku.dk>

Prof Brian Ripley <ripley at stats.ox.ac.uk> writes:

> I haven't seen most of this thread, but this is a classic case of
> passing integers instead of doubles.  And indeed
> 
>      else if(is.numeric(x)) {
>  	storage.mode(x) <- "double"
> 
> has been removed from eigen.R in R-devel in r36952.  So that's the
> culprit.

Thought as much...
 
> [BTW, x86 is not 64-bit, which is x86_64.  I'd take x86 to be ix86,
> what Intel calls ia32 (and AMD does call x86, given what 'i' stands
> for here).]

Actually, I have

> version
               _
platform       x86_64-unknown-linux-gnu
arch           x86_64
os             linux-gnu
system         x86_64, linux-gnu
status         Under development (unstable)
major          2
minor          3.0
year           2006
month          01
day            10
svn rev        37041
language       R
version.string Version 2.3.0 Under development (unstable) (2006-01-10 r37041)

Which is in fact Pentium 4 EMT (not to be confused with ia64 which is
Itanium, AFAIR).

I don't see how anyone could figure that out from the information
given, though...
   
        -p

> On Tue, 10 Jan 2006, Robin Hankin wrote:
> 
> >
> >
> > On 10 Jan 2006, at 14:14, Peter Dalgaard wrote:
> >>
> >
> >
> >>>> Strange and semi-random results on SuSE 9.3 as well:
> >>>>
> >>>>
> >>>>> eigen(matrix(1:100,10,10))$values
> >>>>  [1]  5.4e-311+ 0.0e+00i -2.5e-311+3.7e-311i -2.5e-311-3.7e-311i
> >>>>  [4]  2.5e-312+ 0.0e+00i -2.4e-312+ 0.0e+00i  3.2e-317+ 0.0e+00i
> >>>>  [7]   0.0e+00+ 0.0e+00i   0.0e+00+ 0.0e+00i   0.0e+00+ 0.0e+00i
> >>>> [10]   0.0e+00+ 0.0e+00i
> >>>>
> >>>
> >>> Mine is closer to Robin's, but not the same (EL4 x86).
> >>>
> >>>> eigen(matrix(1:100,10,10))$values
> >>>   [1]  5.208398e+02+0.000000e+00i -1.583980e+01+0.000000e+00i
> >>>   [3]  6.292457e-16+2.785369e-15i  6.292457e-16-2.785369e-15i
> >>>   [5] -1.055022e-15+0.000000e+00i  3.629676e-16+0.000000e+00i
> >>>   [7]  1.356222e-16+2.682405e-16i  1.356222e-16-2.682405e-16i
> >>>   [9]  1.029077e-16+0.000000e+00i -1.269181e-17+0.000000e+00i
> >>>>
> >>>
> >>> But surely, my matrix algebra is a bit rusty, I think this matrix is
> >>> solveable analytically? Most of the eigenvalues shown are almost
> >>> exactly zero, except the first two, actually, which is about 521
> >>> and -16 to the closest integer.
> >>>
> >>> I think the difference between mine and Robin's are rounding errors
> >>> (the matrix is simple enough I expect the solution to be simple
> >>> integers
> >>> or easily expressible analystical expressions, so 8 e-values being
> >>> zero
> >>> is fine). Peter's number seems to be all 10 e-values are zero or one
> >>> being a huge number! So Peter's is odd... and Peter's machine also
> >>> seems
> >>> to be of a different archtecture (64-bit machine)?
> >>>
> >>> HTL
> >>
> >> Notice that Robin got something completely different in _R-devel_
> >> which is where I did my check too.  In R 2.2.1 I get the expected two
> >> non-zero eigenvalues.
> >>
> >> I'm not sure whether (and how) you can work out the eigenvalues
> >> analytically, but since all columns are linear progressions, it is
> >> at least obvious that the matrix must have column rank two.
> >>
> >
> >
> >
> >
> > For everyone's entertainment, here's an example where the analytic
> > solution
> > is known.
> >
> > fact 1:  the first eigenvalue of a magic square is equal to its constant
> > fact 2: the sum of the other eigenvalues of a magic square is zero
> > fact 3: the constant of a magic square of order 10 is 505.
> >
> > R-2.2.0:
> >
> > > library(magic)
> > > round(Re(eigen(magic(10),F,T)$values))
> > [1]  505  170 -170 -105  105   -3    3    0    0    0
> > >
> >
> > answers as expected.
> >
> >
> > R-devel:
> >
> >
> >
> > > a <- structure(c(68, 66, 92, 90, 16, 14, 37, 38, 41, 43, 65, 67, 89,
> > 91, 13, 15, 40, 39, 44, 42, 96, 94, 20, 18, 24, 22, 45, 46, 69,
> > 71, 93, 95, 17, 19, 21, 23, 48, 47, 72, 70, 4, 2, 28, 26, 49,
> > 50, 76, 74, 97, 99, 1, 3, 25, 27, 52, 51, 73, 75, 100, 98, 32,
> > 30, 56, 54, 80, 78, 81, 82, 5, 7, 29, 31, 53, 55, 77, 79, 84,
> > 83, 8, 6, 60, 58, 64, 62, 88, 86, 9, 10, 33, 35, 57, 59, 61,
> > 63, 85, 87, 12, 11, 36, 34), .Dim = c(10, 10))
> >
> > [no magic package!  it fails R CMD check !]
> >
> > > round(Re(eigen(magic(10),F,T)$values))
> > [1] 7.544456e+165  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e
> > +00
> > [6]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e
> > +00
> > >
> >
> >
> > not as expected.
> >
> >
> >
> > --
> > Robin Hankin
> > Uncertainty Analyst
> > National Oceanography Centre, Southampton
> > European Way, Southampton SO14 3ZH, UK
> >  tel  023-8059-7743
> >
> > ______________________________________________
> > 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
> 

-- 
   O__  ---- Peter Dalgaard             ?ster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics     PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark          Ph:  (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)                  FAX: (+45) 35327907



[Rd] eigen()
Peter Dalgaard