R中的慢点积

Question

我试图从331x23152和23152x23152矩阵中获取点积。

在Python和Octave中这是一个简单的操作，但是在R中这似乎非常慢。

N <- 331 
M <- 23152 

mat_1 = matrix(rnorm(N*M,mean=0,sd=1), N, M) 
mat_2 = matrix(rnorm(N*M,mean=0,sd=1), M, M) 
tm3 <- system.time({ 
    mat_3 = mat_1%*%mat_2 
}) 
print(tm3) 

输出是

user system elapsed 
101.95 0.04 101.99 

换句话说，这点积需要超过100秒来执行。

我正在运行64位R-3.4.0，在带有16 GB RAM的i7-4790上运行RStudio v1.0.143。因此，我不希望这项行动花费这么长时间。

我可以俯视吗？我已经开始研究包bigmemory和bigalgebra，但我不禁想到有一个解决方案，而不必诉诸包。

---

编辑

为了让您有时间差的想法，以下是八度的脚本：

n = 331; 
m = 23152; 

mat_1 = rand(n,m); 
mat_2 = rand(m,m); 
tic 
mat_3 = mat_1*mat_2; 
toc 

输出是

Elapsed time is 3.81038 seconds. 

而且在Python：

import numpy as np 
import time 

n = 331 
m = 23152 

mat_1 = np.random.random((n,m)) 
mat_2 = np.random.random((m,m)) 
tm_1 = time.time() 
mat_3 = np.dot(mat_1,mat_2) 
tm_2 = time.time() 
tm_3 = tm_2 - tm_1 
print(tm_3) 

输出是

2.781277894973755 

正如你所看到的，这些数字都没有，即使在同一个球场。

EDIT 2

宋哲元在李的要求，这里是点积玩具的例子。

在R：

mat_1 = matrix(c(1,2,1,2,1,2), nrow = 2, ncol = 3) 
mat_2 = matrix(c(1,1,1,2,2,2,3,3,3), nrow = 3, ncol = 3) 
mat_3 = mat_1 %*% mat_2 
print(mat_3) 

的输出是：

 [,1] [,2] [,3] 
[1,] 3 6 9 
[2,] 6 12 18 

在八度：

mat_1 = [1,1,1;2,2,2]; 
mat_2 = [1,2,3;1,2,3;1,2,3]; 
mat_3 = mat_1*mat_2 

的输出是：

mat_3 = 

    3 6 9 
    6 12 18 

在Python：

import numpy as np 

mat_1 = np.array([[1,1,1],[2,2,2]]) 
mat_2 = np.array([[1,2,3],[1,2,3],[1,2,3]]) 
mat_3 = np.dot(mat_1, mat_2) 
print(mat_3) 

的输出是：

[[ 3 6 9] 
[ 6 12 18]] 

有关矩阵的点产品的更多信息：https://en.wikipedia.org/wiki/Matrix_multiplication

EDIT 3

为sessionInfo()的输出是：

> sessionInfo() 
R version 3.4.0 (2017-04-21) 
Platform: x86_64-w64-mingw32/x64 (64-bit) 
Running under: Windows 7 x64 (build 7601) Service Pack 1 

Matrix products: default 

locale: 
[1] LC_COLLATE=Dutch_Netherlands.1252 LC_CTYPE=Dutch_Netherlands.1252 LC_MONETARY=Dutch_Netherlands.1252 
[4] LC_NUMERIC=C      LC_TIME=Dutch_Netherlands.1252  

attached base packages: 
[1] stats  graphics grDevices utils  datasets methods base  

loaded via a namespace (and not attached): 
[1] compiler_3.4.0 tools_3.4.0 

EDIT 4

我试过bigalgebra包，但这似乎并没有加快速度：

library('bigalgebra') 

N <- 331 
M <- 23152 

mat_1 = matrix(rnorm(N*M,mean=0,sd=1), N, M) 
mat_1 <- as.big.matrix(mat_1) 
mat_2 = matrix(rnorm(N*M,mean=0,sd=1), M, M) 
tm3 <- system.time({ 
    mat_3 = mat_1%*%mat_2 
}) 
print(tm3) 

输出是：

user system elapsed 
101.79 0.00 101.81 

EDIT 5

詹姆斯建议改变我的随机产生的矩阵：

N <- 331 
M <- 23152 

mat_1 = matrix(runif(N*M), N, M) 
mat_2 = matrix(runif(M*M), M, M) 
tm3 <- system.time({ 
    mat_3 = mat_1%*%mat_2 
}) 
print(tm3) 

的输出是：

