The following is at least as much out of intellectual curiosity as for practical reasons. On reviewing some code written by novices to R, I came across:
crossprod(x, y)[1,1]
I thought, "That isn't a very S way of saying that, I wonder what the penalty is for using 'crossprod'." To my surprise the penalty was substantially negative. Handily the client had S-PLUS as well -- there the sign of the penalty was as I had expected, but the order of magnitude was off.
Here are the timings of 1 million computations on vectors of length 1000. This is under Windows, R version 1.9.1 and S-PLUS 6.2 (on the same machine).
Command R S-PLUS sum(x * y) 28.61 97.6 crossprod(x, y)[1,1] 6.77 2256.2
Another example is when computing the sums of the columns of a matrix. For example:
set.seed(1) jjm <- matrix(rnorm(600), 5)
Timings for this under Windows 2000 with R version 2.0.1 (on an old chip running at about 0.7Ghz) for 100,000 computations are:
apply(jjm, 2, sum) 536.59 colSums(jjm) 18.26 rep(1,5) %*% jjm 15.41 crossprod(rep(1,5), jjm) 13.16
(These timings seem to be stable across R versions and on at least one Linux platform.)
Andy Liaw showed another example of 'crossprod' being fast a couple days ago on R-help.
Questions for those with a more global picture of the code:
* Is the speed advantage of 'crossprod' inherent, or is it because more care has been taken with its implementation than the other functions?
* Is 'crossprod' faster than 'sum(x * y)' because 'crossprod' is going to BLAS while 'sum' can't?
For a numeric matrix crossprod ends up calling level 3 BLAS; either dsyrk for the single argument case or dgemm for the two argument case. Especially in accelerated versions of the BLAS like Atlas or Goto's BLAS, those routines are hideously efficient and that's where you are seeing the big gain in speed.
By the way, you didn't mention if you had an accelerated BLAS installed with R. Do you?
______________________________________________ R-devel@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
