velocità fino calcolare la mediana fila-wise di ogni 3-upla di colonne

Question

Se ho un frame di dati come ad esempio:

df = data.frame(matrix(rnorm(100), 5000, 100)) 

posso utilizzare la seguente funzione per ottenere ogni combinazione di fila tre termini mediane -wise:

median_df = t(apply(df, 1, combn, 3, median)) 

Il problema è che questa funzione impiegherà diverse ore per essere eseguita. Il colpevole è mediano(), che richiede circa dieci volte più tempo per essere eseguito rispetto a max() o min().

Come posso accelerare questa funzione, possibilmente scrivendo una versione più veloce di median() o lavorando con i dati originali in modo diverso?

Aggiornamento:

Se faccio funzionare il codice di cui sopra, ma solo per df [, 1: 10] in quanto tale:

median_df = t(apply(df[,1:10], 1, combn, 3, median)) 

prende 29 secondi

fastMedian_df = t(apply(df[,1:10], 1, combn, 3, fastMedian)) 

dal pacchetto ccaPP impiega 6,5 ​​secondi

max_df = t(apply(df[,1:10], 1, combn, 3, max)) 

richiede 2,5 secondi

Quindi vediamo un miglioramento significativo con fastMedian(). Possiamo ancora fare meglio?

Answer 1

Un approccio per accelerare le cose sarebbe notare che la mediana di tre numeri è la loro somma meno il loro massimo meno il loro minimo. Ciò significa che possiamo vettorizzare i nostri calcoli mediani gestendo una volta ogni tripla di colonne (eseguendo la mediana per tutte le righe nello stesso calcolo) invece di gestirla una volta per ogni riga.

set.seed(144) 
# Fully random matrix 
df = matrix(rnorm(50000), 5000, 10) 
original <- function(df) t(apply(df, 1, combn, 3, median)) 
josilber <- function(df) { 
    combos <- combn(seq_len(ncol(df)), 3) 
    apply(combos, 2, function(x) rowSums(df[,x]) - pmin(df[,x[1]], df[,x[2]], df[,x[3]]) - pmax(df[,x[1]], df[,x[2]], df[,x[3]])) 
} 
system.time(res.josilber <- josilber(df)) 
# user system elapsed 
# 0.117 0.009 0.149 
system.time(res.original <- original(df)) 
# user system elapsed 
# 15.107 1.864 16.960 
all.equal(res.josilber, res.original) 
# [1] TRUE 

La vettorizzazione produce un aumento di 110 volte quando ci sono 10 colonne e 5000 righe. Sfortunatamente non ho una macchina con abbastanza memoria per archiviare i 808.5 milioni di numeri nell'output per il tuo esempio completo.

Si potrebbe accelerare ulteriormente implementando una funzione Rcpp che prende come input la rappresentazione vettoriale di una matrice (ovvero il vettore ottenuto leggendo la matrice lungo le colonne) insieme al numero di righe e restituisce la mediana di ciascuna colonna. La funzione si basa pesantemente sulla funzione std::nth_element, che è asintoticamente lineare nel numero di elementi che stai prendendo una mediana di. (Si noti che non faccio la media dei due valori medi quando prendo la mediana di un vettore di lunghezza pari, io invece prendo il più basso dei due).

library(Rcpp) 
cppFunction(
"NumericVector vectorizedMedian(NumericVector x, int chunkSize) { 
const int n = x.size()/chunkSize; 
std::vector<double> input = Rcpp::as<std::vector<double> >(x); 
    NumericVector res(n); 
    for (int i=0; i < n; ++i) { 
    std::nth_element(input.begin()+i*chunkSize, input.begin()+i*chunkSize+chunkSize/2, 
        input.begin()+(i+1)*chunkSize); 
    res[i] = input[i*chunkSize+chunkSize/2]; 
    } 
    return res; 
}") 

Ora dobbiamo solo invocare questa funzione invece di utilizzare rowSums, pmin e pmax:

josilber.rcpp <- function(df) { 
    combos <- combn(seq_len(ncol(df)), 3) 
    apply(combos, 2, function(x) vectorizedMedian(as.vector(t(df[,x])), 3)) 
} 
system.time(josilber.rcpp(df)) 
# user system elapsed 
# 0.049 0.008 0.081 
all.equal(josilber(df), josilber.rcpp(df)) 
# [1] TRUE 

In totale abbiamo quindi ottenere un aumento di velocità 210x; 110x di accelerazione deriva dal passaggio da un'applicazione non vettorizzata di median a un'applicazione vettorizzata e il restante 2x speedup deriva dal passaggio da una combinazione di rowSums, pmin e pmax per il calcolo della mediana in un modo vettorializzato a un sistema basato su Rcpp approccio.