高性能R程序

2015-04-28 | 振导社会 | 程序设计高性能计算 R

程序性能剖析

确定程序运行时间

利用system.time：

system.time(for (i in 1:50) mad(stats::runif(500)))

利用proc.time：

ptm <- proc.time()
for (i in 1:50) mad(stats::runif(500))
proc.time() - ptm

性能监控的函数Rprof¹

可视化性能监控lineprof²

OpenBLAS：加速矩阵运算³

直接通过brew安装支持OpenBLAS的r：

$ brew install openblas --build-from-source
$ brew install r --with-openblas

# 安装可能遇到的问题：
# curl: (7) Failed to connect to rcompletion.googlecode.com port 443: Operation timed out
# Error: Failed to download resource "r--completion"
# Download failed: https://rcompletion.googlecode.com/svn-history/r31/trunk/bash_completion/R

也可以直接下载Revolution R Open (RRO)，默认就支持OpenBLAS。

如果程序已经是多线程，可能会和OpenBLAS发生冲突，可以在环境变量中设置OpenBLAS为单线程：

export OPENBLAS_NUM_THREADS=1

OpenBLAS提升效果：

## 使用了OpenBLAS：
x <- matrix(1:(6000 * 6000), 6000, 6000)
system.time(tmp <- x %*% x)
#   user  system elapsed 
# 13.321   0.323   7.315 

## 没有使用OpenBLAS：
x <- matrix(1:(6000 * 6000), 6000, 6000)
system.time(tmp <- x %*% x)
#    user  system elapsed 
# 206.588   2.216 214.333

parallel：并行计算包⁴⁵

parallel包是从snow包和multicore包合并继承而来，包含了很多非常好用的函数。multicore只能在支持fork的操作系统使用，只能用于单台计算机。snow可以用在Unix系列、Windows或者二者混合的集群上。在单处理器单核上使用multicore和snow没效果。

parallel包可以通过PVM（rpvm包）、MPI（Rmpi包）、NetWorkSpaces（nws包）和raw sockets（如果以上3种都不能使用）平台进行分布计算，支持cluster和多核个人/服务器计算机。原则上，parallel可以通过线程（thread）或轻量级进程（lightweight process）实现并行，但是目前都是依赖于进程（process），实现并行有三种方式：

通过system("Rscript")或类似的方式启动进程。安全机制可能会阻止进程间通过socket通信。按照snow的方式，通过socket监听来自主进程命令的进程池称为节点集群。
通过fork系统调用。fork出的进程副本会共享主进程的内存页，直到其内容发生改变，因此forking方式速度很快。fork的方式最早被multicore采用。由于进程的共享机制，也会共享GUI元素，这回导致havoc⁶。进程间可以通过管道和socket方式通信。
通过系统级机制向其它成员分发任务。snow包利用Rmpi包使用MPI（message passing interface）。这种情况下，通讯过载会增加计算时间，常用于高速内连的网络。在这种工作模式下，CRAN还提供了GridR和Rsge包。

doit <- function(x)(x)^2 + 2*x
system.time(res <- lapply(1:5000000,  doit))

 #   user  system elapsed 
 # 24.624   0.224  25.049 

library(parallel)
cl <- makeCluster(getOption("cl.cores", 4)) # use 4 cores
system.time(res <- parLapply(cl, 1:5000000,  doit))
stopCluster(cl) 

  #  user  system elapsed 
  # 2.405   0.258  10.444 

mc <- getOption("mc.cores", 4)
system.time(res <- mclapply(1:5000000,  doit, mc.cores = mc))

  #  user  system elapsed 
  # 6.023   1.632   5.300

注意：

需要先确定系统处理器核心数目，通常可用detectCores(logical = F)；
注意函数的调用方式是否为Rscript，该方式会复制对象，内存占用大，处理大数据时要当心。

foreach：并行计算包[1]

foreach包是revolution analytics公司贡献给R开源社区的一个包，它能使R中的并行计算更为方便。

