uniknąć dwóch dla pętli w R

Question

Mam kod R, który może zrobić splot dwóch funkcji ...

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 
} 

Czy istnieje sposób, aby usunąć dwa dla pętli i uczynić bieg kodu szybciej?

Dziękuję San

Answer 1

Ponieważ R jest bardzo szybki przy obliczaniu operacje wektorowe, najważniejszą rzeczą, aby pamiętać przy programowaniu dla wydajności jest vectorise jak wiele swoich działań, jak to możliwe.

Oznacza to, że trudno jest zastąpić pętle operacjami wektorowymi. Oto moje rozwiązanie do szybkiego splotu (50 razy szybciej z wejściowych wektorów o długości 1000 każdy):

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 
} 

Można zauważyć, że wewnętrzna pętla (dla j in ...) zniknął. Zamiast tego zastąpiłem go operacją wektorową. j jest teraz zdefiniowany jako wektor (j < - 1: ny). Zauważ też, że odwołuję się do całego wektora y, zamiast go podsegmentować (tj. Y zamiast y [j]).

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

napisałem małą funkcję pomiaru 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") 

} 

Przez dwa wektory o długości 1000 każdy, mam 98% wzrost wydajności:

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. dla wektorów, indeksowania z [], nie [[]], więc używaj xy[ij] itp

2. Splot nie vectorise łatwo, ale jeden wspólny Sztuką jest, aby przełączyć się skompilowany kod. W instrukcji ręcznej wykorzystano splot jako działający przykład i pokazano kilka alternatyw; używamy go również dużo w dokumentacji Rcpp.

Answer 3

Niektórzy mówią, że stosuje się() i sapply (funkcje) są szybsze niż() w pętli R. Można konwertować splot do funkcji i wywołać ją od wewnątrz apply(). Nie ma jednak dowodów przeciwnych http://yusung.blogspot.com/2008/04/speed-issue-in-r-computing-apply-vs.html

Answer 4

Jak o convolve(x, rev(y), type = "open") w 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

Jak mówi Dirk, skompilowany kod może być znacznie szybszy. Musiałem to zrobić w jednym z moich projektów i byłem zaskoczony przyspieszeniem: ~ 40x szybciej niż rozwiązanie 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 

Zrobiłem kilka prób to przyspieszyć w R zanim zdecydowałem, że za pomocą kodu C nie może być tak źle (uwaga: to naprawdę nie było). Wszystkie moje były wolniejsze niż Andrie i były wariantami, które odpowiednio dodawały produkt krzyżowy. Podstawową wersję można zrobić w zaledwie trzech liniach.

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

Ta wersja tylko trochę pomaga.

> 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 

Moja najlepsza wersja była to:

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

To było trochę lepiej, ale nadal nie tak szybko, jak Andrie na

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

Answer 6

Funkcja convolveFast można zoptymalizować trochę przez ostrożne użycie wyłącznie liczby całkowitej i zastąpienie (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 
}