title: Efficiency of matrix rowSums() vs. colSums() in R vs Rcpp vs Armadillo
body: <h2>Background</h2>
  
  <p>Coming from R programming, I'm in the process of expanding to compiled code in the form of C/C++ with <strong>Rcpp</strong>. As a hands on exercise on the effect of loop interchange (and just C/C++ in general), I implemented equivalents to R's <code>rowSums()</code> and <code>colSums()</code> functions for matrices with <strong>Rcpp</strong> (I know these exist as Rcpp sugar and in Armadillo -- this was just an exercise).</p>
  
  <h2>Question</h2>
  
  <p>I have my C++ implementation of <code>rowSums()</code> and <code>colSums()</code> along with <a href="https://github.com/RcppCore/Rcpp/blob/4e50e2bc4405309d3519575ff6cf8e7c57729b3d/inst/include/Rcpp/sugar/functions/rowSums.h#L128" rel="noreferrer">Rcpp sugar</a> and <code>arma::sum()</code> versions in <a href="https://gist.github.com/mikmart/66bf16f0bd329ec468aade6bf81fee96" rel="noreferrer">this <code>matsums.cpp</code> file</a>. Mine are just simple loops like this:</p>
  
  <pre class="lang-cpp prettyprint-override"><code>NumericVector Cpp_colSums(const NumericMatrix&amp; x) {
    int nr = x.nrow(), nc = x.ncol();
    NumericVector ans(nc);
    for (int j = 0; j &lt; nc; j++) {
      double sum = 0.0;
      for (int i = 0; i &lt; nr; i++) {
        sum += x(i, j);
      }
      ans[j] = sum;
    }
    return ans;
  }
  
  NumericVector Cpp_rowSums(const NumericMatrix&amp; x) {
    int nr = x.nrow(), nc = x.ncol();
    NumericVector ans(nr);
    for (int j = 0; j &lt; nc; j++) {
      for (int i = 0; i &lt; nr; i++) {
        ans[i] += x(i, j);
      }
    }
    return ans;
  }
  </code></pre>
  
  <p>(<em>R matrices are stored column-major, so columns in the outer loop should be the more efficient approach. That's what I was testing originally.</em>)</p>
  
  <p>While running benchmarks on these, I ran into something I wasn't expecting: there was a clear performance difference between row sums and col sums (see benchmarks below):</p>
  
  <ol>
  <li>Using the builtin R functions, <code>colSums()</code> is about twice as fast as <code>rowSums()</code>.</li>
  <li>With my own Rcpp and the sugar version, this is reversed: it is <code>rowSums()</code> that is about twice as fast as <code>colSums()</code>.</li>
  <li>And finally, adding the Armadillo implementations, the operations are roughly equal (col sum maybe a bit faster, as I would have expected them to be in R, too).</li>
  </ol>
  
  <p>So my primary question is: <em>why is <code>Cpp_rowSums()</code> significantly faster than <code>Cpp_colSums()</code>?</em></p>
  
  <p>As a secondary interest, I'm also curious why the same difference is reversed in the R implementations. (I skimmed through <a href="https://github.com/wch/r-source/blob/f05c3e60ded95a498dcd6c0ec7e66257f92fc08b/src/main/array.c#L1863" rel="noreferrer">the C source</a>, but could not really make out the significant differences.) (And third, how come Armadillo gets equal performance?)</p>
  
  <h2>Benchmarks</h2>
  
  <p>I tested all 4 implementations of both functions on a <code>10,000 x 10,000</code> symmetric matrix:</p>
  
  <pre class="lang-r prettyprint-override"><code>Rcpp::sourceCpp("matsums.cpp")
  
  set.seed(92136)
  n &lt;- 1e4 # build n x n test matrix
  x &lt;- matrix(rnorm(n), 1, n)
  x &lt;- crossprod(x, x) # symmetric
  
  bench::mark(
    rowSums(x),
    colSums(x),
    Cpp_rowSums(x),
    Cpp_colSums(x),
    Sugar_rowSums(x),
    Sugar_colSums(x),
    Arma_rowSums(x),
    Arma_colSums(x)
  )[, 1:7]
  
  #&gt; # A tibble: 8 x 7
  #&gt;   expression            min     mean   median      max `itr/sec` mem_alloc
  #&gt;   &lt;chr&gt;            &lt;bch:tm&gt; &lt;bch:tm&gt; &lt;bch:tm&gt; &lt;bch:tm&gt;     &lt;dbl&gt; &lt;bch:byt&gt;
  #&gt; 1 rowSums(x)        192.2ms  207.9ms  194.6ms  236.9ms      4.81    78.2KB
  #&gt; 2 colSums(x)         93.4ms   97.2ms   96.5ms  101.3ms     10.3     78.2KB
  #&gt; 3 Cpp_rowSums(x)     73.5ms   76.3ms     76ms   80.4ms     13.1     80.7KB
  #&gt; 4 Cpp_colSums(x)    126.5ms  127.6ms  126.8ms  130.3ms      7.84    80.7KB
  #&gt; 5 Sugar_rowSums(x)   73.9ms   75.6ms   74.3ms   79.4ms     13.2     80.7KB
  #&gt; 6 Sugar_colSums(x)  124.2ms  125.8ms  125.6ms  127.9ms      7.95    80.7KB
  #&gt; 7 Arma_rowSums(x)    73.2ms   74.7ms   73.9ms   79.3ms     13.4     80.7KB
  #&gt; 8 Arma_colSums(x)    62.8ms   64.4ms   63.7ms   69.6ms     15.5     80.7KB
  </code></pre>
  
  <p>(Again, you can find the C++ source file <code>matsums.cpp</code> <a href="https://gist.github.com/mikmart/66bf16f0bd329ec468aade6bf81fee96" rel="noreferrer">here</a>.)</p>
  
  <p><strong>Platform:</strong></p>
  
  <pre><code>&gt; sessioninfo::platform_info()
   setting  value                       
   version  R version 3.5.1 (2018-07-02)
   os       Windows &gt;= 8 x64            
   system   x86_64, mingw32             
   ui       RStudio                     
   language (EN)                        
   collate  English_United States.1252  
   tz       Europe/Helsinki             
   date     2018-08-09  
  </code></pre>
  
  <h1>Update</h1>
  
  <p>Investigating further, I also wrote the same functions using R's traditional C interface: the source is <a href="https://gist.github.com/mikmart/fb191894f29af03f69d3a0705ef1ee48#file-matsums-c" rel="noreferrer">here</a>. I <a href="https://gist.github.com/mikmart/fb191894f29af03f69d3a0705ef1ee48#file-matsums-bat" rel="noreferrer">compiled the functions</a> with <code>R CMD SHLIB</code>, and tested again: the same phenomenon of row sums being faster than col sums persisted (<a href="https://gist.github.com/mikmart/fb191894f29af03f69d3a0705ef1ee48#file-matsums-r-L18-L23" rel="noreferrer">benchmarks</a>). I then also looked at <a href="https://gist.github.com/mikmart/fb191894f29af03f69d3a0705ef1ee48#file-matsums-asm" rel="noreferrer">the disassembly with <code>objdump</code></a>, but it seems to me (with my very limited understanding of asm) that the compiler doesn't really optimize the main loop bodies (<a href="https://gist.github.com/mikmart/fb191894f29af03f69d3a0705ef1ee48#file-matsums-asm-L50-L64" rel="noreferrer">rows</a>, <a href="https://gist.github.com/mikmart/fb191894f29af03f69d3a0705ef1ee48#file-matsums-asm-L127-L141" rel="noreferrer">cols</a>) any further from the C code, either?</p>