<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/tip-optimize-your-rcpp-loops/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://privefl.github.io/blog/Tip-Optimize-your-Rcpp-loops/"> Florian Privé</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<section>
<p>In this post, I will show you how to optimize your <code>Rcpp</code> loops so that they are 2 to 3 times faster than a standard implementation.</p>
<div>
<h2>Context</h2>
<div>
<h3>Real data example</h3>
<p>For this post, I will use a <code>big.matrix</code> which represents genotypes for 15,283 individuals, corresponding to the number of mutations (0, 1 or 2) at 287,155 different loci. Here, I will use only the first 10,000 loci (columns).</p>
<p>What you need to know about the <code>big.matrix</code> format:</p>
<ul>
<li>you can easily and quickly access matrice-like objects stored on disk,</li>
<li>you can use different types of storage (I use type <code>char</code> to store each element on only 1 byte),</li>
<li>it is column-major ordered as standard <code>R</code> matrices,</li>
<li>you can access elements of a <code>big.matrix</code> using <code>X[i, j]</code> in <code>R</code>,</li>
<li>you can access elements of a <code>big.matrix</code> using <code>X[j][i]</code> in <code>Rcpp</code>,</li>
<li>you can get a <code>RcppEigen</code> or <code>RcppArmadillo</code> view of a <code>big.matrix</code> (see Appendix).</li>
<li>for more details, go to <a target="_blank" rel="nofollow" href="https://github.com/kaneplusplus/bigmemory">the GitHub repo</a>.</li>
</ul>
<p>Peek at the data:</p>
<div>
<pre><code><span>print</span>(<span>dim</span>(X))</code></pre>
</div>
<pre><code>## [1] 15283 10000</code></pre>
<div>
<pre><code><span>print</span>(X[<span>1</span>:<span>10</span>, <span>1</span>:<span>12</span>])</code></pre>
</div>
<pre><code>##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
##  [1,]    2    0    2    0    2    2    2    1    2     2     2     2
##  [2,]    2    0    1    2    1    1    1    1    2     1     2     2
##  [3,]    2    0    2    2    2    2    1    1    2     1     2     2
##  [4,]    2    2    0    2    0    0    0    2    2     2     0     2
##  [5,]    2    1    2    2    2    2    2    1    2     2     2     2
##  [6,]    2    1    2    1    2    2    1    1    2     2     2     2
##  [7,]    2    0    2    0    2    2    2    0    2     1     2     2
##  [8,]    2    1    1    2    1    1    1    1    2     1     2     2
##  [9,]    2    1    2    2    2    2    2    2    2     2     2     2
## [10,]    2    0    2    1    2    2    2    0    2     1     2     1</code></pre>
</div>
<div>
<h3>What I needed</h3>
<p>I needed a fast matrix-vector multiplication between a <code>big.matrix</code> and a vector. Moreover, I could not use any <code>RcppEigen</code> or <code>RcppArmadillo</code> multiplication because I needed some options of efficiently subsetting columns or rows in my matrix (see Appendix).</p>
<p>Writing this multiplication in <code>Rcpp</code> is no more than two loops:</p>
<div>
<pre><code><span>// [[Rcpp::depends(RcppEigen, bigmemory, BH)]]</span>
<span>#include &lt;RcppEigen.h></span>
<span>#include &lt;bigmemory/MatrixAccessor.hpp></span>

<span>using</span> <span>namespace</span> Rcpp;

<span>// [[Rcpp::export]]</span>
NumericVector prod1(XPtr&lt;BigMatrix> bMPtr, <span>const</span> NumericVector&amp; x) {
  
  MatrixAccessor&lt;<span>char</span>> macc(*bMPtr);

  <span>int</span> n = bMPtr->nrow();
  <span>int</span> m = bMPtr->ncol();

  NumericVector res(n);
  <span>int</span> i, j;
  
  <span>for</span> (j = <span>0</span>; j &lt; m; j++) {
    <span>for</span> (i = <span>0</span>; i &lt; n; i++) {
      res[i] += macc[j][i] * x[j];
    }
  }

  <span>return</span> res;
}</code></pre>
</div>
<p>One test:</p>
<div>
<pre><code>y &lt;-<span> </span><span>rnorm</span>(<span>ncol</span>(X))

