加快计算每个3元组列的行中位数

Question

如果我有这样的数据框:

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

我可以使用以下函数来逐行获得三项中位数的每个组合:

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

问题是,此功能需要几个小时才能运行.罪魁祸首是中位数(),运行时间比max()或min()大十倍.

如何加快这个功能,可能是通过编写更快版本的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倍的加速.不幸的是,我没有一台机器有足够的内存来存储输出中的808.5百万个数字.

您可以通过实现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

因此,我们总共获得210倍的加速; 加速比的110X是从由非矢量化应用的切换median到矢量化应用程序和剩余的2倍的加速比是从一个组合开关的rowSums,pmin和pmax用于以矢量方式在基于RCPP的方法计算中值.