{"word_count": 539, "character_count": 4381}
My previous post was explaining how mathematically it was possible to parallelize computation to estimate the parameters of a linear regression. More speficially, we have a matrix [latex]\mathbf{X}[/latex] which is [latex]n\times k[/latex] matrix and [latex]\mathbf{y}[/latex] a [latex]n[/latex]-dimensional vector, and we want to compute [latex]\widehat{\mathbf{\beta}}=[\mathbf{X}^T\mathbf{X}]^{-1}\mathbf{X}^T\mathbf{y}[/latex] by spliting the job. Instead of using the [latex]n[/latex] observations, we've seen that it was to possible to compute "something" using the first [latex]n_1[/latex] rows, then the next [latex]n_2[/latex] rows, etc. Then, finally, we "aggregate" the [latex]m[/latex] objects created to get our overall estimate.Parallelizing on multiple coresLet 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 slave, computing what we've seen in the previous post. Here, the data we use arey = cars$distX = data.frame(1,cars$speed)k = ncol(X)On my laptop, I have three cores, so we will split it in [latex]m=3[/latex] chunkslibrary(parallel)library(pbapply)ncl = detectCores()-1cl = makeCluster(ncl)This is more or less what we will do: we have our dataset, and we split the jobs,We can then create lists containing elements that will be sent to each core, as Ewen suggested,chunk = function(x,n) split(x, cut(seq_along(x), n, labels = FALSE))a_parcourir = chunk(seq_len(nrow(X)), ncl)for(i in 1:length(a_parcourir)) a_parcourir[[i]] = rep(i, length(a_parcourir[[i]]))Xlist = split(X, unlist(a_parcourir))ylist = split(y, unlist(a_parcourir))It is also possible to simplify the QR functions we will usecompute_qr = function(x){  list(Q=qr.Q(qr(as.matrix(x))),R=qr.R(qr(as.matrix(x))))}get_Vlist = function(j){  Q3 = QR1[[j]]$Q %*% Q2list[[j]]  t(Q3) %*% ylist[[j]]}clusterExport(cl, c("compute_qr", "get_Vlist"), envir=environment())Then, we can run our functions on each core. The first one is  QR1 = parLapply(cl=cl,Xlist, compute_qr)note that it is also possible to use  QR1 = pblapply(Xlist, compute_qr, cl=cl)which will include a progress bar (that can be nice when the database is rather large). Then use  R1 = pblapply(QR1, function(x) x$R, cl=cl) %&#62;% do.call("rbind", .)  Q1 = qr.Q(qr(as.matrix(R1)))  R2 = qr.R(qr(as.matrix(R1)))  Q2list = split.data.frame(Q1, rep(1:ncl, each=k))  clusterExport(cl, c("QR1", "Q2list", "ylist"), envir=environment())  Vlist = pblapply(1:length(QR1), get_Vlist, cl=cl)  sumV = Reduce('+', Vlist)and finally the ouput issolve(R2) %*% sumV         [,1]X1 -17.579095X2   3.932409which is what we were expecting...Using multiple sourcesIn 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.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 [latex]j[/latex] should run [latex]f_1(\mathbf{X})[/latex] on his part of the data, that function will return [latex]R^{(1)}_j[/latex]. More precisely, to the first provider, sendfunction1 = function(subX){return(qr.R(qr(as.matrix(subX))))}R1 = function1(Xlist[[1]])and actually, send that function to all providers, and aggregate the outputfor(j in 2:m) R1 = rbind(R1,function1(Xlist[[j]]))The create on your side the following objectsQ1 = qr.Q(qr(as.matrix(R1)))R2 = qr.R(qr(as.matrix(R1)))Q2list=list()for(j in 1:m) Q2list[[j]] = Q1[(j-1)*k+1:k,]Finally, contact one last time the providers, and send one of your objectsfunction2=function(subX,suby,Q){Q1=qr.Q(qr(as.matrix(subX)))Q2=Qreturn(t(Q1%*%Q2) %*% suby)}Provider [latex]j[/latex] should then run [latex]f_2(\mathbf{X},\mathbf{y},Q_j^{(2)})[/latex] on his part of the data, using also [latex]Q_j^{(2)}[/latex] as argument (that we obtained on own side) and that function will return [latex](\mathbf{Q}^{(2)}_j\mathbf{Q}^{(1)}_j)^{T}_j\mathbf{y}_j[/latex]. For instance, ask the first provider to runsumV = function2(Xlist[[1]],ylist[[1]], Q2list[[1]])and do the same with all providersfor(j in 2:m) sumV = sumV+ function2(Xlist[[j]],ylist[[j]], Q2list[[j]])solve(R2) %*% sumV         [,1]X1 -17.579095X2   3.932409which is what we were expecting...