密集与稀疏矩阵代数的速度
Speed of dense vs. sparse matrix algebra
我将在 R 中使用相当大 (7 e6 x 4.5 e3) 但非常稀疏的矩阵。所以我试图了解如何有效地处理稀疏矩阵。我有两个相关的问题。
首先:我得到的理解是 Matrix
包 links 自动到 LAPACK 和 SuiteSparse 编译的 dll。 (我在 Windows 工作。)我的印象是,与使用 LAPACK 套件的密集矩阵相比,使用 SuiteSparse 例程会缩短执行时间。但是下面的测试 运行 表明矩阵的稀疏版本的 运行 时间比密集版本慢 很多 。
> library(Matrix)
> sparse <- sparseMatrix(1:4, 1:4, x=rnorm(4))
> dense <- as.matrix(sparse)
> x <- 1:4
> system.time(for (i in 1:10000) sparse %*% x)
user system elapsed
0.23 0.00 0.23
> system.time(for (i in 1:10000) dense %*% x)
user system elapsed
0 0 0
> system.time(for (i in 1:1000) solve(sparse))
user system elapsed
3.94 0.00 3.94
> system.time(for (i in 1:1000) solve(dense))
user system elapsed
0.05 0.00 0.05
a) Matrix
自动与上述两个编译库连接是否正确?如果没有,我如何 link 这些 DLL?
b) 通常使用稀疏矩阵代数实际上比使用密集矩阵代数慢得多吗?
其次:我已经安装了 RcppEigen
和 RcppArmadillo
软件包。我已经能够使用 RcppArmadillo
编译测试程序(使用 Dirk Eddelbuettel 和 Conrad Sanderson 的论文)。但在我的一生中,我还没有找到与 RcppEigen
类似的介绍,它会给我一些我可以用来入门的模型代码。你们中的任何人都可以指出类似于 Eddelbuettel 和 Sanderson 论文的文档可以帮助我开始使用 RcppEigen
吗?
(评论有点太长了。)我将从分析这个更大的矩阵开始;我可以想象,当矩阵很小而不是 非常 稀疏时,稀疏算法处于劣势(例如,在这种情况下,25% 的单元格是非零的)。在下面的示例(1000x1000 矩阵)中,稀疏求解器比密集求解器快 26 倍。您可能会发现 Matrix
例程对于您的目的而言足够快,而无需承担学习 (Rcpp)Eigen
/(Rcpp)Armadillo
...
的额外认知开销
library(rbenchmark)
library(Matrix)
set.seed(101)
sparse <- sparseMatrix(1:1000,1:1000,x=rnorm(1000))
dense <- as.matrix(sparse)
benchmark(solve(sparse),solve(dense),replications=20,
columns = c(
"test", "replications", "elapsed", "relative", "user.self"))
## test replications elapsed relative user.self
## 2 solve(dense) 20 6.932 26.868 6.692
## 1 solve(sparse) 20 0.258 1.000 0.256
我将在 R 中使用相当大 (7 e6 x 4.5 e3) 但非常稀疏的矩阵。所以我试图了解如何有效地处理稀疏矩阵。我有两个相关的问题。
首先:我得到的理解是 Matrix
包 links 自动到 LAPACK 和 SuiteSparse 编译的 dll。 (我在 Windows 工作。)我的印象是,与使用 LAPACK 套件的密集矩阵相比,使用 SuiteSparse 例程会缩短执行时间。但是下面的测试 运行 表明矩阵的稀疏版本的 运行 时间比密集版本慢 很多 。
> library(Matrix)
> sparse <- sparseMatrix(1:4, 1:4, x=rnorm(4))
> dense <- as.matrix(sparse)
> x <- 1:4
> system.time(for (i in 1:10000) sparse %*% x)
user system elapsed
0.23 0.00 0.23
> system.time(for (i in 1:10000) dense %*% x)
user system elapsed
0 0 0
> system.time(for (i in 1:1000) solve(sparse))
user system elapsed
3.94 0.00 3.94
> system.time(for (i in 1:1000) solve(dense))
user system elapsed
0.05 0.00 0.05
a) Matrix
自动与上述两个编译库连接是否正确?如果没有,我如何 link 这些 DLL?
b) 通常使用稀疏矩阵代数实际上比使用密集矩阵代数慢得多吗?
其次:我已经安装了 RcppEigen
和 RcppArmadillo
软件包。我已经能够使用 RcppArmadillo
编译测试程序(使用 Dirk Eddelbuettel 和 Conrad Sanderson 的论文)。但在我的一生中,我还没有找到与 RcppEigen
类似的介绍,它会给我一些我可以用来入门的模型代码。你们中的任何人都可以指出类似于 Eddelbuettel 和 Sanderson 论文的文档可以帮助我开始使用 RcppEigen
吗?
(评论有点太长了。)我将从分析这个更大的矩阵开始;我可以想象,当矩阵很小而不是 非常 稀疏时,稀疏算法处于劣势(例如,在这种情况下,25% 的单元格是非零的)。在下面的示例(1000x1000 矩阵)中,稀疏求解器比密集求解器快 26 倍。您可能会发现 Matrix
例程对于您的目的而言足够快,而无需承担学习 (Rcpp)Eigen
/(Rcpp)Armadillo
...
library(rbenchmark)
library(Matrix)
set.seed(101)
sparse <- sparseMatrix(1:1000,1:1000,x=rnorm(1000))
dense <- as.matrix(sparse)
benchmark(solve(sparse),solve(dense),replications=20,
columns = c(
"test", "replications", "elapsed", "relative", "user.self"))
## test replications elapsed relative user.self
## 2 solve(dense) 20 6.932 26.868 6.692
## 1 solve(sparse) 20 0.258 1.000 0.256