R: Calcul produit cartésien rapide de deux matrices numériques

Question

J'ai deux grandes matrices numériques et je veux calculer leur produit cartésien dans R. Y a-t-il un moyen de le faire avec des performances plus élevées et moins de mémoire qu'avec mon approche actuelle?

EDIT: J'ai ajouté une version Rcpp, qui fonctionne déjà beaucoup mieux que ma première approche-R. Comme je ne suis pas expérimenté avec Rcpp ou RcppArmadillo: Y at-il une manière plus rapide/plus standardisée d'écrire cette fonction Rcpp?

m1 <- matrix(sample(0:9, size=110000, replace = TRUE), ncol = 110) 
m2 <- matrix(sample(0:9, size=110000, replace = TRUE), ncol = 110) 

#Current approach: 
m3 <- apply(m1, 1, function(x) x * t(m2)) 
matrix(m3, ncol = 110, byrow = TRUE) 

#EDIT - Rcpp approach 
library(Rcpp) 
#assuming ncol(m1) == ncol(m2) 
cppFunction('IntegerMatrix cartProd(IntegerMatrix m1, IntegerMatrix m2) { 
    int nrow1 = m1.nrow(), ncol = m1.ncol(), nrow2 = m2.nrow(); 
    int orow = 0; 
    IntegerMatrix out(nrow1 * nrow2, ncol); 

    for (int r1 = 0; r1 < nrow1; r1++) { 
    for (int r2 = 0; r2 < nrow2; r2++) {  
     for (int c = 0; c < ncol; c++){ 
     out(orow, c) = m1(r1, c) * m2(r2, c); 
     } 
     orow++; 
    } 
    } 
    return out; 
}') 
m5 <- cartProd(m1, m2) 

Answer 1

La meilleure approche que vous avez émis l'hypothèse est d'utiliser C++ pour exécuter le produit cartésien que vous désirez. Essayer de porter le code sur Armadillo donnera une légère augmentation par rapport à la version pure de Rcpp, qui est nettement plus rapide que la version R écrite. Pour plus de détails sur l'efficacité de chaque méthode, voir la section sur les points de référence à la fin.

---

La première version est presque un port direct dans tatou et exécute en fait légèrement pire que fonction pure CRPP initiale. La seconde version utilise les fonctions submatrix views et each_row() de l'armadillo pour exploiter l'évaluation en place. Pour obtenir la parité avec la version Rcpp, notez l'utilisation de la référence par renvoi et l'utilisation d'un type entier signé produisant const arma::imat&. Ceci évite une copie en profondeur des deux grandes matrices integer lorsque les types correspondent et une référence est établie.

#include <RcppArmadillo.h> 
// [[Rcpp::depends(RcppArmadillo)] 

// --- Version 1 

// [[Rcpp::export]] 
arma::imat cartProd_arma(const arma::imat& m1, const arma::imat& m2) { 
    int nrow1 = m1.n_rows, ncol = m1.n_cols, nrow2 = m2.n_rows, orow = 0; 
    arma::imat out(nrow1 * nrow2, ncol); 

    for (int r1 = 0; r1 < nrow1; ++r1) { 
    for (int r2 = 0; r2 < nrow2; ++r2) { 
     out.row(orow) = m1.row(r1) % m2.row(r2); 
     orow++; 
    } 
    } 
    return out; 
} 

// --- Version 2 

// [[Rcpp::export]] 
arma::imat cartProd_arma2(const arma::imat& m1, const arma::imat& m2) { 
    int nrow1 = m1.n_rows, ncol = m1.n_cols, nrow2 = m2.n_rows, orow = 0; 
    arma::imat out(nrow1 * nrow2, ncol); 

    for (int r1 = 0; r1 < nrow1; ++r1) { 
    out.submat(orow, 0, orow + nrow2 - 1, ncol - 1) = m1.row(r1) % m2.each_row(); 
    orow += nrow2; 
    } 
    return out; 
} 

---

vérification rapide des détails de mise en œuvre correspondant au produit initial

all.equal(cartProd(m1, m2), cartProd_arma(m1, m2)) 
# [1] TRUE 
all.equal(cartProd(m1, m2), cartProd_arma2(m1, m2)) 
# [1] TRUE 

---

Pour générer les repères, je rangea la fonction initiale légèrement en pré-transposition de la matrice pour éviter plusieurs Transpose les appels chaque fois que l'application a été appelée par ligne. De plus, j'ai inclus l'astuce de la fonction affichée par @ user20650.

# OP's initial R only solution with slight modifications 
op_R = function(m1, m2){ 
    m2 <- t(m2) 
    m3 <- matrix(apply(m1, 1, function(x) x * m2), ncol = ncol(m1), byrow = TRUE) 
} 

# user20650's comment 
so_comment <- function(m1, m2){ 
    m4 <- matrix(rep(t(m1), each=nrow(m1)) * c(m2), ncol=nrow(m1)) 
} 

En conséquence, nous avons la, nous pouvons voir microbenchmark

library("microbenchmark") 
out <- microbenchmark(op_r = op_R(m1, m2), so_comment_r = so_comment(m1, m2), 
       rcpp = cartProd(m1, m2), arma_v1 = cartProd_arma(m1, m2), 
       arma_v2 = cartProd_arma2(m1, m2), 
       times = 50) 
out 
# Unit: milliseconds 
#   expr  min  lq  mean median  uq  max neval 
#   op_r 1615.6572 1693.0526 1793.0515 1771.7353 1886.0988 2053.7050 50 
# so_comment_r 2778.0971 2856.6429 2950.5837 2936.7459 3021.4249 3344.4401 50 
#   rcpp 463.6743 482.3118 565.0525 582.1660 614.3714 699.3516 50 
#  arma_v1 597.9004 620.5888 713.4101 726.7572 783.4225 820.3770 50 
#  arma_v2 384.7205 401.9744 490.5118 503.5007 574.6840 622.9876 50 

Ainsi, de cela, suivant que le cartProd_arma2, la mise en œuvre de tatou sous-matrice, est la meilleure fonction suivie de près par cartProd, la implémentation Rcpp pure.