Pourquoi cette multiplication matricielle naïve est-elle plus rapide que celle de base R?

Question

En R, la multiplication de matrices est très optimisée, c'est-à-dire qu'il s'agit vraiment juste d'un appel à BLAS/LAPACK. Cependant, je suis surpris que ce code C++ très naïf pour la multiplication matrice-vecteur semble systématiquement 30% plus rapide.

 library(Rcpp)

 # Code C++ simple pour la multiplication de matrices
 mm_code = 
 "NumericVector my_mm(NumericMatrix m, NumericVector v){
   int nRow = m.rows();
   int nCol = m.cols();
   NumericVector ans(nRow);
   double v_j;
   for(int j = 0; j < nCol; j++){
     v_j = v[j];
     for(int i = 0; i < nRow; i++){
       ans[i] += m(i,j) * v_j;
     }
   }
   return(ans);
 }
 "
 # Compilation
 my_mm = cppFunction(code = mm_code)

 # Simulation de données à utiliser
 nRow = 10^4
 nCol = 10^4

 m = matrix(rnorm(nRow * nCol), nrow = nRow)
 v = rnorm(nCol)

 system.time(my_ans <- my_mm(m, v))
#>    user  system elapsed 
#>   0.103   0.001   0.103 
 system.time(r_ans <- m %*% v)
#>   user  system elapsed 
#>  0.154   0.001   0.154 

 # Vérification que la réponse est correcte
 max(abs(my_ans - r_ans))
 #> [1] 0

Est-ce que le %*% de base de R effectue une sorte de vérification des données que je suis en train de sauter ?

MODIFICATION:

Après avoir compris ce qui se passait (merci SO!), il convient de noter qu'il s'agit d'un scénario défavorable pour le %*% de R, c'est-à-dire matrice par vecteur. Par exemple, @RalfStubner a souligné qu'utiliser une implémentation RcppArmadillo d'une multiplication matrice-vecteur est encore plus rapide que l'implémentation naïve que j'ai démontrée, impliquant une considérablement plus rapide que R de base, mais est pratiquement identique au %*% de base de R pour la multiplication de matrices (lorsque les deux matrices sont grandes et carrées) :

 arma_code <- 
   "arma::mat arma_mm(const arma::mat& m, const arma::mat& m2) {
 return m * m2;
 };"
 arma_mm = cppFunction(code = arma_code, depends = "RcppArmadillo")

 nRow = 10^3 
 nCol = 10^3

 mat1 = matrix(rnorm(nRow * nCol), 
               nrow = nRow)
 mat2 = matrix(rnorm(nRow * nCol), 
               nrow = nRow)

 system.time(arma_mm(mat1, mat2))
#>   user  system elapsed 
#>   0.798   0.008   0.814 
 system.time(mat1 %*% mat2)
#>   user  system elapsed 
#>   0.807   0.005   0.822  

Donc le %*% actuel de R (v3.5.0) est presque optimal pour les matrices-matrices, mais pourrait être significativement accéléré pour les matrices-vecteurs si vous acceptez de passer outre la vérification.

Answer 1

Un coup d'œil rapide dans names.c (ici en particulier) vous renvoie à la fonction C do_matprod, appelée par %*% et présente dans le fichier array.c. (Il s'avère, de manière intéressante, que les fonctions crossprod et tcrossprod dispatchent également vers cette même fonction). Voici un lien vers le code de do_matprod.

En parcourant la fonction, vous pouvez voir qu'elle gère un certain nombre d'aspects que votre implémentation naïve ne prend pas en compte, notamment :

1. Garde les noms de lignes et de colonnes, lorsque cela a du sens.
2. Permet de dispatcher vers des méthodes S4 alternatives lorsque les deux objets sur lesquels s'applique un appel à %*% sont de classes pour lesquelles de telles méthodes ont été fournies. (C'est ce qui se passe dans cette partie de la fonction.)
3. Prend en charge les matrices réelles et complexes.
4. Met en œuvre une série de règles pour manipuler la multiplication d'une matrice par une matrice, d'un vecteur par une matrice, d'une matrice par un vecteur, et d'un vecteur par un vecteur. (Rappelez-vous qu'en multiplication croisée en R, un vecteur à gauche est traité comme un vecteur ligne, tandis qu'à droite, il est traité comme un vecteur colonne ; c'est le code qui le permet.)

Près de la fin de la fonction, il dispatche vers soit matprod ou ou cmatprod. De manière intéressante (en tout cas pour moi), dans le cas de matrices réelles, si l'une des matrices peut contenir des valeurs NaN ou Inf, alors matprod dispatche (ici) vers une fonction appelée simple_matprod aussi simple et directe que la vôtre. Sinon, il dispatche vers l'une des deux routines BLAS Fortran qui, on le suppose, sont plus rapides, si des éléments de matrice uniformément 'bien comportés' peuvent être garantis.

Answer 2

La réponse de Josh explique pourquoi la multiplication de matrices de R n'est pas aussi rapide que cette approche naïve. J'étais curieux de voir combien on pourrait gagner en utilisant RcppArmadillo. Le code est assez simple :

arma_code <- 
  "arma::vec arma_mm(const arma::mat& m, const arma::vec& v) {
       return m * v;
   };"
