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秒才能执行.

我正在运行R-3.4.0 64位,在带有16 GB RAM的i7-4790上运行RStudio v1.0.143.因此,我没想到这个操作需要这么长时间.

我忽略了什么吗？我已经开始研究bigmemory和bigalgebra的软件包了,但我不禁想到有一个解决方案而不必求助于软件包.

---

编辑

为了让您了解时差,这是Octave的脚本:

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

如你所见,这些数字甚至不在同一个球场.

编辑2

在Zheyuan Li的要求下,这里有点产品的玩具示例.

在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

在Octave:

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

编辑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 

编辑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

编辑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

这是一个微不足道的操作？矩阵乘法在线性代数计算中始终是一项昂贵的操作.

其实我觉得它很快.这个大小的矩阵乘法有

2 * 23.152 * 23.152 * 0.331 = 354.8 GFLOP

在100秒内你的表现是3.5 GFLOP.请注意,在大多数机器上,性能最多为0.8 GLOP - 2 GFLOP,除非您有一个优化的BLAS库.

如果您认为其他地方的实施更快,请检查使用优化BLAS或并行计算的可能性.R使用标准BLAS执行此操作,并且没有并行性.

重要

从R-3.4.0开始,BLAS提供了更多工具.

首先,sessionInfo()现在返回链接的BLAS库的完整路径.是的,这并不是指向符号链接,而是最终的共享对象!这里的另一个答案就是这样:它有OpenBLAS.

时序结果(在另一个答案中)意味着并行计算(通过OpenBLAS中的多线程)到位.我很难说出使用的线程数,但看起来像是超线程,因为"系统"的插槽非常大!

其次,options现在可以通过设置矩阵乘法方法matprod.虽然这是为了处理NA/NaN而引入的,但它也提供了性能测试!

• "internal"是非优化的三重循环嵌套中的实现.这是用C语言编写的,与F77编写的标准(参考)BLAS具有相同的性能;
• "default","blas"和"default.simd"表示使用链接的BLAS进行计算,但检查NA和NaN的方式不同.如果R与标准BLAS相关联,那么如上所述,它与"内部"具有相同的性能; 但除此之外,我们看到了显 另请注意,R团队表示将来可能会删除"default.simd".