为什么R在这个随机排列函数上变慢？

Question

我是R(Revolution Analytics R)的新手,并且已经将一些Matlab函数转换为R.

问题:为什么GRPdur(n)函数这么慢？

GRPdur = function(n){
#
# Durstenfeld's Permute algorithm, CACM 1964
# generates a random permutation of {1,2,...n}
#
p=1:n                           # start with identity p
for (k in seq(n,2,-1)){    
    r    = 1+floor(runif(1)*k); # random integer between 1 and k
    tmp  = p[k];
    p[k] = p[r];                #  Swap(p(r),p(k)).
    p[r] = tmp;                  
} 
return(p)
}

以下是戴尔Precision 690,2xQuadcore Xeon 5345 @ 2.33 GHz,Windows 7 64位的内容:

> system.time(GRPdur(10^6))
   user  system elapsed 
  15.30    0.00   15.32 
> system.time(sample(10^6))
   user  system elapsed 
   0.03    0.00    0.03 

这是我在Matlab 2011b中得到的

>> tic;p = GRPdur(10^6);disp(toc)
    0.1364   

 tic;p = randperm(10^6);disp(toc)
    0.1116

这是我在Matlab 2008a中得到的

>> tic;p=GRPdur(10^6);toc
Elapsed time is 0.124169 seconds.
>> tic;p=randperm(10^6);toc
Elapsed time is 0.211372 seconds.
>> 

---

链接:GRPdur是RPGlab的一部分,RPGlab是我编写的一个Matlab函数包,用于生成和测试各种随机排列生成器.这些注释可以在这里单独查看:关于RPGlab的注释.

最初的Durstenfeld Algol计划就在这里

Answer 1

Matlab和S(后来的R)都是围绕FORTRAN函数的瘦包装开始做数学的.

在S/R中,for循环"总是"很慢,但这一切都很好,因为通常存在表达问题的矢量化方式.此外,R在Fortran或C中有数千个函数可以快速执行更高级别的操作.例如,该sample功能完全与你的for循环相同 - 但更快.

那么为什么MATLAB在执行脚本化的for循环方面要好得多呢？有两个简单的原因:资源和优先事项.

制作MATLAB的MathWorks是一家拥有约2000名员工的大公司.他们多年前决定优先考虑提高脚本的性能.他们聘请了许多编译专家,并花了数年时间开发一个即时编译器(JIT),它接受脚本代码并将其转换为汇编代码.他们也做得很好.感谢他们!

R是开源的,R核心团队在业余时间努力改进R. R core的Luke Tierney努力工作并为R开发了一个编译器包,用于将R脚本编译为字节代码.但它并没有把它变成汇编程序代码,但效果很好.感谢他!

...但是,R编译器与MATLAB编译器的工作量要少得多,因此结果较慢:

system.time(GRPdur(10^6)) # 9.50 secs

# Compile the function...
f <- compiler::cmpfun(GRPdur)
system.time(f(10^6)) # 3.69 secs

如您所见,通过将for循环编译为字节代码,for循环速度提高了3倍.另一个区别是默认情况下不启用R JIT编译器,因为它在MATLAB中.

更新只是为了记录,一个稍微更优化的R版本(基于Knuth的算法),其中随机生成已被矢量化为@joran建议:

f <- function(n) {
  p <- integer(n)
  p[1] <- 1L
  rv <- runif(n, 1, 1:n) # random integer between 1 and k
  for (k in 2:n) {    
    r <- rv[k]
    p[k] = p[r]         #  Swap(p(r),p(k)).
    p[r] = k
  }
  p
}
g <- compiler::cmpfun(f)
system.time(f(1e6)) # 4.84
system.time(g(1e6)) # 0.98

# Compare to Joran's version:
system.time(GRPdur1(10^6)) # 6.43
system.time(GRPdur2(10^6)) # 1.66

...仍然比MATLAB慢一点.但同样,只需使用sample或sample.int显然胜过MATLAB的randperm3倍!

system.time(sample.int(10^6)) # 0.03

Answer 2

因为你在R-skin中编写了一个c程序

n = 10^6L
p = 1:n
system.time( sample(p,n))
0.03    0.00    0.03