6

R のバイナリ ベクトルの大きな行列 (600,000 x 500) で、類似性測定値と呼ばれる Dice 係数を計算する必要があります。速度のために、C / Rcpp を使用します。この関数はうまく動作しますが、私はバックグラウンドでコンピューター科学者ではないため、より高速に実行できるかどうかを知りたい. このコードは並列化に適していますが、C コードを並列化した経験はありません。

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

   1    0
1  a    b
0  c    d

サイコロ係数は次のとおりです。(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
4

2 に答える 2

7

Roland が提案した解決策は、私のユース ケースでは完全に満足できるものではありませんでした。arulesそのため、パッケージのソース コードに基づいて、はるかに高速なバージョンを実装します。このコードは、R の関数arulesを使用する Leisch (2005) のアルゴリズムに依存しています。tcrossproduct()

最初に、2 ～ 3 倍高速な Rcpp / RcppEigen バージョンを作成crossprodしました。これは、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")

次に、Dice 係数を計算する小さな 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)
}

この新しい関数は、古い関数よりも 8 倍高速であり、 の関数よりも 3 倍高速ですarules。

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
于 2013-06-05T15:25:41.137 に答える
5

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

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

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

ここに画像の説明を入力

于 2013-06-05T11:47:08.843 に答える