Perché lo zoo :: rollmean è lento rispetto a una semplice implementazione Rcpp?

Question

zoo::rollmean è una funzione utile che restituisce la media mobile di una serie temporale; per il vettore di lunghezza n e la dimensione della finestra k restituisce il vettore c(mean(x[1:k]), mean(x[2:(k+1)]), ..., mean(x[(n-k+1):n])).

ho notato che sembrava essere in esecuzione lentamente per un po 'di codice che stavo sviluppando, così ho scritto la mia versione con il pacchetto Rcpp e un semplice ciclo for:

library(Rcpp) 
cppFunction("NumericVector rmRcpp(NumericVector dat, const int window) { 
    const int n = dat.size(); 
    NumericVector ret(n-window+1); 
    double summed = 0.0; 
    for (int i=0; i < window; ++i) { 
    summed += dat[i]; 
    } 
    ret[0] = summed/window; 
    for (int i=window; i < n; ++i) { 
    summed += dat[i] - dat[i-window]; 
    ret[i-window+1] = summed/window; 
    } 
    return ret; 
}") 

Con mia grande sorpresa, questa versione di la funzione è molto più veloce rispetto alla funzione zoo::rollmean:

# Time series with 1000 elements 
set.seed(144) 
y <- rnorm(1000) 
x <- 1:1000 
library(zoo) 
zoo.dat <- zoo(y, x) 

# Make sure our function works 
all.equal(as.numeric(rollmean(zoo.dat, 3)), rmRcpp(y, 3)) 
# [1] TRUE 

# Benchmark 
library(microbenchmark) 
microbenchmark(rollmean(zoo.dat, 3), rmRcpp(y, 3)) 
# Unit: microseconds 
#     expr  min  lq  mean median  uq  max neval 
# rollmean(zoo.dat, 3) 685.494 904.7525 1776.88666 1229.2475 1744.0720 15724.321 100 
#   rmRcpp(y, 3) 6.638 12.5865 46.41735 19.7245 27.4715 2418.709 100 

l'aumento di velocità vale anche per i vettori molto più grandi:

# Time series with 5 million elements 
set.seed(144) 
y <- rnorm(5000000) 
x <- 1:5000000 
library(zoo) 
zoo.dat <- zoo(y, x) 

# Make sure our function works 
all.equal(as.numeric(rollmean(zoo.dat, 3)), rmRcpp(y, 3)) 
# [1] TRUE 

# Benchmark 
library(microbenchmark) 
microbenchmark(rollmean(zoo.dat, 3), rmRcpp(y, 3), times=10) 
# Unit: milliseconds 
#     expr  min   lq  mean  median   uq  max 
# rollmean(zoo.dat, 3) 2825.01622 3090.84353 3191.87945 3206.00357 3318.98129 3616.14047 
#   rmRcpp(y, 3) 31.03014 39.13862 42.67216 41.55567 46.35191 53.01875 

Perché un'implementazione semplice di Rcpp è in esecuzione ~ 100 volte più veloce di zoo::rollmean?

Answer 1

Grazie a @DirkEddelbuettel per aver sottolineato che il confronto effettuato nella domanda non era il più giusto perché stavo confrontando il codice C++ con il puro codice R. Quanto segue è una semplice implementazione di base R (senza tutti i controlli del pacchetto zoo); questo è abbastanza simile a come zoo::rollmean implementa il calcolo di base per la media di laminazione:

baseR.rollmean <- function(dat, window) { 
    n <- length(dat) 
    y <- dat[window:n] - dat[c(1, 1:(n-window))] 
    y[1] <- sum(dat[1:window]) 
    return(cumsum(y)/window) 
} 

Rispetto a zoo:rollmean, si vede che questo è ancora un buon affare più veloce:

set.seed(144) 
y <- rnorm(1000000) 
x <- 1:1000000 
library(zoo) 
zoo.dat <- zoo(y, x) 
all.equal(as.numeric(rollmean(zoo.dat, 3)), baseR.rollmean(y, 3), RcppRoll::roll_mean(y, 3), rmRcpp(y, 3)) 
# [1] TRUE 
library(microbenchmark) 
microbenchmark(rollmean(zoo.dat, 3), baseR.rollmean(y, 3), RcppRoll::roll_mean(y, 3), rmRcpp(y, 3), times=10) 
# Unit: milliseconds 
#      expr  min   lq  mean  median   uq  max neval 
#  rollmean(zoo.dat, 3) 507.124679 516.671897 646.813716 563.897005 593.861499 1220.08272 10 
#  baseR.rollmean(y, 3) 46.156480 47.804786 53.923974 49.250144 55.061844 76.47908 10 
# RcppRoll::roll_mean(y, 3) 7.714032 8.513042 9.014886 8.693255 8.885514 11.32817 10 
#    rmRcpp(y, 3) 7.729959 8.045270 8.924030 8.388931 8.996384 12.49042 10 