user system elapsed 
102.46 0.05 103.00 

Answer 1

根据knb和Zheyuan Li的回复，我开始研究优化的BLAS软件包。我遇到了GotoBlas，OpenBLAS和MKL，例如here。

我的结论是，MKL远远超过默认的BLAS。

看来R必须从源码构建，才能合并MKL。相反，我发现R Open。这有MKL（可选）内置，因此安装非常轻松。

用下面的代码：

N <- 331 
M <- 23152 

mat_1 = matrix(rnorm(N*M,mean=0,sd=1), N, M) 
mat_2 = matrix(rnorm(N*M,mean=0,sd=1), M, M) 
tm3 <- system.time({ 
    mat_3 = mat_1%*%mat_2 
}) 
print(tm3) 

的输出是：

user system elapsed 
    10.61 0.10 3.12 

这样，一个解决这个问题的方法是使用MKL而不是默认BLAS。

但是，经过调查，我的真实生活矩阵非常稀疏。通过使用Matrix软件包，我可以充分利用这一优势。在实践中，我使用它如Matrix(x = mat_1, sparse = TRUE)，其中mat_1将是高度稀疏的矩阵。这将执行时间缩短到3秒左右。

Answer 2

这是一个简单的操作??矩阵乘法在线性代数计算中一直是一个昂贵的操作。

其实我觉得它速度很快。在这个大小的矩阵乘法有

2 * 23.152 * 23.152 * 0.331 = 354.8 GFLOP 

用100秒你的表现是3.5 GFLOPs。请注意，在大多数机器上，性能至多为0.8 GLOP - 2 GFLOP，除非您拥有优化的BLAS库。

如果您认为其他地方的实施更快，请检查使用优化的BLAS或并行计算的可能性。 R使用标准的BLAS来做这件事，而且没有并行性。

---

重要

从R-3.4.0，更多的工具可以与BLAS。

首先，sessionInfo()现在返回链接的BLAS库的完整路径。是的，这并不指向符号链接，而是最终的共享对象！这里的其他答案只是表明了这一点：它有OpenBLAS。

时间结果（在另一个答案中）意味着并行计算（通过OpenBLAS中的多线程）已到位。我很难说出所用线程的数量，但看起来像超线程，因为“系统”的插槽相当大！

二，options现在可以通过matprod设置矩阵乘法的方法。尽管这是为了处理NA/NaN而推出的，但它也提供了性能测试！

• “内部”是未优化的三重循环嵌套中的实现。这是用C编写的，并且与F77中编写的标准（参考）BLAS具有相同的性能;
• “default”，“blas”和“default.simd”表示使用链接的BLAS进行计算，但检查NA和NaN的方法不同。如果R与标准BLAS相关联，那么正如所说的那样，它与“内部”具有相同的性能;但否则我们会看到显着的提振。另请注意，R团队表示将来可能会删除“default.simd”。

Answer 3

我有一个类似的机器：Linux的PC，16 GB内存，英特尔4770K，

从sessionInfo()

R version 3.4.0 (2017-04-21) 
Platform: x86_64-pc-linux-gnu (64-bit) 
Running under: Ubuntu 16.04.2 LTS 

Matrix products: default 
BLAS: /usr/lib/openblas-base/libblas.so.3 
LAPACK: /usr/lib/libopenblasp-r0.2.18.so 

locale: 
[1] LC_CTYPE=en_US.UTF-8  LC_NUMERIC=C    LC_TIME=de_DE.UTF-8  LC_COLLATE=en_US.UTF-8  
[5] LC_MONETARY=de_DE.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=de_DE.UTF-8  LC_NAME=C     
[9] LC_ADDRESS=C    LC_TELEPHONE=C    LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C  

attached base packages: 
[1] stats  graphics grDevices utils  datasets methods base  

other attached packages: 
[1] knitr_1.15.1 clipr_0.3.2 tibble_1.3.0 colorout_1.1-2 

loaded via a namespace (and not attached): 
[1] compiler_3.4.0 tools_3.4.0 Rcpp_0.12.10 

在我的机器相关的输出，您的代码段需要约5秒（开始RStudio，创建的空.R文件，跑片断，输出）：

user system elapsed 
27.608 5.524 4.920 

段：

N <- 331 
M <- 23152 

mat_1 = matrix(rnorm(N*M,mean=0,sd=1), N, M) 
mat_2 = matrix(rnorm(N*M,mean=0,sd=1), M, M) 
tm3 <- system.time({ 
     mat_3 = mat_1 %*% mat_2 
}) 
print(tm3)