Избегайте двух петель в R

Question

У меня есть код R, который может свертывать две функции ...

convolveSlow <- function(x, y) {  
nx <- length(x); ny <- length(y)  
xy <- numeric(nx + ny - 1)  
for(i in seq(length = nx)) {  
 xi <- x[[i]]  
        for(j in seq(length = ny)) {  
            ij <- i+j-1  
            xy[[ij]] <- xy[[ij]] + xi * y[[j]]  
        }  
    }  
    xy  
}  

Есть ли способ удалить два цикла for и заставить код работать быстрее?

Спасибо Сан

Answer 1

Поскольку R очень быстр в вычислениях векторных операций, при программировании для повышения производительности важно помнить о том, чтобы векторизовать как можно больше ваших операций.

Это означает, что нужно серьезно подумать о замене циклов навекторные операции.Вот мое решение для быстрой свертки (в 50 раз быстрее с входными векторами длиной 1000 каждый):

convolveFast <- function(x, y) {
    nx <- length(x)
    ny <- length(y)
    xy <- nx + ny - 1
    xy <- rep(0, xy)
    for(i in (1:nx)){
        j <- 1:ny
        ij <- i + j - 1
        xy[i+(1:ny)-1] <- xy[ij] + x[i] * y
    }
    xy
}

Вы заметите, что внутренний цикл (для j in ...) исчез.Вместо этого я заменил его векторной операцией.j теперь определяется как вектор (j <- 1: ny).Также обратите внимание, что я ссылаюсь на весь вектор y, а не на его поднабор (т. Е. Y вместо y [j]). </p>

j <- 1:ny
ij <- i + j - 1
xy[i+(1:ny)-1] <- xy[ij] + x[i] * y

Я написал небольшую функцию для измерения производительности:

measure.time <- function(fun1, fun2, ...){
    ptm <- proc.time()
    x1 <- fun1(...)
    time1 <- proc.time() - ptm

    ptm <- proc.time()
    x2 <- fun2(...)
    time2 <- proc.time() - ptm

    ident <- all(x1==x2)

    cat("Function 1\n")
    cat(time1)
    cat("\n\nFunction 2\n")
    cat(time2)
    if(ident) cat("\n\nFunctions return identical results")

}

Для двух векторов длиной 1000 я получаю улучшение производительности на 98%:

x <- runif(1000)
y <- runif(1000)
measure.time(convolveSlow, convolveFast, x, y)

Function 1
7.07 0 7.59 NA NA

Function 2
0.14 0 0.16 NA NA

Functions return identical results

Answer 2

1. Для векторов вы индексируете с помощью [], а не [[]], поэтому используйте xy[ij] и т. Д.

2. Свертка не легко векторизуется, кроме одногоОбычный трюк - переключиться на скомпилированный код.Руководство Writing R Extensions использует свертку в качестве рабочего примера и показывает несколько альтернатив;мы также часто используем его в документации Rcpp .

Answer 3

Функция convolveFast можно немного оптимизировать, осторожно используя только целочисленную математику и заменяя (1: ny) -1L на seq.int (0L, ny-1L):

convolveFaster <- function(x, y) {
    nx <- length(x)
    ny <- length(y)
    xy <- nx + ny - 1L
    xy <- rep(0L, xy)
    for(i in seq_len(nx)){
        j <- seq_len(ny)
        ij <- i + j - 1L
        xy[i+seq.int(0L, ny-1L)] <- xy[ij] + x[i] * y
    }
    xy
}

Answer 4

Как говорит Дирк, скомпилированный код может быть намного быстрее.Я должен был сделать это для одного из моих проектов и был удивлен ускорением: ~ в 40 раз быстрее, чем решение Andrie.

> a <- runif(10000)
> b <- runif(10000)
> system.time(convolveFast(a, b))
   user  system elapsed 
  7.814   0.001   7.818 
> system.time(convolveC(a, b))
   user  system elapsed 
  0.188   0.000   0.188 

Я сделал несколько попыток ускорить это в R, прежде чем решил, что с помощью кода на Cне может быть так плохо (примечание: это действительно не было).Все мои были медленнее, чем у Андри, и были вариантами правильного сложения перекрестного продукта.Элементарная версия может быть сделана всего за три строки.

convolveNotAsSlow <- function(x, y) {
  xyt <- x %*% t(y)
  ds <- row(xyt)+col(xyt)-1
  tapply(xyt, ds, sum)
}

Эта версия только немного помогает.

> a <- runif(1000)
> b <- runif(1000)
> system.time(convolveSlow(a, b))
   user  system elapsed 
  6.167   0.000   6.170 
> system.time(convolveNotAsSlow(a, b))
   user  system elapsed 
  5.800   0.018   5.820 

Моя лучшая версия была такой:

convolveFaster <- function(x,y) {
  foo <- if (length(x)<length(y)) {y %*% t(x)} else { x %*% t(y) }
  foo.d <- dim(foo)
  bar <- matrix(0, sum(foo.d)-1, foo.d[2])
  bar.rc <- row(bar)-col(bar)
  bar[bar.rc>=0 & bar.rc<foo.d[1]]<-foo
  rowSums(bar)
}

Это было немного лучше, но все же не так быстро, как у Андри

> system.time(convolveFaster(a, b))
   user  system elapsed 
  0.280   0.038   0.319 

Answer 5

Как насчет convolve(x, rev(y), type = "open") в stats?

> x <- runif(1000)
> y <- runif(1000)
> system.time(a <- convolve(x, rev(y), type = "o"))
   user  system elapsed 
  0.032   0.000   0.032 
> system.time(b <- convolveSlow(x, y))
   user  system elapsed 
 11.417   0.060  11.443 
> identical(a,b)
[1] FALSE
> all.equal(a,b)
[1] TRUE

Answer 6

Некоторые говорят, что функции apply () и sapply () работают быстрее, чем циклы for () в R. Вы можете преобразовать свертку в функцию и вызывать ее из apply (). Однако есть доказательства обратного http://yusung.blogspot.com/2008/04/speed-issue-in-r-computing-apply-vs.html