R中的慢点积
Slow dot product in R
我正在尝试从 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
应李哲元的要求，这里有点积的玩具示例。
在 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
编辑 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
James 建议更改我随机生成的矩阵：
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
这是一个微不足道的操作？？矩阵乘法在线性代数计算中一直是一个昂贵的运算。
其实我觉得挺快的。这种大小的矩阵乘法有
2 * 23.152 * 23.152 * 0.331 = 354.8 GFLOP
在 100 秒内，您的性能为 3.5 GFLOPs。请注意，在大多数机器上，性能最多为 0.8 GLOPs - 2 GFLOPs，除非您有优化的 BLAS 库。
如果您认为其他地方的实施速度更快，请检查使用优化的 BLAS 或并行计算的可能性。 R 使用标准 BLAS 执行此操作，没有并行性。
重要
从 R-3.4.0 开始，BLAS 提供了更多工具。
首先，sessionInfo() 现在 returns linked BLAS 库的 完整路径 。是的，这不是指向符号link，而是最终的共享对象！这里的另一个答案只是说明了这一点：它有 OpenBLAS。
计时结果（在另一个答案中）暗示并行计算（通过 OpenBLAS 中的多线程）已经到位。我很难说出使用的线程数，但看起来超线程已打开，因为 "system" 的插槽相当大！
其次，options 现在可以通过 matprod 设置矩阵乘法方法。虽然这是为了处理 NA / NaN 而引入的，但它也提供了性能测试！
- "internal" 是非优化三重循环嵌套中的实现。这是用 C 编写的，与用 F77 编写的标准（参考）BLAS 具有同等性能；
- "default"、"blas"和"default.simd"表示使用linked BLAS进行计算，但检查NA和NaN的方式不同。如果将 R linked 为标准 BLAS，则如前所述，它与 "internal" 具有相同的性能；但除此之外，我们看到了显着的提升。另请注意，R 团队表示 "default.simd" 将来可能会被删除。
我有一台类似的机器：Linux PC，16 GB RAM，intel 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)
根据knb和李哲元的回复，我开始研究优化的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 秒。
我正在尝试从 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
应李哲元的要求，这里有点积的玩具示例。
在 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
编辑 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
James 建议更改我随机生成的矩阵：
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
这是一个微不足道的操作？？矩阵乘法在线性代数计算中一直是一个昂贵的运算。
其实我觉得挺快的。这种大小的矩阵乘法有
2 * 23.152 * 23.152 * 0.331 = 354.8 GFLOP
在 100 秒内，您的性能为 3.5 GFLOPs。请注意，在大多数机器上，性能最多为 0.8 GLOPs - 2 GFLOPs，除非您有优化的 BLAS 库。
如果您认为其他地方的实施速度更快，请检查使用优化的 BLAS 或并行计算的可能性。 R 使用标准 BLAS 执行此操作，没有并行性。
重要
从 R-3.4.0 开始，BLAS 提供了更多工具。
首先，sessionInfo() 现在 returns linked BLAS 库的 完整路径 。是的，这不是指向符号link，而是最终的共享对象！这里的另一个答案只是说明了这一点：它有 OpenBLAS。
计时结果（在另一个答案中）暗示并行计算（通过 OpenBLAS 中的多线程）已经到位。我很难说出使用的线程数，但看起来超线程已打开，因为 "system" 的插槽相当大！
其次，options 现在可以通过 matprod 设置矩阵乘法方法。虽然这是为了处理 NA / NaN 而引入的，但它也提供了性能测试！
- "internal" 是非优化三重循环嵌套中的实现。这是用 C 编写的，与用 F77 编写的标准（参考）BLAS 具有同等性能；
- "default"、"blas"和"default.simd"表示使用linked BLAS进行计算，但检查NA和NaN的方式不同。如果将 R linked 为标准 BLAS，则如前所述，它与 "internal" 具有相同的性能；但除此之外，我们看到了显着的提升。另请注意，R 团队表示 "default.simd" 将来可能会被删除。
我有一台类似的机器：Linux PC，16 GB RAM，intel 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)
根据knb和李哲元的回复，我开始研究优化的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 秒。