<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/parallelizing-linear-regression-or-using-multiple-sources/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://freakonometrics.hypotheses.org/53283"> R-english – Freakonometrics</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>My <a rel="nofollow" href="https://freakonometrics.hypotheses.org/53269" target="_blank">previous post</a> was explaining how mathematically it was possible to parallelize computation to estimate the parameters of a linear regression. More speficially, we have a matrix <img alt="\mathbf{X}" title="\mathbf{X}" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7BX%7D&bg=ffffff&fg=000&s=0" class="latex"> which is <img alt="n\times k" title="n\times k" src="https://s0.wp.com/latex.php?latex=n%5Ctimes+k&bg=ffffff&fg=000&s=0" class="latex"> matrix and <img alt="\mathbf{y}" title="\mathbf{y}" src="https://s0.wp.com/latex.php?latex=%5Cmathbf%7By%7D&bg=ffffff&fg=000&s=0" class="latex"> a <img alt="n" title="n" src="https://s0.wp.com/latex.php?latex=n&bg=ffffff&fg=000&s=0" class="latex">-dimensional vector, and we want to compute <img alt="\widehat{\mathbf{\beta}}=[\mathbf{X}^T\mathbf{X}]^{-1}\mathbf{X}^T\mathbf{y}" title="\widehat{\mathbf{\beta}}=[\mathbf{X}^T\mathbf{X}]^{-1}\mathbf{X}^T\mathbf{y}" src="https://s0.wp.com/latex.php?latex=%5Cwidehat%7B%5Cmathbf%7B%5Cbeta%7D%7D%3D%5B%5Cmathbf%7BX%7D%5ET%5Cmathbf%7BX%7D%5D%5E%7B-1%7D%5Cmathbf%7BX%7D%5ET%5Cmathbf%7By%7D&bg=ffffff&fg=000&s=0" class="latex"> by spliting the job. Instead of using the <img alt="n" title="n" src="https://s0.wp.com/latex.php?latex=n&bg=ffffff&fg=000&s=0" class="latex"> observations, we’ve seen that it was to possible to compute “something” using the first <img alt="n_1" title="n_1" src="https://s0.wp.com/latex.php?latex=n_1&bg=ffffff&fg=000&s=0" class="latex"> rows, then the next <img alt="n_2" title="n_2" src="https://s0.wp.com/latex.php?latex=n_2&bg=ffffff&fg=000&s=0" class="latex"> rows, etc. Then, finally, we “aggregate” the <img alt="m" title="m" src="https://s0.wp.com/latex.php?latex=m&bg=ffffff&fg=000&s=0" class="latex"> objects created to get our overall estimate.</p>
<h2>Parallelizing on multiple cores</h2>
<p>Let us see how it works from a computational point of view, to run each computation on a different core of the machine. Each core will see a <a rel="nofollow" href="https://en.wikipedia.org/wiki/Master/slave_(technology)" target="_blank">slave</a>, computing what we’ve seen in the <a rel="nofollow" href="https://freakonometrics.hypotheses.org/53269" target="_blank">previous post</a>. Here, the data we use are</p>
<div><table><tr><td><pre>1
2
3
</pre></td><td><pre>y <span>=</span> <span>cars</span>$dist
X <span>=</span> <span>data.<span>frame</span></span><span>(</span><span>1</span>,<span>cars</span>$speed<span>)</span>
k <span>=</span> <span>ncol</span><span>(</span>X<span>)</span></pre></td></tr></table></div>
<p>On my laptop, I have three cores, so we will split it in <img alt="m=3" title="m=3" src="https://s0.wp.com/latex.php?latex=m%3D3&bg=ffffff&fg=000&s=0" class="latex"> chunks</p>
<div><table><tr><td><pre>1
2
3
4
</pre></td><td><pre><span>library</span><span>(</span>parallel<span>)</span>
<span>library</span><span>(</span>pbapply<span>)</span>
ncl <span>=</span> detectCores<span>(</span><span>)</span><span>-</span><span>1</span>
cl <span>=</span> makeCluster<span>(</span>ncl<span>)</span></pre></td></tr></table></div>
<p>This is more or less what we will do: we have our dataset, and we split the jobs,</p>
<p><img src="https://i2.wp.com/freakonometrics.hypotheses.org/files/2018/06/Capture-d%E2%80%99e%CC%81cran-2018-06-21-a%CC%80-18.14.05.png?resize=348%2C193&ssl=1" alt="" width="348" data-recalc-dims="1" height="193"></p>
<p>We can then create lists containing elements that will be sent to each core, as <a rel="nofollow" href="https://twitter.com/3wen" target="_blank">Ewen</a> suggested,</p>
<div><table><tr><td><pre>1
2
3
4
5
</pre></td><td><pre>chunk <span>=</span> <span>function</span><span>(</span>x,n<span>)</span> <span>split</span><span>(</span>x, <span>cut</span><span>(</span><span>seq_along</span><span>(</span>x<span>)</span>, n, <span>labels</span> <span>=</span> FALSE<span>)</span><span>)</span>
a_parcourir <span>=</span> chunk<span>(</span><span>seq_len</span><span>(</span><span>nrow</span><span>(</span>X<span>)</span><span>)</span>, ncl<span>)</span>
<span>for</span><span>(</span>i <span>in</span> <span>1</span><span>:</span><span>length</span><span>(</span>a_parcourir<span>)</span><span>)</span> a_parcourir<span>[</span><span>[</span>i<span>]</span><span>]</span> <span>=</span> <span>rep</span><span>(</span>i, <span>length</span><span>(</span>a_parcourir<span>[</span><span>[</span>i<span>]</span><span>]</span><span>)</span><span>)</span>
Xlist <span>=</span> <span>split</span><span>(</span>X, <span>unlist</span><span>(</span>a_parcourir<span>)</span><span>)</span>
ylist <span>=</span> <span>split</span><span>(</span>y, <span>unlist</span><span>(</span>a_parcourir<span>)</span><span>)</span></pre></td></tr></table></div>
<p>It is also possible to simplify the QR functions we will use</p>
<div><table><tr><td><pre>1
2
3
4
5
6
7
8
</pre></td><td><pre>compute_qr <span>=</span> <span>function</span><span>(</span>x<span>)</span><span>{</span>
  <span>list</span><span>(</span>Q<span>=</span><span>qr.<span>Q</span></span><span>(</span><span>qr</span><span>(</span><span>as.<span>matrix</span></span><span>(</span>x<span>)</span><span>)</span><span>)</span>,R<span>=</span><span>qr.<span>R</span></span><span>(</span><span>qr</span><span>(</span><span>as.<span>matrix</span></span><span>(</span>x<span>)</span><span>)</span><span>)</span><span>)</span>
<span>}</span>
get_Vlist <span>=</span> <span>function</span><span>(</span>j<span>)</span><span>{</span>
  Q3 <span>=</span> QR1<span>[</span><span>[</span>j<span>]</span><span>]</span>$Q <span>%*%</span> Q2list<span>[</span><span>[</span>j<span>]</span><span>]</span>
  <span>t</span><span>(</span>Q3<span>)</span> <span>%*%</span> ylist<span>[</span><span>[</span>j<span>]</span><span>]</span>
<span>}</span>
clusterExport<span>(</span>cl, <span>c</span><span>(</span><span>&quot;compute_qr&quot;</span>, <span>&quot;get_Vlist&quot;</span><span>)</span>, envir<span>=</span><span>environment</span><span>(</span><span>)</span><span>)</span></pre></td></tr></table></div>
<p>Then, we can run our functions on each core. The first one is</p>
<div><table><tr><td><pre>1
</pre></td><td><pre>  QR1 <span>=</span> parLapply<span>(</span>cl<span>=</span>cl,Xlist, compute_qr<span>)</span></pre></td></tr></table></div>
<p>note that it is also possible to use</p>
<div><table><tr><td><pre>1
</pre></td><td><pre>  QR1 <span>=</span> pblapply<span>(</span>Xlist, compute_qr, cl<span>=</span>cl<span>)</span></pre></td></tr></table></div>
<p>which will include a progress bar (that can be nice when the database is rather large). Then use</p>
<div><table><tr><td><pre>1
2
3
4
5
6
7
</pre></td><td><pre>  R1 <span>=</span> pblapply<span>(</span>QR1, <span>function</span><span>(</span>x<span>)</span> x$R, cl<span>=</span>cl<span>)</span> <span>%&amp;</span>gt<span>;%</span> <span>do.<span>call</span></span><span>(</span><span>&quot;rbind&quot;</span>, .<span>)</span>
  Q1 <span>=</span> <span>qr.<span>Q</span></span><span>(</span><span>qr</span><span>(</span><span>as.<span>matrix</span></span><span>(</span>R1<span>)</span><span>)</span><span>)</span>
  R2 <span>=</span> <span>qr.<span>R</span></span><span>(</span><span>qr</span><span>(</span><span>as.<span>matrix</span></span><span>(</span>R1<span>)</span><span>)</span><span>)</span>
  Q2list <span>=</span> <span>split.<span>data</span>.<span>frame</span></span><span>(</span>Q1, <span>rep</span><span>(</span><span>1</span><span>:</span>ncl, each<span>=</span>k<span>)</span><span>)</span>
  clusterExport<span>(</span>cl, <span>c</span><span>(</span><span>&quot;QR1&quot;</span>, <span>&quot;Q2list&quot;</span>, <span>&quot;ylist&quot;</span><span>)</span>, envir<span>=</span><span>environment</span><span>(</span><span>)</span><span>)</span>
  Vlist <span>=</span> pblapply<span>(</span><span>1</span><span>:</span><span>length</span><span>(</span>QR1<span>)</span>, get_Vlist, cl<span>=</span>cl<span>)</span>
  sumV <span>=</span> <span>Reduce</span><span>(</span><span>'+'</span>, Vlist<span>)</span></pre></td></tr></table></div>
