也谈提高R语言的运算效率

用过底层语言做计算的人转入R语言的时候一般都会觉得R语言的运算太慢,这是一个常见的对R的误解或者对R的设计的不理解。在二三十年前Chambers等一群人在贝尔实验室设计S语言之前,统计计算主要是通过那些底层语言实现的,典型的如Fortran。当时有一个基于Fortran的统计分析库,用它的麻烦就在于无论做什么样的数据分析,都涉及到一大摞繁琐的底层代码,这让数据分析变得很没劲,统计学家有统计学家的事情,天天陷在计算机程序代码中也不是个办法,要摆脱底层语言,那就只能设计高层语言了。有所得必有所失,众所周知,高层语言一般来说比底层语言低效,但对用户来说更友好。举个简单的例子,用C语言计算均值时,我们得对一个向量(数组)从头循环到尾把每个值累加起来,得到总和,然后除以向量的长度,而均值在统计当中简直是再家常便饭不过了,要是每次都要来这么个循环,大家也都甭干活儿了,天天写循环好了。

前两天COS论坛上有个帖子提到“R语言的执行效率问题”,大意如下:

有120000行数据,用perl写成12万行R命令做t.test,然后执行,大概几分钟就算完了。如果只用R语言,把所有数据先读入,然后用循环,每一行做t.test,花了几个小时,到现在还没有算完。这说明一个问题,在R里执行单行命令要比用循环快,涉及到循环的问题,最好写成单行命令执行。

我不清楚作者是如何写这“12万行”R命令的,本文假设是把t.test(dat[i, ]), i = 1, 2, ..., 120000写入一个文件,然后用source()执行之。面对这样一个问题,我们有若干种改进计算的可能性。首先我们看“硬”写入程序代码的方法:

## 为了使计算可重复,设定随机数种子
set.seed(123)
## 生成数据,随机数,120000行 x 100列矩阵
dat = matrix(rnorm(120000 * 100), 120000)
nr = nrow(dat)
nc = ncol(dat)
## 六种方法的P值向量
p1 = p2 = p3 = p4 = p5 = p6 = numeric(nr)

## via source()
f = file("test.t")
writeLines(sprintf("p1[%d] = t.test(dat[%d, ])$p.value",
    seq(nr), seq(nr)), f)
system.time({
    source(f)
})
#   user  system elapsed
#  95.36    0.19   95.86
close(f)
unlink("test.t")

1、向量式计算:apply()以及*apply()

当我们需要对矩阵的行或者列逐一计算时,apply()系列函数可能会提高效率。本例是对矩阵的行做t检验,那么可以这样计算:

## via apply()
system.time({
    p2 = apply(dat, 1, function(x) {
        t.test(x)$p.value
    })
})
#   user  system elapsed
#  63.12    0.06   63.50
identical(p1, p2)
# [1] TRUE

结果比第一种方法快了大约半分钟,而且计算结果完全一样。apply()本质上仍然是循环,但它在某些情况下比直接用显式循环要快:

## via for-loop
system.time({
    for (i in seq(nr)) p3[i] = t.test(dat[i, ])$p.value
})
#   user  system elapsed
#  69.88    0.03   70.05
identical(p2, p3)
# [1] TRUE

不过apply()系列函数在提高运算速度上优势并不会太明显,提倡用它的原因是它和统计中的矩阵运算相似,可以简化代码,相比起$\sum_{i=1}^n x_i/n$,我们可能更愿意看$\bar{x}$这样的表达式。很多R内置函数都是用底层语言写的,我们需要做的就是把一个对象作为整体传给函数去做计算,而不要自行把对象分解为一个个小部分计算,这个例子可能更能体现向量式计算的思想:

system.time(apply(dat, 1, mean))
#   user  system elapsed
#   5.28    0.04    5.25
system.time({
    x = numeric(nr)
    for (i in 1:nr) x[i] = mean(dat[i, ])
})
#   user  system elapsed
#   4.44    0.02    4.42
system.time(rowMeans(dat))
#   user  system elapsed
#   0.11    0.00    0.13

