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

当循环数很大时(1万以上),`foreach`会变得非常慢。

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

1. 测试

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

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

  • SqrtR:用循环非并行操作每个元素。这种方法在R语言编程中不推荐,而应该尽量“向量化”操作。

  • SqrtRforeachforeach并行版本。

  • SqrtRParSapply: parSapply并行版本。

  • SqrtCppC++非并行版本。

  • SqrtCppParaRcppParallel包的C++并行版本。

  • sqrt:R内置的向量化方法,C非并行版本。

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

5 versions
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
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胜出。

top 3 versions
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
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对象。测试代码如下:

Synchronization
1
2
3
4
5
6
7
8
9
10
11
12
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代码中的第2731行,可以发现测试结果不正确。

参考网址

更新记录

2018年10月5日

Comments