加快计算每列三元组的行中值

Question

如果我有这样的数据框：

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

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

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

问题是，这个函数需要几个小时才能运行。罪魁祸首是 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 倍的加速。不幸的是，我的机器没有足够的内存来存储完整示例输出中的 8.085 亿个数字。

您可以通过实现一个 Rcpp 函数来进一步加快速度，该函数将矩阵的向量表示（也就是通过向下读取矩阵的列获得的向量）以及行数和 returns 每列的中位数。该函数在很大程度上依赖于 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 倍的加速； 110 倍的加速来自从 median 的非矢量化应用程序切换到矢量化应用程序，其余 2 倍的加速来自从 rowSums、pmin 和 [= 的组合切换16=] 用于以矢量化方式计算基于 Rcpp 的方法的中位数。

加快计算每列三元组的行中值

Speed up computing the row-wise median of every 3-tuple of columns

r

performance

function

dataframe

median