<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/performance-avoid-coercing-indices-to-doubles/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://www.jottr.org/2018/04/02/coercion-of-indices/"> JottR on R</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p><img src="https://i2.wp.com/www.jottr.org/post/1or1L.png?w=456&ssl=1" alt="" data-recalc-dims="1"></p>
<p><code>x[idxs + 1]</code> or <code>x[idxs + 1L]</code>? That is the question.</p>
<p>Assume that we have a vector $x$ of $n = 100,000$ random values, e.g.</p>
<pre><code>> n &lt;- 100000
> x &lt;- rnorm(n)
</code></pre>
<p>and that we wish to calculate the $n-1$ first-order differences $y=(y_1, y_2, …, y_{n-1})$ where $y_i=x_{i+1} – x_i$. In R, we can calculate this using the following vectorized form:</p>
<pre><code>> idxs &lt;- seq_len(n - 1)
> y &lt;- x[idxs + 1] - x[idxs]
</code></pre>
<p>We can certainly do better if we turn to native code, but is there a more efficient way to implement this using plain R code? It turns out there is (*). The following <strong>calculation is ~15-20% faster</strong>:</p>
<pre><code>> y &lt;- x[idxs + 1L] - x[idxs]
</code></pre>
<p>The reason for this is because the index calculation:</p>
<pre><code>idxs + 1
</code></pre>
<p>is <strong>inefficient due to a coercion of integers to doubles</strong>. We have that <code>idxs</code> is an integer vector but <code>idxs + 1</code> becomes a double vector because <code>1</code> is a double:</p>
<pre><code>> typeof(idxs)
[1] &quot;integer&quot;
> typeof(idxs + 1)
[1] &quot;double&quot;
> typeof(1)
[1] &quot;double&quot;
</code></pre>
<p>Note also that doubles (aka “numerics” in R) take up <strong>twice the amount of memory</strong>:</p>
<pre><code>> object.size(idxs)
400040 bytes
> object.size(idxs + 1)
800032 bytes
</code></pre>
<p>which is because integers are stored as 4 bytes and doubles as 8 bytes.</p>
<p>By using <code>1L</code> instead, we can avoid this coercion from integers to doubles:</p>
<pre><code>> typeof(idxs)
[1] &quot;integer&quot;
> typeof(idxs + 1L)
[1] &quot;integer&quot;
> typeof(1L)
[1] &quot;integer&quot;
</code></pre>
<p>and we save some, otherwise wasted, memory;</p>
<pre><code>> object.size(idxs + 1L)
400040 bytes
</code></pre>
<p><strong>Does it really matter for the overall performance?</strong> It should because <strong>less memory is allocated</strong> which always comes with some overhead. Possibly more importantly, by using objects that are smaller in memory, the more likely it is that elements can be found in the memory cache rather than in the RAM itself, i.e. the <strong>chance for <em>cache hits</em> increases</strong>. Accessing data in the cache is orders of magnitute faster than in RAM. Furthermore, we also <strong>avoid coercion/casting</strong> of doubles to integers when R adds one to each element, which may add some extra CPU overhead.</p>
<p>The performance gain is confirmed by running <strong><a rel="nofollow" href="https://cran.r-project.org/package=microbenchmark" target="_blank">microbenchmark</a></strong> on the two alternatives:</p>
<pre><code>> microbenchmark::microbenchmark(
+   y &lt;- x[idxs + 1 ] - x[idxs],
+   y &lt;- x[idxs + 1L] - x[idxs]
+ )
Unit: milliseconds
                        expr  min   lq mean median   uq  max neval cld
  y &lt;- x[idxs + 1] - x[idxs] 1.27 1.58 3.71   2.27 2.62 80.6   100   a
 y &lt;- x[idxs + 1L] - x[idxs] 1.04 1.25 2.38   1.34 2.20 76.5   100   a
</code></pre>
<p>From the median (which is the most informative here), we see that using <code>idxs + 1L</code> is ~15-20% faster than <code>idxs + 1</code> in this case (it depends on $n$ and the overall calculation performed).</p>
<p><strong>Is it worth it?</strong> Although it is “only” an absolute difference of ~1 ms, it adds up if we do these calculations a large number times, e.g. in a bootstrap algorithm. And if there are many places in the code that result in coercions from index calculations like these, that also adds up. Some may argue it’s not worth it, but at least now you know it does indeed improve the performance a bit if you specify index constants as integers, i.e. by appending an <code>L</code>.</p>
<p>To wrap it up, here is look at the cost of subsetting all of the $1,000,000$ elements in a vector using various types of integer and double index vectors:</p>
<pre><code>> n &lt;- 1000000
> x &lt;- rnorm(n)
> idxs &lt;- seq_len(n)          ## integer indices
> idxs_dbl &lt;- as.double(idxs) ## double indices

