2013-06-05 33 views
5

Rのバイナリベクトルの大きな行列(600,000 x 500)にわたってダイス係数を呼び出す類似度を計算する必要があります。速度に関してはC/Rcppを使用します。この機能は素晴らしく実行されますが、私はコンピュータ科学者ではないので、より速く実行できるかどうかを知りたいと思います。このコードは並列化に適していますが、Cコードを並列化する経験はありません。C/Rcppのダイス係数の計算速度を向上させる

ダイス係数は、類似性/非類似性の単純な尺度です(どのようにそれをとるかによって異なります)。これは、非対称バイナリベクトルを比較することを目的としており、組み合わせ(通常0-0)は重要ではなく、一致(1-1ペア)は不一致(1-0または0-1ペア)よりも重みが大きいことを意味します。次分割表を想像:

1 0 
1 a b 
0 c d 

サイコロCOEFは次のとおりです。(2 * A)/(2 * A + B + C)ここで

は私のRcpp実装です:

library(Rcpp) 
cppFunction(' 
    NumericMatrix dice(NumericMatrix binaryMat){ 
     int nrows = binaryMat.nrow(), ncols = binaryMat.ncol(); 
     NumericMatrix results(ncols, ncols); 
     for(int i=0; i < ncols-1; i++){ // columns fixed 
      for(int j=i+1; j < ncols; j++){ // columns moving 
       double a = 0; 
       double d = 0; 
       for (int l = 0; l < nrows; l++) { 
        if(binaryMat(l, i)>0){ 
         if(binaryMat(l, j)>0){ 
          a++; 
         } 
        }else{ 
         if(binaryMat(l, j)<1){ 
          d++; 
         } 
        } 
       } 
       // compute Dice coefficient   
       double abc = nrows - d; 
       double bc = abc - a; 
       results(j,i) = (2*a)/(2*a + bc);   
      } 
     } 
     return wrap(results); 
    } 
') 

x <- rbinom(1:200000, 1, 0.5) 
X <- matrix(x, nrow = 200, ncol = 1000) 
system.time(dice(X)) 
    user system elapsed 
    0.814 0.000 0.814 

答えて

6

ローランドによって提案された解決策は、私のユースケースのために完全に満足しませんでした。したがって、arulesパッケージのソースコードに基づいて、はるかに高速なバージョンを実装します。 arulesのコードは、私は2〜3時間も高速ですcrossprodのRcpp/RcppEigenバージョンを書いた、

まずR.でtcrossproduct()機能を使用して(2005)Leischから、アルゴリズムに依存しています。これは、RcppEigenビネットのサンプルコードに基づいています。

library(Rcpp) 
library(RcppEigen) 
library(inline) 
crossprodCpp <- ' 
using Eigen::Map; 
using Eigen::MatrixXi; 
using Eigen::Lower; 

const Map<MatrixXi> A(as<Map<MatrixXi> >(AA)); 

const int m(A.rows()), n(A.cols()); 

MatrixXi AtA(MatrixXi(n, n).setZero().selfadjointView<Lower>().rankUpdate(A.adjoint())); 

return wrap(AtA); 
' 

fcprd <- cxxfunction(signature(AA = "matrix"), crossprodCpp, "RcppEigen") 

次に、ダイス係数を計算するための小さなR関数を書きました。

diceR <- function(X){ 
    a <- fcprd(X) 

nx <- ncol(X) 
rsx <- colSums(X) 

c <- matrix(rsx, nrow = nx, ncol = nx) - a 
# b <- matrix(rsx, nrow = nx, ncol = nx, byrow = TRUE) - a 
b <- t(c) 

m <- (2 * a)/(2*a + b + c) 
return(m) 
} 

この新機能は、より速くarulesで1より〜8時間より速く、古いものと〜3時間未満です。

m <- microbenchmark(dice(X), diceR(X), dissimilarity(t(X), method="dice"), times=100) 
m 
# Unit: milliseconds 
#         expr  min  lq median  uq  max neval 
#        dice(X) 791.34558 809.8396 812.19480 814.6735 910.1635 100 
#        diceR(X) 62.98642 76.5510 92.02528 159.2557 507.1662 100 
# dissimilarity(t(X), method = "dice") 264.07997 342.0484 352.59870 357.4632 520.0492 100 
+0

ニース。時間がある場合は、多分それを少しきれいにして、[Rcpp Gallery](http://gallery.rcpp.org)の投稿にしてください。 –

+0

ありがとう!しましょう。私はgithubに投稿する予定のパッケージを作ります。 –

+0

良い解決策を見つけてくれてうれしいです。あなたの答えを受け入れることを忘れないでください。 – Roland

4

私は仕事であなたの機能を実行することはできませんが、結果はこれと同じですか?

library(arules) 
plot(dissimilarity(X,method="dice")) 

system.time(dissimilarity(X,method="dice")) 
#user system elapsed 
#0.04 0.00 0.04 

enter image description here

+0

結果は同じではありません。 m < - dissimilarity(X、method = "dice")の後に:abs(as.matrix(m) - 1)を指定する必要があります。しかし、それはより速いです。 –

+0

これはほとんど同じタイミングを与えます。 – Roland

+0

私は、m < - 相違度(t(X)、method = "dice")を意味しました。彼らはLeisch(2005)のアルゴリズムを使用します。 –

関連する問題