我希望能够在Rcpp中实现一个简单的“分割-应用-组合”程序,其中数据集(矩阵)被分成组,并返回每组列的总和。这是一个在R中容易实现但常常需要相当长时间的过程。我已经成功地实现了一个Rcpp解决方案,它超越了R的性能,但我想知道是否可以进一步改进。为了说明,这里有一些代码,首先是使用R的代码:
然而,根据这个问题的回答,我了解到在
对这些解决方案进行基准测试,包括@Roland提供的
n <- 50000
k <- 50
set.seed(42)
X <- matrix(rnorm(n*k), nrow=n)
g=rep(1:8,length.out=n )
use.for <- function(mat, ind){
sums <- matrix(NA, nrow=length(unique(ind)), ncol=ncol(mat))
for(i in seq_along(unique(ind))){
sums[i,] <- colSums(mat[ind==i,])
}
return(sums)
}
use.apply <- function(mat, ind){
apply(mat,2, function(x) tapply(x, ind, sum))
}
use.dt <- function(mat, ind){ # based on Roland's answer
dt <- as.data.table(mat)
dt[, cvar := ind]
dt2 <- dt[,lapply(.SD, sum), by=cvar]
as.matrix(dt2[,cvar:=NULL])
}
事实证明,for
循环实际上非常快,并且使用Rcpp
实现最容易(对我来说)。它的工作原理是为每个组创建一个子矩阵,然后在矩阵上调用colSums
。这是使用RcppArmadillo
实现的:
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
using namespace arma;
// [[Rcpp::export]]
arma::mat use_arma(arma::mat X, arma::colvec G){
arma::colvec gr = arma::unique(G);
int gr_n = gr.n_rows;
int ncol = X.n_cols;
arma::mat out = zeros(gr_n, ncol);
for(int g=0; g<gr_n; g++){
int g_id = gr(g);
arma::uvec subvec = find(G==g_id);
arma::mat submat = X.rows(subvec);
arma::rowvec res = sum(submat,0);
out.row(g) = res;
}
return out;
}
然而,根据这个问题的回答,我了解到在
C++
(就像在R
中一样),创建副本是代价昂贵的,但循环不像在R
中那么糟糕。由于arma
方案依赖于为每个组创建矩阵(代码中的submat
),因此我猜想避免这种情况将进一步加速该过程。因此,这里是第二个仅使用循环基于Rcpp
的实现:#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericMatrix use_Rcpp(NumericMatrix X, IntegerVector G){
IntegerVector gr = unique(G);
std::sort(gr.begin(), gr.end());
int gr_n = gr.size();
int nrow = X.nrow(), ncol = X.ncol();
NumericMatrix out(gr_n, ncol);
for(int g=0; g<gr_n; g++){
int g_id = gr(g);
for (int j = 0; j < ncol; j++) {
double total = 0;
for (int i = 0; i < nrow; i++) {
if (G(i) != g_id) continue; // not sure how else to do this
total += X(i, j);
}
out(g,j) = total;
}
}
return out;
}
对这些解决方案进行基准测试,包括@Roland提供的
use_dt
版本(我的先前版本不公平地歧视了data.table
),以及@beginneR建议的dplyr
解决方案,得到以下结果: library(rbenchmark)
benchmark(use.for(X,g), use.apply(X,g), use.dt(X,g), use.dplyr(X,g), use_arma(X,g), use_Rcpp(X,g),
+ columns = c("test", "replications", "elapsed", "relative"), order = "relative", replications = 1000)
test replications elapsed relative
# 5 use_arma(X, g) 1000 29.65 1.000
# 4 use.dplyr(X, g) 1000 42.05 1.418
# 3 use.dt(X, g) 1000 56.94 1.920
# 1 use.for(X, g) 1000 60.97 2.056
# 6 use_Rcpp(X, g) 1000 113.96 3.844
# 2 use.apply(X, g) 1000 301.14 10.156
我的直觉是(use_Rcpp
比use_arma
更好),但事实证明我错了。话虽如此,我猜测在我的use_Rcpp
函数中的这行代码if (G(i) != g_id) continue;
减慢了整个过程。如果有其他可替代方案,我很愿意学习。
我很高兴自己完成了与R
相同的任务所需时间的一半,但也许若干Rcpp比R快得多的
例子误导了我的期望值,我想知道是否还能进一步提速。请问有人有什么想法吗?因为我对Rcpp
和C++
都相对陌生,所以我也欢迎任何编程或代码的评论。
setDT
避免这种情况。 - RolandX
和求和输出进行多次矩阵乘法。出于这个原因,一开始就将所有东西设置为矩阵。虽然我可以将一些内容更改为data.frame形式,但是除非我将所有内容转换为矩阵,否则我不知道如何继续操作。data.table是否允许对整个表进行类似矩阵的操作,例如取逆等? - coffeinjunkycvar := ind
。您可以直接在聚合中使用ind
,如下所示:dt[,lapply(.SD, sum), by=ind]
。 - Arun