doParallel包是foreach包并行计算的后端，它提供了并行执行foreach循环的机制。foreach必须采用doParallel这样的包才能实现并行计算。用户在使用时必须注册并行计算后端，否则即使用了%dopar%程序也串行执行。doParallel包起着foreach包和parallel包之间接口的作用。默认情况，doParallel包在Unix系列操作系统使用multicore功能，在Windows系统使用snow功能。[2]

# snow-like
library(doParallel)
cl <- makeCluster(2)
registerDoParallel(cl)
foreach(i=1:3) %dopar% sqrt(i)

# multicore-like
library(doParallel)
registerDoParallel(cores=2)
foreach(i=1:3) %dopar% sqrt(i)

## 该环境的后续程序都按multicore模式进行。

并行的boostrap：

## 已经注册了并行方式，不需要再注册……

x <- iris[which(iris[,5] != "setosa"), c(1,5)]
trials <- 10000

# 并行方案
ptime <- system.time({
    r <- foreach(icount(trials), .combine=cbind) %dopar% {
      ind <- sample(100, 100, replace=TRUE)
      result1 <- glm(x[ind,2]~x[ind,1], family=binomial(logit))
      coefficients(result1)
  }
})[3] 

# 串行方案：
stime <- system.time({
    r <- foreach(icount(trials), .combine=cbind) %do% {
      ind <- sample(100, 100, replace=TRUE)
      result1 <- glm(x[ind,2]~x[ind,1], family=binomial(logit))
      coefficients(result1)
  }
})[3] 

c(ptime, stime)

# elapsed elapsed 
#  26.372  37.126

memoise：本地缓存包⁷

memoise是一个简单的缓存包，主要用来减少重复计算，从而提升CPU性能。当你用相同的参数执行计算的时候，你会得到之前计算过的结果，而不是重算一遍。缓存技术对于有并发访问的应用来说，是性价比最高的性能提升方案。memoise包只有2个函数：forget重置缓存函数，memoize定义缓存函数。

#定义缓存函数
fun <- memoise(function(x) { Sys.sleep(1); runif(1) })

#第一次执行fun函数
system.time(print(fun()))
# [1] 0.4342335
#    user  system elapsed 
#   0.002   0.002   1.004 

#第二次执行fun函数 
system.time(print(fun()))
# [1] 0.4342335
#    user  system elapsed 
#   0.001   0.000   0.000 

#重置缓存函数
forget(fun)

#第三次执行fun函数 
system.time(print(fun()))
# [1] 0.786522
#    user  system elapsed 
#   0.003   0.003   1.002

compiler：编译功能包⁸

执行函数之前，把它编译成二进制程序。

library(compiler)

myFunction<-function() {for(i in 1:1e7) {1*(1+1)}}
myCompiledFunction <- cmpfun(myFunction) # 编译函数

system.time(myFunction())
  #  user  system elapsed 
  # 3.448   0.024   3.486 

system.time(myCompiledFunction())
  #  user  system elapsed 
  # 0.611   0.017   0.637

Rcpp：R中融合C++[3]

library(Rcpp)

cppFunction(
    'int fib_cpp_0(int n){
       if(n == 1 || n == 2) return 1;
       return(fib_cpp_0(n - 1) + fib_cpp_0( n - 2));
   }'
   )

fib_r <- function(n){
    if(n == 1 || n == 2) return(1)
        return(fib_r(n - 1) + fib_r(n - 2))
}

system.time(fib_cpp_0(30))

  #  user  system elapsed 
  # 0.002   0.000   0.002 

system.time(fib_r(30))

  #  user  system elapsed 
  # 1.697   0.021   1.739

Rcpp简化了在R中集成C++代码，它将各种R对象映射为特定的C++类，使得C++和R之间的对象管理变得简单、灵活，并提供了对STL等的广泛支持。C++代码可以被编译、链接并动态加载，或者通过包加载。

