r - 为什么 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)
}

以下是我在 Dell 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 的一部分，我编写的一个 Matlab 函数包，用于生成和测试各种随机排列生成器。这些注释可以在这里单独查看:Notes on RPGlab .

原始的 Durstenfeld Algol 程序是 here

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 脚本编译为字节码。然而，它不会将其转换为汇编代码，但效果很好。向他致敬!

...但是与 MATLAB 编译器相比，R 编译器投入的精力要少得多，因此结果更慢:

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

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

如您所见，通过将 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 的 randperm 快了 3 倍!

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