<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/package-bigstatsr-statistics-with-matrices-on-disk-user-2017/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://privefl.github.io/blog/package-bigstatsr-statistics-with-matrices-on-disk-user-2017/"> blog</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 talk about my package <strong>bigstatsr</strong>, which I’ve just presented in a lightning talk of 5 minutes at useR!2017. You can listen to me in action <a rel="nofollow" href="https://t.co/aYt0q8MeXJ" target="_blank">there</a>. I should have chosen a longer talk to explain more about this package, maybe next time. I will use this post to give you a more detailed version of the talk I gave in Brussels.</p>
<div>
<h2>Motivation behind bigstatsr</h2>
<p>I’m a PhD student in predictive human genetics. I’m basically trying to predict someone’s risk of disease based on their DNA mutations. These DNA mutations are in the form of large matrices so that I’m currently working with a matrix of 15K rows and 300K columns. This matrix would take approximately 32GB of RAM if stored as a standard R matrix.</p>
<p>When I began studying this dataset, I had only 8GB of RAM on my computer. I now have 64GB of RAM but it would take only copying this matrix once to make my computer begin swapping and therefore slowing down. I found a convenient solution by using the object <code>big.matrix</code> provided by the R package <strong>bigmemory</strong> <span>(Kane, Emerson, and Weston 2013)</span>. With this solution, you can access a matrix that is stored on disk almost as if it were a standard R matrix in memory.</p>
<p><img width="70%" src="https://cdn.rawgit.com/privefl/useR-2017/f89f5928/memory-solution.svg"></p>
<p>Yet, for these special matrices, some useful statistical functions were missing or not fast enough. So, I implemented these. It was a good experience about programming optimized and parallelized algorithms. I’m aware that there are other packages that come with <strong>bigmemory</strong>, such as <strong>biganalytics</strong> and <strong>bigalgebra</strong>, that already implement some of the features I will talk about in this post. I will discuss why I don’t use these packages. However, I use the work of other packages such as <strong>biglasso</strong> and <strong>RSpectra</strong>.</p>
</div>
<div>
<h2>Introduction to bigmemory</h2>
<div>
<pre><code><span># loading package bigstatsr (and bigmemory)</span>
<span>library</span>(bigstatsr)
<span># initializing some matrix on disk: wrapper to bigmemory::big.matrix()</span>
mat &lt;-<span> </span><span>FBM</span>(<span>backingroot =</span> <span>&quot;matrix-on-disk&quot;</span>, <span>descriptor =</span> <span>FALSE</span>)(<span>5e3</span>, <span>10e3</span>)</code></pre>
</div>
<pre><code>## Creating directory &quot;backingfiles&quot; which didn't exist..</code></pre>
<div>
<pre><code><span>dim</span>(mat)</code></pre>
</div>
<pre><code>## [1]  5000 10000</code></pre>
<div>
<pre><code>mat[<span>1</span>:<span>5</span>, <span>1</span>:<span>5</span>]</code></pre>
</div>
<pre><code>##      [,1] [,2] [,3] [,4] [,5]
## [1,]    0    0    0    0    0
## [2,]    0    0    0    0    0
## [3,]    0    0    0    0    0
## [4,]    0    0    0    0    0
## [5,]    0    0    0    0    0</code></pre>
<div>
<pre><code>mat[<span>1</span>, <span>1</span>] &lt;-<span> </span><span>2</span>
mat[<span>1</span>:<span>5</span>, <span>1</span>:<span>5</span>]</code></pre>
</div>
<pre><code>##      [,1] [,2] [,3] [,4] [,5]
## [1,]    2    0    0    0    0
## [2,]    0    0    0    0    0
## [3,]    0    0    0    0    0
## [4,]    0    0    0    0    0
## [5,]    0    0    0    0    0</code></pre>
<div>
<pre><code>mat[<span>2</span>:<span>4</span>] &lt;-<span> </span><span>3</span>
mat[<span>1</span>:<span>5</span>, <span>1</span>:<span>5</span>]</code></pre>
</div>
<pre><code>##      [,1] [,2] [,3] [,4] [,5]
## [1,]    2    0    0    0    0
## [2,]    3    0    0    0    0
## [3,]    3    0    0    0    0
## [4,]    3    0    0    0    0
## [5,]    0    0    0    0    0</code></pre>
<div>
<pre><code>mat[, <span>2</span>:<span>3</span>] &lt;-<span> </span><span>rnorm</span>(<span>2</span> *<span> </span><span>nrow</span>(mat))
mat[<span>1</span>:<span>5</span>, <span>1</span>:<span>5</span>]</code></pre>
</div>
<pre><code>##      [,1]       [,2]       [,3] [,4] [,5]
## [1,]    2 -0.1611891  3.1218981    0    0
## [2,]    3 -0.6026680 -1.3935111    0    0
## [3,]    3  0.4743703 -0.6566549    0    0
## [4,]    3  0.1514252  0.8247594    0    0
## [5,]    0 -1.5167186  0.7329260    0    0</code></pre>
<p>What we can see is that big matrices (<code>big.matrix</code> objects) can be accessed (read/write) almost as if they were standard R matrices, but you have to be cautious. For example, doing <code>mat[1, ]</code> isn’t recommended. Indeed, big matrices, as standard R matrices, are stored by column so that it is in fact a big vector with columns stored one after the other, contiguously. So, accessing the first row would access elements that are not stored contiguously in memory, which is slow. One should always access columns rather than rows.</p>
</div>
<div>
<h2>Apply an R function to a big matrix</h2>
<p>An easy strategy to apply an R function to a big matrix would be the split-apply-combine strategy <span>(Wickham 2011)</span>. For example, you could access only a block of columns at a time, apply a (vectorized) function to this block, and then combine the results of all blocks. This is implemented in function <code>big_apply()</code>.</p>
<p><img width="70%" src="https://cdn.rawgit.com/privefl/useR-2017/f89f5928/split-apply-combine.svg"></p>
<div>
<pre><code><span># Compute the sums of the first 1000 columns</span>
colsums_1 &lt;-<span> </span><span>colSums</span>(mat[, <span>1</span>:<span>1000</span>])
<span># Compute the sums of the second block of 1000 columns</span>
colsums_2 &lt;-<span> </span><span>colSums</span>(mat[, <span>1001</span>:<span>2000</span>])
<span># Combine the results</span>
colsums_1_2 &lt;-<span> </span><span>c</span>(colsums_1, colsums_2)
<span># Do this automatically with big_apply()</span>
colsums_all &lt;-<span> </span><span>big_apply</span>(mat, <span>a.FUN =</span> function(X, ind) <span>colSums</span>(X[, ind]), 
                         <span>a.combine =</span> <span>'c'</span>)</code></pre>
