那些在R中进行并行计算的坑

最近在疯狂做模拟数据分析,需要大量在R中进行并行运算。在这过程中被R无情的虐的体无完肤,所以整理一下踩过的坑,希望能帮助到有相同困难的盆友。

doParallelforeach

这两个包一般合并在一起使用,可以实现类似与for循环这样的结构。总的来说,foreach包还是非常有效率和简洁的。并行计算的模板是foreach(args) %dopar% {expression}

然而,我发现一些问题有:

  • 很难计算单个运算单元的运行时间。总共所有条件$\times$重复用的时间可以用system.time()来记录,但是对于单个循环的时间很难记录。这样的话在一个很多条件的模拟研究中很难事前估记大概要花多久。只能先进行一轮计算,然后大概估计。
  • 对于内存的占用率很高。一般来说这是fork并行运算的通病,当然如果是socket方式的话io之间的信息传递速度就是瓶颈了。当内存过高后,foreach会提示no enough vmem。一个解决方法是使用sequential和parallel共同使用的方法,这样的话内存使用率就会限制在conditions之间而不是iterations $\times$ conditions。比如:
for(i in 1:iterations){
  foreach(n in 1:conditions) %dopar% {
    expression...
  }
}

multidplyr

这个包是Hadley用来填补tidyverse没有并行处理从而开发的。这个包的原理是将数据用partition()来分布到不同的clusters,然后用group_bysummarisedplyr传统方法来进行运算。然而,在尝试之后,我觉得非常不好用:

  • 使用范围非常局限。必须将所有的conditions的数据用nested方法打包成一个tibble才好使用。这样的一个超大文件会大量占用内存空间,如果你想把这样的文件保存成一个RDS文件的话,会超大无比。比如一个10000 conditions的实验,算上原始数据,设定参数和运算结果会占用1~10GB。你不会想把这些都一股脑的放进一个R session的。
  • 必须对所有的数据进行统一的运算。用summarise方法进行呈现的话,只能提供一个或多个方法的结果。而不能实现诸如 调取数据-运算-得出结果 这样三段式方法。如果不同条件要进行不同的处理的话就会有问题。
  • 并没有觉得很快。可能是这种long tibble拆分的方式局限在数据层,而没有在方法层进行并行处理。
  • 前面的foreach包有的缺点都有。

HPC物理方法

除了软件方法,还可以用物理方法。通过合理使用高性能计算机集群(High Performance Cluster),可以实现更加强大的并行处理。我们爱荷华大学的HPC使用Argon进行任务管理。在实现上可以使用多job提交的方式来实现并行处理,比如不同的condition提供不同的jobID。学生好像可以一次性并行提交5个job tasks。

以下是一个模板。makeArgoJobFile方法可以批量制造hpc的job文件。使用system()来让R运行命令行。使用qsub -t 1-[nreps] ****.job的bash命令来进行同一个条件重复模拟(nreps替换成重复数量,比如10000)并提交给HPC集群。

# loop through conditions and submit job file with replications as array
index = 0
items = n_items_per_factor[1]
for (items in n_items_per_factor){

  obs = n_obs[1]
  for (obs in n_obs){

    factors = n_Factors[1]
    for (factors in n_Factors){

      index = index+1
      job = makeArgonJobFile(n_items_per_factor = items,
                             n_obs = obs,
                             n_Factors = factors,
                             argon_Rscript_path = "~/PPMCsim/CFA/Script/argon_irt_onerep.R",
                             argon_job_path = paste0("~/PPMCsim/CFA/Results", index, ".sh"),
                             seed = initial.seed+index*100000,
                             n_reps = 1)

      write.table(job, file = paste0("~/PPMCsim/CFA/Results/argonjob_", obs, "_", factors,
                                     "_", items, ".job"), quote=FALSE, row.names=FALSE, col.names=FALSE)

      system(paste0("qsub -t 1-", n_reps, " ~/PPMCsim/CFA/Results/argonjob_", obs, "_", factors, "_", items, ".job"))

    }
  }
}

用这个方法的问题是这个方法局限于对同一个条件(condition)进行重复(repetition)。而有些项目不是重复同一个条件,而是无数条件中抽样而不重复,这样的话就无法享受qsub -t这个命令了。

Bonus:使用OpenBLAS对矩阵运算进行优化

下载openblas,替换掉R自带对blas库。

brew install openblas
ln -sf /usr/local/opt/openblas/lib/libblas.dylib /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.dylib
~/Downloads                                                                                                                             
▶ Rscript 1-r-benchmark-25.R         
Loading required package: Matrix


   R Benchmark 2.5
   ===============
Number of times each test is run__________________________:  3

   I. Matrix calculation
   ---------------------
