我在 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 中是否有任何内置函数可以使这种乘法更快？
【问题讨论】：
根据您的代码、工作量和硬件，有很多方法可以解决这个问题。
最简单的方法是使用crossprod，它与t(a)%*% b 相同（注意——这只会稍微提高速度）
crossprod(a,b)
C++ 可能会更大程度地提高您的代码速度。使用线性代数库也可能会进一步帮助这一点（因此 Eigen 和 Armadillo）。然而，这假设您愿意编写一些 C++。
在此之后，您将查看 BLAS 后端，例如 OpenBLAS、Atlas 等。将这些后端连接到 R 取决于您的操作系统。如果您使用的是像 Ubuntu 这样的 Debian 系统，这很容易。你可以找到一个演示here。这些有时可以被 Armadillo 和 Eigen 等库进一步利用。
如果您有 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);
}
测试.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
【讨论】：
crossprod 进行比较，不幸的是，速度并没有太大提升。
添加到 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
【讨论】：
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)
为什么 crossprod() 是我服务器中最快的方法？
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。我猜微软优化了矩阵运算。
> 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