Yulong Niu

个人博客

使用RcppParallel并行计算

Posted at — Jul 19, 2018

在之前的博文中，我详细讨论了使用多种R包实现并行计算。其中，提到一个非常重要的问题：

当循环数很大时（1万以上），`foreach`会变得非常慢。

这个问题在Florian Privé的A guide to parallelism in R中也提到，解释是foreach每次只合并100个循环结果。

1. 测试

我尝试使用RcppParallel包调用C++的并行方法。结论是：在循环数很大时，RcppParallel包提供的并行方法优于foreach。

一个简单的测试场景：对一个数值向量的每个元素做平方根运算，结果按原始顺序返回。在Gist1和Gist2中，分别实现了：

首先，比较5种实现效率，并行计算调用8个线程（Intel i7-4790 CPU@3.6GHz）。测试结果显示SqrtRPara（使用foreach）和非向量化的R版本SqrtR效率较低。

tmp1 <- runif(10e3)

all.equal(SqrtCpp(tmp1),
          sqrt(tmp1),
          SqrtR(tmp1),
          SqrtRforeach(tmp1),
          SqrtRParSapply(tmp1),
          SqrtCppPara(tmp1))

## TRUE

microbenchmark(
    SqrtCpp(tmp1),
    sqrt(tmp1),
    SqrtR(tmp1),
    SqrtRforeach(tmp1),
    SqrtRParSapply(tmp1),
    SqrtCppPara(tmp1))

## Unit: microseconds
##                 expr         min          lq         mean       median
##        SqrtCpp(tmp1)      56.295      72.648 9.338755e+01      82.0335
##           sqrt(tmp1)      36.216      46.074 4.865115e+01      48.3090
##          SqrtR(tmp1)    3030.682    3116.718 4.229971e+03    3947.9380
##   SqrtRforeach(tmp1) 1488851.181 1532937.096 1.561865e+06 1547849.9610
## SqrtRParSapply(tmp1)  954757.348  963478.755 9.701841e+05  969925.9090
##    SqrtCppPara(tmp1)      23.837      79.314 1.069003e+02     104.5975
##           uq         max neval
##      89.0800    1183.279   100
##      52.2995      66.875   100
##    4560.0930   10391.379   100
## 1584297.5760 1750382.995   100
##  974233.5690 1012400.281   100
##     111.9160    1331.442   100

然后，增加循环数，比较效率较高的前三种方法。测试结果显示调用RcppParallel包的C++并行版本SqrtCppPara胜出。

tmp1 <- runif(10e6)

all.equal(SqrtCpp(tmp1),
          sqrt(tmp1),
          SqrtCppPara(tmp1))

## TRUE

microbenchmark(
  SqrtCpp(tmp1),
  sqrt(tmp1),
  SqrtCppPara(tmp1))
  
## Unit: milliseconds
##               expr      min       lq     mean   median       uq       max neval
##      SqrtCpp(tmp1) 76.68263 78.55146 82.51442 79.48709 87.03026 100.56873   100
##         sqrt(tmp1) 52.19705 53.67441 58.16940 54.60512 66.67642  70.94672   100
##  SqrtCppPara(tmp1) 37.10116 38.34199 42.23896 39.17889 42.98785  61.94529   100

2. 使用vector代替List

在使用RcppParallel并行计算时，不能在并行循环中调用Rcpp::List对象。一个解决办法是：使用std::vector替代Rcpp:List。例如，List中都是数值向量，那么可以建立std::vector<Rcpp::NumericVector>对象替代。 Gist3中提供了一个例子。这种方法的局限在于List中每一个元素的类型需要相同。

3. 同步

如果多个线程同时操作某一个共享内存对象，需要在RcppParallel包中使用“锁”。如Gist4所示，多个线程都需要操作estcount对象。测试代码如下：

library('Rcpp')

sourceCpp('TestSynchron.cpp')

n <- 10
g <- 1000
ecin <- sample(0:9, g*n, replace = TRUE) %>%
  split(1:g)

ecin %>% unlist %>% table

TestShare(ecin, 10)

如果去掉Gist4代码中的第27和31行，可以发现测试结果不正确。

参考网址

更新记录

2018年10月5日