</div>
<p>When the split-apply-combine strategy can be used for a given function, you could use <code>big_apply()</code> to get the results, while accessing only small blocks of columns (or rows) at a time. Package <strong>biganalytics</strong>, by the creators of <strong>bigmemory</strong>, provides a way to apply an R function to margins of a <code>big.matrix</code>. Yet, for example, if the <code>big.matrix</code> has a lot of columns, it would be much slower to loop through all columns rather that applying a vectorized function to blocks of columns. You can find more examples of applications for <code>big_apply()</code> <a rel="nofollow" href="https://privefl.github.io/bigstatsr/reference/big_apply.html" target="_blank">there</a>.</p>
<div>
<pre><code>colsums_all2 &lt;-<span> </span>biganalytics::<span>apply</span>(mat, <span>2</span>, sum)
<span>all.equal</span>(colsums_all2, colsums_all)</code></pre>
</div>
<pre><code>## [1] TRUE</code></pre>
</div>
<div>
<h2>Use Rcpp with a big matrix</h2>
<p>Using Rcpp with a <code>big.matrix</code> is super easy. Let’s use the previous example, i.e. the computation of the colsums of a <code>big.matrix</code>. We will do it in 3 different ways.</p>
<div>
<h3>1. Using the bigmemory way</h3>
<div>
<pre><code><span>// [[Rcpp::depends(bigmemory, BH)]]</span>
<span>#include </span></code></pre></div></div></div></section>