Przyspieszenie obliczania row-mądry medianę każdym 3-krotki kolumn

Question

Jeśli mam ramki danych jako takie:

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

mogę korzystać z następujących funkcji, aby każde połączenie trzech perspektywie median rzędu -wise:

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

problem polega na tym, funkcja ta będzie trwać kilka godzin, aby uruchomić. Winowajcą jest mediana(), która trwa około dziesięć razy dłużej niż max() lub min().

W jaki sposób mogę przyspieszyć tę funkcję, prawdopodobnie poprzez pisanie szybszej wersji median() lub pracę z oryginalnymi danymi w inny sposób?

Update:

Jeżeli uruchomić powyższy kod ale tylko df [1 10] w następujący sposób:

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

trwa 29 sekund

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

z pakiet ccaPP trwa 6,5 ​​sekundy

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

trwa 2,5 sekundy

Widzimy więc znaczącą poprawę dzięki fastMedian(). Czy możemy jeszcze lepiej?

Answer 1

Jednym ze sposobów przyspieszenia tego procesu byłoby odnotowanie, że mediana trzech liczb jest ich sumą minus ich maks. Minus min. Oznacza to, że możemy wektoryzować nasze medianowe obliczenia, traktując każdą potrójną kolumnę tylko raz (wykonując medianę dla wszystkich wierszy z tym samym obliczeniem) zamiast obchodzić się z nią raz dla każdego rzędu.

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 

Wektoryzacja daje przyspieszenie 110x, gdy jest 10 kolumn i 5000 rzędów. Niestety nie mam maszyny z wystarczającą pamięcią do przechowywania 808,5 miliona numerów w wynikach dla pełnego przykładu.

Można przyspieszyć to dalej poprzez zaimplementowanie funkcji Rcpp, która pobiera jako dane wejściowe reprezentację wektorową macierzy (podobnie jak wektor uzyskany przez odczytanie macierzy w dół kolumn) wraz z liczbą wierszy i zwraca medianę każdego z nich. kolumna. Funkcja ta w dużej mierze opiera się na funkcji std::nth_element, która jest asymptotycznie liniowa w liczbie elementów, z których bierzesz medianę. (Zwróć uwagę, że nie uśredniam środkowych dwóch wartości, gdy przyjmuję medianę wektora o parzystej długości, ale zamiast tego biorę niższą z dwóch wartości).

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

Teraz tylko wywołać tę funkcję zamiast korzystania rowSums, pmin i 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 

W sumie zatem uzyskać 210x przyspieszenie; 110x przyspieszenia polega na przełączaniu z niewizułowanej aplikacji median na wektorową aplikację, a pozostałe 2x przyspieszenie polega na przełączaniu z kombinacji rowSums, pmin i pmax w celu obliczenia mediany w wektoryzowanym sposobie na podstawie Rcpp. podejście.