R中大型稀疏矩阵的特征值

提问者：小点点

R中大型稀疏矩阵的特征值

我有一些大的（对称实）稀疏矩阵，我想计算一些最大和最小的（按大小）特征值。我研究了Armadillo，但它链接到用Fortran编写的ARPACK库。作为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编写的第一个也是唯一一个东西。

根据要求，我可以将矩阵上传到某个地方并分享链接。我可以在Matlab中计算相同矩阵的特征值，但这主要是一个手动过程，我必须将矩阵转移到另一台机器上（也是16 GB的RAM，但Matlab是32位的，所以理论上它的工作空间要有限得多），而这台机器恰好在欧洲，而我在澳大利亚，所以需要很长时间。如果可能的话，我更喜欢在本地机器上使用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

这也可能意味着这个问题与最小特征值使用的shift和inversion模式有关。它真的在这个过程中反转矩阵吗？如果是这样，我可以看到为什么它运行内存溢出，这种矩阵的逆将是完全密集的。

共1个答案

匿名用户

在与rARPACK包的开发人员讨论之后，很明显，问题不在于矩阵被转换为密集，而在于LU分解必须显着改变稀疏矩阵的顺序以避免数值问题，从而大大填充矩阵。

新版本的包将使用LDLT因式分解，看起来没有这个问题，快速计算特征值。

在新版本发布之前，可以获取C源代码的副本并将它们用于：

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
};