正确地将 for 循环与 R 并行化

Question

我一直在为我的主管开发一个简单的函数集合，这些函数将执行一些简单的初始基因组规模统计，这很容易让我的团队快速了解可能需要更多时间的未来分析 - 因为例如 RDP4 或 BioC（只是为了解释为什么我没有直接去 BioConductor）。我想加快一些速度以允许更大的 contig 大小，所以我决定使用 doParallel 和 foreach 来编辑一些 for 循环以允许这样做。下面是一个简单的函数，它可以识别某些序列（存储为矩阵）中相同的碱基并删除它们。

strip.invar <- function(x) {
  cat("
          Now removing invariant sites from DNA data matrix, this may take some time...
      ")
  prog <- txtProgressBar(min=0, max=ncol(x), style=3)
  removals<-c()
  for(i in 1:ncol(x)){
    setTxtProgressBar(prog, i)
    if(length(unique(x[,i])) == 1) { 
      removals <- append(removals, i)
    }
  }
  newDnaMatrix <- x[,-removals]
  return(newDnaMatrix)
}

在阅读了 doParallel 和 foreach 的介绍后，我尝试制作一个版本以容纳更多内核 - 在我的 Mac 上这是 8 个 - 每个内核有两个线程的四核 - 8 个虚拟内核：

strip.invar <- function(x, coresnum=detectCores()){
  cat("
          Now removing invariant sites from DNA data matrix, this may take some time...
          ")
  prog <- txtProgressBar(min=0, max=ncol(x), style=3)
  removals<-c()
  if(coresnum > 1) {
    cl <- makeCluster(coresnum)
    registerDoParallel(cl)
    foreach(i=1:ncol(x)) %dopar% {
      setTxtProgressBar(prog, i)
      if(all(x[,i] == x[[1,i]])){
        removals <- append(removals, i)
      }
    }
  } else {
    for(i in 1:ncol(x)){
      setTxtProgressBar(prog, i)
      if(length(unique(x[,i])) == 1) { 
        removals <- append(removals, i)
      }
    }
  }
  newDnaMatrix <- x[,-removals]
  return(newDnaMatrix)
}

但是，如果我运行它并将核心数设置为 8，我不完全确定它是否有效 - 我看不到进度条做任何事情，但后来我听说打印到屏幕和涉及图形的东西设备在 R 中的并行计算很棘手。但这似乎仍然需要一些时间，而且我的笔记本电脑变得“非常”热，所以我不确定我是否正确地做到了这一点，我在看到几个例子后尝试过（我成功地在小插图中运行了一个很好的引导示例），但我一定会遇到学习障碍。顺便说一句，我想我也想问问人们的意见，对于涉及循环或应用的 R 代码瓶颈的最佳加速是什么 - 并行化它，还是 Rcpp？

谢谢。

Answer 1

我的另一个答案不正确，因为 colmean 等于第一个值不足以作为唯一值数量的测试。所以这里是另一个答案：

您可以使用apply 优化循环。

set.seed(42)
dat <- matrix(sample(1:5,1e5,replace=TRUE),nrow=2)
res1 <- strip.invar(dat)


strip.invar2 <- function(dat) {
  ix <- apply(dat,2,function(x) length(unique(x))>1)
  dat[,ix]}

res2 <- strip.invar2(dat)

all.equal(res1,res2)
#TRUE
library(microbenchmark)
microbenchmark(strip.invar(dat),strip.invar2(dat),times=10)
#Unit: milliseconds
#             expr       min        lq    median       uq      max neval
#strip.invar(dat)  2514.7995 2529.2827 2547.6751 2678.464 2792.405    10
#strip.invar2(dat)  933.3701  945.5689  974.7564 1008.589 1018.400    10

这大大提高了性能，但不如矢量化可能实现的那么多。

并行化在这里不会提供更好的性能，因为每次迭代本身不需要太多性能，因此并行化开销实际上会增加所需的时间。但是，您可以并行拆分数据和处理块。

Answer 2

首先，尝试运行cl <- makeCluster( coresnum-1 )。主 R 进程已经在使用您的一个核心，并用于分派和接收从属作业的结果，因此您有 7 个空闲核心用于从属作业。我认为您将有效地排队您的 foreach 循环之一，以等待之前的循环之一完成，因此这项工作将需要更长的时间才能完成。

其次，您通常在非并行环境中运行此功能的控制台上看到的内容仍会打印到控制台，只是每个作业输出都会打印到从属进程控制台，因此您不会看到它。但是，您可以将不同 foreach 循环的输出保存到文本文件中以检查它们。 Here is an example 如何保存控制台输出。将代码粘贴到您的 foreach 语句中。

您的笔记本电脑会变得非常热，因为在您运行此作业时，您的所有内核都在以 100% 的容量工作。

我发现foreach 包提供了一组出色的函数来提供简单的并行处理。 Rcpp 可能（会？！）给您带来更好的性能，但是您在编写 C++ 代码方面的表现如何？这个函数的运行时间是多少，它的使用频率如何？我总是首先考虑这些事情。