r - 对 Rcpp NumericMatrix 的列进行排序以进行中值计算

Question

我一直在测试 Rcpp 和 RcppArmadillo 来计算大矩阵的汇总统计数据。这比基础 R colMeans 或 Armadillo 在大约 400 万行、45 列上要快得多(快 5 或 10 倍)。

colMeansRcpp <- cxxfunction(signature(X_="integer"), 
                            plugin='Rcpp',
                            body='
                            Rcpp::IntegerMatrix X = X_;
                            int ncol = X.ncol(); int nrow = X.nrow();                      
                            Rcpp::NumericVector out(ncol);
                            for(int col = 0; col < ncol; col++){
                              out[col]=Rcpp::sum(X(_, col));
                            }                             
                            return wrap(out/nrow);
                          ')

我真的很想计算中位数，也许还可以计算其他分位数来绘图 - 因为它需要某种排序，所以更需要 C++ 外包。 Armadillo 似乎有点慢，所以我想对类似于上面的代码进行就地排序，但我只是无法获得正确的语法......这就是我正在尝试的......

# OK I'm aware this floor(nrow/2) is not **absolutely** correct 
# I'm simplifying here
    colMedianRcpp <- cxxfunction(signature(X_="integer"), 
                          plugin='Rcpp',
                          body='
                          Rcpp::IntegerMatrix X = clone(X_);
                          int ncol = X.ncol(); int nrow = X.nrow();                           
                          Rcpp::NumericVector out(ncol);
                          for(int col = 0; col < ncol; col++){
                          X(_,col)= std::sort((X_,col).begin, (X_,col).end));
                          out[col]=X(floor(nrow/2), col));
                          }
                        return wrap(out);
                        ')

基本上就是这条线

X(_,col)= std::sort((X_,col).begin, (X_,col).end));

我不知道如何用 Rcpp 糖和 std C++ 的混合物来表达“就地排序列”。抱歉，我知道我正在做的事情是错误的，但是有关正确语法的提示会很可爱。

ps 我对吗？我需要执行此clone()，这样我就不会更改 R 对象？

编辑我添加了 RcppArmadillo 代码和基准比较来解决下面的答案/评论。为了快速回复，基准测试仅针对 50k 行，但我记得它与更多行类似。我知道您是 Rcpp 作者......非常感谢您的宝贵时间!

我想到也许我对 RcppArmadillo 代码做了一些愚蠢的事情，使其运行速度比基本 colMeans 或 Rcpp 版本慢得多？

colMeansRcppArmadillo <- cxxfunction(signature(X_="integer"), 
                                     plugin="RcppArmadillo",
                                      body='
                                      arma::mat X = Rcpp::as<arma::mat > (X_);
                                      arma::rowvec MD= arma::mean(X, 0);
                                      return wrap(MD);
                                    ')

基准是......

(mb = microbenchmark(
+                     colMeans(fqSmallMatrix), 
+                     colMeansRcpp(fqSmallMatrix), 
+                     colMeansRcppArmadillo(fqSmallMatrix),
+                     times=50))
Unit: milliseconds
                                 expr       min       lq    median        uq        max neval
              colMeans(fqSmallMatrix) 10.620919 10.63289 10.640819 10.648882  10.907145    50
          colMeansRcpp(fqSmallMatrix)  2.649038  2.66832  2.676709  2.700839   2.841012    50
 colMeansRcppArmadillo(fqSmallMatrix) 25.687067 26.23488 33.168589 33.792489 113.832495    50

Answer 1

您可以使用

将列复制到新向量中

NumericVector y = x(_,j);

完整示例:

library(Rcpp)
cppFunction('
  NumericVector colMedianRcpp(NumericMatrix x) {
    int nrow = x.nrow();
    int ncol = x.ncol();
    int position = nrow / 2; // Euclidian division
    NumericVector out(ncol);
    for (int j = 0; j < ncol; j++) {
      NumericVector y = x(_,j); // Copy the column -- the original will not be modified
      std::nth_element(y.begin(), y.begin() + position, y.end());
      out[j] = y[position];
    }
    return out;
  }
')
x <- matrix( sample(1:12), 3, 4 )
x
colMedianRcpp(x)
x   # Unchanged