[Rcpp-devel] Most efficient way to apply function over rows of a matrix with Rcpp?

[Romain Francois]

Le 18/01/13 03:09, Kevin Ushey a écrit :
> Hi guys,
>
> I'm trying to use Rcpp to implement an 'apply'-type function for
> matrices, whereby I apply a function to each row or column of a matrix.
> I'm testing it here with the 'mean' just to see how well we can do vs.
> the highly optimized rowMeans, colMeans functions, but of course any
> other function taking a vector and returning a scalar might fit here.
>
> I'm wondering if I can do better. Please see the gist in the link here;
> you can sourceCpp it in R.
>
> https://gist.github.com/4561281
>
> I try to limit the amount of copying as much as possible with
> NumericMatrix::Column and NumericMatrix::Row.
>
> Note that the Rcpp solution over columns is just as fast as colMeans,
> but over rows it's a fair bit slower relative to rowMeans. Is there any
> way I could improve this?
>
> I plan to submit an expanded version of this to the Rcpp gallery so any
> advice is appreciated!
>
> Thanks,
> -Kevin

The use of the NumericMatrix::Row has some cost, each time we access 
tmp[i], we must calculate the right offset :

inline int offset( int i, int j) const { return i + nrows * j ; }

You could implement it yourself taking into account the structure of the 
matrix:

// [[Rcpp::export]]
NumericVector row_means_2( NumericMatrix& X ) {

   int nRows = X.nrow(), nCols = X.ncol() ;
   NumericVector out = no_init(nRows);

   for( int j=0, k=0; j < nCols; j++ ) {
     for( int i=0; i < nRows; i++, k++ ) {
       out[i] += X[k] ;
     }
   }
   for( int i =0; i<nRows; i++) out[i] /= nCols ;

   return out;

}

With this I get:

 > benchmark(replications = 5, order = NULL, apply(x,
+     1, mean), rowMeans(x), row_means(x), row_means_2(x))
                test replications elapsed  relative user.self sys.self
1 apply(x, 1, mean)            5   0.266 24.181818     0.253    0.013
2       rowMeans(x)            5   0.012  1.090909     0.011    0.000
3      row_means(x)            5   0.041  3.727273     0.041    0.000
4    row_means_2(x)            5   0.011  1.000000     0.011    0.000
   user.child sys.child
1          0         0
2          0         0
3          0         0
4          0         0


Maybe there is something we can do at our end to improve this by making 
a more efficient iterator over a matrix row.

Romain

-- 
Romain Francois
Professional R Enthusiast
+33(0) 6 28 91 30 30

R Graph Gallery: http://gallery.r-enthusiasts.com

blog:            http://romainfrancois.blog.free.fr
|- http://bit.ly/RE6sYH : OOP with Rcpp modules
`- http://bit.ly/Thw7IK : Rcpp modules more flexible