Question

我一直在測試Rcpp和RcppArmadillo以計算大矩陣的彙總統計信息。這比基礎的R colMeans或者Arm​​adillo在約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++外包。犰狳似乎有點慢，所以我想在那種地方做一個類似的代碼上面，但我只是不能得到正確的語法...這裏就是我想..

# 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); 
         ') 

基本上它是line

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

我不知道如何用Rcpp糖和std C++的混合物來表達「就地排序」。對不起，我可以看到我在做什麼是錯誤的，但在正確的語法提示將是可愛的。

ps我是否正確我需要做這個clone（），所以我不改變R對象？

編輯 我添加RcppArmadillo代碼和基準比較，以解決下面的答案/評論。該基準僅在5萬行獲得了快速回復，但我記得它與其他許多類似。我意識到你是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

你實際上並沒有顯示RcppArmadillo代碼 - 我一直用RcppArmadillo代碼，我需要行/列列子集的表現相當滿意。

您可以通過Rcpp實例化Armadillo矩陣，就像高效（不復制，重新使用R對象內存）一樣，所以我會試試。

而你：你希望clone()爲不同的副本，我想你會得到免費的，如果你使用默認的RcppArmadillo ctor（而不是更有效的兩步）。

編輯了幾個小時後

你已經離開了，爲什麼你的犰狳是緩慢的一個懸而未決的問題。與此同時，文森特爲您解決了這個問題，但這是一個使用您的代碼以及文森特的重新清潔的解決方案。

現在，它如何實例化Armadillo矩陣而無需複製 - 因此速度更快。它還避免了混合整數和數字矩陣。代碼首先：

#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> 

Answer 2

您可以將列複製到一個新的載體與

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