online - ¿Por qué esta multiplicación de matrices ingenua es más rápida que la base R?

multiplicacion de matrices por un escalar

Question

En R, la multiplicación de matrices está muy optimizada, es decir, es solo una llamada a BLAS / LAPACK. Sin embargo, me sorprende que este código muy ingenuo de C ++ para la multiplicación de matrices vectoriales sea confiable 30% más rápido.

 library(Rcpp)

 # Simple C++ code for matrix multiplication
 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);
 }
 "
 # Compiling
 my_mm = cppFunction(code = mm_code)

 # Simulating data to use
 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 

 # Double checking answer is correct
 max(abs(my_ans - r_ans))
 #> [1] 0

¿El %*% base R realiza algún tipo de verificación de datos que estoy saltando?

EDITAR: Después de entender lo que está pasando (¡gracias SO!), Vale la pena señalar que este es el peor escenario para el %*% R, es decir, matriz por vector. Por ejemplo, @RalfStubner señaló que el uso de RcppArmadillo es incluso más rápido que la implementación ingenua, pero es virtualmente idéntico para la matriz-matriz multiplicada (cuando ambas matrices son grandes y cuadradas):

 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  

Por lo tanto, el %*% actual de R (v3.5.0) es casi óptimo para la matriz-matriz, pero podría acelerarse significativamente para el vector-matriz si está bien omitiendo la comprobación.

Answer 1

La respuesta de Josh explica por qué la multiplicación de la matriz de R no es tan rápida como este enfoque ingenuo. Tenía curiosidad por ver cuánto se podía ganar con RcppArmadillo. El código es bastante 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")

Punto de referencia:

> microbenchmark::microbenchmark(my_mm(m,v), m %*% v, arma_mm(m,v), times = 10)
Unit: milliseconds
          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

Entonces RcppArmadillo nos brinda una sintaxis más agradable y un mejor rendimiento.

La curiosidad sacó lo mejor de mí. Aquí una solución para usar BLAS directamente:

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 <R_ext/BLAS.h>")

Punto de referencia:

Unit: milliseconds
          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 y BLAS (OpenBLAS en mi caso) son casi lo mismo. Y el código BLAS es lo que R hace al final también. Entonces, 2/3 de lo que R hace es la comprobación de errores, etc.

Answer 2

Una rápida ojeada en names.c ( aquí en particular ) lo señala a do_matprod , la función C a la que llama %*% y que se encuentra en el archivo array.c . (Resulta interesante que resulta que tanto crossprod como tcrossprod envían a la misma función). Aquí hay un enlace al código de do_matprod .

Al desplazarse por la función, puede ver que se encarga de una serie de cosas que su implementación ingenua no hace, entre ellas:

1. Mantiene los nombres de las filas y columnas, donde tiene sentido.
2. Permite el envío a métodos S4 alternativos cuando los dos objetos operados por una llamada a %*% son clases para las cuales se han proporcionado dichos métodos. (Eso es lo que está sucediendo en esta parte de la función).
3. Maneja matrices reales y complejas.
4. Implementa una serie de reglas sobre cómo manejar la multiplicación de una matriz y una matriz, un vector y una matriz, una matriz y un vector, y un vector y un vector. (Recuerde que bajo la multiplicación cruzada en R, un vector en el LHS se trata como un vector de fila, mientras que en el RHS, se trata como un vector de columna; este es el código que lo hace así).

Cerca del final de la función , se envía a matprod o matprod . De manera interesante (al menos para mí), en el caso de matrices reales, si una matriz puede contener valores NaN o Inf , entonces matprod envía ( here ) una función llamada simple_matprod que es tan simple y directa como la suya. De lo contrario, se envía a una de un par de rutinas BLAS Fortran que, presumiblemente, son más rápidas, si se pueden garantizar elementos de matriz uniformemente "bien educados".