arma_mm = cppFunction(code = arma_code, depends = "RcppArmadillo")

Comparaison :

> microbenchmark::microbenchmark(my_mm(m,v), m %*% v, arma_mm(m,v), times = 10)
Unité: millisecondes
          expr      min       lq      mean    median        uq       max neval
   my_mm(m, v) 71.23347 75.22364  90.13766  96.88279  98.07348  98.50182    10
       m %*% v 92.86398 95.58153 106.00601 111.61335 113.66167 116.09751    10
 arma_mm(m, v) 41.13348 41.42314  41.89311  41.81979  42.39311  42.78396    10

Ainsi, RcppArmadillo nous offre une syntaxe plus agréable et de meilleures performances.

La curiosité l'a emporté sur moi. Voici une solution pour utiliser BLAS directement :

blas_code = "
NumericVector blas_mm(NumericMatrix m, NumericVector v){
  int nRow = m.rows();
  int nCol = m.cols();
  NumericVector ans(nRow);
  char trans = 'N';
  double one = 1.0, zero = 0.0;
  int ione = 1;
  F77_CALL(dgemv)(&trans, &nRow, &nCol, &one, m.begin(), &nRow, v.begin(),
           &ione, &zero, ans.begin(), &ione);
  return ans;
}"
blas_mm <- cppFunction(code = blas_code, includes = "#include ")

Comparaison :

Unité: millisecondes
          expr      min       lq      mean    median        uq       max neval
   my_mm(m, v) 72.61298 75.40050  89.75529  96.04413  96.59283  98.29938    10
       m %*% v 95.08793 98.53650 109.52715 111.93729 112.89662 128.69572    10
 arma_mm(m, v) 41.06718 41.70331  42.62366  42.47320  43.22625  45.19704    10
 blas_mm(m, v) 41.58618 42.14718  42.89853  42.68584  43.39182  44.46577    10

Armadillo et BLAS (OpenBLAS dans mon cas) sont presque identiques. Et le code BLAS est ce que R fait à la fin aussi. Ainsi, 2/3 de ce que fait R est la vérification des erreurs, etc.

Answer 3

Pour ajouter un autre point à la solution de Ralf Stubner, vous pouvez utiliser la version C++ suivante pour

1. faire plusieurs colonnes en même temps pour éviter de relire le vecteur de sortie plusieurs fois.

