以下所有内容摘自《Prallel Computing for Data Science with Examples in R, C++ and CUDA》（数据科学中的并行计算：以R、C++和CUDA为例）这本书，作者是美国加州大学戴维斯分校的统计系创始人之一，统计学和计算机双教授，还精通汉语（是的，你没看错，这老头就是这么牛逼。。。），后面花点时间准备再摘几个经典并实用的，对这块感兴趣的强烈推荐这本书！
library(parallel)
doichunk <- function(ichunk) {
tot <- 0
nr <- nrow(lnks) # lnks在worker处是全局变量
for (i in ichunk) {
tmp <- lnks[(i+1):nr, ] %*% lnks[i, ]
tot <- tot + sum(tmp)
}
tot
}
mutoutpar <- function(cls, lnks) {
nr <- nrow(lnks) # lnks在manager处是全局变量
clusterExport(cls, "lnks")
ichunks <- clusterSplit(cls, 1:(nr - 1)) #把块１:125发送给第一个worker......375:499发送给第四个worker
tots <- clusterApply(cls, ichunks, doichunk) # 返回的列表包含四个元素，而不是之前的499个
Reduce(sum, tots) / nr
}
initcls <- function(workers) {
makeCluster(spec = workers)
}
cls <- initcls(4)
clusterSplit(cls, 1:50) # 返回四个列表
snowsim <- function(nr,nc, cls) {
lnks <<-
matrix(sample(0:1, (nr*nc), replace= TRUE),
nrow = nr)
system.time(mutoutpar(cls, lnks))
}
initmc <- function(nworkers) {
makeCluster(nworkers)
}
cl2 <- initmc(2)
snowsim(2000, 2000, cl2)
可以看到，相对之前没有分块，elapsed时间明显减少
主要的问题是预测变量的选择：一方面，要在回归方程中包含所有相关的预测变量。另一方面，我们必须避免过拟合。
假设我们有ｎ个观测值，ｐ个预测变量。通过无偏估计$adjusted$ $R^2$来选择模型
一共有$2^p$种可能的模型，因此计算量会相当大————这非常适合使用并行计算。在这里有两种可能：
选项(a)存在问题。对一个给定的包含m个预测变量的集合，我们必须先计算各种平方和与乘积和。每个加法都有ｎ个被加数，一共有$O(m^2)$个加法计算，这使得计算复杂度为$O(nm^2)$。之后必须进行矩阵求逆（或者其它等价计算，比如QR分解)，其复杂度为$O(m^３)^2$。
矩阵求逆迄今仍然不是一个易并行的计算，所以选项(b)简单很多，是易并行的，事实上它包含一个循环。
下面并行由snow包实现，它在所有满足条件的模型中计算校正决定系数值，这些模型预测变量集合的大小最大为k。用户可以选择静态或者动态调度，或者反向调度，用户可以指定一个（不变的）块的大小。
算法的整体策略：
# 回归响应变量列Ｙ以及Xi预测变量的所有可能子集
# 子集规模最大为Ｋ
# 返回每个子集的校正决定系数
# 调度参数：
# static(clusterApply())
# dynamic(clusterApply())
# 颠倒任务的顺序
# 块大小（动态情形下）
# 参数：
# cls: Snow集群
# x: 预测变量矩阵，每列一个
# y: 响应变量向量
# k: 预测变量集合大小的最大值
# reverse: TRUE表示对迭代顺序进行颠倒
# dyn: TRUE表示动态调度
# chunksize: 调度块的规模
# 返回值：
# R矩阵，显示校正决定系数
# 使用预测集合来进行索引
snowapr <- function(cls, x, y, k, reverse = F, dyn = F, chunksize = 1) {
require(parallel)
p <- ncol(x)
# 生成预测变量子集，一个R的list，1个元素代表一个预测变量子集
allcombs <- genallcombs(p, k)
ncombs <- length(allcombs)
clusterExport(cls, "dolpset")
# 设定任务索引
tasks <- if(!reverse)
seq(1, ncombs, chunksize)
else
seq(ncombs, 1, -chunksize)
if(!dyn) {
out <- clusterApply(cls, tasks, dochunk, x, y, allcombs, chunksize)
}
# out的每个元素都由校正决定系数(adjusted R^2)和产生这个值的预测变量集合的索引构成，然后把这些变量集合到一个矩阵中
Reduce(rbind, out)
}
# 生成１..p的所有大小 <= k的非空集合；
# 返回一个索引向量形式的R list，每个元素代表一个预测变量集合
genallcombs <- function(p, k) {
allcombs <- list()
for (i in 1:k) {
tmp <- combn(1:p, i)
allcombs <- c(allcombs, matrixtolist(tmp, rc = 2))
}
allcombs
}
# 从矩阵中提取行(rc = 1)或列(rc = 2), 生成一个list
matrixtolist <- function(rc, m) {
if (rc == 1) {
Map(function(rownum) m[rownum, ], 1:nrow(m))
} else
Map(function(colnum) m[, colnum], 1:ncol(m))
}
# 处理allcombs块中的所有预测变量集合
# 这个分块的第一个索引是psetsstart
dochunk <- function(psetsstart, x, y, allcombs, chunksize) {
ncombs <- length(allcombs)
lasttask <- min(psetsstart + chunksize - 1, ncombs)
t(sapply(allcombs[psetsstart:lasttask], dolpset, x, y))
}
# 找到指定预测变量集合onepset的校正决定系数(adjusted R^2)；
# 返回值是校正决定系数，紧跟着用０来填充空位的预测变量集合的索引
# 为了方便，对dolpset()的调用所返回的向量的长度都是k+1
# 例如，对k=4,(0.28, 1, 3, 0, 0)意味着预测变量集合是由x的第１和第３列，Ｒ平方值为0.28
dolpset <- function(onepset, x, y) {
slm <- summary(lm(y ~ x[, onepset]))
nOs <- ncol(x) - length(onepset)
c(slm$adj.r.squared, onepset, rep(0, nOs))
}
# 看起来最好的预测变量集合
snowtest <- function(cls, n, p, k, chunksize = 1, dyn = F, rvrs = F) {
gendata(n, p)
snowapr(cls, x, y, k, rvrs, dyn, chunksize)
}
gendata <- function(n, p) {
x <<- matrix(rnorm(n*p), ncol = p)
y <<- x %*% c(rep(0.5, p)) + rnorm(n)
}
initmc <- function(nworkers) {
makeCluster(nworkers)
}
cl8 <- initmc(8)
# 在这里我们生成了大小n = 100的模拟数据，有p = 4个预测变量，预测变量集合的最大k =2，对使用预测变量１和３的模型，
# 即ｘ的第2列和第4列，校正决定系数的最大值大约是0.398
snowtest(cl8, 100, 4, 2)