R 的高性能计算

王诗翔 · 2018-08-29

分类: r
标签: r HPC

《R语言编程指南》

代码性能

简单看函数运行时间，用system.time()函数。

高级点，使用microbenchmark包的microbenchmark()函数。

例如：

library(microbenchmark)
x = rnorm(10000)
microbenchmark(x)
#> Unit: nanoseconds
#>  expr min lq mean median uq  max neval
#>     x  32 33 76.5     34 35 4029   100

给出了类似summary()函数结果的样子，默认计算100次。

代码性能分析

R提供了内置函数Rprof()对代码性能进行分析。

使用方法：调用Rprof()开始分析，运行想分析的代码，然后调用Rprof(NULL)停止分析，最后调用summaryRprof()查看分析结果：

x = rnorm(10000)
tmp = tempfile(fileext = ".out")
Rprof(tmp)
for (i in 1:1000) {
    cumsum(x)
}
Rprof(NULL)
summaryRprof(tmp)
#> $by.self
#>        self.time self.pct total.time total.pct
#> "eval"      0.04      100       0.04       100
#> 
#> $by.total
#>                       total.time total.pct self.time self.pct
#> "eval"                      0.04       100      0.04      100
#> "block_exec"                0.04       100      0.00        0
#> "call_block"                0.04       100      0.00        0
#> "eval.parent"               0.04       100      0.00        0
#> "evaluate"                  0.04       100      0.00        0
#> "evaluate::evaluate"        0.04       100      0.00        0
#> "evaluate_call"             0.04       100      0.00        0
#> "handle"                    0.04       100      0.00        0
#> "in_dir"                    0.04       100      0.00        0
#> "knitr::knit"               0.04       100      0.00        0
#> "local"                     0.04       100      0.00        0
#> "process_file"              0.04       100      0.00        0
#> "process_group"             0.04       100      0.00        0
#> "process_group.block"       0.04       100      0.00        0
#> "rmarkdown::render"         0.04       100      0.00        0
#> "timing_fn"                 0.04       100      0.00        0
#> "withCallingHandlers"       0.04       100      0.00        0
#> "withVisible"               0.04       100      0.00        0
#> 
#> $sample.interval
#> [1] 0.02
#> 
#> $sampling.time
#> [1] 0.04

当Rprof()设定line.profiling = TRUE时，将改为分析代码行，而不是函数，这有利于分析复杂的计算过程。

高级点，使用RStudio提供的profvis()，它是一个交互式可视化的工具。

install.packages("profvis")

提供代码性能

使用内置函数
使用向量化
使用字节码编译器，再写完函数后用cmpfun()函数编译它，然后再使用
使用微软的多线程计算MRO
使用并行计算
使用Rcpp

前两点不细讲了。

使用字节码编译器

如果不是R内置的函数，一般通过字节码进行编译后再使用会提升速度，这里说下流程：

#载入包
library(compiler)
#调用函数进行编译
# 假设我们写了一个自定义的sum函数
sum_cmp = cmpfun(sum)
# 然后再使用sum_cmp函数

使用并行计算

由于不同的操作系统有不同的线程和线程模型的实现，因此某些功能在Linux和MacOS系统与Windows系统不同。

我们先创建一个伪任务，以方便进行并行计算：

simulate = function(x) {
    s = 0
    for (i in x){
        while(i >= 1){
            s = i * (i - 1)
            i = i - 1
        }
    }
    s
}

system.time(simulate(1:11000))
#>    user  system elapsed 
#>   4.443   0.028   4.489

在Windows下使用

Windows下需要先创建多个R会话的本地集群：

library(parallel)
cl = makeCluster(detectCores())

detectCores()返回我们计算机有的CPU核数目，然后调用parLapply()函数，它是并行版本的lapply()函数：

system.time(parLapply(cl, 1:11000, simulate))
#>    user  system elapsed 
#>   0.007   0.001   1.531

当我们不需要这个集群时，调用stopCluster()终止刚才创建的R会话。

注意，当我们使用这种方式进行并行计算时，创建了新的R会话，因此环境里没有用户定义的变量。如果有我们自己的数据，需要提前导入。

使用clusterExport()可以将数据导入到集群，使用clusterEvalQ()可以在每个集群节点上计算表达式。

或者使用clusterCall()和<<-在每个节点创建全局变量，而<-创建局部变量。

stopCluster()
#> Error in defaultCluster(cl): no cluster 'cl' supplied and none is registered

Linux和MacOS使用并行计算

在这两个系统上使用并行计算要比Windows容易的多，mclapply()可直接将当期R会话分配到多个会话中，保留所有内容，并为每个子会话安排任务，并行运行：

system.time( mclapply(1:1100, simulate, mc.cores = detectCores()))
#>    user  system elapsed 
#>   0.106   0.089   0.053

我们无需导出变量到集群中，因为在每个分配的进程中，它们是可以直接使用的。

我们还可以用非常灵活的方式创建并行作业，例如创建一个生成10个随机数的作业：

job1 = mcparallel(rnorm(10), "job1")

然后调用mccollect()函数收集作业结果：

mccollect(job1)
#> $`9584`
#>  [1]  0.479  0.558 -0.953  1.489  2.162  0.457 -1.489  3.329  1.414  0.432

还可以通过编程创建并运行多项作业，例如我们创建8项作业，每个随机休眠一段时间：

jobs = lapply(1:8, function(i){
    mcparallel({
        t = rbinom(1, 5, 0.6)
        Sys.sleep(t)
        t
    }, paste0("job", i))
    
})
system.time(res <- mccollect(jobs))
#>    user  system elapsed 
#>   0.010   0.016   4.926

使用Rcpp

并行计算只有在每次迭代都是独立的情况下才可行，这样最终结果才不会依赖运行顺序，然后并非所有的任务都像这样理想。另一种让算法更快的方式是使用Rcpp（http://www.rcpp.org）。

Rcpp是一个扩展包，它使我们能够利用R和C++的无缝整合来编写C++代码。使用Rcpp可以编写C++代码，并且代码中还可以调用R函数，利用R数据结构的优势。

使用Rcpp，先要使用正确的工具链，在Windows系统下使用Rtools，在Linux和MacOS下也需要安装正确的C/C++工具链。然后我们就安装Rcpp包：

install.packages("Rcpp")

先在Rcpp目录下创建一个C++源文件：

#include <Rcpp.h>
usingnamespace Rcpp;
// [[Rcpp::export]]
NumericVector timesTwo (NumericVector x){
    return x * 2;
}

如果你不熟悉C++语法，可以通过http://www.learncpp.com学习最简单的部分。

C++是强类型语言，需要指定函数参数的类型和函数返回的类型。使用[[Rcpp::export]]注释的函数会被Rcpp捕获，当在RStudio中执行一个脚本文件，或者直接使用Rcpp::sourceCpp，这些C++函数将被自动编译并移植到R的工作环境中。

Rcpp::sourceCpp("../../R/Rcpp/rcpp-demo.cpp")
#> Warning in normalizePath(dirname(file)): path[1]="../../R/Rcpp": No such file or
#> directory
#> Warning in normalizePath(file, winslash = "/"): path[1]="../../R/Rcpp/rcpp-
#> demo.cpp": No such file or directory
#> Error in Rcpp::sourceCpp("../../R/Rcpp/rcpp-demo.cpp"): file not found: '../../R/Rcpp/rcpp-demo.cpp'
timesTwo
#> Error in eval(expr, envir, enclos): object 'timesTwo' not found

Rcpp官网：http://www.rcpp.org/ 一些相关技术与包：

OpenMP
RcppParallel