Rcpp: mon programme de masortingce de distance est plus lent que la fonction du package

Je voudrais calculer la masortingce de distance euclidienne par paires. J’ai écrit les programmes Rcpp par la suggestion de Dirk Eddelbuettel comme suit

NumericMasortingx calcPWD1 (NumericMasortingx x){ int outrows = x.nrow(); double d; NumericMasortingx 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) ; } 

Mais je trouve que mon programme est plus lent que la fonction dist .

 > benchmark(as.masortingx(dist(b)),calcPWD1(b)) test replications elapsed relative user.self sys.self user.child sys.child 1 as.masortingx(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 

Avez-vous des suggestions? Ma masortingce est très simple. Il n’y a pas de nom de colonne ou de nom de ligne, mais simplement une masortingce simple (par exemple, comme b=masortingx(c(rnorm(1000*10)),1000,10) ). Voici le programme de 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) }   

Je m’attends à ce que mon programme soit plus rapide que dist car dans dist , il y a trop de choses à vérifier (comme method , diag ).

Fonctions Rcpp vs. R internes (C / Fortran)

Tout d’abord, ce n’est pas parce que vous écrivez l’algorithme en utilisant Rcpp que l’équivalent R vaincra, en particulier si la fonction R appelle une routine C ou Fortran pour effectuer la majeure partie des calculs. Dans les autres cas où la fonction est écrite uniquement en R, il existe une forte probabilité que sa transformation en Rcpp produise le gain de vitesse souhaité.

N’oubliez pas que lors de la réécriture des fonctions internes, l’un d’eux se heurte à l’équipe de R Core composée de programmeurs C complètement fous qui seront probablement gagnants.

Implémentation de base de dist()

Deuxièmement, le calcul de distance utilisé par R est effectué en C comme indiqué par:

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

, qui est la dernière ligne de la source R de la fonction dist() . Cela lui confère un léger avantage par rapport au C ++, car il est plus granulaire que modélisé.

De plus, l’ implémentation C utilise OpenMP lorsqu’il est disponible pour paralléliser le calcul.

Modification proposée

Troisièmement, en modifiant légèrement l’ordre des sous-ensembles et en évitant de créer une variable supplémentaire, les délais entre les versions diminuent.

 #include  // [[Rcpp::export]] Rcpp::NumericMasortingx calcPWD1 (const Rcpp::NumericMasortingx & x){ unsigned int outrows = x.nrow(), i = 0, j = 0; double d; Rcpp::NumericMasortingx 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; } 

Vous étiez presque là. Mais votre corps de boucle interne a essayé de faire trop dans une ligne. La programmation par modèle est déjà assez difficile et il est parfois préférable d’étaler un peu les instructions pour donner une meilleure chance au compilateur. Donc, je viens de faire cinq déclarations, et construit immédiatement.

Nouveau code:

 #include  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) */ 

Le lancer:

 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> 

Vous voudrez peut-être vérifier votre logique d’indexation. On dirait que vous avez manqué plus de comparaisons.

Edit : Pour les coups de pied, voici une version plus compacte de votre fonction de distance:

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