<p>and finally the ouput is</p>
<div><table><tr><td><pre>1
2
3
4
</pre></td><td><pre><span>solve</span><span>(</span>R2<span>)</span> <span>%*%</span> sumV
         <span>[</span>,<span>1</span><span>]</span>
X1 <span>-</span><span>17.579095</span>
X2   <span>3.932409</span></pre></td></tr></table></div>
<p>which is what we were expecting…</p>
<h2>Using multiple sources</h2>
<p>In practice, it might also happen that various “servers” have the data, but we cannot get a copy. But it is possible to run some functions on their server, and get some output, that we can use afterwards.</p>
<p><img src="https://i2.wp.com/freakonometrics.hypotheses.org/files/2018/06/Capture-d%E2%80%99e%CC%81cran-2018-06-21-a%CC%80-18.14.13.png?resize=348%2C193&ssl=1" alt="" width="348" data-recalc-dims="1" height="193"></p>
<p>Datasets are supposed to be available somewhere. We can send a request, and get a matrix. Then we we aggregate all of them, and send another request. That’s what we will do here. Provider <img alt="j" title="j" src="https://s0.wp.com/latex.php?latex=j&bg=ffffff&fg=000&s=0" class="latex"> should run <img alt="f_1(\mathbf{X})" title="f_1(\mathbf{X})" src="https://s0.wp.com/latex.php?latex=f_1%28%5Cmathbf%7BX%7D%29&bg=ffffff&fg=000&s=0" class="latex"> on his part of the data, that function will return <img alt="R^{(1)}_j" title="R^{(1)}_j" src="https://s0.wp.com/latex.php?latex=R%5E%7B%281%29%7D_j&bg=ffffff&fg=000&s=0" class="latex">. More precisely, to the first provider, send</p>
<div><table><tr><td><pre>1
2
3
</pre></td><td><pre>function1 <span>=</span> <span>function</span><span>(</span>subX<span>)</span><span>{</span>
<span>return</span><span>(</span><span>qr.<span>R</span></span><span>(</span><span>qr</span><span>(</span><span>as.<span>matrix</span></span><span>(</span>subX<span>)</span><span>)</span><span>)</span><span>)</span><span>}</span>
R1 <span>=</span> function1<span>(</span>Xlist<span>[</span><span>[</span><span>1</span><span>]</span><span>]</span><span>)</span></pre></td></tr></table></div>
<p>and actually, send that function to all providers, and aggregate the output</p>
<div><table><tr><td><pre>1
</pre></td><td><pre><span>for</span><span>(</span>j <span>in</span> <span>2</span><span>:</span>m<span>)</span> R1 <span>=</span> <span>rbind</span><span>(</span>R1,function1<span>(</span>Xlist<span>[</span><span>[</span>j<span>]</span><span>]</span><span>)</span><span>)</span></pre></td></tr></table></div>
<p>The create on your side the following objects</p>
<div><table><tr><td><pre>1
2
3
4
</pre></td><td><pre>Q1 <span>=</span> <span>qr.<span>Q</span></span><span>(</span><span>qr</span><span>(</span><span>as.<span>matrix</span></span><span>(</span>R1<span>)</span><span>)</span><span>)</span>
R2 <span>=</span> <span>qr.<span>R</span></span><span>(</span><span>qr</span><span>(</span><span>as.<span>matrix</span></span><span>(</span>R1<span>)</span><span>)</span><span>)</span>
Q2list<span>=</span><span>list</span><span>(</span><span>)</span>
<span>for</span><span>(</span>j <span>in</span> <span>1</span><span>:</span>m<span>)</span> Q2list<span>[</span><span>[</span>j<span>]</span><span>]</span> <span>=</span> Q1<span>[</span><span>(</span>j<span>-</span><span>1</span><span>)</span><span>*</span>k<span>+</span><span>1</span><span>:</span>k,<span>]</span></pre></td></tr></table></div>
<p>Finally, contact one last time the providers, and send one of your objects</p>
<div><table><tr><td><pre>1
2
3
4
</pre></td><td><pre>function2<span>=</span><span>function</span><span>(</span>subX,suby,Q<span>)</span><span>{</span>
Q1<span>=</span><span>qr.<span>Q</span></span><span>(</span><span>qr</span><span>(</span><span>as.<span>matrix</span></span><span>(</span>subX<span>)</span><span>)</span><span>)</span>
Q2<span>=</span>Q
<span>return</span><span>(</span><span>t</span><span>(</span>Q1<span>%*%</span>Q2<span>)</span> <span>%*%</span> suby<span>)</span><span>}</span></pre></td></tr></table></div>
<p>Provider <img alt="j" title="j" src="https://s0.wp.com/latex.php?latex=j&bg=ffffff&fg=000&s=0" class="latex"> should then run <img alt="f_2(\mathbf{X},\mathbf{y},Q_j^{(2)})" title="f_2(\mathbf{X},\mathbf{y},Q_j^{(2)})" src="https://s0.wp.com/latex.php?latex=f_2%28%5Cmathbf%7BX%7D%2C%5Cmathbf%7By%7D%2CQ_j%5E%7B%282%29%7D%29&bg=ffffff&fg=000&s=0" class="latex"> on his part of the data, using also <img alt="Q_j^{(2)}" title="Q_j^{(2)}" src="https://s0.wp.com/latex.php?latex=Q_j%5E%7B%282%29%7D&bg=ffffff&fg=000&s=0" class="latex"> as argument (that we obtained on own side) and that function will return <img alt="(\mathbf{Q}^{(2)}_j\mathbf{Q}^{(1)}_j)^{T}_j\mathbf{y}_j" title="(\mathbf{Q}^{(2)}_j\mathbf{Q}^{(1)}_j)^{T}_j\mathbf{y}_j" src="https://s0.wp.com/latex.php?latex=%28%5Cmathbf%7BQ%7D%5E%7B%282%29%7D_j%5Cmathbf%7BQ%7D%5E%7B%281%29%7D_j%29%5E%7BT%7D_j%5Cmathbf%7By%7D_j&bg=ffffff&fg=000&s=0" class="latex">. For instance, ask the first provider to run</p>
<div><table><tr><td><pre>1
</pre></td><td><pre>sumV <span>=</span> function2<span>(</span>Xlist<span>[</span><span>[</span><span>1</span><span>]</span><span>]</span>,ylist<span>[</span><span>[</span><span>1</span><span>]</span><span>]</span>, Q2list<span>[</span><span>[</span><span>1</span><span>]</span><span>]</span><span>)</span></pre></td></tr></table></div>
<p>and do the same with all providers</p>
<div><table><tr><td><pre>1
</pre></td><td><pre><span>for</span><span>(</span>j <span>in</span> <span>2</span><span>:</span>m<span>)</span> sumV <span>=</span> sumV<span>+</span> function2<span>(</span>Xlist<span>[</span><span>[</span>j<span>]</span><span>]</span>,ylist<span>[</span><span>[</span>j<span>]</span><span>]</span>, Q2list<span>[</span><span>[</span>j<span>]</span><span>]</span><span>)</span></pre></td></tr></table></div>
<div><table><tr><td><pre>1
2
3
4
</pre></td><td><pre><span>solve</span><span>(</span>R2<span>)</span> <span>%*%</span> sumV
         <span>[</span>,<span>1</span><span>]</span>
X1 <span>-</span><span>17.579095</span>
X2   <span>3.932409</span></pre></td></tr></table></div>
<p>which is what we were expecting…</p>
	<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://freakonometrics.hypotheses.org/53283"> R-english – Freakonometrics</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/Z73_T9V0IDw" height="1"></a></div>