2、明确计算的目的

很多情况下,R函数返回的不仅仅是一个数字作为结果,而是会得到一系列诸如统计量、P值、各种系数等对象,在我们调用R函数之前如果能想清楚我们究竟需要什么,也许对计算的速度提升有帮助。比如本例中,也许我们仅需要12万个双边P值,其它数字我们都用不着,那么可以考虑“手工”计算P值:

## "hand" calculation in R
system.time({
    p4 = 2 * pt(apply(dat, 1, function(x, mu = 0) -abs((mean(x) -
        mu)/sqrt(var(x)/nc))), nc - 1)
})
#   user  system elapsed
#  15.97    0.07   16.08
identical(p3, p4)
# [1] TRUE

上面的计算更“纯净”,所以计算速度有了本质的提升,而且计算的结果和前面毫无差异。在做计算之前,人的脑子多思考一分钟,也许计算机的“脑子”会少转一个小时。

3、把四则运算交给底层语言

R是高层语言,把它拿来做简单的四则运算是很不划算的,而且容易导致程序低效。加加减减的活儿是C和Fortran等底层语言的强项,所以可以交给它们去做。以下我们用一段C程序来计算t统计量,然后用R CMD SHLIB将它编译为dll(Windows)或so(Linux)文件,并加载到R中,用.C()调用,最终用R函数pt()计算P值:

