Ускорение вычисления коэффициента кости в C/Rcpp

Question

Мне нужно вычислить меру подобия, чтобы коэффициент Кости на больших матрицах (600 000 x 500) двоичных векторов в R. Для скорости я использую C/Rcpp. Функция работает отлично, но поскольку я не компьютерный ученый, я бы хотел знать, сможет ли он работать быстрее. Этот код подходит для параллелизации, но у меня нет опыта, параллельного C-коду.

Коэффициент кости - это простая мера сходства/несходства (в зависимости от того, как вы его принимаете). Предназначен для сравнения асимметричных двоичных векторов, то есть одна из комбинаций (обычно 0-0) не имеет значения, а соглашение (1-1 пары) имеет больший вес, чем несогласие (1-0 или 0-1 пары). Представьте себе следующую таблицу сопряженности:

1 0 
1 a b 
0 c d 

коэффици- Dice является: (2 * а)/(2 * а + Ь + с)

Вот моя реализация 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 

Answer 1

Решение, предложенное Roland не был полностью удовлетворительным для моего случая использования. Поэтому на основе исходного кода из пакета arules я реализую гораздо более быструю версию. Код в arules полагаться на алгоритме из Leisch (2005) с помощью функции tcrossproduct() в R.

Во-первых, я написал версию Rcpp/RcppEigen из crossprod, которое в 2-3 раза быстрее. Это основано на примере кода в виньетке 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) 
} 

Эта новая функция ~ 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 

Answer 2

Я не могу запустить вашу функцию на работе, но результат такой же, как этот?

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

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