<section>
<p>In this post, I compare different approaches to get first principal components of large matrices in R.</p>
<div>
<h2>Comparison</h2>
<div>
<pre><code><span>library</span>(bigstatsr)
<span>library</span>(tidyverse)</code></pre>
</div>
<div>
<h3>Data</h3>
<div>
<pre><code><span># Create two matrices, one with some structure, one without</span>
n &lt;-<span> </span><span>20e3</span>
seq_m &lt;-<span> </span><span>c</span>(<span>1e3</span>, <span>3e3</span>, <span>10e3</span>)
sizes &lt;-<span> </span><span>seq_along</span>(seq_m)
X &lt;-<span> </span>E &lt;-<span> </span><span>list</span>()
<span>for</span> (i <span>in</span> sizes) {
  m &lt;-<span> </span>seq_m[i]
  U &lt;-<span> </span><span>matrix</span>(<span>0</span>, n, <span>10</span>); U[] &lt;-<span> </span><span>rnorm</span>(<span>length</span>(U))
  V &lt;-<span> </span><span>matrix</span>(<span>0</span>, m, <span>10</span>); V[] &lt;-<span> </span><span>rnorm</span>(<span>length</span>(V))
  E[[i]] &lt;-<span> </span><span>matrix</span>(<span>rnorm</span>(n <span>*</span><span> </span>m), n, m)
  X[[i]] &lt;-<span> </span><span>tcrossprod</span>(U, V) <span>+</span><span> </span>E[[i]]
}</code></pre>
</div>
<p>I use matrices of different sizes. Some are structured with 10 hidden components, and some with only random data.</p>
</div>
<div>
<h3>Optimized math library</h3>
<p>I linked my R installation with OpenBLAS, an optimized parallel matrix library.</p>
<div>
<pre><code>(NCORES &lt;-<span> </span>RhpcBLASctl<span>::</span><span>get_num_cores</span>())</code></pre>
</div>
<pre><code>## [1] 6</code></pre>
<div>
<pre><code>RhpcBLASctl<span>::</span><span>blas_set_num_threads</span>(NCORES)</code></pre>
</div>
<pre><code>## detected function goto_set_num_threads</code></pre>
</div>
<div>
<h3>Compared methods</h3>
<div>
<pre><code>models &lt;-<span> </span><span>tribble</span>(
  <span>~</span>method,                <span>~</span>fun,                      <span>~</span>params,
  <span>&quot;bigstatsr - 1 core&quot;</span>,   bigstatsr<span>::</span>big_randomSVD,  <span>list</span>(<span>k =</span> <span>10</span>),
  <span>&quot;bigstatsr - 6 cores&quot;</span>,  bigstatsr<span>::</span>big_randomSVD,  <span>list</span>(<span>k =</span> <span>10</span>, <span>ncores =</span> NCORES),
  <span>&quot;Rspectra&quot;</span>,             RSpectra<span>::</span>svds,            <span>list</span>(<span>k =</span> <span>10</span>),
  <span>&quot;irlba&quot;</span>,                irlba<span>::</span>irlba,              <span>list</span>(<span>nv =</span> <span>10</span>, <span>nu =</span> <span>10</span>),
  <span>&quot;svd&quot;</span>,                  svd<span>::</span>propack.svd,          <span>list</span>(<span>neig =</span> <span>10</span>),
  <span>&quot;rsvd&quot;</span>,                 rsvd<span>::</span>rsvd,                <span>list</span>(<span>k =</span> <span>10</span>)
) <span>%>%</span>
<span> </span><span>mutate</span>(<span>size =</span> <span>list</span>(sizes), <span>structured =</span> <span>list</span>(<span>c</span>(<span>TRUE</span>, <span>FALSE</span>))) <span>%>%</span>
<span> </span><span>unnest</span>(size, <span>.drop =</span> <span>FALSE</span>) <span>%>%</span>
<span> </span><span>unnest</span>(structured, <span>.drop =</span> <span>FALSE</span>) <span>%>%</span>
<span> </span><span>mutate</span>(<span>user_time =</span> <span>NA</span>, <span>real_time =</span> <span>NA</span>, <span>pcs =</span> <span>list</span>(<span>NA</span>))</code></pre>
</div>
</div>
<div>
<h3>Computing</h3>
<div>
<pre><code><span># Filling this data frame with times and PC scores for each method and dataset</span>
<span>for</span> (i <span>in</span> <span>rows_along</span>(models)) {

  mat &lt;-<span> `</span><span>if</span><span>`</span>(models<span>$</span>structured[[i]], X, E)[[models<span>$</span>size[[i]]]]

  time &lt;-<span> </span><span>system.time</span>({
    <span>if</span> (<span>grepl</span>(<span>&quot;bigstatsr&quot;</span>, models<span>$</span>method[[i]])) mat &lt;-<span> </span><span>as_FBM</span>(mat)
    res &lt;-<span> </span><span>do.call</span>(models<span>$</span>fun[[i]], <span>args =</span> <span>c</span>(<span>list</span>(mat), models<span>$</span>params[[i]]))
  })

  models[[<span>&quot;user_time&quot;</span>]][[i]] &lt;-<span> </span>time[<span>1</span>]
  models[[<span>&quot;real_time&quot;</span>]][[i]] &lt;-<span> </span>time[<span>3</span>]
  models[[<span>&quot;pcs&quot;</span>]][[i]]  &lt;-<span> </span>res
}</code></pre>
</div>
<pre><code>## Warning in (function (X, neig = min(m, n), opts = list()) : Only 4 singular triplets converged
## within 50 iterations.</code></pre>
<pre><code>## Warning in (function (X, neig = min(m, n), opts = list()) : Only 5 singular triplets converged
## within 50 iterations.</code></pre>
<pre><code>## Warning in (function (X, neig = min(m, n), opts = list()) : Only 5 singular triplets converged
## within 50 iterations.</code></pre>
<div>
<pre><code>models &lt;-<span> </span><span>mutate</span>(models, <span>size =</span> seq_m[size])</code></pre>
</div>
</div>
<div>
<h3>Timings</h3>
<div>
<pre><code>models <span>%>%</span>
<span> </span><span>ggplot</span>(<span>aes</span>(size <span>/</span><span> </span><span>1000</span>, real_time, <span>color =</span> method)) <span>+</span>
<span> </span><span>theme_bigstatsr</span>() <span>+</span>
<span> </span><span>geom_point</span>(<span>cex =</span> <span>6</span>) <span>+</span>
<span> </span><span>geom_line</span>(<span>aes</span>(<span>linetype =</span> method), <span>lwd =</span> <span>2</span>) <span>+</span>
<span> </span><span>facet_grid</span>(structured <span>~</span><span> </span>., <span>scales =</span> <span>&quot;free&quot;</span>) <span>+</span>
<span> </span><span>theme</span>(<span>legend.position =</span> <span>c</span>(<span>0.25</span>, <span>0.87</span>),
        <span>legend.key.width =</span> <span>unit</span>(<span>6</span>, <span>&quot;line&quot;</span>)) <span>+</span>
<span> </span><span>labs</span>(<span>x =</span> <span>sprintf</span>(<span>&quot;ncol (x1000) (nrow = %d)&quot;</span>, n), <span>y =</span> <span>&quot;Time (in seconds)&quot;</span>,
       <span>color =</span> <span>&quot;Methods:&quot;</span>, <span>linetype =</span> <span>&quot;Methods:&quot;</span>)</code></pre>
</div>
<p><img width="70%" src="https://privefl.github.io/blog/knitr_files/post-bench-svd_files/figure-html/unnamed-chunk-6-1.svg"></p>
<div>
<pre><code>models <span>%>%</span>
<span> </span><span>filter</span>(size <span>==</span><span> </span><span>max</span>(seq_m)) <span>%>%</span>
<span> </span><span>select</span>(method, structured, user_time, real_time)</code></pre>
</div>
<pre><code>## # A tibble: 12 x 4
##    method              structured user_time real_time
##    </code></pre></div></div></section>