Accelerare il calcolo del coefficiente di Dice in C/Rcpp

Question

Ho bisogno di calcolare una misura di similitudine chiamare il coefficiente Dice su grandi matrici (600.000 x 500) di vettori binari in R. Per velocità uso C/Rcpp. La funzione funziona alla grande ma, poiché non sono un informatico di base, vorrei sapere se potrebbe funzionare più velocemente. Questo codice è adatto per la parallelizzazione ma non ho esperienza parallelizzando il codice C.

Il coefficiente dei dadi è una semplice misura di somiglianza/dissomiglianza (a seconda di come la si prende). Si intende confrontare i vettori binari asimmetrici, ovvero una delle combinazioni (solitamente 0-0) non è importante e l'accordo (1-1 coppie) ha più peso che disaccordo (coppie 1-0 o 0-1). Immaginate la seguente tabella di contingenza:

1 0 
1 a b 
0 c d 

Il coef dadi è: (2 * a)/(2 * a + b + c)

Qui è la mia realizzazione 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); 
    } 
') 

ed ecco un esempio in esecuzione:

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 

Answer 1

La soluzione proposta da Roland non era del tutto soddisfacente per il mio caso d'uso. Quindi, in base al codice sorgente del pacchetto arules implemento una versione molto più veloce. Il codice in arules si basano su un algoritmo da Leisch (2005) utilizzando la funzione tcrossproduct() in R.

In primo luogo, ho scritto una versione Rcpp/RcppEigen di crossprod che è di 2-3 il tempo più veloce. Questo è basato sul codice di esempio nella vignetta 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") 

Quindi ho scritto una piccola funzione R per calcolare il coefficiente Dice.

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) 
} 

Questa nuova funzione è ~ 8 tempo più veloce di quello vecchio e ~ 3 volte più rapida di quella di 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 

Answer 2

non è possibile eseguire la vostra funzione sul posto di lavoro, ma è il risultato lo stesso di questo?

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

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