> microbenchmark::microbenchmark(unit = &quot;ms&quot;,
+   x[],
+   x[idxs],
+   x[idxs + 0L],
+   x[idxs_dbl],
+   x[idxs_dbl + 0],
+   x[idxs_dbl + 0L],
+   x[idxs + 0]
+ )
Unit: milliseconds
             expr    min     lq   mean median     uq    max neval  cld
              x[] 0.7056 0.7481 1.6563 0.7632 0.8351 74.682   100 a   
          x[idxs] 3.9647 4.0638 5.1735 4.2020 4.7311 78.038   100  b  
     x[idxs + 0L] 5.7553 5.8724 6.2694 6.0810 6.6447  7.845   100  bc 
      x[idxs_dbl] 6.6355 6.7799 7.9916 7.1305 7.6349 77.696   100   cd
  x[idxs_dbl + 0] 7.7081 7.9441 8.6044 8.3321 8.9432 12.171   100    d
 x[idxs_dbl + 0L] 8.0770 8.3050 8.8973 8.7669 9.1682 12.578   100    d
      x[idxs + 0] 7.9980 8.2586 8.8544 8.8924 9.2197 12.345   100    d
</code></pre>
<p>(I ordered the entries by their ‘median’ processing times.)</p>
<p>In all cases, we are extracting the complete vector of <code>x</code>. We see that</p>
<ol>
<li>subsetting using an integer vector is faster than using a double vector,</li>
<li><code>x[idxs + 0L]</code> is faster than <code>x[idxs + 0]</code> (as seen previously),</li>
<li><code>x[idxs + 0L]</code> is still faster than <code>x[idxs_dbl]</code> despite also involving an addition, and</li>
<li><code>x[]</code> is whoppingly fast (probably because it does not have to iterate over an index vector) and serves as a lower-bound reference for the best we can hope for.</li>
</ol>
<p>(*): There already exists a highly efficient implementation for calculating the first-order differences, namely <code>y &lt;- diff(x)</code>. But for the sake of the take-home message of this blog post, let’s ignore that.</p>
<p><strong>Bonus</strong>: Did you know that <code>sd(y) / sqrt(2)</code> is an estimator of the standard deviation of the above <code>x</code>:s (von Neumann et al., 1941)? It’s actually not too hard to derive this – give it a try by deriving the variance when <code>x</code> is independent, identically distributed Gaussian random variables. This property is useful in cases where we are interested in the noise level of <code>x</code> and <code>x</code> has a piecewise constant mean level which changes at a small number of locations, e.g. a DNA copy-number profile of a tumor. In such cases we cannot use <code>sd(x)</code>, because the estimate would be biased due to the different mean levels. Instead, by taking the first-order differences <code>y</code>, changes in mean levels of <code>x</code> become sporadic outliers in <code>y</code>. If we could trim off these outliers, <code>sd(y) / sqrt(2)</code> would be a good estimate of the standard deviation of <code>x</code> after subtracting the mean levels. Even better, by using a robust estimator, such as the median absolute deviation (MAD) – <code>mad(y) / sqrt(2)</code> – we do not have to worry about have to identify the outliers. Efficient implementations of <code>sd(diff(x)) / sqrt(2))</code> and <code>mad(diff(x)) / sqrt(2))</code> are <code>sdDiff(x)</code> and <code>madDiff(x)</code> of the <strong><a rel="nofollow" href="https://cran.r-project.org/package=matrixStats" target="_blank">matrixStats</a></strong> package.</p>
<h1>References</h1>
<p>J. von Neumann et al., The mean square successive difference. <em>Annals of Mathematical Statistics</em>, 1941, 12, 153-162.</p>
<h1>Session information</h1>
<details>
<pre><code>> sessionInfo()
R version 3.4.4 (2018-03-15)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.4 LTS

Matrix products: default
BLAS: /usr/lib/atlas-base/atlas/libblas.so.3.0
LAPACK: /usr/lib/atlas-base/atlas/liblapack.so.3.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
[1] compiler_3.4.4
</code></pre>
</details>
	<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://www.jottr.org/2018/04/02/coercion-of-indices/"> JottR on R</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 rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/OQ5no1YsY7U" height="1"></a></div>