加快計算列

Question

的每3元組的在行中位數。如果我有一個數據幀，例如：

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

我可以使用下面的函數獲取三個學期中位數排的每個組合-wise：

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

問題是，這個函數需要幾個小時才能運行。罪魁禍首是median（），比max（）或min（）運行時間要長10倍。

如何通過寫入更高版本的median（）或使用原始數據以不同方式加速此功能？

更新：

如果我運行上面的代碼但僅針對DF [，1：10]，例如：

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

需要29秒

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

從包ccaPP需要6.5秒

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

需要2.5秒

所以我們看到fastMedian（）的顯着改進。我們還可以做得更好嗎？

Answer 1

加快速度的一種方法是注意三個數字的中位數是它們的總和減去它們的最大值減去它們的最小值。這意味着我們可以通過處理每個三列的列（向同一計算中的所有行執行中間值）而不是爲每行處理一次來對矢量化中值計算。

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 

當有10列和5000行時，矢量化會產生110倍加速。不幸的是，我沒有一臺具有足夠內存的機器來存儲輸出中的8.085億個數字。

您可以通過實現一個Rcpp函數來進一步提高速度，該函數將矩陣的向量表示（也就是通過讀取矩陣向下的列獲得的向量）與行數一起作爲輸入，並返回每個行的中值柱。該函數在很大程度上依賴於std::nth_element函數，該函數在中位數的元素數量上漸近線性。 （請注意，當我取一個長度爲偶數的向量的中值時，我不會平均中間的兩個值;而是取兩個中的較低者）。

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

現在我們只需調用，而不是使用rowSums，pmin和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 

總共所以我們得到了一個210X加速;加速的110倍是從median的非矢量化應用程序切換到矢量化應用程序，剩餘的2倍加速是通過從rowSums,pmin和pmax的組合中切換來計算向量化方式中的中值到基於Rcpp的方式做法。