eigenvalue calculations

[Douglas Bates]

I should have remembered that there was a problem with eigen() in
0.64.0.  In the patched versions of R-release (available under
src/devel at the CRAN sites) that bug has been fixed.

In case anyone else is interested, I redid the determinant
calculations in 
 Version 0.64.0 Patched (unreleased snapshot) (April 19, 1999)
using the method from Stephan Steinhaus's script (det0), the method
based on calculating the eigenvalues alone (det1) and the method based
on the QR decomposition (det2).  The results from the three methods
are consistent except that the eigenvalue methods both produce a
non-trivial imaginary component.

I enclose the timing results below.  Notice that the determinant
becomes extremely large as the matrix size increases and that the
error introduced by the complex arithmetic grows with it.

R> det0 <- function(x) prod(eigen(x)$values)
R> det1 <- function(x) prod(eigen(x, only = TRUE)$values)
R> det2 <- function(x) prod(diag(qr(x)$qr))*(-1)^(ncol(x)-1)
R> x1 <- array(runif(100*100), c(100,100))
R> x2 <- array(runif(200*200), c(200,200))
R> x3 <- array(runif(300*300), c(300,300))
R> system.time(print(det0(x1)))
[1] -9.267e+24+912457728i
[1] 1.36 0.00 2.00 0.00 0.00
R> system.time(print(det1(x1)))
[1] -9.267e+24+1041334272i
[1] 0.32 0.00 0.00 0.00 0.00
R> system.time(print(det2(x1)))
[1] -9.267e+24
[1] 0.04 0.00 0.00 0.00 0.00
R> system.time(print(det0(x2)))
[1] -9.4162e+79-1.458e+64i
[1] 10.28  0.00 11.00  0.00  0.00
R> system.time(print(det1(x2)))
[1] -9.4162e+79-2.9602e+64i
[1] 2.7 0.0 3.0 0.0 0.0
R> system.time(print(det2(x2)))
[1] -9.4162e+79
[1] 0.36 0.00 1.00 0.00 0.00
R> system.time(print(det0(x3)))
[1] -5.6519e+145+2.9732e+130i
[1] 34.86  0.01 34.00  0.00  0.00
R> system.time(print(det1(x3)))
[1] -5.6519e+145+9.494e+129i
[1] 10.49  0.00 10.00  0.00  0.00
R> system.time(print(det2(x3)))
[1] -5.6519e+145
[1] 1.49 0.00 1.00 0.00 0.00
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._