Somma cumulativa più veloce

Question

Sto provando a prendere somme cumulative per ogni colonna di una matrice. Ecco il mio codice R:

testMatrix = matrix(1:65536, ncol=256); 
microbenchmark(apply(testMatrix, 2, cumsum), times=100L); 

Unit: milliseconds 
         expr  min  lq  mean median  uq  max neval 
apply(testMatrix, 2, cumsum) 1.599051 1.766112 2.329932 2.15326 2.221538 93.84911 10000 

ho usato Rcpp per il confronto:

cppFunction('NumericMatrix apply_cumsum_col(NumericMatrix m) { 
    for (int j = 0; j < m.ncol(); ++j) { 
     for (int i = 1; i < m.nrow(); ++i) { 
      m(i, j) += m(i - 1, j); 
     } 
    } 
    return m; 
}'); 
microbenchmark(apply_cumsum_col(testMatrix), times=10000L); 

Unit: microseconds 
         expr  min  lq  mean median  uq  max neval 
apply_cumsum_col(testMatrix) 205.833 257.719 309.9949 265.986 276.534 96398.93 10000 

Così il codice C++ è 7,5 volte più veloce. È possibile fare meglio di apply(testMatrix, 2, cumsum) in puro R? Mi sembra di avere un ordine di grandezza sopra la testa senza motivo.

Answer 1

L'utilizzo di un ciclo per byte compilato da byte è leggermente più veloce della chiamata apply sul sistema. Mi aspettavo che fosse più veloce perché fa meno lavoro di apply. Come previsto, il ciclo R è ancora più lento della semplice funzione C++ che hai scritto.

colCumsum <- compiler::cmpfun(function(x) { 
    for (i in 1:ncol(x)) 
    x[,i] <- cumsum(x[,i]) 
    x 
}) 

testMatrix <- matrix(1:65536, ncol=256) 
m <- testMatrix 
require(microbenchmark) 
microbenchmark(colCumsum(m), apply_cumsum_col(m), apply(m, 2, cumsum), times=100L) 
# Unit: microseconds 
#     expr  min  lq median  uq  max neval 
#  matrixCumsum(m) 1478.671 1540.5945 1586.1185 2199.9530 37377.114 100 
# apply_cumsum_col(m) 178.214 192.4375 204.3905 234.8245 1616.030 100 
# apply(m, 2, cumsum) 1879.850 1940.1615 1991.3125 2745.8975 4346.802 100 
all.equal(colCumsum(m), apply(m, 2, cumsum)) 
# [1] TRUE 

Answer 2

È difficile battere C++ con solo codice R. Il modo più veloce che posso pensare di farlo è se sei disposto a dividere la tua matrice in una lista. In questo modo, R sta usando le funzioni primitive e non copia l'oggetto con ogni iterazione (apply è essenzialmente un ciclo piuttosto). Puoi vedere che C++ vince ancora, ma c'è un notevole aumento di velocità con l'approccio list se vuoi davvero usare il codice R.

fun1 <- function(){ 
    apply(testMatrix, 2, cumsum) 
} 

testList <- split(testMatrix, col(testMatrix)) 

fun2 <- function(){ 
    lapply(testList, cumsum) 
} 

microbenchmark(fun1(), 
       fun2(), 
       apply_cumsum_col(testMatrix), 
       times=100L) 


Unit: microseconds 
         expr  min  lq  mean median  uq  max neval 
         fun1() 3298.534 3411.9910 4376.4544 3477.608 3699.2485 9249.919 100 
         fun2() 558.800 596.0605 766.2377 630.841 659.3015 5153.100 100 
apply_cumsum_col(testMatrix) 219.651 282.8570 576.9958 311.562 339.5680 4915.290 100 

EDIT Si prega di notare che questo metodo è più lento di fun1 se si include il tempo di dividere la matrice in un elenco.