Creation, transp., deformation of a 2500x2500 matrix (sec):  0.722333333333334 
2400x2400 normal distributed random matrix ^1000____ (sec):  0.183333333333333 
Sorting of 7,000,000 random values__________________ (sec):  0.824 
2800x2800 cross-product matrix (b = a' * a)_________ (sec):  0.221 
Linear regr. over a 3000x3000 matrix (c = a \ b')___ (sec):  0.200333333333333 
                      --------------------------------------------
                 Trimmed geom. mean (2 extremes eliminated):  0.317415197050083 

   II. Matrix functions
   --------------------
FFT over 2,400,000 random values____________________ (sec):  0.362666666666669 
Eigenvalues of a 640x640 random matrix______________ (sec):  0.589 
Determinant of a 2500x2500 random matrix____________ (sec):  0.242666666666667 
Cholesky decomposition of a 3000x3000 matrix________ (sec):  0.261666666666668 
Inverse of a 1600x1600 random matrix________________ (sec):  0.286666666666666 
                      --------------------------------------------
                Trimmed geom. mean (2 extremes eliminated):  0.300753769812312 

   III. Programmation
   ------------------
3,500,000 Fibonacci numbers calculation (vector calc)(sec):  0.244333333333332 
Creation of a 3000x3000 Hilbert matrix (matrix calc) (sec):  0.277666666666666 
Grand common divisors of 400,000 pairs (recursion)__ (sec):  0.259 
Creation of a 500x500 Toeplitz matrix (loops)_______ (sec):  0.0336666666666664 
Escoufier's method on a 45x45 matrix (mixed)________ (sec):  0.255000000000003 
                      --------------------------------------------
                Trimmed geom. mean (2 extremes eliminated):  0.252701345790627 


Total time for all 15 tests_________________________ (sec):  4.96333333333334 
Overall mean (sum of I, II and III trimmed means/3)_ (sec):  0.288945177268809 
                      --- End of test ---


~/Downloads                                                                                                                             
▶ ln -sf /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.dylib 

~/Downloads                                                                                                                             
▶ Rscript 1-r-benchmark-25.R
Loading required package: Matrix


   R Benchmark 2.5
   ===============
Number of times each test is run__________________________:  3

   I. Matrix calculation
   ---------------------
Creation, transp., deformation of a 2500x2500 matrix (sec):  0.747666666666667 
2400x2400 normal distributed random matrix ^1000____ (sec):  0.328 
Sorting of 7,000,000 random values__________________ (sec):  0.790999999999999 
2800x2800 cross-product matrix (b = a' * a)_________ (sec):  14.37 
Linear regr. over a 3000x3000 matrix (c = a \ b')___ (sec):  6.509 
                      --------------------------------------------
                 Trimmed geom. mean (2 extremes eliminated):  1.5672306794424 

   II. Matrix functions
   --------------------
FFT over 2,400,000 random values____________________ (sec):  0.255333333333335 
Eigenvalues of a 640x640 random matrix______________ (sec):  0.912333333333336 
Determinant of a 2500x2500 random matrix____________ (sec):  3.26666666666666 
Cholesky decomposition of a 3000x3000 matrix________ (sec):  5.97066666666668 
Inverse of a 1600x1600 random matrix________________ (sec):  3.55833333333334 
                      --------------------------------------------
                Trimmed geom. mean (2 extremes eliminated):  2.1970249936814 

   III. Programmation
   ------------------
3,500,000 Fibonacci numbers calculation (vector calc)(sec):  0.250333333333335 
Creation of a 3000x3000 Hilbert matrix (matrix calc) (sec):  0.294999999999997 
Grand common divisors of 400,000 pairs (recursion)__ (sec):  0.305333333333323 
Creation of a 500x500 Toeplitz matrix (loops)_______ (sec):  0.0336666666666569 
Escoufier's method on a 45x45 matrix (mixed)________ (sec):  0.370000000000005 
                      --------------------------------------------
                Trimmed geom. mean (2 extremes eliminated):  0.282512914658612 


Total time for all 15 tests_________________________ (sec):  37.9633333333333 
Overall mean (sum of I, II and III trimmed means/3)_ (sec):  0.99083668118365 
                      --- End of test ---

结语

这篇小日志只是简单了列举了一些我最近在使用的工具。但是还有很多方法可以进行并行处理,比如用snow包和parallel包。或者使用更底层的mpi进行不同运算单元的数据交换。如果我以后学习到更多并行处理的知识,我也会继续更新的。Keep learning!

张吉鸿
张吉鸿
自由的灵魂

我的研究兴趣是贝叶斯分析和潜变量分析