Why does R slow down this random permutation function?

Question

I am new to R (Revolution Analytics R) and have translated some Matlab functions to R.

Question : Why is the function GRPdur (n) so slow?

GRPdur = function(n){
#
# Durstenfeld 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)
}

Here is what I get from the 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 

Here is what I get in Matlab 2011b

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

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

Here is what I get in 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.
>> 

---

REFERENCES : GRPdur is part of RPGlab , the Matlab function package I wrote that generates and tests various random permutation generators. Notes can be viewed here separately: Notes on RPGlab .

The original Durstenfeld Algol program is here

Answer 1

Matlab S ( R) FORTRAN .

S/R for-loops "" , , . , R Fortran C, . , sample, , for-loop, .

MATLAB ? : .

MathWorks, MATLAB, - , 2000 . . - Just-In-Time (JIT), script . . !

R , R R . R core 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-loop 3 , . , R JIT , MATLAB.

UPDATE - R ( ), , @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 version:
system.time(GRPdur1(10^6)) # 6.43
system.time(GRPdur2(10^6)) # 1.66

... , MATLAB. sample sample.int, , -, MATLAB randperm 3x!

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

Answer 2

c- R-skin

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

Answer 3

OP , , :

#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

10x, , 3x Tommy . :

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

So the 10x comment (perhaps not surprisingly based on a single run system.time) is optimistic, but vectorization gives you a higher speed than what the byte compiler does.