evitare due cicli for in R

Question

Ho un codice R che può fare convoluzione di due funzioni ...

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 
} 

C'è un modo per rimuovere i due cicli for e rendere il codice di correre più veloce?

Grazie San

Answer 1

Poiché R è molto veloce a calcolare le operazioni vettoriali, la cosa più importante da tenere a mente durante la programmazione per le prestazioni è quello di vectorise come molte delle vostre operazioni possibili.

Ciò significa riflettere seriamente sulla sostituzione di loop con operazioni vettoriali. Ecco la mia soluzione per convoluzione veloce (50 volte più veloce con vettori di ingresso di lunghezza 1000 ciascuno):

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 
} 

si noterà che il ciclo interno (per j in ...) è scomparso. Invece, l'ho sostituito con un'operazione vettoriale. j è ora definito come un vettore (j < - 1: ny). Si noti inoltre che mi riferisco all'intero vettore y, piuttosto che al suo sottotitolo (cioè y invece di y [j]).

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

ho scritto una piccola funzione per misurare 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") 

} 

Per due vettori di lunghezza 1000 ciascuna, ottengo un miglioramento delle prestazioni del 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. per i vettori, è indice con [], non [[]], in modo da utilizzare xy[ij] ecc

2. Convoluzione non vectorise facilmente, ma un trucco comune è quello di passare alla compilazione codice. Le estensioni per scrittura R utilizzano la convoluzione come esempio in esecuzione e mostrano diverse alternative; lo usiamo anche molto nella documentazione Rcpp.

Answer 3

Alcuni dicono che la applicano() e sapply (funzioni) sono più veloci rispetto a() loop in R. Si potrebbe convertire il convoluzione a una funzione e chiamare da applicare all'interno di(). Tuttavia, non v'è prova del contrario http://yusung.blogspot.com/2008/04/speed-issue-in-r-computing-apply-vs.html

Answer 4

Come su convolve(x, rev(y), type = "open") in 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

Come dice Dirk, il codice compilato può essere molto più veloce. Ho dovuto farlo per uno dei miei progetti e sono rimasto sorpreso dall'accelerazione: ~ 40 volte più veloce della soluzione di 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 

ho fatto diversi tentativi per accelerare questo in R prima ho deciso che utilizzando il codice C non poteva essere così male (nota: in realtà non era). Tutti i miei erano più lenti di quelli di Andrie, e c'erano varianti per aggiungere il prodotto incrociato in modo appropriato. Una versione rudimentale può essere fatta in sole tre righe.

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

Questa versione aiuta solo un po '.

> 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 

La mia migliore versione è stata questa:

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

Questo è stato un po 'meglio, ma ancora non così veloce come di

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

Answer 6

La funzione convolveFast può essere ottimizzato un po Andrie usando con attenzione solo matematica intera e sostituendo (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 
}