50-R茶话会 （十：R编程效率提升指北）
参考：https://www.math.pku.edu.cn/teachers/lidf/docs/Rbook/html/_Rbook/prog-prof.html
1. R 的运行效率
R是解释型语言，在执行单个运算时， 效率与编译代码相近；在执行迭代循环时， 效率较低， 与编译代码的速度可能相差几十倍。在循环中对变量进行修改尤其低效， 因为R在修改某些数据类型的子集时会复制整个数据对象。（这个在前面提到过） R以向量、矩阵为基础运算单元， 在进行向量、矩阵运算时效率很高， 应尽量采用向量化编程。
为了提高R程序的运行效率，需要尽可能利用R的向量化特点， 尽可能使用已有的高效函数， 还可以把运行速度瓶颈部分改用C++、FORTRAN等编译语言实现， 可以用R的profiler工具查找运行瓶颈。对于大量数据的长时间计算， 可以借助于现代的并行计算工具。
对已有的程序，仅在运行速度不满意时才需要进行改进，否则没必要花费宝贵的时间用来节省几秒钟的计算机运行时间。
另外，在改进已有程序的效率时， 第一要注意的就是不要把原来的正确算法改成一个速度更快但是结果错误的算法。这个问题可以通过建立试验套装， 用原算法与新算法同时试验看结果是否一致来避免。多种解决方案的正确性都可以这样保证，也可以比较多种解决方案的效率。（改错还不如不改）
如果要实现一个比较单纯的不需要利用R已有功能的算法， 发现用R计算速度很慢的时候， 也可以考虑先用Julia语言实现。Julia语言设计比R更先进，运算速度快得多，运算速度经常能与编译代码相比， 缺点是刚刚诞生几年的时间，可用的软件包还比较少。（之前学python 的时候也有人推荐过Julia，可其社区还是不够成熟）
2. 提高R 运行效率的几个策略
2.1 尽量使用已有函数及向量化
在计算总和、元素乘积或者每个向量元素的函数变换时， 应使用相应的函数，如sum, prod, sqrt, log等。
其中有的内建函数， 如sum, prod, cumsum, cumprod, mean, var, sd等。这些函数以编译程序的速度运行， 不存在效率损失。
而有的函数则是向量化的， 可以直接对输入向量的每个元素进行变换。这个我们先前已经提到过了。
参见：www.cnblogs.com/hyb221512/p/10783242.html 这里有一个prod 的例子。
prod 用来计算矩阵或向量元素的乘积。如果是矩阵，还可以指定参数dim 选定行或列：
> test
[1] 1.185311 1.318823 1.379850 1.362257 1.174529 1.373174 1.098781 1.486921
[9] 1.026443 1.646543 1.249221 1.430974 1.578645 1.488934 1.393388 1.674874
[17] 1.503912 1.842919 1.095980 1.890939
> prod(test)
[1] 737.0558
cumprod 则可以得到累积的乘积，比如向量中的a,b,c 三个数，prod 是直接获得a*b*c 的结果，而cumprod 则是分别得到a, a*b, a*b*c 的结果。
可以通过下面的例子来比较一下速度：
f1 <- function(){
ny <- 365
x <- numeric(ny)
for(n in 1:ny){
s <- 1
for(j in 0:(n-1)){
s <- s * (365-j)/365
}
x[n] <- 1 - s
}
x
}
f2 <- function(){
ny <- 365
x <- numeric(ny)
for(n in 1:ny){
x[n] <- 1 - prod((365:(365-n+1))/365)
}
x
}
f3 <- function(){
ny <- 365
x <- 1 - cumprod((ny:1)/ny)
x
}
microbenchmark::microbenchmark(
f1(),
f2(),
f3()
)
## Unit: microseconds
## expr min lq mean median uq max neval
## f1() 2534.807 2577.730 2679.48244 2615.256 2712.9270 3408.187 100
## f2() 323.855 333.365 414.81692 344.160 369.3485 5868.456 100
## f3() 1.028 1.542 2.52415 2.056 2.5700 25.189 100
2.2 避免显式循环
显式循环是R运行速度较慢的部分， 有循环的程序也比较冗长， 与R的向量化简洁风格不太匹配。
所谓显式循环，也就是在代码中不直接调用for 或while 这些循环函数。
使用apply 族函数
apply 族包括apply, sapply, lapply, tapply, vapply，比如可以使用它们对向量、数据框、列表、矩阵等多种类型数据进行处理。
expr min lq mean median
for (i in 1:10) print(mean(my_mat[i, ])) 2201.132 2471.456 3349.322 2650.3690
apply(my_mat, 1, mean) 58.058 71.520 120.866 115.9945
uq max neval
3828.0745 9859.091 100
137.0115 459.627 100
可见apply 函数快了15倍之多。
replicate() 函数
其用法比for 和apply 都要简单，类似于for()循环但是没有计数变量。其用法为replicate(n, fun)，n 表示重复的次数，fun 表示方法。
该函数比较常用在随机模拟上（获得多个模拟结果），比如生成随机数：
> replicate(10, {x = runif(10,1,100); mean(x)})
[1] 43.46471 56.03731 67.77589 48.84824 50.86738 41.22307 32.98558 49.75672
[9] 42.35581 60.92928
2.3 事先分配合理长度数据结构
如果事先清楚需要创建的数据结构和其长度，可以事先声明，这样的程序结构更清晰， 效率更高， 而且循环次数越多， 比x <- c(x, ...)这样的做法的优势越大。
ps:如果是列表的话，可以使用vector(n, mode = "list") 创建长度为n的空列表:
set.seed(101)
system.time({
M <- 1E5
x <- numeric(M)
for(i in seq(M)){
x[[i]] <- diff(range(runif(10)))
}
mean(x)
})
## 用户 系统 流逝
## 1.59 0.01 1.61
2.4 避免制作副本
在循环内修改数据子集，例如数据框子集， 可能会先制作副本再修改， 这当然会损失很多效率。R 3.1.0版本以后列表元素在修改时不制作副本， 但数据框还会制作副本。
因此，对于重复较多且大的数据框对象，我们可以先将其用列表处理，最后再转换成数据框：
set.seed(101)
m <- 2E4; n <- 100
x <- as.data.frame(matrix(
runif(n*M), nrow=n, ncol=m))
system.time({
for(j in seq(m)){
x[[j]] <- x[[j]] + 1
}
})
## 用户 系统 流逝
## 14.92 0.02 14.94
set.seed(101)
m <- 2E4; n <- 100
x <- replicate(m,
runif(n),
simplify=FALSE)
system.time({
for(j in seq(m)){
x[[j]] <- x[[j]] + 1
}
})
## 用户 系统 流逝
## 0.01 0.01 0.03
x <- as.data.frame(x)
replicate() 函数中用simplify=FALSE 使结果总是返回列表。要注意的是， 上面第二个程序中的as.data.frame(x)也是效率较差的。将数据保存在列表中比保存在数据框中访问效率高， 数据框提供的功能更丰富。
2.5 R 的并行运算
R 提供了parallel 及snowfall 进行apply 族函数的并行运算，foreach 提供了 for 函数的并行。另外，WIN OS 下还提供了特别的R 版本，可以实现更加方便的R 的并行运算。不过在使用R 的并行时需要注意合理分配线程及内存释放的管理。
2.7 通过Rprof 进行代码性能分析
参见：http://blog.fens.me/r-perform-rprof-profr/
要改善运行速度， 首先要找到运行的瓶颈， 这可以用专门的性能分析（profiling）功能实现。R软件中的Rprof()函数可以执行性能分析的数据收集工作， 收集到的性能数据用summaryRprof()函数可以显示运行最慢的函数。如果使用RStudio软件，可以用Profile菜单执行性能数据收集与分析， 可以在图形界面中显示程序中哪些部分运行花费时间最多。
我们可以直接使用Rstudio 的profile 工具，对选定的代码进行分析：
可这个结果有点难看懂啊，比如：
猜测是分别比较上下两段代码的空间和时间的差异吧。
仔细学习函数及参数
- system.time 这个函数比较直接，可以直接以秒/s 为单位，显示代码的运行时间：比如：
> system.time(for (i in 1:10000){
+ my_df[i,4] <- my_df[i,1] * my_df[i,3]
+ })
用户 系统 流逝
1.050 0.372 1.444
> my_df <- data.frame(a = 1:10000, b = sample(LETTERS, 10000, replace = T),
+ c = rnorm(10000, 1))
> system.time(my_df[,4] <- my_df[,1] * my_df[,3])
用户 系统 流逝
0.000 0.000 0.001
看来还是向量化编程完胜啊。
- Rprof() 这个函数用起来还是挺奇怪的：
file<-"test.out" # 指定一个输出文件夹
Rprof(file) # 指定输出文件夹
### 下面部分就可以执行需要测试的代码了
fun1()
###
Rprof(NULL) # 执行完毕，关闭Rprof
输出的文件长这个样子：
这个文件是没有任何意义哒，我们需要使用其他命令来分析它。ps：prof 指令运行完毕后不会自动停止（不过思考一下，也确实，因为如果不关闭执行，Rprof 还是要继续等待输入命令的嘛）
可以看一下我的一些代码的运行结果：使用summaryRprof 命令查看。
> summaryRprof(file)
$by.self
self.time self.pct total.time total.pct
"[<-.data.frame" 1.44 92.31 1.48 94.87
"paste" 0.04 2.56 0.04 2.56
"[.data.frame" 0.02 1.28 0.04 2.56
"%in%" 0.02 1.28 0.02 1.28
"anyNA" 0.02 1.28 0.02 1.28
"max" 0.02 1.28 0.02 1.28
$by.total
total.time total.pct self.time self.pct
"[<-.data.frame" 1.48 94.87 1.44 92.31
"[<-" 1.48 94.87 0.00 0.00
"paste" 0.04 2.56 0.04 2.56
"[.data.frame" 0.04 2.56 0.02 1.28
"[" 0.04 2.56 0.00 0.00
"head" 0.04 2.56 0.00 0.00
"%in%" 0.02 1.28 0.02 1.28
"anyNA" 0.02 1.28 0.02 1.28
"max" 0.02 1.28 0.02 1.28
$sample.interval
[1] 0.02
$sampling.time
[1] 1.56
by.self: 当前函数的耗时情况 by.total:整体函数调用的耗时情况
这么看还是不够直观，尝试使用profr 包画图（直接install 安装一下）。
该包中的函数parse_rprof 同样可以对Rprof 后的文件进行分析：
> profr::parse_rprof(file)
level g_id t_id f start end n leaf time source
1 1 1 1 [<- 0.00 0.12 1 FALSE 0.12 base
2 1 2 1 [ 0.12 0.14 1 FALSE 0.02 base
3 1 3 1 [<- 0.14 0.88 1 FALSE 0.74 base
4 1 4 1 [ 0.88 0.90 1 FALSE 0.02 base
5 1 5 1 [<- 0.90 1.52 1 FALSE 0.62 base
6 1 6 1 head 1.52 1.56 1 FALSE 0.04 utils
7 2 1 1 [<-.data.frame 0.00 0.12 1 TRUE 0.12 base
8 2 2 1 [.data.frame 0.12 0.14 1 TRUE 0.02 base
9 2 3 1 [<-.data.frame 0.14 0.88 1 TRUE 0.74 base
10 2 4 1 [.data.frame 0.88 0.90 1 FALSE 0.02 base
11 2 5 1 [<-.data.frame 0.90 1.52 1 TRUE 0.62 base
12 2 6 1 paste 1.52 1.56 1 TRUE 0.04 base
13 3 1 1 anyNA 0.76 0.78 1 TRUE 0.02 base
14 3 2 2 %in% 0.88 0.90 1 TRUE 0.02 base
15 3 3 3 max 1.50 1.52 1 TRUE 0.02 base
我们可以直接对这个结果画图：
ggplot2::ggplot(profr::parse_rprof(file))
可还是觉得基础包的绘图直观一些：
纵坐标显示层级关系，横坐标显示其执行时间。
Rprof 命令行使用
我们还可以直接在终端中调用R：
R CMD BATCH file # 范例
## 参数含义
-h, –help: 打印帮助信息
-v, –version: 打印版本信息
–lines: 打印显示多行
–total: 只显示总数
–self: 只显示自己
–linesonly: 只显示单行(配合–lines使用)
–min%total=: 显示total的不低于X的百分比
–min%self=: 显示self的不低于X的百分比
## 例子
R CMD Rprof --total --min%total=50 fun1_rprof.out
Each sample represents 0.02 seconds.
Total run time: 0.14 seconds.
Total seconds: time spent in function and callees.
Self seconds: time spent in function alone.
% total % self
total seconds self seconds name
100.0 0.14 0.0 0.00 "fun1"
71.4 0.10 0.0 0.00 "ddply"
57.1 0.08 0.0 0.00 "ldply"
其他的命令运行检查工具microbenchmark
直接install 安装：
f1 <- function(x){
n <- length(x)
mhat <- median(x)
s <- 0.0
for(i in 1:n){
s <- s + abs(x[i] - mhat)
}
s <- s/n
return(s)
}
# 循环计算
f2 <- function(x) mean( abs(x - median(x)) ) # 向量化计算
x <- runif(10000) # 获得随机数，类似rnrom
microbenchmark::microbenchmark(
f1(x),
f2(x)
)
## Unit: microseconds
## expr min lq mean median uq max neval cld
## f1(x) 2301.9 2386.75 3064.035 2466.85 3273 11499.9 100 b
## f2(x) 400.6 421.25 548.897 436.55 559 4940.8 100 a
该结果的单位为毫秒，可见，使用向量化计算比循环平均快了整整六倍。
腾讯云开发者
