R 中大型稀疏矩阵的特征值
Eigenvalues of a large sparse matrix in R
我有一些大的（对称实数）稀疏矩阵，我想计算一些最大和最小（按量级）的特征值。我查看了犰狳，但是 links 到 ARPACK 库，它是用 Fortran 编写的。作为 Widnows 用户，我手边没有 Fortran 编译器，所以我决定使用为 Windows.
预编译的 ARPACK R 包
我的矩阵是以MatrixMarket格式存储的，在R中，我是这样读的：
library(Matrix)
f <- file("file://./lambda.mtx", "r")
str(m <- readMM(f))
close(f)
# read sparse matrix
矩阵是对称的，只有下半部分实际存储在mtx文件中。我试过使用对称矩阵，但我一无所获。所以我写了以下（对我来说相当可怕）代码来形成稀疏矩阵的两半：
str(A <- sparseMatrix(c(m@i, m@j), c(m@j, m@i), index1 = FALSE, x = c(m@x, m@x), giveCsparse = TRUE, use.last.ij = TRUE))
# do not treat it as symmetric, this doubles the memory but avoids
# conversion to a full matrix in eigs() so it is a good deal ...
str(nnzM <- length(m@x))
str(nnzA <- A@p[A@Dim[2] + 1])
if(nnzM * 2 - m@Dim[2] != nnzA) {
stop("failed to form full matrix A")
}
if(A@x[1] != m@x[1]) {
if(A@x[1] == 2 * m@x[1])
warning("failed to form full matrix A: diagonal elements added")
else
warning("failed to form full matrix A: bad diagonal value")
}
# make sure the repeated values of the diagonal were eliminated instead of added (use.last.ij = TRUE)
最后，我使用 rARPACK 计算特征值：
library(rARPACK)
str(le <- eigs(A, k = 5, which = "LM", opts = list(retvec = FALSE))) # or dsyMatrix
str(se <- eigs(A, k = 5, sigma = 0, opts = list(retvec = FALSE))) # or dsyMatrix
le$values[1]
se$values[se$nconv]
这适用于小矩阵。但是在大矩阵上，我在 eigs() 中得到 bad_alloc 个错误。请注意，我是 运行 64 位版本的 R，盒子里有 16 GB 的 RAM（有问题的稀疏矩阵是 174515 x 174515，有 9363966 个非零值，mtx 中有超过 300 MB（文本格式)).
我以前在使用对称稀疏矩阵时遇到过内存问题。我的猜测是矩阵的格式不正确，仍然会在某处转换为密集矩阵。否则我不明白它怎么会需要这么多内存。任何帮助将不胜感激，这是迄今为止我用 R 编写的第一个也是唯一一个东西。
根据要求，我可以将矩阵上传到某处并分享 link。我可以在 Matlab 中计算相同矩阵的特征值，但这主要是一个手动过程，我必须将矩阵传输到另一台机器（也是 16 GB 的 RAM，但 Matlab 是 32 位的，所以理论上它有更多的限制working space), 而这台机器恰好在欧洲，而我在澳大利亚，所以它需要很长时间。如果可能的话，我更喜欢在我的本地机器上使用 R。
编辑
复制问题所需的所有数据和脚本都在 https://sourceforge.net/projects/slamfrontend/files/data/eigenvalues/。此外，该进程实际上在它死亡之前确实分配了大约 12 GB 的内存。
编辑 2
这是脚本的输出。居然在计算最小特征值的时候崩溃，最大的计算速度很快，没有崩溃：
E:\my-projects\Robotics\MonaschGit\wdir>"C:\Program Files\R\R-3.2.2\bin\x64\R" --no-save -q --slave 0<get_lambda_eigenvalues.R
Formal class 'dsTMatrix' [package "Matrix"] with 7 slots
..@ i : int [1:9363966] 0 1 2 3 4 5 6 1 2 3 ...
..@ j : int [1:9363966] 0 0 0 0 0 0 0 1 1 1 ...
..@ Dim : int [1:2] 174515 174515
..@ Dimnames:List of 2
.. ..$ : NULL
.. ..$ : NULL
..@ x : num [1:9363966] 1.62e+08 -2.49e+01 -2.49e+06 2.19e+07 1.79e+09 ...
..@ uplo : chr "L"
..@ factors : list()
Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
..@ i : int [1:18553417] 0 1 2 3 4 5 6 7 8 9 ...
..@ p : int [1:174516] 0 14362 28724 43086 57448 71810 86172 100534 115525 130516 ...
..@ Dim : int [1:2] 174515 174515
..@ Dimnames:List of 2
.. ..$ : NULL
.. ..$ : NULL
..@ x : num [1:18553417] 1.62e+08 -2.49e+01 -2.49e+06 2.19e+07 1.79e+09 ...
..@ factors : list()
int 9363966
int 18553417
List of 4
$ values : num [1:5] 1.38e+11 1.31e+11 1.30e+11 1.22e+11 1.16e+11
$ vectors: NULL
$ nconv : int 5
$ niter : int 60
Error: std::bad_alloc
Execution halted
这也可能意味着问题与用于最小特征值的移位和反转模式有关。它实际上是否在此过程中反转了矩阵？如果是这样，我明白为什么它会耗尽内存，这样的矩阵的逆将是完全密集的。
在使用了 rARPACK 包的 discussion with the developer 之后，很明显问题不在于矩阵被转换为密集矩阵，而是在于 LU 分解必须显着改变稀疏矩阵的顺序以避免数值问题，从而大量填充矩阵。
新版包会使用LDLT分解，貌似没有这个问题，而且计算特征值很快
在发布新版本之前，可以获取 C++ sources 的副本并将它们用于：
class CSymmetricSparseMatrix_InvProduct {
protected:
typedef Eigen::SparseMatrix<double, Eigen::ColMajor> SpMat;
typedef Eigen::SimplicialLDLT<SpMat> SpLDLTSolver;
const size_t n;
SpMat m_mat;
SpLDLTSolver solver;
public:
CSymmetricSparseMatrix_InvProduct(const cs *p_matrix)
:m_mat(int(p_matrix->m), int(p_matrix->n)), n(p_matrix->n)
{
// could maybe use Eigen::internal::viewAsEigen(p_matrix);
_ASSERTE(p_matrix->m == p_matrix->n);
// must be square
std::vector<Eigen::Triplet<double> > triplets;
triplets.reserve(p_matrix->p[p_matrix->n]);
for(size_t i = 0; i < n; ++ i) {
for(size_t p = p_matrix->p[i], e = p_matrix->p[i + 1]; p < e; ++ p) {
double f_val = p_matrix->x[p];
size_t n_row = p_matrix->i[p], n_col = i;
triplets.push_back(Eigen::Triplet<double>(int(n_row), int(n_col), f_val));
}
}
// expand the CSparse matrix to triplets
m_mat.setFromTriplets(triplets.begin(), triplets.end());
// fill the Eigen matrix
m_mat.makeCompressed();
// compress the Eigen matrix
}
size_t rows() const
{
return n;
}
size_t cols() const
{
return n;
}
void set_shift(double sigma)
{
SpMat I((int)n, (int)n);
I.setIdentity();
solver.compute(m_mat - I * sigma); // also better use solver.setShift(-sigma, 1.0);
/*size_t n_fac_nnz = ((SpMat)solver.matrixL()).nonZeros(); // type cast required, otherwise it loops forever
fprintf(stderr, "debug: the LDLT factorization has " PRIsize " nonzeros\n", n_fac_nnz);*/
}
// y_out = inv(A - sigma * I) * x_in
void perform_op(const double *x_in, double *y_out) const
{
Eigen::Map<const Eigen::VectorXd, Eigen::DontAlign> x(x_in, n);
Eigen::Map<Eigen::VectorXd, Eigen::DontAlign> y(y_out, n);
y.noalias() = solver.solve(x);
}
private:
CSymmetricSparseMatrix_InvProduct(const CSymmetricSparseMatrix_InvProduct &r_other); // no-copy
CSymmetricSparseMatrix_InvProduct &operator =(const CSymmetricSparseMatrix_InvProduct &r_other); // no-copy
};
我有一些大的（对称实数）稀疏矩阵，我想计算一些最大和最小（按量级）的特征值。我查看了犰狳，但是 links 到 ARPACK 库，它是用 Fortran 编写的。作为 Widnows 用户，我手边没有 Fortran 编译器，所以我决定使用为 Windows.
预编译的 ARPACK R 包我的矩阵是以MatrixMarket格式存储的，在R中，我是这样读的：
library(Matrix)
f <- file("file://./lambda.mtx", "r")
str(m <- readMM(f))
close(f)
# read sparse matrix
矩阵是对称的，只有下半部分实际存储在mtx文件中。我试过使用对称矩阵，但我一无所获。所以我写了以下（对我来说相当可怕）代码来形成稀疏矩阵的两半：
str(A <- sparseMatrix(c(m@i, m@j), c(m@j, m@i), index1 = FALSE, x = c(m@x, m@x), giveCsparse = TRUE, use.last.ij = TRUE))
# do not treat it as symmetric, this doubles the memory but avoids
# conversion to a full matrix in eigs() so it is a good deal ...
str(nnzM <- length(m@x))
str(nnzA <- A@p[A@Dim[2] + 1])
if(nnzM * 2 - m@Dim[2] != nnzA) {
stop("failed to form full matrix A")
}
if(A@x[1] != m@x[1]) {
if(A@x[1] == 2 * m@x[1])
warning("failed to form full matrix A: diagonal elements added")
else
warning("failed to form full matrix A: bad diagonal value")
}
# make sure the repeated values of the diagonal were eliminated instead of added (use.last.ij = TRUE)
最后，我使用 rARPACK 计算特征值：
library(rARPACK)
str(le <- eigs(A, k = 5, which = "LM", opts = list(retvec = FALSE))) # or dsyMatrix
str(se <- eigs(A, k = 5, sigma = 0, opts = list(retvec = FALSE))) # or dsyMatrix
le$values[1]
se$values[se$nconv]
这适用于小矩阵。但是在大矩阵上，我在 eigs() 中得到 bad_alloc 个错误。请注意，我是 运行 64 位版本的 R，盒子里有 16 GB 的 RAM（有问题的稀疏矩阵是 174515 x 174515，有 9363966 个非零值，mtx 中有超过 300 MB（文本格式)).
我以前在使用对称稀疏矩阵时遇到过内存问题。我的猜测是矩阵的格式不正确，仍然会在某处转换为密集矩阵。否则我不明白它怎么会需要这么多内存。任何帮助将不胜感激，这是迄今为止我用 R 编写的第一个也是唯一一个东西。
根据要求，我可以将矩阵上传到某处并分享 link。我可以在 Matlab 中计算相同矩阵的特征值，但这主要是一个手动过程，我必须将矩阵传输到另一台机器（也是 16 GB 的 RAM，但 Matlab 是 32 位的，所以理论上它有更多的限制working space), 而这台机器恰好在欧洲，而我在澳大利亚，所以它需要很长时间。如果可能的话，我更喜欢在我的本地机器上使用 R。
编辑
复制问题所需的所有数据和脚本都在 https://sourceforge.net/projects/slamfrontend/files/data/eigenvalues/。此外，该进程实际上在它死亡之前确实分配了大约 12 GB 的内存。
编辑 2
这是脚本的输出。居然在计算最小特征值的时候崩溃，最大的计算速度很快，没有崩溃：
E:\my-projects\Robotics\MonaschGit\wdir>"C:\Program Files\R\R-3.2.2\bin\x64\R" --no-save -q --slave 0<get_lambda_eigenvalues.R
Formal class 'dsTMatrix' [package "Matrix"] with 7 slots
..@ i : int [1:9363966] 0 1 2 3 4 5 6 1 2 3 ...
..@ j : int [1:9363966] 0 0 0 0 0 0 0 1 1 1 ...
..@ Dim : int [1:2] 174515 174515
..@ Dimnames:List of 2
.. ..$ : NULL
.. ..$ : NULL
..@ x : num [1:9363966] 1.62e+08 -2.49e+01 -2.49e+06 2.19e+07 1.79e+09 ...
..@ uplo : chr "L"
..@ factors : list()
Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
..@ i : int [1:18553417] 0 1 2 3 4 5 6 7 8 9 ...
..@ p : int [1:174516] 0 14362 28724 43086 57448 71810 86172 100534 115525 130516 ...
..@ Dim : int [1:2] 174515 174515
..@ Dimnames:List of 2
.. ..$ : NULL
.. ..$ : NULL
..@ x : num [1:18553417] 1.62e+08 -2.49e+01 -2.49e+06 2.19e+07 1.79e+09 ...
..@ factors : list()
int 9363966
int 18553417
List of 4
$ values : num [1:5] 1.38e+11 1.31e+11 1.30e+11 1.22e+11 1.16e+11
$ vectors: NULL
$ nconv : int 5
$ niter : int 60
Error: std::bad_alloc
Execution halted
这也可能意味着问题与用于最小特征值的移位和反转模式有关。它实际上是否在此过程中反转了矩阵？如果是这样，我明白为什么它会耗尽内存，这样的矩阵的逆将是完全密集的。
在使用了 rARPACK 包的 discussion with the developer 之后，很明显问题不在于矩阵被转换为密集矩阵，而是在于 LU 分解必须显着改变稀疏矩阵的顺序以避免数值问题，从而大量填充矩阵。
新版包会使用LDLT分解，貌似没有这个问题，而且计算特征值很快
在发布新版本之前，可以获取 C++ sources 的副本并将它们用于：
class CSymmetricSparseMatrix_InvProduct {
protected:
typedef Eigen::SparseMatrix<double, Eigen::ColMajor> SpMat;
typedef Eigen::SimplicialLDLT<SpMat> SpLDLTSolver;
const size_t n;
SpMat m_mat;
SpLDLTSolver solver;
public:
CSymmetricSparseMatrix_InvProduct(const cs *p_matrix)
:m_mat(int(p_matrix->m), int(p_matrix->n)), n(p_matrix->n)
{
// could maybe use Eigen::internal::viewAsEigen(p_matrix);
_ASSERTE(p_matrix->m == p_matrix->n);
// must be square
std::vector<Eigen::Triplet<double> > triplets;
triplets.reserve(p_matrix->p[p_matrix->n]);
for(size_t i = 0; i < n; ++ i) {
for(size_t p = p_matrix->p[i], e = p_matrix->p[i + 1]; p < e; ++ p) {
double f_val = p_matrix->x[p];
size_t n_row = p_matrix->i[p], n_col = i;
triplets.push_back(Eigen::Triplet<double>(int(n_row), int(n_col), f_val));
}
}
// expand the CSparse matrix to triplets
m_mat.setFromTriplets(triplets.begin(), triplets.end());
// fill the Eigen matrix
m_mat.makeCompressed();
// compress the Eigen matrix
}
size_t rows() const
{
return n;
}
size_t cols() const
{
return n;
}
void set_shift(double sigma)
{
SpMat I((int)n, (int)n);
I.setIdentity();
solver.compute(m_mat - I * sigma); // also better use solver.setShift(-sigma, 1.0);
/*size_t n_fac_nnz = ((SpMat)solver.matrixL()).nonZeros(); // type cast required, otherwise it loops forever
fprintf(stderr, "debug: the LDLT factorization has " PRIsize " nonzeros\n", n_fac_nnz);*/
}
// y_out = inv(A - sigma * I) * x_in
void perform_op(const double *x_in, double *y_out) const
{
Eigen::Map<const Eigen::VectorXd, Eigen::DontAlign> x(x_in, n);
Eigen::Map<Eigen::VectorXd, Eigen::DontAlign> y(y_out, n);
y.noalias() = solver.solve(x);
}
private:
CSymmetricSparseMatrix_InvProduct(const CSymmetricSparseMatrix_InvProduct &r_other); // no-copy
CSymmetricSparseMatrix_InvProduct &operator =(const CSymmetricSparseMatrix_InvProduct &r_other); // no-copy
};