3
$\begingroup$

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 ?

$\endgroup$

1 Answer 1

3
$\begingroup$

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.

$\endgroup$
6
  • $\begingroup$ thanks: it really didn't occur to me that a transpose could make such a large difference. $\endgroup$
    – user603
    Commented Feb 8, 2011 at 8:46
  • $\begingroup$ @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. $\endgroup$
    – mpiktas
    Commented Feb 8, 2011 at 8:49
  • $\begingroup$ @mpiktas (+1) Nicely done. (oh, and good to know that you're on a Mac :-) $\endgroup$
    – chl
    Commented Feb 8, 2011 at 11:15
  • $\begingroup$ @chl, only because of hardware :) Mac OS X's woeful swap treatment is killing me. Better than Windows, but way worse than Linux. $\endgroup$
    – mpiktas
    Commented Feb 8, 2011 at 11:41
  • $\begingroup$ @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. $\endgroup$
    – mpiktas
    Commented Feb 8, 2011 at 12:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.