# Function timings in the Matrix package

In order to test the Matrix package in R, I ran the following piece of code:

library(Matrix)

p <- 50; n <- 500
a0 <- sample(1:(n*p), n*p*0.1)
d1 <- (a0-1) %% n
d2 <- floor((a0-d1)/n)
d0 <- cbind(d1+1, d2+1)
x0 <- Matrix(0, n, p, forceCheck=TRUE)
x0[d0] <- rnorm(n*p*0.1, 10, 1)
x1 <- matrix(rnorm(n*p), n, p)

fx01 <- function(ll, x0) tcrossprod(x0, rnorm(ncol(x0)))
fx02 <- function(ll, x1) crossprod(t(x1), rnorm(ncol(x1)))

system.time(lapply(1:1000, fx01, x0=x0))
system.time(lapply(1:1000, fx02, x1=x1))


Is it normal that the first function takes 4-5 times more time to run than the second when the second does not use a sparse matrix ?

There is an error in definition of fx01, the matrices are not compatible the way you defined. If we define

fx01 <- function(ll, x0) tcrossprod(rnorm(ncol(x0)),x0)


Then on my Macbook Pro with 2.53 GHz Intel Core 2 Duo processor with 4GB RAM and R 2.11.1, I get

> system.time(lapply(1:1000, fx02, x1=x1))
user  system elapsed
0.459   0.092   0.546
> system.time(lapply(1:1000, fx02, x1=x0))
user  system elapsed
0.632   0.005   0.633
> system.time(lapply(1:1000, fx01, x0=x1))
user  system elapsed
0.078   0.001   0.079
> system.time(lapply(1:1000, fx01, x0=x0))
user  system elapsed
3.810   0.014   3.805


Now this might give a hint (I cropped the output to make it clearer):

> showMethods("crossprod")
x="dgCMatrix", y="numeric"
(inherited from: x="CsparseMatrix", y="numeric")

> showMethods("tcrossprod")
x="numeric", y="dgCMatrix"
(inherited from: x="numeric", y="Matrix")


So I suspect that in tcrossprod there is a conversion of dgCMatrix (class of x0) to Matrix, meaning we convert sparse matrix to usual matrix. There is no such conversion in case of crossproduct. This might explain the difference, since sometimes unnecessary conversions are costly.

Further inspection confirms this, some sparse matrix code is used:

> getMethod("crossprod",signature=c(x="CsparseMatrix", y="numeric"))
Method Definition:

function (x, y = NULL)
.Call(Csparse_dense_crossprod, x, y)
<environment: namespace:Matrix>

Signatures:
x               y
target  "CsparseMatrix" "numeric"
defined "CsparseMatrix" "numeric"


No sparse matrix code is used:

> getMethod("tcrossprod",signature=c(x="numeric", y="Matrix"))
Method Definition:

function (x, y = NULL)
{
dim(x) <- c(1L, length(x))
callGeneric(x, y)
}
<environment: namespace:Matrix>

Signatures:
x         y
target  "numeric" "Matrix"
defined "numeric" "Matrix"


Side note: This is a first time I managed to get code from S4 methods (thanks for the question which forced me to figure it out), so I still may be wrong.

• thanks: it really didn't occur to me that a transpose could make such a large difference. Commented Feb 8, 2011 at 8:46
• @user603, you might want to file a bug for Matrix package developers. It seems that this is a simple matter of forgetting to add appropriate method definitions. Commented Feb 8, 2011 at 8:49
• @mpiktas (+1) Nicely done. (oh, and good to know that you're on a Mac :-)
– chl
Commented Feb 8, 2011 at 11:15
• @chl, only because of hardware :) Mac OS X's woeful swap treatment is killing me. Better than Windows, but way worse than Linux. Commented Feb 8, 2011 at 11:41
• @chl, if your swap partition is larger than 5GB sure :) If you know that you will need more swap than your current swap partition size, you can mount additional space as swap. This needs manual intervention of some kind, but on other hand it is doable via R. Commented Feb 8, 2011 at 12:08