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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
> 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 matterThis might be related to the problem (bad values from eigen) I reported a couple of weeks ago. I think it was tracked down to a problem in the header files. As I understand it, the bug has been around for awhile, but may not always bite - unless you're using a snapshot where it might be fixed. (I've reverted to 0.63.3 as it didn't seem to cause me a problem there.) Paul Gilbert -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
On Tue, 20 Apr 1999, Douglas Bates wrote:> 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) >That's the problem. On R : Copyright 1999, The R Development Core Team Version 0.64.0 Patched (unreleased snapshot) (April 19, 1999)> all(eigen(x,only.values=TRUE)$values == eigen(x)$values)[1] TRUE As the NEWS says: BUG FIXES o eigen() now should work (again) in all cases. There is quite a long list of bug fixes in 0.64.0, so I recommend the R-release version. Brian -- Brian D. Ripley, ripley@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 272860 (secr) Oxford OX1 3TG, UK Fax: +44 1865 272595 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._