Как заставить цикл работать быстрее в R?

Question

Я хочу использовать arms(), чтобы каждый раз получать по одному сэмплу и делать цикл, подобный следующему в моей функции.Работает очень медленно.Как я мог заставить его работать быстрее?Спасибо.

library(HI)    
dmat <- matrix(0, nrow=100,ncol=30)
system.time(
    for (d in 1:100){
        for (j in 1:30){
            y <- rep(0, 101)
            for (i in 2:100){

                y[i] <- arms(0.3, function(x) (3.5+0.000001*d*j*y[i-1])*log(x)-x,
                    function(x) (x>1e-4)*(x<20), 1)       
            }
        dmat[d, j] <- sum(y)
        }
    }
) 

Answer 1

Это версия, основанная на ответе Томми, но избегающая всех петель:

library(multicore) # or library(parallel) in 2.14.x
set.seed(42)
m = 100
n = 30
system.time({
    arms.C <- getNativeSymbolInfo("arms")$address
    bounds <- 0.3 + convex.bounds(0.3, dir = 1, function(x) (x>1e-4)*(x<20))
    if (diff(bounds) < 1e-07) stop("pointless!")
    # create the vector of z values
    zval <- 0.00001 * rep(seq.int(n), m) * rep(seq.int(m), each = n)
    # apply the inner function to each grid point and return the matrix
    dmat <- matrix(unlist(mclapply(zval, function(z)
            sum(unlist(lapply(seq.int(100), function(i)
                .Call(arms.C, bounds, function(x) (3.5 + z * i) * log(x) - x,
                      0.3, 1L, parent.frame())
            )))
        )), m, byrow=TRUE)
}) 

На многоядерной машине это будет очень быстро, поскольку распределяет нагрузки по ядрам.На одноядерном компьютере (или для слабых пользователей Windows) вы можете заменить mclapply выше на lapply и получить лишь небольшое ускорение по сравнению с ответом Томми.Но обратите внимание, что результат будет отличаться для параллельных версий, поскольку он будет использовать разные последовательности RNG.

Обратите внимание, что любой код C, который должен оценивать функции R, будет по своей сути медленным (поскольку интерпретируемый код медленный).Я добавил arms.C только для того, чтобы удалить все накладные расходы R-> C, чтобы осчастливить moli;), но это не имеет значения.

Вы можете выжать еще несколько миллисекунд, используя column-основная обработка (код вопроса был главным по строке, который требует повторного копирования, поскольку матрицы R всегда являются основными по столбцам).

Редактировать: Я заметил, что Моли слегка изменила вопрос, так как Томми ответил- поэтому вместо части sum(...) вы должны использовать цикл, поскольку y[i] являются зависимыми, поэтому function(z) будет выглядеть как

function(z) { y <- 0
    for (i in seq.int(99))
         y <- y + .Call(arms.C, bounds, function(x) (3.5 + z * y) * log(x) - x,
                        0.3, 1L, parent.frame())
    y }

Answer 2

Ну, один эффективный способ - избавиться от накладных расходов внутри arms.Он делает некоторые проверки и вызывает indFunc каждый раз, хотя результат всегда одинаков в вашем случае.Некоторые другие оценки также могут быть выполнены вне цикла.Эти оптимизации сокращают время с 54 секунд до 6,3 секунд на моей машине.... и ответ идентичен.

set.seed(42)
#dmat2 <- ##RUN ORIGINAL CODE HERE##

# Now try this:
set.seed(42)
dmat <- matrix(0, nrow=100,ncol=30)
system.time({
    e <- new.env()
    bounds <- 0.3 + convex.bounds(0.3, dir = 1, function(x) (x>1e-4)*(x<20))
    f <- function(x) (3.5+z*i)*log(x)-x
    if (diff(bounds) < 1e-07) stop("pointless!")
    for (d in seq_len(nrow(dmat))) {
        for (j in seq_len(ncol(dmat))) {
            y <- 0
            z <- 0.00001*d*j
            for (i in 1:100) {
                y <- y + .Call("arms", bounds, f, 0.3, 1L, e)
            }
            dmat[d, j] <- y
        }
    }
}) 

all.equal(dmat, dmat2) # TRUE

Answer 3

почему бы не так?

dat <- expand.grid(d=1:10, j=1:3, i=1:10)

arms.func <- function(vec) {
  require(HI)
  dji <- vec[1]*vec[2]*vec[3]
  arms.out <- arms(0.3, 
                   function(x,params) (3.5 + 0.00001*params)*log(x) - x,
                   function(x,params) (x>1e-4)*(x<20),
                   n.sample=1,
                   params=dji)

  return(arms.out)
}

dat$arms <- apply(dat,1,arms.func)

library(plyr)
out <- ddply(dat,.(d,j),summarise, arms=sum(arms))

matrix(out$arms,nrow=length(unique(out$d)),ncol=length(unique(out$j)))

Тем не менее, он все еще одноядерный и требует много времени.Но это не медленный R, это функция рук.