di approfondire il motivo per cui di vedere un aumento di velocità 10x durante l'utilizzo di base R, ho utilizzato strumento lineprof di Hadley, afferrando codice sorgente dalla sorgente zoo pacchetto dove necessario:

lineprof(rollmean.zoo(zoo.dat, 3)) 
#  time alloc release dups ref     src 
# 1 0.001 0.954  0 26 #27 rollmean.zoo/unclass 
# 2 0.001 0.954  0 0 #28 rollmean.zoo/:  
# 3 0.002 0.954  0 1 #28 rollmean.zoo   
# 4 0.001 1.431  0 0 #28 rollmean.zoo/seq_len 
# 5 0.001 0.000  0 0 #28 rollmean.zoo/c  
# 6 0.006 2.386  0 1 #28 rollmean.zoo   
# 7 0.002 0.954  0 2 #31 rollmean.zoo/cumsum 
# 8 0.001 0.000  0 0 #31 rollmean.zoo//  
# 9 0.005 1.912  0 1 #33 rollmean.zoo   
# 10 0.013 2.898  0 14 #33 rollmean.zoo/[<-  
# 11 0.299 28.941  0 127 #34 rollmean.zoo/na.fill 

Quasi tutto il tempo viene speso nella funzione na.fill, che in realtà viene chiamata dopo che i valori della media mobile sono già stati calcolati.

lineprof(na.fill.zoo(zoo.dat, fill=NA, 2:999999)) 
#  time alloc release dups ref     src 
# 1 0.004 1.913  0 39 #26 na.fill.zoo/seq  
# 2 0.002 1.921  0 9 #33 na.fill.zoo/coredata 
# 3 0.002 1.921  0 6 #37 na.fill.zoo/[<-  
# 4 0.001 0.955  0 10 #46 na.fill.zoo   
# 5 0.008 3.838  0 19 #46 na.fill.zoo/[<-  
# 6 0.003 0.959  0 2 #52 na.fill.zoo   
# 7 0.006 0.972  0 21 #52 na.fill.zoo/[<-  
# 8 0.001 0.486  0 0 #57 na.fill.zoo/seq_len 
# 9 0.005 0.959  0 6 #66 na.fill.zoo   
# 10 0.124 11.573  0 34 #66 na.fill.zoo/[ 

Quasi sempre vengono spesi sottoinsiemi dell'oggetto zoo:

lineprof("[.zoo"(zoo.dat, 2:999999)) 
# time alloc release dups   ref   src 
# 1 0.004 0.004  0 0 character(0)    
# 2 0.002 1.922  0 4   #4 [.zoo/coredata 
# 3 0.038 11.082  0 29   #19 [.zoo/zoo  
# 4 0.004 0.000  0 1   #28 [.zoo 

Quasi tutti i sottoinsiemi di tempo è trascorso costruzione di un nuovo oggetto zoo con la funzione zoo:

lineprof(zoo(y[2:999999], 2:999999)) 
# time alloc release dups    ref  src 
# 1 0.021 4.395  0 8 c("zoo", "unique") zoo/unique 
# 2 0.012 0.477  0 8 c("zoo", "ORDER") zoo/ORDER 
# 3 0.001 0.477  0 1    "zoo" zoo  
# 4 0.001 0.954  0 0  c("zoo", ":") zoo/:  
# 5 0.015 3.341  0 5    "zoo" zoo  

Varie operazioni necessarie per impostare un nuovo oggetto zoo (ad es. Determinare punti temporali unici e ordinarli).

In conclusione, il pacchetto zoo sembra aver aggiunto un sacco di spese generali per le sue operazioni di rotazione media costruendo un nuovo oggetto zoo invece di utilizzare le parti interne dell'oggetto attuale zoo; questo crea un rallentamento di 10 volte rispetto a un'implementazione di base R e un rallentamento di 100 volte rispetto a un'implementazione di Rcpp.

Answer 2

Rovistando in zoo sembrare che i metodi rollmean.* sono tutti in implementato in R.

Mentre si implementato uno in C++. Probabilmente il codice R confezionato fa anche qualche altro controllo ecc. Pp quindi forse non è così sorprendente che lo battete?