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".