eigenvalue/eigenvector calculations

[Douglas Bates]

Some of you may have seen a message on s-news by Stefan Steinhaus
regarding his paper on "Comparison of mathematical programs for data
analysis".  He compares S-PLUS 4.5 with several other programs.  He
does not include R in the comparisons.  On p. 28 of his report he
gives the URL the Auckland site along with URL's for two other systems
but comments that "I didn't received (sic) any support for my work by
the software developers".  I'm not sure what that means.

Anyway, one of his comparisons for timing is for the calculation of the
determinant of a 1000x1000 matrix of uniform (0,1) random numbers.  It
is not the wisest computation because it is usually Infinity but we
may still want to look at some of the properties.  He used
 determinant <- function(x) prod(eigen(x)$values)
There are many problems with this calculation and I address some of
them in a note to him with a copy to S-news.

I was cross-checking my results in R and came up with the following
peculiar behaviour.  If you use only.values = TRUE in eigen, you get
noticeably different eigenvalues and I don't think it is just a matter
of the ordering of the results.  Their products are noticeably
different.

R : Copyright 1999, The R Development Core Team
Version 0.64.0  (April 8, 1999)

> x <- array( runif( 10000 ), c(100, 100))
> determinant <- function(x) prod(eigen(x)$values)
> det1 <- function(x) prod(eigen(x, only.values = TRUE)$values)
> det2 <- function(x) prod(diag(qr(x)$qr)) * (-1)^(ncol(x) - 1)
> system.time(print(determinant(x)))
[1] -27320631-8.050401e-10i
[1] 0.06 0.04 0.00 0.00 0.00
> system.time(print(det1(x)))
[1] 36072865748-5.155196e-06i
[1] 0.06 0.00 0.00 0.00 0.00
> system.time(print(det2(x)))
[1] -9.854994e+25
[1] 0.03 0.00 0.00 0.00 0.00
> all(eigen(x,only.values=TRUE)$values == eigen(x)$values)
[1] FALSE
> rel.diff <- function(x,y) Mod(x-y)/Mod(x)
> rel.diff(eigen(x,only.values=TRUE)$values, eigen(x)$values)
  [1] 1.395689e-07 1.793713e-04 1.793713e-04 7.604683e-06 7.604683e-06
  [6] 6.027408e-05 6.027408e-05 2.001094e+00 1.887877e+00 2.088785e+00
 [11] 1.889478e+00 1.263914e+00 4.822391e-04 4.822391e-04 1.033824e-13
 [16] 1.582459e-05 1.590463e+00 1.684538e+00 1.861574e+00 1.891079e+00
 [21] 4.260913e-01 1.762078e+00 1.796430e+00 1.483012e+00 5.403557e-01
 [26] 5.403557e-01 1.217432e+00 1.217432e+00 5.188953e-01 5.188953e-01
 [31] 1.621965e+00 1.942157e+00 1.815673e+00 1.625482e+00 1.914144e+00
 [36] 1.914177e+00 1.728506e+00 1.724235e+00 7.271222e-01 7.271222e-01
 [41] 9.928653e-01 1.971190e+00 1.931997e+00 1.920628e+00 1.816204e+00
 [46] 1.946860e+00 1.565290e+00 8.955887e-01 1.924923e+00 1.388633e+00
 [51] 1.960728e+00 1.822560e+00 1.059483e+00 1.585604e+00 1.105094e+00
 [56] 5.622539e-02 5.622539e-02 1.795988e+00 1.846025e+00 1.333497e-01
 [61] 1.418881e+00 1.418881e+00 9.874141e-02 2.524210e-01 1.615536e+00
 [66] 1.195681e+00 1.793995e+00 6.952189e-01 1.268999e-01 1.472427e+00
 [71] 1.446002e+00 1.627224e-01 1.870035e+00 1.888633e+00 1.145509e-01
 [76] 1.625031e-01 1.860097e+00 1.469557e-01 1.438975e+00 1.505856e+00
 [81] 1.955270e+00 1.260539e+00 1.440042e+00 1.575450e-01 1.469263e-01
 [86] 1.769835e+00 1.780567e+00 2.119072e-01 3.139617e-01 2.748911e-01
 [91] 3.652874e-01 2.611222e-01 1.746714e+00 1.345816e-01 1.560194e+00
 [96] 1.225182e+00 2.020903e-01 1.836618e+00 1.962674e+00 9.724950e-16

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._