Rcpp包提供了在C++层次无缝访问、扩展和修改R对象的API。R的API基于SEXP上的函数与宏操作，SEXP是R对象的内部表示。这些API的关键功能包括：C++类对R对象的轻量级封装、自动垃圾回收策略、代码内连、R与C++的数据交换，以及错误处理。

Rcpp包API的两个典型应用场景：

用C++代码替代R代码以提升程序性能；
方便调用其它库提供的函数。

以下代码是采用Rcpp计算卷积：

#include <Rcpp.h>

RcppExport SEXP convolve3cpp(SEXP a, SEXP b) {
    Rcpp::NumericVector xa(a);
    Rcpp::NumericVector xb(b);
    int n_xa = xa.size(), n_xb = xb.size();
    int nab = n_xa + n_xb - 1;
    Rcpp::NumericVector xab(nab);
    for (int i = 0; i < n_xa; i++)
        for (int j = 0; j < n_xb; j++)
            xab[i + j] += xa[i] * xb[j];
    return xab; 
}

以上程序展示了使用Rcpp的几个重要方法：

使用Rcpp的API只需要一个头文件Rcpp.h；
RcppExport是方便从C调用C++的宏；
两个SEXP类型的输入变量，输出变量类型通过R的API的.Call()定义；
Rcpp将两个输入变量转换成了C++的向量类型；
通过成员函数size()查看对象大小，通过[]索引向量元素；
内存管理仍然由R完成；
返回值自动实现从NumericVector到SEXP的转换。

Rcpp sugar能在C++中使用类似R的语法。它不仅提供了漂亮的语法，而且使程序运行更高效。

其它技术

提升读入数据效率

read.table()、read.csv()适合读取小规模数据框，有效地读取大数值矩阵要使用更为底层的read.delim()，甚至scan()函数。⁹

在读取大型数据时，设定comment.char=""，以读取目标的原子向量类型（逻辑型，整型，数值型，复数型，字符型等）。事先设置好每列的colClasses，给定需要读入的行数nrows（适当地高估一点比不设置这个参数还要快）等措施会部分地提高效率。如果需要试探数据，可以把nrows设置为10或者更小，这样就可以只读取并查看数据的前几行。⁹¹⁰

参考资料

[1]R. Analytics and S. Weston, doParallel: Foreach parallel adaptor for the parallel package. 2014. [Online]
[2]S. Weston and R. Calaway, “Getting Started with doParallel and foreach,” 2014.
[3]D. Eddelbuettel and R. François, “Rcpp: Seamless R and C++ Integration,” Journal of Statistical Software, vol. 40, no. 8, pp. 1–18, 2011. [Online]

R语言性能监控工具Rprof ↩
R语言性能可视化lineprof ↩
用 OpenBLAS 加速 R 的矩阵运算 ↩
并行化你的运算-初识parallel包 ↩
R使用parallel包并行计算 ↩
Some precautions are taken on Mac OS X: for example the event loops for R.app and the quartz device are inhibited in the child. This information is available at C level in the Rboolean variable R_isForkedChild. ↩
R语言本地缓存memoise ↩
简单加速R ↩
R数据读取的几点tips ↩ ↩²
R语言读取数据-性能优化 ↩

打赏作者

2016-10-24 » NNML（03）：BP 学习
2016-10-16 » NNML（02）：感知器学习
2016-10-09 » 鲁棒及自适应控制（2）：模型
2016-10-09 » NNML（01）：引言
2016-09-24 » 家用监控设备用于电网的可行性分析
2016-09-19 » 无人机电缆隧道巡检可行性调研报告
2016-09-18 » 鲁棒及自适应控制（1）：概论
2016-09-13 » Gradient-Based Learning Applied to Document Recognition

上一篇：R包开发下一篇：arules：频繁项集与关联规则的挖掘

程序性能剖析

确定程序运行时间

性能监控的函数Rprof1

可视化性能监控lineprof2

OpenBLAS：加速矩阵运算3

parallel：并行计算包45