2. ajouter __restrict__ pour potentiellement permettre des opérations vectorielles (probablement pas important ici car il s'agit seulement de lectures, je suppose).

   include

   using namespace Rcpp;

   inline void mat_vec_mult_vanilla (double const restrict m, double const restrict v, double restrict const res, size_t const dn, size_t const dm) noexcept { for(size_t j = 0; j < dm; ++j, ++v){ double r = res; for(size_t i = 0; i < dn; ++i, ++r, ++m) r += m v; } }

   inline void mat_vec_mult (double const restrict const m, double const restrict const v, double restrict const res, size_t const dn, size_t const dm) noexcept { size_t j(0L); double const vj = v,

   • mi = m; constexpr size_t const ncl(8L); { double const mvals[ncl]; size_t const end_j = dm - (dm % ncl), inc = ncl dn; for(; j < end_j; j += ncl, vj += ncl, mi += inc){ double r = res; mvals[0] = mi; for(size_t i = 1; i < ncl; ++i) mvals[i] = mvals[i - 1L] + dn; for(size_t i = 0; i < dn; ++i, ++r) for(size_t ii = 0; ii < ncl; ++ii) r += (vj + ii) *mvals[ii]++; } }

     mat_vec_mult_vanilla(mi, vj, res, dn, dm - j); }

   // [[Rcpp::export("mat_vec_mult", rng = false)]] NumericVector mat_vec_mult_cpp(NumericMatrix m, NumericVector v){ size_t const dn = m.nrow(), dm = m.ncol(); NumericVector res(dn); mat_vec_mult(&m[0], &v[0], &res[0], dn, dm); return res; }

   // [[Rcpp::export("mat_vec_mult_vanilla", rng = false)]] NumericVector mat_vec_mult_vanilla_cpp(NumericMatrix m, NumericVector v){ size_t const dn = m.nrow(), dm = m.ncol(); NumericVector res(dn); mat_vec_mult_vanilla(&m[0], &v[0], &res[0], dn, dm); return res; }

Le résultat avec -O3 dans mon fichier Makevars et gcc-8.3 est

set.seed(1)
dn <- 10001L
dm <- 10001L
m <- matrix(rnorm(dn * dm), dn, dm)
lv <- rnorm(dm)

all.equal(drop(m %*% lv), mat_vec_mult(m = m, v = lv))
#R> [1] TRUE
all.equal(drop(m %*% lv), mat_vec_mult_vanilla(m = m, v = lv))
#R> [1] TRUE

bench::mark(
  R              = m %*% lv, 
  `Version de l'OP` = my_mm(m = m, v = lv), 
  `BLAS`         = blas_mm(m = m, v = lv),
  `C++ vanilla`  = mat_vec_mult_vanilla(m = m, v = lv), 
  `C++`          = mat_vec_mult(m = m, v = lv), check = FALSE)
#R> # Un tibble: 5 x 13
#R>   expression        min   median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc total_time result memory                 time          gc               
#R>                                                            
#R> 1 R             147.9ms    151ms      6.57    78.2KB        0     4     0      609ms      
#R> 2 Version de l'OP   56.9ms   57.1ms     17.4     78.2KB        0     9     0      516ms      
#R> 3 BLAS           90.1ms   90.7ms     11.0     78.2KB        0     6     0      545ms      
#R> 4 C++ vanilla    57.2ms   57.4ms     17.4     78.2KB        0     9     0      518ms      
#R> 5 C++              51ms   51.4ms     19.3     78.2KB        0    10     0      519ms    

Une légère amélioration. Le résultat peut toutefois dépendre beaucoup de la version de BLAS. La version que j'ai utilisée est

sessionInfo()
#R> #...
#R> Matrix products: default
#R> BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
#R> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
#R> ...

Le fichier entier que j'ai Rcpp::sourceCpp()é est

#include 
#include 
using namespace Rcpp;

inline void mat_vec_mult_vanilla
(double const * __restrict__ m, 
 double const * __restrict__ v, 
 double * __restrict__ const res, 
 size_t const dn, size_t const dm) noexcept {
  for(size_t j = 0; j < dm; ++j, ++v){
    double * r = res;
    for(size_t i = 0; i < dn; ++i, ++r, ++m)
      *r += *m * *v;
  }
}

inline void mat_vec_mult
(double const * __restrict__ const m, 
 double const * __restrict__ const v, 
 double * __restrict__ const res, 
 size_t const dn, size_t const dm) noexcept {
  size_t j(0L);
  double const * vj = v,
               * mi = m;
  constexpr size_t const ncl(8L);
  {
    double const * mvals[ncl];
    size_t const end_j = dm - (dm % ncl),
                   inc = ncl * dn;
    for(; j < end_j; j += ncl, vj += ncl, mi += inc){
      double *r = res;
      mvals[0] = mi;
      for(size_t i = 1; i < ncl; ++i)
        mvals[i] = mvals[i - 1L] + dn;
      for(size_t i = 0; i < dn; ++i, ++r)
        for(size_t ii = 0; ii < ncl; ++ii)
          *r += *(vj + ii) * *mvals[ii]++;
    }
  }

  mat_vec_mult_vanilla(mi, vj, res, dn, dm - j);
}

// [[Rcpp::export("mat_vec_mult", rng = false)]]
NumericVector mat_vec_mult_cpp(NumericMatrix m, NumericVector v){
  size_t const dn = m.nrow(), 
               dm = m.ncol();
  NumericVector res(dn);
  mat_vec_mult(&m[0], &v[0], &res[0], dn, dm);
  return res;
}

// [[Rcpp::export("mat_vec_mult_vanilla", rng = false)]]
NumericVector mat_vec_mult_vanilla_cpp(NumericMatrix m, NumericVector v){
  size_t const dn = m.nrow(), 
               dm = m.ncol();
  NumericVector res(dn);
  mat_vec_mult_vanilla(&m[0], &v[0], &res[0], dn, dm);
  return res;
}

// [[Rcpp::export(rng = false)]]
NumericVector my_mm(NumericMatrix m, NumericVector v){
  int nRow = m.rows();
  int nCol = m.cols();
  NumericVector ans(nRow);
  double v_j;
  for(int j = 0; j < nCol; j++){
    v_j = v[j];
    for(int i = 0; i < nRow; i++){
      ans[i] += m(i,j) * v_j;
    }
  }
  return(ans);
}

// [[Rcpp::export(rng = false)]]
NumericVector blas_mm(NumericMatrix m, NumericVector v){
  int nRow = m.rows();
  int nCol = m.cols();
  NumericVector ans(nRow);
  char trans = 'N';
  double one = 1.0, zero = 0.0;
  int ione = 1;
  F77_CALL(dgemv)(&trans, &nRow, &nCol, &one, m.begin(), &nRow, v.begin(),
           &ione, &zero, ans.begin(), &ione);
  return ans;
}

/*** R
set.seed(1)
dn <- 10001L
dm <- 10001L
m <- matrix(rnorm(dn * dm), dn, dm)
lv <- rnorm(dm)

all.equal(drop(m %*% lv), mat_vec_mult(m = m, v = lv))
all.equal(drop(m %*% lv), mat_vec_mult_vanilla(m = m, v = lv))

bench::mark(
  R              = m %*% lv, 
  `Version de l'OP` = my_mm(m = m, v = lv), 
  `BLAS`         = blas_mm(m = m, v = lv),
  `C++ vanilla`  = mat_vec_mult_vanilla(m = m, v = lv), 
  `C++`          = mat_vec_mult(m = m, v = lv), check = FALSE)
*/