¿Por qué R es lento en esta función de permutación aleatoria?

Question

Soy nuevo a R (Revolution Analytics R) y he estado traduciendo algunas funciones de Matlab en R.

Pregunta: ¿Por qué es la función GRPdur (n) es tan lenta?

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

Aquí es lo que me pasa en un Dell Precision 690, 2xQuadcore Xeon 5345 @ 2,33 GHz, Windows 7 de 64 bits:

> 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 

Aquí es lo que me pasa en Matlab 2011b

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

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

Aquí es lo que me pasa en 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. 
>> 

---

ENLACES: GRPdur es parte de RPGlab, un paquete de funciones de Matlab que escribí que genera y prueba varios generadores de permutación aleatoria. Las notas se pueden ver por separado aquí: Notes on RPGlab.

El programa original Durstenfeld Algol es here

Answer 1

Debido a que escribió un programa en C en un R-piel

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

Answer 2

Tanto Matlab y S (más tarde R) comenzaron envolturas delgadas alrededor funciones FORTRAN para hacer cosas de matemáticas

En S/R los bucles for-loops "siempre" han sido lentos, pero eso ha estado bien porque generalmente hay formas vectorizadas de expresar el problema. Además, R tiene miles de funciones en Fortran o C que hacen cosas de alto nivel rápidamente. Por ejemplo, la función sample que hace exactamente lo que hace for-loop, pero mucho más rápido.

Entonces, ¿por qué entonces MATLAB es mucho mejor al ejecutar for-loop con scripts? Dos razones simples: RECURSOS y PRIORIDADES.

MathWorks que hacen MATLAB es una empresa bastante grande con alrededor de 2000 empleados. Decidieron hace años priorizar la mejora del rendimiento de los scripts. Contrataron a un grupo de expertos en compilación y pasaron años desarrollando un compilador Just-In-Time (JIT) que toma el código del script y lo convierte en código de ensamblador. Hicieron un muy buen trabajo también. Felicitaciones a ellos!

R es de código abierto, y el equipo central R trabaja para mejorar R en su tiempo libre. Luke Tierney de R core ha trabajado duro y ha desarrollado un paquete compilador para R que compila los guiones R con el código de bytes. NO lo convierte en código ensamblador, pero funciona bastante bien. Felicitaciones a él!

... Pero la cantidad de esfuerzo puesto en el compilador R vs el compilador MATLAB es simplemente mucho menos, y por lo tanto el resultado es más lento:

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

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

Como se puede ver, el bucle de se volvió 3 veces más rápido compilando el código de bytes. Otra diferencia es que el compilador R JIT no está habilitado por defecto, ya que está en MATLAB.

ACTUALIZACIÓN Sólo para que conste, una versión optimizada R poco más (basado en el algoritmo de Knuth), donde la generación aleatoria se ha vectorizado como @joran sugirió:

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 

... todavía una magnitud más lento que MATLAB. Pero nuevamente, solo use sample o sample.int que aparentemente supera el randperm de MATLAB por 3x.

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

Answer 3

En respuesta a la petición de la OP era demasiado largo para caber en un comentario, así que aquí es lo que me refería:

#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 

No exactamente 10 veces, pero más que el 3x Tommy mostró simplemente utilizando la compilador. Durante un tiempo un poco más preciso:

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 

Así que el comentario 10x era (quizás no es sorprendente que, basándose en una sola system.time plazo) optimista, pero la vectorización le permite ganar un poco justo más velocidad sobre lo que hace el byte-compilador .