Почему это наивное матричное умножение быстрее, чем базовое R?

Question

В R умножение матриц очень оптимизировано, то есть на самом деле это просто вызов BLAS / LAPACK. Тем не менее, я удивлен, что этот очень наивный C ++ код для умножения матрицы на вектор кажется надежно на 30% быстрее.

 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

Базис R %*% выполняет какую-либо проверку данных, которую я пропускаю?

EDIT:

После понимания того, что происходит (спасибо ТАК!), Стоит отметить, что это худший вариант для R %*%, то есть матрица за вектором. Например, @RalfStubner указал, что использование RcppArmadillo реализации умножения матрицы на вектор даже быстрее, чем наивная реализация, которую я продемонстрировал, подразумевая значительно быстрее, чем базовая R, но практически идентична базовой R * %*% для матрицы-матрицы. умножить (когда обе матрицы большие и квадратные):

 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  

Таким образом, ток R (v3.5.0) %*% почти оптимален для матрицы-матрицы, но может быть значительно ускорен для матрицы-вектора, если вы в порядке, пропуская проверку.

Answer 1

Ответ Джоша объясняет, почему умножение матрицы R не так быстро, как этот наивный подход. Мне было любопытно посмотреть, сколько можно получить, используя RcppArmadillo. Код достаточно прост:

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

Benchmark:

> 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

Таким образом, RcppArmadillo дает нам более хороший синтаксис и лучшую производительность.

Любопытство одолело меня. Вот решение для использования BLAS напрямую:

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

Benchmark:

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 и BLAS (в моем случае OpenBLAS) почти одинаковы. И код BLAS - это то, что R делает в конце концов. Итак, 2/3 того, что делает R, это проверка ошибок и т. Д.

Answer 2

Быстрый взгляд на names.c ( здесь, в частности ) указывает на do_matprod, функцию C, которая вызывается %*% и находится в файле array.c. (Интересно, что оказывается, что и crossprod, и tcrossprod отправляют и той же функции). Вот ссылка на код do_matprod.

Прокручивая функцию, вы можете увидеть, что она заботится о многих вещах, которых не делает ваша наивная реализация, включая:

1. Сохраняет имена строк и столбцов, где это имеет смысл.
2. Позволяет отправлять альтернативные методы S4, когда два объекта, управляемые вызовом %*%, относятся к классам, для которых были предоставлены такие методы. (Это то, что происходит в этой части функции.)
3. Обрабатывает как действительные, так и сложные матрицы.
4. Реализует ряд правил для обработки умножения матрицы и матрицы, вектора и матрицы, матрицы и вектора, а также вектора и вектора. (Напомним, что при перекрестном умножении в R вектор на LHS обрабатывается как вектор строки, тогда как на RHS он рассматривается как вектор столбца; это код, который делает это так.)

Ближе к концу функции отправляется либо в matprod, либо в cmatprod. Интересно (по крайней мере мне), в случае реальных матриц, , если любая матрица может содержать NaN или Inf значений, то matprod отправляет ( здесь ) в функция с именем simple_matprod примерно такая же простая и понятная, как и ваша. В противном случае он отправляет одну из пары подпрограмм BLAS Fortran, которые, предположительно, работают быстрее, если можно гарантировать равномерно «хорошо себя ведущие» матричные элементы.