我有两个 R 语言中的矩阵需要相乘：
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
考虑到维度较大时，矩阵乘法需要很长时间，有没有特定的方法可以使计算更快？并且在R中是否有内置函数可以加速这样的乘法？
3个回答
36
根据您的代码、努力和硬件，有许多方法可以解决这个问题。
- 使用最适合该任务的函数
最简单的方法是使用 crossprod，它与 t(a)%*% b 相同（注意 - 这只会带来少量速度提升）。
crossprod(a,b)
- 使用Rcpp（以及可能的RcppEigen/RcppArmadillo）。
C++语言可以显著提高代码的速度。使用线性代数库也有助于进一步提高速度（因此选择Eigen和Armadillo）。但这假设您愿意编写一些C++代码。
- 使用更好的后端
在此之后，您将查看BLAS后端，例如OpenBLAS、Atlas等。连接这些到R取决于您的操作系统。如果您使用像Ubuntu这样的Debian系统，则非常容易。您可以在此处找到演示here。这些有时可以通过类似Armadillo和Eigen的库进一步利用。
- GPU计算
如果您有GPU（例如AMD、NVIDIA等），则可以利用其中的许多核心大大加快计算速度。其中一些可能会有用，包括gpuR、gputools和gmatrix
编辑 - 为了回应@jenesaiquoi的评论，解释了使用Rcpp的好处
test.cpp
// [[Rcpp::depends(RcppArmadillo, RcppEigen)]]
#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){
Eigen::MatrixXd C = A * B;
return Rcpp::wrap(C);
}
// [[Rcpp::export]]
SEXP eigenMapMatMult(const Eigen::Map<Eigen::MatrixXd> A, Eigen::Map<Eigen::MatrixXd> B){
Eigen::MatrixXd C = A * B;
return Rcpp::wrap(C);
}
test.R
library(Rcpp)
A <- matrix(rnorm(10000), 100, 100)
B <- matrix(rnorm(10000), 100, 100)
library(microbenchmark)
sourceCpp("test.cpp")
microbenchmark(A%*%B, armaMatMult(A, B), eigenMatMult(A, B), eigenMapMatMult(A, B))
Unit: microseconds
expr min lq mean median uq max neval
A %*% B 885.846 892.1035 933.7457 901.1010 938.9255 1411.647 100
armaMatMult(A, B) 846.688 857.6320 915.0717 866.2265 893.7790 1421.557 100
eigenMatMult(A, B) 205.978 208.1295 233.1882 217.0310 229.4730 369.369 100
eigenMapMatMult(A, B) 192.366 194.9835 207.1035 197.5405 205.2550 366.945 100
- cdeterman
16
3我使用了microbenchmark来与
crossprod进行比较，但不幸的是并没有太大的速度提升。 - RaadeigenMapMatMult 总是更好，还是要根据具体情况而定？在我的情况下，速度提高了约8倍，这非常好。 （谢谢！） - generic_user对于提供的方法，我相信它总是更快的，因为它避免了复制。唯一的例外可能是使用不同的线性代数后端，比如OpenBLAS。那么armadillo可能会节省一些时间，但这需要在您特定的机器上进行更多的基准测试。 - cdeterman
@Mooks 我对MKL获胜并不感到意外。这是一个高度优化的数学库，Microsoft（Revolution R）努力将其与R良好地集成在一起。我不熟悉细节，但我认为它很可能利用了类似于“映射”特征结构以及并行化的方面。 - cdeterman
显示剩余11条评论
6
为cdeterman的回答做一个补充：
你可以利用eigen内置的并行化功能来进行稠密矩阵乘法。为了实现这个功能，你需要使用已激活open mp的编译器。
// [[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);
}
以下是一些基准测试结果:
请注意，如果N = k = 100，并行化不一定会提高性能。
当矩阵的尺寸变大时， 并行化开始产生影响（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
- A.Fischer
-1
Rcpp 代码：
// [[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 代码：
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)
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
我正在使用 Microsoft R 4.0。 我猜 Microsoft 优化了矩阵运算。
> 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
- Sheng Luan
3
请添加更多细节以扩展您的答案，例如工作代码或文档引用。 - Community
如果您有新的问题，请通过点击“提问”按钮来提出。如果有必要，可以包含此问题的链接以提供上下文。 - Rosalie Bruel
网页内容由stack overflow 提供, 点击上面的可以查看英文原文，