如果我有这样的数据框:
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()有了显着的改进.我们还能做得更好吗?
加快速度的一种方法是注意三个数字的中位数是它们的总和减去最大值减去它们的最小值.这意味着我们可以通过处理每个三次列(在同一计算中执行所有行的中位数)来矢量化我们的中值计算,而不是每行处理一次.
set.seed(144)
# Fully random matrix
df = matrix(rnorm(50000), 5000, 10)
original
josilber
combos
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
# user system elapsed
# 0.117 0.009 0.149
system.time(res.original
# 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 input = Rcpp::as<:vector> >(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
combos
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的加速是从非矢量化中值应用切换到矢量化应用,剩余的2x加速来自rowSums,pmin和pmax的组合切换,用于以矢量化方式计算中值到基于Rcpp的方式做法.