Rcpp：私の距離行列プログラムがパッケージ内の関数よりも遅い

Question

私はペアワイドユークリッド距離行列を計算したいと思います。

NumericMatrix calcPWD1 (NumericMatrix x){ 
    int outrows = x.nrow(); 
    double d; 
    NumericMatrix out(outrows,outrows); 

    for (int i = 0 ; i < outrows - 1; i++){ 
    for (int j = i + 1 ; j < outrows ; j ++){ 
     NumericVector v1= x.row(i); 
     NumericVector v2= x.row(j); 
     NumericVector v3=v1-v2; 
     d = sqrt(sum(pow(v3,2))); 
     out(j,i)=d; 
     out(i,j)=d; 
    } 
    } 
    return (out) ; 
} 

を次のように私はダークEddelbuettelの提案によりRcppプログラムを書いたしかし、私は私のプログラムがdist関数よりも遅いです見つけます。

> benchmark(as.matrix(dist(b)),calcPWD1(b)) 
       test replications elapsed relative user.self sys.self user.child sys.child 
1 as.matrix(dist(b))   100 24.831 1.000 24.679 0.010   0   0 
2  calcPWD1(b)   100 27.362 1.102 27.346 0.007   0   0 

皆さんお気軽にお問い合わせください。私の行列はとてもシンプルです。列名や行名はありません。単純な行列です（たとえば、b=matrix(c(rnorm(1000*10)),1000,10)など）。ここ は、私がdistで、（method、diagのように）チェックする必要があるためにあまりにも多くのものがありますので、私のプログラムがdistよりも高速である期待dist

> dist 
function (x, method = "euclidean", diag = FALSE, upper = FALSE, 
    p = 2) 
{ 
    if (!is.na(pmatch(method, "euclidian"))) 
     method <- "euclidean" 
    METHODS <- c("euclidean", "maximum", "manhattan", "canberra", 
     "binary", "minkowski") 
    method <- pmatch(method, METHODS) 
    if (is.na(method)) 
     stop("invalid distance method") 
    if (method == -1) 
     stop("ambiguous distance method") 
    x <- as.matrix(x) 
    N <- nrow(x) 
    attrs <- if (method == 6L) 
     list(Size = N, Labels = dimnames(x)[[1L]], Diag = diag, 
      Upper = upper, method = METHODS[method], p = p, call = match.call(), 
      class = "dist") 
    else list(Size = N, Labels = dimnames(x)[[1L]], Diag = diag, 
     Upper = upper, method = METHODS[method], call = match.call(), 
     class = "dist") 
    .Call(C_Cdist, x, method, attrs, p) 
} 
<bytecode: 0x56b0d40> 
<environment: namespace:stats> 

のプログラムです。

Answer 1

Rcppという理由だけで、内部R機能（C/Fortranの）すべての

まず対：キックのために、ここにあなたの距離関数のよりコンパクトなバージョンがありますRcppを使ってアルゴリズムを書いているとはいえ、特にR関数が大量の計算を実行するCまたはFortranルーチンを呼び出す場合、必ずしもR同等を打ち消すことを意味しません。関数が純粋にRで書かれている他のケースでは、Rcppで変換することで所望の速度利得が得られる可能性が高い。

内部機能を書き直すと、非常に狂ったCプログラマーのRコアチームが勝つ可能性が高いことに注意してください。 dist()関数のRソースの最後の行である、

.Call(C_Cdist, x, method, attrs, p) 

：によって示されるように第二に、算出Rが使用する距離はCで行われdist()

の

基本実装。これはテンプレート化されているのではなく、より細分化されているため、C++に比べてわずかな利点があります。

また、C implementation uses OpenMPが利用可能な場合は、計算を並列化します。

提案変形

第三に、わずかに部分集合の順序を変更して追加の変数を作成回避することにより、バージョン間のタイミングは減少。

#include <Rcpp.h> 

// [[Rcpp::export]] 
Rcpp::NumericMatrix calcPWD1 (const Rcpp::NumericMatrix & x){ 
    unsigned int outrows = x.nrow(), i = 0, j = 0; 
    double d; 
    Rcpp::NumericMatrix out(outrows,outrows); 

    for (i = 0; i < outrows - 1; i++){ 
    Rcpp::NumericVector v1 = x.row(i); 
    for (j = i + 1; j < outrows ; j ++){ 
     d = sqrt(sum(pow(v1-x.row(j), 2.0))); 
     out(j,i)=d; 
     out(i,j)=d; 
    } 
    } 

    return out; 
} 

Answer 2

ほとんどでした。しかし、あなたの内側のループボディは、1行であまりにも多くをやろうとしました。テンプレートのプログラミングはそれほど難しくありません。また、コンパイラにより良い機会を与えるために、少しだけ命令を広げる方が良い場合もあります。だから私はちょうど5つのステートメントを作って、immediateltを作りました。

新コード：それを実行

#include <Rcpp.h> 

using namespace Rcpp; 

double dist1 (NumericVector x, NumericVector y){ 
    int n = y.length(); 
    double total = 0; 
    for (int i = 0; i < n ; ++i) { 
    total += pow(x(i)-y(i),2.0); 
    } 
    total = sqrt(total); 
    return total; 
} 

// [[Rcpp::export]] 
NumericMatrix calcPWD (NumericMatrix x){ 
    int outrows = x.nrow(); 
    int outcols = x.nrow(); 
    NumericMatrix out(outrows,outcols); 

    for (int i = 0 ; i < outrows - 1; i++){ 
    for (int j = i + 1 ; j < outcols ; j ++) { 
     NumericVector v1 = x.row(i); 
     NumericVector v2 = x.row(j-1); 
     double d = dist1(v1, v2); 
     out(j-1,i) = d; 
     out(i,j-1)= d; 
    } 
    } 
    return (out) ; 
} 

/*** R 
M <- matrix(log(1:9), 3, 3) 
calcPWD(M) 
*/ 

：

R> sourceCpp("/tmp/mikebrown.cpp") 

R> M <- matrix(log(1:9), 3, 3) 

R> calcPWD(M) 
     [,1]  [,2] [,3] 
[1,] 0.000000 0.740322 0 
[2,] 0.740322 0.000000 0 
[3,] 0.000000 0.000000 0 
R> 

あなたはしかし、あなたのインデックスのロジックをチェックすることもできます。あなたがより多くの比較を見逃したように見えます。

編集：

// [[Rcpp::export]] 
double dist2(NumericVector x, NumericVector y){ 
    double d = sqrt(sum(pow(x - y, 2))); 
    return d; 
}