[Rcpp-devel] Significant difference in performance when computing Ax

[Søren Højsgaard]

Dear all,
When forming the matrix-vector-product one can form the inner products between As rows and x (mvprod1 below) or form a linear combination of As columns (mvprod2 below). The difference in computing time is striking (at least for me) for large matrices and it really illustrates the "caching issue" (I assume):

> microbenchmark(mvprod1(A,x), mvprod2(A,x), A%*%x)
Unit: milliseconds
          expr               min            lq                  median         uq          max              neval
mvprod1(A, x) 35.783812 36.310956 36.600186 36.991113 38.922889   100
mvprod2(A, x)   4.408892    4.436882   4.475368   4.568668   5.477176   100
       A %*% x         8.514091    8.592230   8.734978   8.811951   9.738653   100

Just wanted to share this!
Cheers
Søren



#include <Rcpp.h>
using namespace Rcpp;

//[[Rcpp::export]]
NumericVector mvprod1 (NumericMatrix A, NumericVector x){
  int nrA=A.nrow(), ncA=A.ncol(), i, j;
  NumericVector y(nrA);
  for (i=0; i<nrA; ++i){
    for (j=0; j<ncA; ++j){
      y[i] += A(i,j)*x[j];
    }
  }
  return y;
}

//[[Rcpp::export]]
NumericVector mvprod2 (NumericMatrix A, NumericVector x){
  int nrA=A.nrow(), ncA=A.ncol(), i, j;
  NumericVector y(nrA);
  for (j=0; j<ncA; ++j){
    for (i=0; i<nrA; ++i){
      y[i] += A(i,j)*x[j];
    }
  }
  return y;
}



/*** R
d <- 5
A <- matrix(1.0*(1:d^2),nrow=d)
x <- 1.0*(1:d)

A%*%x
mvprod1(A,x)
mvprod2(A,x)

d <- 2000
A <- matrix(1.0*(1:d^2),nrow=d)
x <- 1.0*(1:d)

library(microbenchmark)
microbenchmark(mvprod1(A,x), mvprod2(A,x), A%*%x)
*/
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20140425/932ec89d/attachment.html>