[Rcpp-devel] Optimising 2d convolution

[Douglas Bates]

On Fri, Jan 6, 2012 at 2:58 PM, Douglas Bates <bates at stat.wisc.edu> wrote:
> On Fri, Jan 6, 2012 at 2:02 PM, Hadley Wickham <hadley at rice.edu> wrote:
>>> What are you doing the timing on?  On a modest desktop (2.6 GHz Athlon
>>> X4) I get less than a second for this
>> ...
>>>> k <- matrix(runif(20* 20), ncol = 20)
>>>> system.time(convolve_2d(x, k))
>>>   user  system elapsed
>>>  0.864   0.000   0.862
>>
>> Ooops, I must have been using a old k - I now get
>>
>>> system.time(convolve_2d(x, k))
>>   user  system elapsed
>>  2.265   0.006   2.298
>>
>> I'm on a 2nd gen macbook air, so the processor isn't the fastest (2.13
>> GHz). Maybe my intuition is wrong but it seems like you should be able
>> to do it faster than that.
>>
>> But let me think - I'm doing 1000 * 1000 * 20 * 20 = 4e08 multiplies +
>> adds in ~2s.  So that's 5e-9 s = 5 ns per operation.  Looking at it
>> that way it does seem pretty fast already.
>
> You can make this about 40% faster by using block operations in Eigen
> through RcppEigen.  In this code fragment I did remember to use a
> typedef and a const modifier.  I allocated the result as an
> Rcpp::NumericMatrix then mapped it to avoid an extra allocation and
> copy but that is a minor speedup for matrices of this size.

I just saw a superfluous statement after I changed to allocating an
Rcpp::NumericMatrix.  The .setZero() is no longer needed because Rcpp
zeros the storage.
-------------- next part --------------

R version 2.14.1 (2011-12-22)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(inline)
> # 2d convolution -------------------------------------------------------------
> 
> convolve_2d <- cxxfunction(signature(sampleS = "numeric", kernelS =
+ "numeric"), plugin = "Rcpp", '
+    Rcpp::NumericMatrix sample(sampleS), kernel(kernelS);
+    int x_s = sample.nrow(), x_k = kernel.nrow();
+    int y_s = sample.ncol(), y_k = kernel.ncol();
+ 
+    Rcpp::NumericMatrix output(x_s + x_k - 1, y_s + y_k - 1);
+    for (int row = 0; row < x_s; row++) {
+      for (int col = 0; col < y_s; col++) {
+        for (int i = 0; i < x_k; i++) {
+          for (int j = 0; j < y_k; j++) {
+            output(row + i, col + j) += sample(row, col) * kernel(i, j);
+          }
+        }
+      }
+    }
+    return output;
+ ')
> 
> convolve_2d_Eigen <- cxxfunction(signature(sampleS = "numeric", kernelS =
+ "numeric"), plugin = "RcppEigen", '
+    typedef Eigen::Map<Eigen::MatrixXd>  MMat;
+ 
+    const MMat sample(as<MMat>(sampleS)), kernel(as<MMat>(kernelS));
+    int x_s = sample.rows(), x_k = kernel.rows();
+    int y_s = sample.cols(), y_k = kernel.cols();
+    int x_o = x_s + x_k - 1, y_o = y_s + y_k - 1;
+ 
+    Rcpp::NumericMatrix out(x_o, y_o);
+    MMat  output(out.begin(), x_o, y_o);
+ 
+    for (int row = 0; row < x_s; row++) {
+      for (int col = 0; col < y_s; col++) {
+        output.block(row, col, x_k, y_k) += sample(row, col) * kernel;
+      }
+    }
+    return out;
+ ')
Loading required package: lattice

Attaching package: ?Matrix?

The following object(s) are masked from ?package:base?:

    det

> 
> x <- diag(1000)
> set.seed(1)
> k <- matrix(runif(20* 20), ncol = 20)
> system.time(loop <- convolve_2d(x, k))
   user  system elapsed 
  0.840   0.004   0.847 
> system.time(block <- convolve_2d_Eigen(x, k))
   user  system elapsed 
  0.580   0.004   0.584 
> all.equal(loop, block)
[1] TRUE
> 
> library(microbenchmark)
> microbenchmark(loop  = convolve_2d(x, k),
+                block = convolve_2d_Eigen(x, k), times=30)
Unit: milliseconds
   expr      min       lq   median       uq      max
1 block 473.0885 476.9951 478.4295 481.5343 492.4735
2  loop 830.6481 833.8049 835.1170 841.9029 854.5755
> 
> proc.time()
   user  system elapsed 
 51.726   1.180  52.486