数据表和并行计算

Question

在这篇文章之后：multicore and data.table in R，我想知道在使用 data.table 时是否有办法使用所有内核，通常按组进行计算可以并行化。 plyr 似乎是有意允许此类操作的。

Answer 1

首先要检查的是data.table FAQ 3.1 第 2 点已经陷入：

只为最大的组分配一个内存，然后 内存被其他组重用。垃圾很少 收集。

这就是 data.table 分组快速的原因之一。但是这种方法不适合并行化。并行化意味着将数据复制到其他线程，而不是花费时间。但是，我的理解是data.table 分组通常比plyr 更快，无论如何.parallel。这取决于每组任务的计算时间，以及是否可以轻松减少该计算时间。移动数据通常占主导地位（在对大数据任务运行 1 或 3 次进行基准测试时）。

到目前为止，更多时候，实际上是一些问题在[.data.table 的j 表达式中引起了人们的注意。例如，最近我们看到 data.table 分组表现不佳，但罪魁祸首竟然是 min(POSIXct) (Aggregating in R over 80K unique ID's)。避免这个问题会产生超过 50 倍的加速。

所以口头禅是：Rprof、Rprof、Rprof。

此外，同一常见问题解答中的第 1 点可能很重要：

只有该列被分组，其他 19 列被忽略，因为 data.table 检查 j 表达式并意识到它没有使用 其他列。

所以，data.table 根本不遵循拆分-应用-组合范式。它的工作方式不同。 split-apply-combine 适合并行化，但它确实不能扩展到大数据。

另请参阅 data.table 介绍小插图中的脚注 3：

我们想知道有多少人将并行技术部署到代码中 也就是矢量扫描

这是在说“当然，并行的速度要快得多，但是使用高效的算法真正需要多长时间？”。

但是，如果您已经进行了分析（使用 Rprof），并且每个组的任务确实是计算密集型的，那么 datatable-help 上的 3 个帖子（包括“多核”一词）可能会有所帮助：

multicore posts on datatable-help

当然，有很多任务在 data.table 中并行化会很好，并且有一种方法可以做到这一点。但它还没有完成，因为通常其他因素会咬人，所以它的优先级很低。如果您可以发布具有基准和 Rprof 结果的可重现虚拟数据，那将有助于提高优先级。

Answer 2

我已经按照@matt dowle 之前的 Rprof、Rprof、Rprof 口头禅做了一些测试。

我发现并行化的决定取决于上下文；但可能很重要。根据测试操作（例如下面的foo，可以自定义）和使用的核心数量（我尝试了8个和24个），我得到了不同的结果。

以下结果：

1. 使用 8 个内核，我看到在这个并行化示例中 21% 的改进
2. 使用 24 个内核，我看到 14% 的改进。

我还查看了一些真实世界（不可共享）的数据/操作，这些数据/操作显示了更大的（33% 或25%，两个不同的测试）改进并行与 24 个内核。 2018 年 5 月编辑一组新的真实示例案例显示，与 1000 个组的并行操作相比，改进了接近 85%。

R> sessionInfo() # 24 core machine:
R version 3.3.2 (2016-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods
[8] base

other attached packages:
[1] microbenchmark_1.4-2.1 stringi_1.1.2          data.table_1.10.4

R> sessionInfo() # 8 core machine:
R version 3.3.2 (2016-10-31)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: macOS Sierra 10.12.4

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] microbenchmark_1.4-2.1 stringi_1.1.5          data.table_1.10.4     

以下示例：

library(data.table)
library(stringi)
library(microbenchmark)

set.seed(7623452L)
my_grps <- stringi::stri_rand_strings(n= 5000, length= 10)

my_mat <- matrix(rnorm(1e5), ncol= 20)
dt <- data.table(grps= rep(my_grps, each= 20), my_mat)

foo <- function(dt) {
  dt2 <- dt ## needed for .SD lock
  nr <- nrow(dt2)

  idx <- sample.int(nr, 1, replace=FALSE)

  dt2[idx,][, `:=` (
    new_var1= V1 / V2,
    new_var2= V4 * V3 / V10,
    new_var3= sum(V12),
    new_var4= ifelse(V10 > 0, V11 / V13, 1),
    new_var5= ifelse(V9 < 0, V8 / V18, 1)
  )]


  return(dt2[idx,])
}

split_df <- function(d, var) {
  base::split(d, get(var, as.environment(d)))
}

foo2 <- function(dt) {
  dt2 <- split_df(dt, "grps")

  require(parallel)
  cl <- parallel::makeCluster(min(nrow(dt), parallel::detectCores()))
  clusterExport(cl, varlist= "foo")
  clusterExport(cl, varlist= "dt2", envir = environment())
  clusterEvalQ(cl, library("data.table"))

  dt2 <- parallel::parLapply(cl, X= dt2, fun= foo)

  parallel::stopCluster(cl)
  return(rbindlist(dt2))
}

print(parallel::detectCores()) # 8

microbenchmark(
  serial= dt[,foo(.SD), by= "grps"],
  parallel= foo2(dt),
  times= 10L
)

Unit: seconds
     expr      min       lq     mean   median       uq      max neval cld
   serial 6.962188 7.312666 8.433159 8.758493 9.287294 9.605387    10   b
 parallel 6.563674 6.648749 6.976669 6.937556 7.102689 7.654257    10  a 

print(parallel::detectCores()) # 24

Unit: seconds
     expr       min        lq     mean   median       uq      max neval cld
   serial  9.014247  9.804112 12.17843 13.17508 13.56914 14.13133    10   a
 parallel 10.732106 10.957608 11.17652 11.06654 11.30386 12.28353    10   a

分析：

我们可以使用this answer 来更直接地回应@matt dowle 对分析的原始评论。

因此，我们确实看到大部分计算时间由base 处理，而不是data.table。正如预期的那样，data.table 操作本身非常快。虽然有些人可能会争辩说这证明data.table 内不需要并行性，但我认为这个工作流/操作集不是非典型的。也就是说，我强烈怀疑大多数大型data.table 聚合涉及大量非data.table 代码；并且这与交互式使用与开发/生产使用相关。因此，我得出结论，对于大型聚合，data.table 中的并行性将是有价值的。

library(profr)

prof_list <- replicate(100, profr::profr(dt[,foo(.SD), by= "grps"], interval = 0.002),
                       simplify = FALSE)

pkg_timing <- fun_timing <- vector("list", length= 100)
for (i in 1:100) {
  fun_timing[[i]] <- tapply(prof_list[[i]]$time, paste(prof_list[[i]]$source, prof_list[[i]]$f, sep= "::"), sum)
  pkg_timing[[i]] <- tapply(prof_list[[i]]$time, prof_list[[i]]$source, sum)
}

sort(sapply(fun_timing, sum)) #  no large outliers

fun_timing2 <- rbindlist(lapply(fun_timing, function(x) {
  ret <- data.table(fun= names(x), time= x)
  ret[, pct_time := time / sum(time)]
  return(ret)
}))

pkg_timing2 <- rbindlist(lapply(pkg_timing, function(x) {
  ret <- data.table(pkg= names(x), time= x)
  ret[, pct_time := time / sum(time)]
  return(ret)
}))

fun_timing2[, .(total_time= sum(time),
                avg_time= mean(time),
                avg_pct= round(mean(pct_time), 4)), by= "fun"][
  order(avg_time, decreasing = TRUE),][1:10,]

pkg_timing2[, .(total_time= sum(time),
                avg_time= mean(time),
                avg_pct= round(mean(pct_time), 4)), by= "pkg"][
  order(avg_time, decreasing = TRUE),]

结果：

                      fun total_time avg_time avg_pct
 1:               base::[    670.362  6.70362  0.2694
 2:      NA::[.data.table    667.350  6.67350  0.2682
 3:       .GlobalEnv::foo    335.784  3.35784  0.1349
 4:              base::[[    163.044  1.63044  0.0655
 5:   base::[[.data.frame    133.790  1.33790  0.0537
 6:            base::%in%    120.512  1.20512  0.0484
 7:        base::sys.call     86.846  0.86846  0.0348
 8: NA::replace_dot_alias     27.824  0.27824  0.0112
 9:           base::which     23.536  0.23536  0.0095
10:          base::sapply     22.080  0.22080  0.0089

          pkg total_time avg_time avg_pct
1:       base   1397.770 13.97770  0.7938
2: .GlobalEnv    335.784  3.35784  0.1908
3: data.table     27.262  0.27262  0.0155

交叉发布在github/data.table

Answer 3

是的（不过，这可能不值得，正如@Alex W 所指出的那样）。

以下提供了一个简单的模式来做到这一点。为简单起见，我使用了一个不值得的示例（使用mean 函数），但它显示了模式。

例子：

假设您想在 iris 数据集中按物种计算 Petal.Length 的平均值。

您可以直接使用 data.table 作为：

as.data.table(iris)[by=Species,,.(MPL=mean(Petal.Length))]
      Species   MPL
1:     setosa 1.462
2: versicolor 4.260
3:  virginica 5.552

但是，如果 mean 是一个运行时间足够长且成本高昂的计算（也许是通过分析确定的，尽管有时它只是“显而易见的”），您可能希望使用 parallel::mclapply。由于最小化与 mclapply 生成的所有子进程的通信可以大大减少整体计算，而不是将选择从 data.table 传递到每个子进程，您只想传递选择的索引。此外，通过首先对 data.table 进行排序，您可以只传递这些索引的范围（最大值和最小值）。像这样：

> o.dt<-as.data.table(iris)[order(Species)] # note: iris happens already to be ordered
> i.dt<-o.dt[,by=Species,.(irange=.(range(.I)))]
> i.dt
      Species  irange
1:     setosa    1,50
2: versicolor  51,100
3:  virginica 101,150


> result<-mclapply(seq(nrow(i.dt)),function(r) o.dt[do.call(seq,as.list(i.dt[r,irange][[1]])),.(MPL=mean(Petal.Length))])
> result
[[1]]
     MPL
1: 1.462

[[2]]
    MPL
1: 4.26

[[3]]
     MPL
1: 5.552

> result.dt<-cbind(i.dt,rbindlist(result))[,-2]
> result.dt
      Species   MPL
1:     setosa 1.462
2: versicolor 4.260
3:  virginica 5.552

回顾模式：

• 对输入进行排序。
• 计算每个组的索引范围。
• 定义匿名function 以提取包含组成员的行，并执行所需的计算（在本例中为均值）。
• 在索引范围的行索引上使用 mclapply 将该函数应用于每个组。
• 使用rbindlist 将结果作为data.table 获取，cbind 将其添加到输入中，然后将其删除索引列（除非您出于其他原因需要保留它们）。

注意事项：

• 最终的rbindlist 通常很昂贵，可以根据您的应用跳过）。

待办事项：

• 让 data.table 团队相信此模式足够通用且足够有用，以至于其他 data.table 索引选项应该调用它。想象一下，传递 mc=TRUE 将调用此模式，并在 ...
中支持其他并行选项

iris.dt[by=Species,,.(MPL=mean(Petal.Length)), mc=TRUE, mc.preschedule=FALSE, mc.set.seed=TRUE,...]