2

我一直在测试 Rcpp 和 RcppArmadillo 来计算大矩阵的汇总统计数据。这比大约 400 万行、45 列的基本 R colMeans 或犰狳要快得多（快 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++ 外包。犰狳似乎有点慢，所以我想对与上面类似的代码进行就地排序，但我只是无法正确理解语法......这就是我正在尝试的......

# 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 糖和标准 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
4

2 回答 2

5

您可以将列复制到一个新的向量中

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
于 2013-04-04T21:54:49.030 回答
2

您实际上并没有显示 RcppArmadillo 代码——我对 RcppArmadillo 代码的性能非常满意，我需要行/列子集。

您可以通过 Rcpp 以同样有效的方式实例化 Armadillo 矩阵（无需复制，重新使用 R 对象内存），所以我会尝试这样做。

而你：你想要clone()一个独特的副本，如果你使用默认的 RcppArmadillo ctor（而不是更有效的两步），我认为你会免费获得它。

几个小时后编辑

你留下了一个悬而未决的问题，为什么你的犰狳很慢。与此同时，Vincent 为您解决了这个问题，但这里使用您的代码和 Vincent 的代码重新审视了一个更简洁的解决方案。

现在它如何在没有副本的情况下实例化犰狳矩阵——所以它更快。它还避免了混合整数和数字矩阵。先上代码：

#include <RcppArmadillo.h> 

using namespace Rcpp;

// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
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 column -- original will not be mod
        std::nth_element(y.begin(), y.begin() + position, y.end()); 
        out[j] = y[position];  
    }
    return out;
}

// [[Rcpp::export]]
arma::rowvec colMeansRcppArmadillo(NumericMatrix x){
    arma::mat X = arma::mat(x.begin(), x.nrow(), x.ncol(), false); 
    return arma::mean(X, 0); 
}

// [[Rcpp::export]]
NumericVector colMeansRcpp(NumericMatrix 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);
} 

/*** R
set.seed(42)
X <- matrix(rnorm(100*10), 100, 10)
library(microbenchmark)

mb <- microbenchmark(colMeans(X), colMeansRcpp(X), colMeansRcppArmadillo(X),
                     colMedianRcpp(X), times=50)  
print(mb)
*/

这是我机器上的结果，简明的犰狳版本和你的一样快，中值稍微慢一点，因为它必须做更多的工作：

R> sourceCpp("/tmp/stephen.cpp") 
R> set.seed(42)
R> X <- matrix(rnorm(100*10), 100, 10)
R> library(microbenchmark)
R> mb <- microbenchmark(colMeans(X), colMeansRcpp(X), colMeansRcppArmadillo(X),
+                      colMedianRcpp(X), times=50) 
R> print(mb)
Unit: microseconds
                     expr    min     lq  median     uq    max neval
              colMeans(X)  9.469 10.422 11.5810 12.421 30.597    50 
          colMeansRcpp(X)  3.922  4.281  4.5245  5.306 18.020    50 
 colMeansRcppArmadillo(X)  4.196  4.549  4.9295  5.927 11.159    50 
         colMedianRcpp(X) 15.615 16.291 16.7290 17.971 27.026    50 
R>
于 2013-04-04T20:16:26.500 回答