Почему R медленно работает с этой функцией случайной перестановки? - PullRequest
Новая
       39

Почему R медленно работает с этой функцией случайной перестановки?

6 голосов
/ 14 января 2012

Я новичок в 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 ГГц, 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.
>>

LINKS : GRPdur является частью RPGlab , пакетаФункции Matlab, которые я написал, генерируют и тестируют различные генераторы случайных перестановок.Примечания можно посмотреть отдельно здесь: Примечания к RPGlab .

Оригинальная программа Durstenfeld Algol - здесь

Ответы [ 3 ]

12 голосов
/ 15 января 2012

И Matlab, и S (позже R) начинали как тонкие обертки вокруг функций FORTRAN для выполнения математических задач.

В S / R циклы for "всегда" были медленными, но это было нормально, потому что обычно существуют векторизованные способы выражения проблемы.Кроме того, R имеет тысячи функций в Fortran или C, которые быстро выполняют высокоуровневые вещи.Например, функция sample, которая делает именно то, что делает ваш цикл for, но гораздо быстрее.

Так почему же тогда MATLAB намного лучше выполняет скриптовые циклы for?Две простые причины: РЕСУРСЫ и ПРИОРИТЕТЫ.

MathWorks, которые делают MATLAB, - довольно большая компания, в которой работает около 2000 сотрудников.Несколько лет назад они решили сделать приоритетным улучшение производительности скриптов.Они наняли группу экспертов по компиляции и потратили годы на разработку компилятора Just-In-Time (JIT), который берет код скрипта и превращает его в код ассемблера.Они тоже очень хорошо поработали.Престижность им!

R с открытым исходным кодом, и основная команда R работает над улучшением R в свободное время.Люк Тирни из ядра R усердно работал и разработал пакет компилятора для 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 стал в 3 раза быстрее, скомпилировав его в байт-код.Другое отличие состоит в том, что компилятор R JIT по умолчанию не включен, как в MATLAB.

ОБНОВЛЕНИЕ Просто для записи, немного более оптимизированная версия 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'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
6 голосов
/ 14 января 2012

Потому что вы написали c-программу на R-skin 

n = 10^6L
p = 1:n
system.time( sample(p,n))
0.03    0.00    0.03
2 голосов
/ 15 января 2012

Ответ на запрос ОП был слишком длинным, чтобы уместиться в комментарии, поэтому вот что я имел в виду: 

#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

Не совсем 10х, но больше, чем 3х Томми показал только при использовании компилятора.Для более точной синхронизации: 

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

Таким образом, 10-кратный комментарий был (возможно, неудивительно, основываясь на одном system.time прогоне) оптимистичным, но векторизация дает вам чуть больше скорости, чемэто делает байтовый компилятор.

...