Tại sao R chậm trên chức năng hoán vị ngẫu nhiên này?

Question

Tôi mới vào R (Revolution Analytics R) và đã được dịch một số chức năng Matlab vào R.

Câu hỏi: Tại sao là hàm GRPdur (n) để làm chậm?

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) 
} 

Dưới đây là những gì tôi nhận được trên một Dell Precision 690, 2xQuadcore Xeon 5345 @ 2.33 GHz, Windows 7 64-bit:

> 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 

Dưới đây là những gì tôi nhận được trong Matlab 2011b

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

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

Dưới đây là những gì tôi nhận được trong 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. 
>> 

---

LIÊN KẾT: GRPdur là một phần của RPGlab, một gói chức năng Matlab mà tôi đã viết tạo và kiểm tra các trình tạo hoán vị ngẫu nhiên khác nhau. Các ghi chú có thể được xem riêng tại đây: Notes on RPGlab.

Bản gốc chương trình Durstenfeld Algol là here

Answer 1

Bởi vì bạn đã viết một c-chương trình trong một R-da

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

Answer 2

Cả Matlab và S (sau này là R) bắt đầu ra giấy gói mỏng xung quanh chức năng FORTRAN để làm công cụ toán học.

Trong S/R, các vòng lặp "luôn" bị chậm, nhưng điều đó đã được chấp nhận vì thường có cách vector hóa để diễn tả vấn đề. Ngoài ra, R có hàng ngàn chức năng trong Fortran hoặc C làm những việc cấp cao hơn một cách nhanh chóng. Ví dụ: chức năng sample thực hiện chính xác những gì mà vòng lặp của bạn thực hiện - nhưng nhanh hơn nhiều.

Vậy tại sao MATLAB tốt hơn nhiều khi thực thi kịch bản cho các vòng lặp? Hai lý do đơn giản: TÀI NGUYÊN và ƯU TIÊN.

MathWorks làm cho MATLAB là một công ty khá lớn với khoảng 2000 nhân viên. Họ đã quyết định cách đây nhiều năm để ưu tiên cải thiện hiệu suất của các tập lệnh. Họ thuê một loạt các chuyên gia biên dịch và dành nhiều năm để phát triển một trình biên dịch Just-In-Time (JIT) lấy mã kịch bản và biến nó thành mã assembly. Họ cũng đã làm rất tốt. Kudos cho họ!

R là nguồn mở và nhóm lõi R hoạt động để cải thiện R trong thời gian rảnh rỗi. Luke Tierney của R lõi đã làm việc chăm chỉ và phát triển một gói trình biên dịch cho R biên dịch R script để mã byte. Nó không biến nó thành mã assembler tuy nhiên, nhưng hoạt động khá tốt. Kudo với anh ấy!

... Nhưng những nỗ lực đưa vào trình biên dịch R vs trình biên dịch MATLAB chỉ đơn giản là ít hơn nhiều, và do đó kết quả là chậm:

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

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

Như bạn có thể thấy, cho vòng lặp đã trở thành 3x nhanh hơn bằng cách biên dịch nó thành mã byte. Một khác biệt nữa là trình biên dịch R JIT không được kích hoạt mặc định vì nó nằm trong MATLAB.

CẬP NHẬT Chỉ cần cho các hồ sơ, một phiên bản R hơi tối ưu hơn (dựa trên thuật toán Knuth), nơi mà các thế hệ ngẫu nhiên đã được vector hóa như @joran gợi ý:

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 

... vẫn còn là một cường độ chậm hơn MATLAB. Nhưng một lần nữa, chỉ cần sử dụng sample hoặc sample.int mà dường như đánh bại randperm MATLAB của 3x!

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

Answer 3

Hưởng ứng yêu cầu của OP là quá dài để phù hợp trong một chú thích, vì vậy đây là những gì tôi đã đề cập đến:

#Create r outside for loop 
GRPdur1 <- function(n){ 
    p <- 1:n       
    k <- seq(n,2,-1) 
    r <- 1 + floor(runif(length(k)) * k) 
    for (i in 1:length(k)){  
     tmp <- p[k[i]]; 
     p[k[i]] <- p[r[i]];     
     p[r[i]] <- tmp;     
    } 
return(p) 
} 

library(compiler) 
GRPdur2 <- cmpfun(GRPdur1) 

set.seed(1) 
out1 <- GRPdur(100) 

set.seed(1) 
out2 <- GRPdur1(100) 

#Check the GRPdur1 is generating the identical output 
identical(out1,out2) 

system.time(GRPdur(10^6)) 
    user system elapsed 
12.948 0.389 13.232 
system.time(GRPdur2(10^6)) 
    user system elapsed 
1.908 0.018 1.910 

Không hẳn 10x, nhưng nhiều hơn 3x Tommy cho thấy chỉ bằng cách sử dụng trình biên dịch. Đối với một thời gian khá chính xác hơn:

library(rbenchmark) 
benchmark(GRPdur(10^6),GRPdur2(10^6),replications = 10) 
      test replications elapsed relative user.self sys.self 
1 GRPdur(10^6)   10 127.315 6.670946 124.358 3.656   
2 GRPdur2(10^6)   10 19.085 1.000000 19.040 0.222 

Vì vậy, những nhận xét 10x là (có lẽ không đáng ngạc nhiên, được dựa trên một system.time chạy đơn) lạc quan, nhưng các vector giành bạn tốc độ nhanh hơn một chút công bằng hơn những gì trình biên dịch byte không .