## "hand" calculation in C for t-statistic
writeLines("#include <math.h>

void calc_tstat(double *x, int *nr, int *nc, double *mu, double *tstat)
{
    int i, j;
    double sum = 0.0, sum2 = 0.0, mean, var;
    for (i = 0; i < *nr; i++) {
        for (j = 0; j < *nc; j++) {
            sum += x[i + j * *nr];
        }
        mean = sum / (double) *nc;
        sum = 0.0;
        for (j = 0; j < *nc; j++) {
            sum2 += (x[i + j * *nr] - mean) * (x[i + j * *nr] - mean);
        }
        var = sum2 / (double) (*nc - 1);
        sum2 = 0.0;
        tstat[i] = (mean - *mu) / sqrt(var / (*nc - 1));
    }
}", "calc_tstat.c")
system("R CMD SHLIB calc_tstat.c")
dyn.load(sprintf("calc_tstat%s", .Platform$dynlib.ext))
system.time({
    p5 = 2 * pt(-abs(.C("calc_tstat", as.double(dat), nr, nc,
        0, double(nrow(dat)))[[5]]), nc - 1)
})
#   user  system elapsed
#   0.86    0.06    0.92
dyn.unload(sprintf("calc_tstat%s", .Platform$dynlib.ext))

因为R可以用system()调用系统命令,所以整个过程全都可以用R完成,Windows用户需要安装Rtools并设置系统环境变量PATH才能使用R CMD SHLIB

更进一步,因为R自身的一些C程序也是可供用户的C程序调用的,所以我们可以把整个P值的计算过程全都扔进C代码中,一步完成:

## "hand" calculation in C calling Rmath.h
writeLines("#include <Rmath.h>
void calc_pvalue(double *x, int *nr, int *nc, double *mu, double *pval)
{
    int i, j;
    double sum = 0.0, sum2 = 0.0, mean, var;
    for (i = 0; i < *nr; i++) {
        for (j = 0; j < *nc; j++) {
            sum += x[i + j * *nr];
        }
        mean = sum / (double) *nc;
        sum = 0.0;
        for (j = 0; j < *nc; j++) {
            sum2 += (x[i + j * *nr] - mean) * (x[i + j * *nr] - mean);
        }
        var = sum2 / (double) (*nc - 1);
        sum2 = 0.0;
        pval[i] = 2 * pt(-fabs((mean - *mu) / sqrt(var / (*nc - 1))),
                      (double) (*nc - 1), 1, 0);
    }
}", "calc_pvalue.c")
system("R CMD SHLIB calc_pvalue.c")
dyn.load(sprintf("calc_pvalue%s", .Platform$dynlib.ext))
system.time({
    p6 = .C("calc_pvalue", as.double(dat), nr, nc, as.double(0),
        double(nrow(dat)))[[5]]
})
#   user  system elapsed
#   0.83    0.07    0.91
dyn.unload(sprintf("calc_pvalue%s", .Platform$dynlib.ext))

头文件Rmath.h的引入使得我们可以调用很多基于C程序的R函数,详情参考手册Writing R Extensions。通过C计算出来的P值和前面用R算的略有差异,下面画出p6 - p1 vs p1以及p6 - p5 vs p5的图:

P值的差异
P值的差异

导致差异的原因此处不细究,感兴趣的读者可以帮忙检查一下。

小结

  1. 若你熟悉底层语言,计算又不太复杂,那么可用底层语言写,然后用R调之;
  2. 否则把R对象当做整体去计算,能做x + 1就不要做for (i in length(x)) x[i] + 1
  3. 不要低估R core们的编程水平,他们已经做了很多工作让用户脱离底层编程

注:本文中的运算时间可能不可重复,这与计算机所处的状态有关,但大体来说,运算速度的快慢是可以肯定的。本文仅仅是关于统计计算的一个微小的例子,以后若有更好的例子,可能会更新本文;也欢迎各位提供更多示例。

关于谢益辉

RStudio码了个工,Iowa State University统计系博了个士。统计之都网站创办者;研究兴趣为统计图形及数据可视化,对统计模型方法的发展感兴趣但不喜欢纯粹抽象的数学理论,以直观、实用为学习标准;偏好以R语言为工具;Email:[email protected];个人主页:http://yihui.name

也谈提高R语言的运算效率》有18个想法

  1. 感谢师兄的这篇文章!也希望大家能把自己编程中提高效率的一些经验与大家进行分享。
    此外我记得R中有一个包叫做inline,它其实就是简化了上面的writeLines()和system()的操作,通过写临时文件和删临时文件并调用系统命令来编译C程序。


  2. 用C写 加加减减 的程序, 尤其是矩阵计算里要用的加加减减, 速度确实提高很多.

  3. 再有R 的核心还应该有改进的余地,我粗粗的看了一下, src/main/names.c 中

    attribute_hidden FUNTAB R_FunTab[] 这个结构在使用的时候,是通过

    for 循环进行比较 R_FunTab.name 的, 这个循环方式不经济,

    可以通过 hash R_FunTab.name 的方式,快速定位,这样就快了点点,

    不过这个改进还只是治标的小改动。

    1. 我觉得当R的影响力足够大的时候,肯定会有一批人专门来研究怎么提高核心代码的效率的,就像是S语言一代一代的更新。

      另:这位就是多次跟师兄提及的牟先生哦,呵呵。

      1. 当时我要提自己事情主要是担心只报了上海而不能参加北京活动(不规则的参与者的行动或给组织者造成了麻烦)。我在统计学领域还只是幼儿园级别,特别感谢各位能够提供机会让我能够学习。

      2. 哦,谢谢相告,我之前还搜了好半天邮件,人名和ID对不起来了。也谢谢lyxmoo的热心支持!

  4. 请问师兄为什么我直接运行第一段代码显示的是

    用户 系统 流逝 
       0    0    0 

    需要做什么改动吗?另外test.t是什么格式?

    1. That’s an arbitrary file name; you can name it as anything you want, e.g. ttest.txt. It is a pure text file.

      I’ve tested the code under both Windows and Linux, so it should not be a problem with the computation. I guess the problem is whether you have enough privilege to write a file to your working directory getwd(). (AFAIK, in Win7 or Vista there can be problems to access your C disk.)

  5. The version of R in my computer is REvolution R 1.3.0. Now I have changed to the newest R version 2.10.1 to test the same code. The problem dispeared! Maybe there is something in relation to the function source(). When using REvovlution, it create a file n my hard disk. However when sourcing it, there is nothing in computation, at the same time, the file size become to be 0 kb…
    PS:My OS is Windows XP SP3.

    1. 这样看来可能是REvolution对连接(?connection)的处理还有所欠缺,R里面很多函数都支持各种连接,通过文件名访问文件只是最基础的一种,其它方式包括socket、URL、gzip等。我猜的话,source('test.t')应该是可以运行的。

  6. 用R做像EM之类的处理通常很痛苦。问题就在循环、四则运算、拷贝等问题上面。因为R语言是变量和LISP语言一样,是一张表(当然也可以描述成树)。
    普通的一个
    A<-matrix(runif(100),nrow=10);

    B<-A

    for (i in 1:10) {for (j in 1:10) B[i,j]=A[i,j]}
    可是差得太多的。虽然,如果写成Fortran或者C++,那是一样的东西。问题出在R是解释语言,不是编译后的东西(Python的速度可能要例外)。
    可以想见的结果就是通常的代数运算会很慢。通常的速度问题出现在矩阵操作上。

    至于使用较为底层的C/C++,F/F99,Pascal等,问题在于代码有的时候非常难写。比如说高维矩阵求逆,或者稀疏矩阵求广义逆等,这些是都是非常麻烦的。要牵涉到很复杂的算法、精度控制等(比方说双double精度)。

    R提供了一个很好的接口。R本身也是通过调用Lapack/BLAS等的功能完成上述矩阵的。但是这一点是操作系统相关的。上述的库在Win平台上是很麻烦的,甚至很多功能不兼容。比如(c)BLAS至今有很多功能也不能和VS2005等Win上主流C/C++编译器兼容。所以,打开R的etc,我们发现实际上用的是GCC完成代码编译的(编译选项是 -O2, 会的同学可以试图改成 -O2 -march=[视CPU而定]。运算速度有时能提高不少。 )

    所以我的建议是,在进行这种复杂操作时:
    1. 底层语言和R语言混合编程,事实上,同样的问题存在于Octave/Matlab(R)上面;
    2. 重要项目用LAPACK等开源库完成代码,写成R的包;
    3. 如果能换平台,用Xuinx类系统,这样执行效率要高很多,也能突破内存限制。

  7. 在R2.10里面貌似循环速度大大加快,只比apply慢了300%,几分钟时间完全可以接受了。
    我机器,笔记本酷睿2 t7300 主频2.0GHz 内存ddr2 667 2G
    也搞了18万行的循环做个小测试

    > b=rep(0,1800)#这里当时打错了,不过R能自动扩展向量长度的嘛
    > a=matrix(rnorm(900000,10,1),nc=5)
    > system.time(for(i in 1:180000){b[i]=t.test(a[i,])})
    用户 系统 流逝
    331.36 2.49 339.83
    警告多于50个(用warnings()来显示第一个到第五十个)
    > system.time(b<-apply(a,1,t.test))
    用户 系统 流逝
    102.81 0.06 103.20

    现在看来已经没有几个小时那么夸张了嘛,多等两分钟问题不是很大吧。这时候程序可读性貌似更重要,因为自我感觉向量化会让语句看起来比较累,然后这样可以把精力放在解决具体问题上而不是如何搞向量化……我是菜鸟了,意见仅供参考。

    1. 补充一点,系统是xp sp3。 类似测试在fedora 12下也进行过,结果依然类似。

  8. 简单来说,R底层基本上都是C做的,调用的库也是久经考验的。
    程序运行的效率核心在于对语言提供的功能十分有足够的了解和体会,如果把base包中的内容认真仔细的通读下,再掌握常见运算的核心所在,相信就不会抱怨工具的问题——相对来说R对计算方面考虑的已经非常好了。

  9. 嗯 尽量避免通过R来对大量数据进行预处理 数据的一些预处理可以交给C或者Java 而用R只进行实验计算等

  10. 请问dependent iteration(原谅我中英夹杂,因为我不知道怎么正确翻译)有什么好的解决方法了吗?实在是不想写for循环啊

发表评论

电子邮件地址不会被公开。 必填项已用*标注