Machen kumulative Summe schneller

Question

Ich versuche kumulative Summen für jede Spalte einer Matrix zu nehmen. Hier ist mein Code in 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 

I verwendet RCPP zum Vergleich:

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 

So ist der C++ Code 7,5 mal so schnell. Ist es möglich, besser als apply(testMatrix, 2, cumsum) in reinem R zu tun? Es fühlt sich an, als hätte ich eine Größenordnung ohne jeden Grund.

Answer 1

Die Verwendung einer Byte-kompilierten For-Schleife ist etwas schneller als der apply Aufruf auf meinem System. Ich habe erwartet, dass es schneller ist, weil es weniger Arbeit als apply leistet. Wie erwartet, ist die R-Schleife immer noch langsamer als die einfache C++ - Funktion, die Sie geschrieben haben.

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

Es ist schwierig, C++ nur mit R-Code zu schlagen. Der schnellste Weg, an den ich denken kann, ist, wenn Sie bereit sind, Ihre Matrix in eine Liste aufzuteilen. Auf diese Weise verwendet R primitive Funktionen und kopiert das Objekt nicht mit jeder Iteration (apply ist im Wesentlichen eine hübsche Schleife). Sie können sehen, dass C++ immer noch gewinnt, aber es gibt eine deutliche Beschleunigung mit dem list Ansatz, wenn Sie wirklich nur R-Code verwenden möchten.

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 Bitte beachten Sie, dass diese Methode ist langsamer als fun1 wenn Sie die Zeit umfassen die Matrix in einer Liste zu teilen.