2

大きな行列の要約統計を計算するために、Rcpp と RcppArmadillo をテストしてきました。これは、400 万行、45 列で、ベースの R colMeans または Armadillo よりもはるかに高速 (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);
                          ')

私は本当に、プロットのために中央値とおそらく他の分位数を計算したいと思っています. アルマジロは少し遅いように見えるので、上記のようなコードでインプレースソートを実行したかったのですが、正しい構文を取得できません...これが私が試みていることです..

# 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 R オブジェクトを変更しないように、この clone() を実行する必要がありますか?

EDIT RcppArmadilloコードとベンチマーク比較を追加して、以下の回答/コメントに対処します。ベンチマークは、迅速な返信のために 50k 行のみでしたが、それ以上の行でも同様だったことを思い出します。私はあなたが Rcpp の作者であることを理解しています.. お時間をいただきありがとうございます!

ベースの colMeans または Rcpp バージョンよりもはるかに遅く実行するために、おそらく RcppArmadillo コードで何か馬鹿なことをしているのだろうか?

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()個別のコピーが必要で、(より効率的な 2 ステップではなく) デフォルトの RcppArmadillo ctor を使用すれば、無料で入手できると思います。

数時間後に編集

アルマジロが遅い理由について未解決の質問を残しました。それまでの間、Vincent が問題を解決してくれましたが、Vincent のコードと同様にあなたのコードを使用した、再検討されたよりクリーンなソリューションを次に示します。

コピーなしで 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)
*/

そして、これが私のマシンでの結果です。簡潔なArmadilloバージョンはあなたのものとほぼ同じくらい速く、中央値はより多くの作業を行う必要があるため少し遅くなります:

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 に答える