[R] Sparse matrix calculation

[Douglas Bates]

On 11/15/06, YONGWAN CHUN <chun.49 at osu.edu> wrote:
> I work on large matrices and found something interesting. For multiplication of matrices, the order has a huge influence on computing time when one of them is a sparse matrix. In the below example, M is a full matrix and A is a sparse matrix in a regular matrix class. A %*% M takes much more time than M %*% A; moreover, t(t(M) %*% t(A)) is much faster than A %*% M with same result. I would like to know how it is possible. Even though I do not have an exact reason, this fact may be helpful to others.


> > n <- 1000
> > M <- diag(n) - matrix(1/n,n,n)
> > A <- matrix(rnorm(n*n)>2,n,n)
> > system.time(M %*% A)
> [1] 0.10 0.03 0.12   NA   NA
> > system.time(A %*% M)
> [1] 3.47 0.03 3.50   NA   NA
> > system.time(t(t(M) %*% t(A)))
> [1] 0.23 0.00 0.23   NA   NA

It is difficult to give a general answer to a question like this
because the behavior could change with different compilers or levels
of optimization or, more importantly, with the use of an accelerated
BLAS (Basic Linear Algebra Subroutines) library.  From Tamas' response
that he was unable to reproduce this behavior it appears that you may
be using an accelerated BLAS that creates a matrix product by using
the elements of a column of the right operand to define a linear
combination of the columns of the left operand.  If an element of the
right operand is zero then an entire 'AXPY' (y <- a * x + y) operation
can be skipped so it is an advantage to have the sparse matrix on the
right.

operand.