<span>print</span>(<span>system.time</span>(
  test &lt;-<span> </span><span>prod1</span>(X@address, y)
))</code></pre>
</div>
<pre><code>##    user  system elapsed 
##   0.664   0.004   0.668</code></pre>
<p><strong>What comes next should be transposable to other applications and other types of data.</strong></p>
</div>
</div>
<div>
<h2>Unrolling optimization</h2>
<p>While searching for optimizing my multiplication, I came across <a target="_blank" rel="nofollow" href="http://stackoverflow.com/a/12289513/6103040">this Stack Overflow answer</a>.</p>
<p>Unrolling in action:</p>
<div>
<pre><code><span>// [[Rcpp::depends(RcppEigen, bigmemory, BH)]]</span>
<span>#include &lt;RcppEigen.h></span>
<span>#include &lt;bigmemory/MatrixAccessor.hpp></span>

<span>using</span> <span>namespace</span> Rcpp;

<span>// [[Rcpp::export]]</span>
NumericVector prod4(XPtr&lt;BigMatrix> bMPtr, <span>const</span> NumericVector&amp; x) {
  
  MatrixAccessor&lt;<span>char</span>> macc(*bMPtr);
  
  <span>int</span> n = bMPtr->nrow();
  <span>int</span> m = bMPtr->ncol();
  
  NumericVector res(n);
  <span>int</span> i, j;
  
  <span>for</span> (j = <span>0</span>; j &lt;= m - <span>4</span>; j += <span>4</span>) {
    <span>for</span> (i = <span>0</span>; i &lt; n; i++) { <span>// unrolling optimization</span>
      res[i] += (x[j] * macc[j][i] + x[j<span>+1</span>] * macc[j<span>+1</span>][i]) +
        (x[j<span>+2</span>] * macc[j<span>+2</span>][i] + x[j<span>+3</span>] * macc[j<span>+3</span>][i]);
    } <span>// The parentheses are somehow important. Try without.</span>
  }
  <span>for</span> (; j &lt; m; j++) {
    <span>for</span> (i = <span>0</span>; i &lt; n; i++) {
      res[i] += x[j] * macc[j][i];
    }
  }
  
  <span>return</span> res;
}</code></pre>
</div>
<div>
<pre><code><span>require</span>(microbenchmark)

<span>print</span>(<span>microbenchmark</span>(
  <span>PROD1 =</span> test1 &lt;-<span> </span><span>prod1</span>(X@address, y),
  <span>PROD4 =</span> test2 &lt;-<span> </span><span>prod4</span>(X@address, y),
  <span>times =</span> <span>5</span>
))</code></pre>
</div>
<pre><code>## Unit: milliseconds
##   expr      min       lq     mean   median       uq      max neval
##  PROD1 609.0916 612.6428 613.7418 613.3740 616.4907 617.1096     5
##  PROD4 262.2658 267.7352 267.0268 268.0026 268.0785 269.0521     5</code></pre>
<div>
<pre><code><span>print</span>(<span>all.equal</span>(test1, test2))</code></pre>
</div>
<pre><code>## [1] TRUE</code></pre>
<p>Nice! Let’s try more. Why not using 8 or 16 rather than 4?</p>
<div>
<pre><code>Rcpp::<span>sourceCpp</span>(<span>'https://privefl.github.io/blog/code/prods.cpp'</span>)

<span>print</span>(bench &lt;-<span> </span><span>microbenchmark</span>(
  <span>PROD1 =</span> <span>prod1</span>(X@address, y),
  <span>PROD2 =</span> <span>prod2</span>(X@address, y),
  <span>PROD4 =</span> <span>prod4</span>(X@address, y),
  <span>PROD8 =</span> <span>prod8</span>(X@address, y),
  <span>PROD16 =</span> <span>prod16</span>(X@address, y),
  <span>times =</span> <span>5</span>
))</code></pre>
</div>
<pre><code>## Unit: milliseconds
##    expr      min       lq     mean   median       uq      max neval
##   PROD1 620.9375 627.9209 640.6087 631.1818 659.4236 663.5798     5
##   PROD2 407.6275 418.1752 417.1746 418.4589 419.0665 422.5451     5
##   PROD4 267.1687 271.4726 283.1928 271.9553 279.6698 325.6979     5
##   PROD8 241.5542 242.9120 255.4974 246.5218 267.7683 278.7307     5
##  PROD16 212.4335 213.5228 217.4781 217.1801 221.5119 222.7423     5</code></pre>
<div>
<pre><code>time &lt;-<span> </span><span>summary</span>(bench)[, <span>&quot;median&quot;</span>]
step &lt;-<span> </span><span>2</span>^(<span>0</span>:<span>4</span>)
<span>plot</span>(step, time, <span>type =</span> <span>&quot;b&quot;</span>, <span>xaxt =</span> <span>&quot;n&quot;</span>, <span>yaxt =</span> <span>&quot;n&quot;</span>, 
     <span>xlab =</span> <span>&quot;size of each step&quot;</span>)
