evitar dos bucles en I

Question

Tengo un código R que puede hacer convolución de dos funciones ...

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 
} 

¿Hay una manera de quitar los dos bucles y hacer que el código se ejecute más rápido?

Gracias San

Answer 1

Puesto que R es muy rápido en el cálculo de operaciones vectoriales, lo más importante a tener en cuenta durante la programación para el rendimiento es Vectorise ya que muchos de sus operaciones como sea posible.

Esto significa pensar mucho sobre la sustitución de bucles con operaciones vectoriales. Aquí está mi solución para convolución rápida (50 veces más rápido con los vectores de entrada de longitud 1.000 cada uno):

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 
} 

Usted notará que el bucle interno (para j en ...) ha desaparecido. En cambio, lo reemplacé con una operación vectorial. j ahora se define como un vector (j < - 1: ny). Observe también que me refiero a todo el vector y, en lugar de subconjunto (es decir, y en lugar de y [j]).

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

escribí una pequeña función para medir peformance:

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

} 

Para dos vectores de longitud 1.000 cada uno, tengo un 98% de mejora del rendimiento:

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. Para los vectores, se indexa con [], no [[]], a fin de utilizar xy[ij] etc

2. convolución no Vectorise con facilidad, pero un truco común es para cambiar al compilado código. El manual Writing R Extensions usa convolución como un ejemplo en ejecución y muestra varias alternativas; también lo usamos mucho en la documentación Rcpp.

Answer 3

Algunos dicen que la aplican() y sapply() funciones son más rápidas que para() bucles en R. Se puede convertir la convolución de una función y llamaremos a partir de aplicar dentro de(). Sin embargo, no existe evidencia de lo contrario http://yusung.blogspot.com/2008/04/speed-issue-in-r-computing-apply-vs.html

Answer 4

¿Qué tal convolve(x, rev(y), type = "open") en 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 5

Como dice Dirk, el código compilado puede ser mucho más rápido. Tuve que hacer esto para uno de mis proyectos y me sorprendió la aceleración: ~ 40 veces más rápido que la solución de 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 

Me hicieron varios intentos para acelerar este proceso en I antes de que decidiera que el uso de código C no podía ser tan malo (nota: en realidad no lo era). Todos los míos fueron más lentos que los de Andrie, y fueron variantes al sumar el producto cruzado de manera apropiada. Una versión rudimentaria se puede hacer en solo tres líneas.

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

Esta versión solo ayuda un poco.

> 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 

Mi mejor versión era la siguiente:

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

Este fue un poco mejor, pero todavía no es tan rápido como de

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

Answer 6

La función convolveFast se puede optimizar un poco Andrie utilizando cuidadosamente matemáticas enteras solamente y reemplazando (1: ny) -1L con 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 
}