I have two matrices in R that I want to multiply:
a = matrix(rnorm(20*10000, mean=0, sd=5), 20, 10000)
b = matrix(rnorm(20*10000, mean=0, sd=5), 20, 10000)
t(a)%*%b
Given that the dimension in larger this matrix multiplication takes a lot of time, is there a specific way to make computations faster ? And are there any build in functions in R to make such multiplications faster?
To add to cdeterman's answer: You can use eigen's build in parallelization for dense matrix products. In order to do so, you need to compile with open mp activated.
// [[Rcpp::depends(RcppArmadillo, RcppEigen)]]
// [[Rcpp::plugins(openmp)]]
#include <omp.h>
#include <RcppArmadillo.h>
#include <RcppEigen.h>
// [[Rcpp::export]]
SEXP armaMatMult(arma::mat A, arma::mat B){
arma::mat C = A * B;
return Rcpp::wrap(C);
}
// [[Rcpp::export]]
SEXP eigenMatMult(Eigen::MatrixXd A,
Eigen::MatrixXd B,
int n_cores){
Eigen::setNbThreads(n_cores);
//qDebug() << Eigen::nbThreads( );
Eigen::MatrixXd C = A * B;
return Rcpp::wrap(C);
}
// [[Rcpp::export]]
SEXP eigenMapMatMult2(const Eigen::Map<Eigen::MatrixXd> A,
Eigen::Map<Eigen::MatrixXd> B,
int n_cores){
Eigen::setNbThreads(n_cores);
Eigen::MatrixXd C = A * B;
return Rcpp::wrap(C);
}
Here are some benchmarks:
Note that if N = k = 100, parallelization does not necessarily improve performance.
If the matrix dimensions get larger, parallelization starts to have an impact (N = k = 1000):
library(microbenchmark)
# Benchmark 1: N = k = 100
N <- 100
k <- 100
A <- matrix(rnorm(N*k), N, k)
B <- matrix(rnorm(N*k), k, N)
microbenchmark(A%*%B,
armaMatMult2(A, B),
eigenMatMult2(A, B, n_cores = 1),
eigenMatMult2(A, B, n_cores = 2),
eigenMatMult2(A, B, n_cores = 4),
eigenMapMatMult2(A, B, n_cores = 1),
eigenMapMatMult2(A, B, n_cores = 2),
eigenMapMatMult2(A, B, n_cores = 4),
times = 100
# Unit: microseconds
# expr min lq mean median uq max neval
# A %*% B 535.6 540.75 552.594 551.25 554.50 650.2 100
# armaMatMult2(A, B) 542.0 549.10 560.975 556.35 560.25 738.1 100
# eigenMatMult2(A, B, n_cores = 1) 147.1 152.65 159.165 159.65 162.90 180.5 100
# eigenMatMult2(A, B, n_cores = 2) 97.1 109.90 124.496 119.60 127.50 391.8 100
# eigenMatMult2(A, B, n_cores = 4) 71.7 88.15 155.220 115.55 216.95 507.3 100
# eigenMapMatMult2(A, B, n_cores = 1) 139.1 150.10 154.889 154.20 158.35 244.3 100
# eigenMapMatMult2(A, B, n_cores = 2) 93.4 105.70 116.808 113.55 120.40 323.7 100
# eigenMapMatMult2(A, B, n_cores = 4) 66.8 82.60 161.516 196.25 210.40 598.9 100
)
# Benchmark 2: N = k = 1000
N <- 1000
k <- 1000
A <- matrix(rnorm(N*k), N, k)
B <- matrix(rnorm(N*k), k, N)
microbenchmark(A%*%B,
armaMatMult2(A, B),
eigenMatMult2(A, B, n_cores = 1),
eigenMatMult2(A, B, n_cores = 2),
eigenMatMult2(A, B, n_cores = 4),
eigenMapMatMult2(A, B, n_cores = 1),
eigenMapMatMult2(A, B, n_cores = 2),
eigenMapMatMult2(A, B, n_cores = 4),
times = 100
)
Unit: milliseconds
expr min lq mean median uq
A %*% B 597.1293 605.56840 814.52389 665.86650 1025.5896
armaMatMult2(A, B) 603.3894 620.25675 830.98947 693.22355 1078.4853
eigenMatMult2(A, B, n_cores = 1) 131.4696 135.22475 186.69826 193.37870 219.8727
eigenMatMult2(A, B, n_cores = 2) 67.8948 71.71355 114.52759 74.17380 173.3060
eigenMatMult2(A, B, n_cores = 4) 41.8564 48.87075 79.55535 72.00705 106.8572
eigenMapMatMult2(A, B, n_cores = 1) 125.3890 129.26125 175.09933 177.23655 213.0536
eigenMapMatMult2(A, B, n_cores = 2) 62.2866 65.78785 115.74248 79.92470 167.0217
eigenMapMatMult2(A, B, n_cores = 4) 35.2977 40.42480 68.21669 63.13655 97.2571
max neval
1217.6475 100
1446.5127 100
419.2043 100
217.9513 100
139.9629 100
298.2859 100
230.6307 100
118.2553 100
Rcpp codes:
// [[Rcpp::depends(RcppArmadillo, RcppEigen)]]
// [[Rcpp::plugins(openmp)]]
#include <omp.h>
#include <RcppArmadillo.h>
#include <RcppEigen.h>
// [[Rcpp::export]]
SEXP armaMatMult(arma::mat A, arma::mat B){
arma::mat C = A * B;
return Rcpp::wrap(C);
}
// [[Rcpp::export]]
SEXP eigenMatMult(Eigen::MatrixXd A,
Eigen::MatrixXd B,
int n_cores){
Eigen::setNbThreads(n_cores);
//qDebug() << Eigen::nbThreads( );
Eigen::MatrixXd C = A * B;
return Rcpp::wrap(C);
}
// [[Rcpp::export]]
SEXP eigenMapMatMult2(const Eigen::Map<Eigen::MatrixXd> A,
Eigen::Map<Eigen::MatrixXd> B,
int n_cores){
Eigen::setNbThreads(n_cores);
Eigen::MatrixXd C = A * B;
return Rcpp::wrap(C);
}
R codes:
library(microbenchmark)
# Benchmark 1: N = k = 100
N <- 1000
k <- 1000
A <- matrix(rnorm(N*k), N, k)
B <- matrix(rnorm(N*k), k, N)
microbenchmark(A%*%B,
crossprod(A,B),
armaMatMult(A, B),
eigenMatMult(A, B, n_cores = 1),
eigenMatMult(A, B, n_cores = 2),
eigenMatMult(A, B, n_cores = 4),
eigenMapMatMult2(A, B, n_cores = 1),
eigenMapMatMult2(A, B, n_cores = 2),
eigenMapMatMult2(A, B, n_cores = 4),
times = 100)
why crossprod() is the fastest way in my server?
Unit: milliseconds
expr min lq mean median uq max neval
A %*% B 9.630260 9.943458 11.54563 10.10754 11.04614 118.87243 100
crossprod(A, B) 9.668597 9.846900 10.54206 10.07100 11.05562 18.93269 100
armaMatMult(A, B) 12.460832 13.250503 19.07578 20.17903 23.33987 30.40554 100
eigenMatMult(A, B, n_cores = 1) 57.248345 58.737071 60.26937 60.08885 61.47355 71.58064 100
eigenMatMult(A, B, n_cores = 2) 31.465149 32.908495 34.27428 34.33102 35.28012 38.14678 100
eigenMatMult(A, B, n_cores = 4) 18.724495 19.509954 21.49996 20.44093 22.34533 38.44640 100
eigenMapMatMult2(A, B, n_cores = 1) 54.040737 55.628182 57.16926 57.13763 58.41168 67.12156 100
eigenMapMatMult2(A, B, n_cores = 2) 28.965028 30.021924 31.26321 30.87678 32.45385 35.07308 100
eigenMapMatMult2(A, B, n_cores = 4) 16.208926 16.830164 18.32303 17.30694 19.64176 33.19359 100
I am using Microsoft R 4.0. I guess Microsoft optimized the matrix operation.
> sessionInfo()
R version 4.0.2 (2020-06-22)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.5 LTS
Matrix products: default
BLAS: /opt/microsoft/ropen/4.0.2/lib64/R/lib/libRblas.so
LAPACK: /opt/microsoft/ropen/4.0.2/lib64/R/lib/libRlapack.so
locale:
[1] LC_CTYPE=zh_CN.UTF-8 LC_NUMERIC=C LC_TIME=zh_CN.UTF-8 LC_COLLATE=zh_CN.UTF-8
[5] LC_MONETARY=zh_CN.UTF-8 LC_MESSAGES=zh_CN.UTF-8 LC_PAPER=zh_CN.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=zh_CN.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] microbenchmark_1.4-7 miraculix_1.0.5 RandomFieldsUtils_1.0.6 RevoUtils_11.0.2 RevoUtilsMath_11.0.0
loaded via a namespace (and not attached):
[1] Rcpp_1.0.7 lattice_0.20-44 digest_0.6.27 grid_4.0.2
[5] evaluate_0.14 rlang_0.4.11 RcppArmadillo_0.10.6.0.0 Matrix_1.3-4
[9] rmarkdown_2.9 tools_4.0.2 RcppEigen_0.3.3.7.0 tinytex_0.32
[13] nadiv_2.17.1 parallel_4.0.2 xfun_0.24 yaml_2.2.1
[17] compiler_4.0.2 htmltools_0.5.1.1 knitr_1.33