<span>axis</span>(<span>side =</span> <span>1</span>, <span>at =</span> step)
<span>axis</span>(<span>side =</span> <span>2</span>, <span>at =</span> <span>round</span>(time))</code></pre>
</div>
<p><img src="https://i1.wp.com/privefl.github.io/blog/knitr_files/post-optimize-rcpp-loops_files/figure-html/unnamed-chunk-7-1.png?w=456&ssl=1" data-recalc-dims="1"></p>
</div>
<div>
<h2>Conclusion</h2>
<p>We have seen that unrolling can dramatically improve performances on loops. Steps of size 8 or 16 are of relatively little extra gain compared to 2 or 4.</p>
<p>As pointed out in the SO answer, it can behave rather differently between systems. So, if it is for your personal use, use the maximum gain (try 32!), but as I want my function to be used by others in a package, I think it’s safer to choose a step of 4.</p>
</div>
<div>
<h2>Appendix</h2>
<p>You can do a <code>big.matrix</code>-vector multiplication easily with <code>RcppEigen</code> or <code>RcppArmadillo</code> (see <a target="_blank" rel="nofollow" href="https://privefl.github.io/blog/code/prods2.cpp">this code</a>) but it lacks of efficient subsetting option.</p>
<p>Indeed, you still can’t use subsetting in <code>Eigen</code>, but this will come as said in <a target="_blank" rel="nofollow" href="http://eigen.tuxfamily.org/bz/show_bug.cgi?id=329">this feature request</a>. For <code>Armadillo</code>, you can but it is rather slow:</p>
<div>
<pre><code>Rcpp::<span>sourceCpp</span>(<span>'https://privefl.github.io/blog/code/prods2.cpp'</span>)

n &lt;-<span> </span><span>nrow</span>(X)
ind &lt;-<span> </span><span>sort</span>(<span>sample</span>(n, <span>size =</span> n/<span>2</span>))

<span>print</span>(<span>microbenchmark</span>(
  <span>EIGEN =</span> test3 &lt;-<span> </span><span>prodEigen</span>(X@address, y),
  <span>ARMA =</span> test4 &lt;-<span> </span><span>prodArma</span>(X@address, y),
  <span>ARMA_SUB =</span> test5 &lt;-<span> </span><span>prodArmaSub</span>(X@address, y, ind -<span> </span><span>1</span>),
  <span>times =</span> <span>5</span>
))</code></pre>
</div>
<pre><code>## Unit: milliseconds
##      expr       min        lq      mean    median        uq       max neval
##     EIGEN  567.5607  570.1843  717.2433  572.9402  576.2028 1299.3285     5
##      ARMA 1242.3581 1263.8803 1329.1212 1264.7070 1284.5612 1590.0993     5
##  ARMA_SUB  455.1174  457.5862  466.3982  461.5883  465.9056  491.7935     5</code></pre>
<div>
<pre><code><span>print</span>(<span>all</span>(
  <span>all.equal</span>(test3, test), 
  <span>all.equal</span>(<span>as.numeric</span>(test4), test),
  <span>all.equal</span>(<span>as.numeric</span>(test5), test[ind])
))</code></pre>
</div>
<pre><code>## [1] TRUE</code></pre>
</div>
</section>
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://privefl.github.io/blog/Tip-Optimize-your-Rcpp-loops/"> Florian Privé</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/NmDl3nFTRCs"></a></div>