Skip to content

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

4 messages · Dirk Eddelbuettel, Kevin Ushey, Romain Francois

Original
Kevin Ushey
# 
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20130117/6a3d2b6a/attachment.html>
Dirk Eddelbuettel
# 
Hi Kevin,

First off, congrats regarding RcppRool. Nice to see another package on
CRAN. I had some follow-up questions as to whether you planned to push this
further: the mean at n+1 is really just adding 1/n of the difference between
the newest and the oldest... but getting these things right is hard. Anyway,
I didn't really have time poke around so consider these questions not asked :)
On 17 January 2013 at 18:09, Kevin Ushey wrote:
| 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 never measured. I always stopped at the simple sapply/lapply
implementations I showed in the talks/workshops.

One problem is column-wise storage of vectors. For kicks, can you see what
happens when for rowMeans you do a transpose and then compute colMeans ?

| I plan to submit an expanded version of this to the Rcpp gallery so any advice
| is appreciated!

That;s the spirit -- let's grow http://gallery.Rcpp.org which is meant to be
open for contributions.

Cheers,  Dirk
Kevin Ushey
#
On Thu, Jan 17, 2013 at 7:25 PM, Dirk Eddelbuettel <edd at debian.org> wrote:

            

      
      
      
    

        
Thanks! I do plan on extending the package further, but I need to iron out
some of the bugs and makes the pseudo-API available within more accessible
first. I'm definitely open to new ideas / additions though.

            

      
      
      
    

        
It's still faster to just compute the rowMeans compared to column means on
the transpose, for different sizes of matrices... Interestingly though, the
Rcpp rowMeans actually runs faster (!) than base-R rowMeans for a matrix w/
1E4 rows, 1E2 columns, so the relative speed does depend on the size /
shape of the matrix. I'm always a little surprised to learn these things.

Regardless, an Rcpp 'apply' solution is almost always faster than the
corresponding R 'apply' solution.

            

      
      
      
    

        
I will have something nifty to submit tomorrow. :)

            

      
      
      
    

        
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20130117/eb491f5a/attachment-0001.html>

  
  


  


      

    

    
12 days later

    

      
      

        


  

        

      Romain Francois
    

    

      
      #
    

  


  

      
Le 18/01/13 03:09, Kevin Ushey a ?crit :

            

      
      
      
    

        
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