Running out of memory in inverse and shift mode on sparse matrices

[swajnautcz]

Hello,

I am trying to get the Eigenvalues of a large sparse matrix in R. For the large eigenvalues, I have no problems but for the small Eigenvalues, the "SM" does not converge even if I increase the number of iterations and tolerance. If I use sigma = 0, then I get a bad_alloc in a while. The matrix in question is real symmetric 174515 x 174515, with 9363966 nonzeroes and has just over 300 MB in the mtx (text format). I'm running 64-bit version of R and there is 16 GB of RAM in the box, I see little reason why this should happen. Obviously I'm doing something wrong. My R code is as follows:

library(Matrix)
f <- file("file://./lambda.mtx", "r")
str(m <- readMM(f))
close(f)
# read sparse matrix

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)

library(rARPACK)
str(le <- eigs(A, k = 5, which = "LM", opts = list(retvec = FALSE))) # or dsyMatrix # runs ok
str(se <- eigs(A, k = 5, sigma = 0, opts = list(retvec = FALSE))) # or dsyMatrix # bad_alloc
le$values[1]
se$values[se$nconv]

I have never written anything in R, maybe that's the issue. Any insights?

[swajnautcz]

I was also unable to use eigs_sym(), I was getting "matrix too large" errors in cholmod_dense.c.

[yixuan]

Hi @swajnautcz
Are you using the version on CRAN or the development version here?

[swajnautcz]

I'm using the CRAN pre-build stable version for Windows (rARPACK_0.7-0.zip on https://cran.r-project.org/web/packages/rARPACK/index.html). Is there a difference between those?

[yixuan]

The development version is 0.8-0, and you may want to give it a try. A pre-built package can be found at http://win-builder.r-project.org/r0kJrdtI7NyO/

However, the problem might still happen, given the difficulty of the computation. To use the shift-and-invert mode, we need to factorize a large matrix, whose storage is generally not sparse and hence may be out of memory. Another aspect is that if the matrix is singular, then it will for sure run into error.

Usually for large matrices we are only interested in calculating the largest eigenvalues. Is there any specific reason why you also want the smallest ones?

[swajnautcz]

The matrix is positive definite (and hence invertible). I'm doing some analysis of least squares parameterizations, and I'm looking for the condition numbers of the matrices. One way to get a condition number is by taking a ratio of the eigenvalues. I will give the 0.8-0 a try and let you know.

[swajnautcz]

The matrix in question can be downloaded at https://sourceforge.net/projects/slamfrontend/files/data/eigenvalues/ in case you want to try.

[swajnautcz]

Ok, so admittedly, the 0.8-0 did calculate the solution:

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 5
 $ 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 4
 $ nops   : int 60
List of 5
 $ values : num [1:5] 1.32e-08 1.19e-08 1.18e-08 1.09e-08 4.85e-09
 $ vectors: NULL
 $ nconv  : int 5
 $ niter  : int 3
 $ nops   : int 45
[1] 138045534935
[1] 4.850883e-09

and the results do match the Matlab results:

big =
  1.3805e+011

small =
  4.8458e-009

Elapsed time is 161.076702 seconds.

However, to do so, it needed over 15 GB of RAM and about 2 hours runtime. What kind of matrix factorization is being done in this case? My code factorizes the matrix in a second, solving is instant.

I'm quite sure that either I'm calling the code from R inappropriately and somewhere a dense decomposition is called, or the decomposition is being calculated without fill-reducing ordering.

[yixuan]

Thanks for the test.

Since this package is using Eigen (http://eigen.tuxfamily.org/) for linear algebra computation, the factorization I'm using for sparse matrix is an implementation of SuperLU (http://eigen.tuxfamily.org/dox/classEigen_1_1SparseLU.html). This can be improved since we can use sparse LDL to optimize for symmetric matrices.

Eigen also provides other linear solvers including iterative methods (http://eigen.tuxfamily.org/dox/group__TopicSparseSystems.html). I have little experience on large sparse linear solver, so could you give me some hints on which solver I should use? Seeing that you could factorize the matrix at this scale in a second, I'm eager to learn what technique you are using. Thank you!

[swajnautcz]

Oh, I see. The Eigen's sparse LU should be fine, unless something fishy is going on. I'm currently using sparse LLT factorization, as my matrix is positive definite, but that should be only about 1/3 faster than the LU and use about 1/2 of memory, so the LU should still be very reasonable, yet versatile. For non-square problems I would use sparse QR factorization. All these methods require a good ordering in order to run well (in Eigen that's the analyzePattern() function).

Iterative solving on the other hand is slightly less precise, and I believe it is slightly faster only if you solve once. If you solve many times, then the costly factorization is amortized by much simpler solving which degrades to backsubstitution (similar cost as matrix vector multiply, usually two of these are needed).

Could you point me to the code where you calculate the inverse of a sparse matrix? I started using the C++ code directly, rather than through R, and was able to obtain more than reasonable results with this:

#include <csparse/cs.hpp>

class CPosDefSparseMatrix_InvProduct {
    css *m_p_symbolic_decomposition;
    csn *m_p_factor;
    mutable std::vector<double> m_temp;

public:
    CPosDefSparseMatrix_InvProduct(const cs *p_matrix)
        :m_temp(p_matrix->n)
    {
        m_p_symbolic_decomposition = cs_schol(1, p_matrix);
        m_p_factor = cs_chol(p_matrix, m_p_symbolic_decomposition);
        // do sparse cholesky
    }

    ~CPosDefSparseMatrix_InvProduct()
    {
        cs_sfree(m_p_symbolic_decomposition);
        cs_nfree(m_p_factor);
    }

    void set_shift(double f_shift)
    {
        assert(!f_shift);
        // can't solve for other shifts (could solve for positive shifts, but not for negative ones)
    }

    size_t rows() const
    {
        return m_p_factor->L->m;
    }

    size_t cols() const
    {
        return m_p_factor->L->n;
    }

    void operator ()(const double *p_x, double *p_y) const
    {
        csi n_col_num = m_p_factor->L->n;
        cs_ipvec(m_p_symbolic_decomposition->pinv, p_x, &m_temp[0], n_col_num); // temp = P*x
        cs_lsolve(m_p_factor->L, &m_temp[0]); // x = L\x
        cs_ltsolve(m_p_factor->L, &m_temp[0]); // x = L'\x
        cs_pvec(m_p_symbolic_decomposition->pinv, &m_temp[0], p_y, n_col_num); // y = P'*temp
    }

private:
    CPosDefSparseMatrix_InvProduct(const CPosDefSparseMatrix_InvProduct &r_other); // no-copy
    CPosDefSparseMatrix_InvProduct &operator =(const CPosDefSparseMatrix_InvProduct &r_other); // no-copy
};

[yixuan]

The code is at https://github.com/yixuan/rARPACK/blob/master/inst/include/RMatOp/RealShift_sparseMatrix.h

I tried to test your data file, but it seems that the link is broken. Is there any other access to the file? I do want to investigate the weird behavior of sparse LU.

Also, if you prefer C++ code, you may also take a look at the eigen solver (https://github.com/yixuan/arpack-eigen) that the R interface is based on.

[swajnautcz]

Yes, so the problem appears to be that you are not calling solver.analyzePattern() before factorizing the matrix. The improved code would be:

void set_shift(double sigma)
{
    // Create a sparse idendity matrix
    SpMat I(n, n);
    I.setIdentity();
    SpMat M = mat - sigma * I; // or better "mat - SpMat::Identity(n, n) * I;" if SpMat supports it
    solver.analyzePattern(M);
    solver.compute(M);
}

If you do not call analyzePattern(), it will calculate the factorization with the natural ordering, which may introduce considerable fill-in (more nonzero entries in the sparse matrix) which causes both more memory use and also longer runtimes for both the factorization and solving. See e.g. the example section in http://au.mathworks.com/help/matlab/ref/colamd.html or http://au.mathworks.com/help/matlab/ref/symamd.html. Notice how dense the LU gets without the ordering.

As for the link, http://sourceforge.net/projects/slamfrontend/files/data/eigenvalues/ should now work. Appologies for that, there was a permission problem. You can verify that there is a large difference in runtime with / without calling the .analyzePattern().

[yixuan]

It's very nice of you for the careful explanation!

Actually I notice that the compute() method is simply a combination of analyzePattern() and factorize() (http://eigen.tuxfamily.org/dox/SparseLU_8h_source.html, line 112-118), and SparseLU will automatically do an implicit COLAMD ordering. I suspect that the difference is due to the different implementations between Eigen and CSparse.

Later I'll try the sparse Choleskey decomposition, to see if it could give better performance.

[swajnautcz]

Ah, I see, I have never used that one before so I did not know that. I'll try to use it in my code then and will let you know what the result is.

[yixuan]

I'd also like to test the Choleskey solver on this data set, but I saw errors when I tried to extract it. Could you double check the validity of the zip file? Thanks.

[swajnautcz]

I see, the problems never end :). Sorry, I'm on a vacation in Australia and the connection here is just rubbish, the upload probably dropped the first time around without me noticing. I re-uploaded the file. It should indeed have been 115 MB instead of 50 MB.

I have tried using your code and indeed the Eigen sparse LU wrapper seems to be eating loads of memory. Investigating further, I have used CSparse to get the LU decomposition as well (since Eigen does not allow access to the L and U matrices in the sparse LU class). While the Cholesky factorization is quite sparse (10 M nonzeroes), the LU factorization sparsity also depends on the pivoting. Apparently, my matrix is not very well conditioned and the pivoting changes the ordering quite a lot (220 M nonzeroes), introducing many nonzeroes. I'm not sure if using sparse LDLT would make a difference here.

[swajnautcz]

I don't have the complete results yet, but as it stands, Cholesky and LDLT are both fast on the given matrix, while LU and QR are slow and need a lot of memory. Not a fault of the implementation, rather the fault of the matrix.

[yixuan]

You are correct.
I've updated the R package so that a sparse LDLT factorization will be used for symmetric sparse matrices (a symmetry test will be conducted before calculating eigenvalues). Below is the timing result:

> system.time(le <- eigs(A, k = 5, which = "LM", opts = list(retvec = FALSE)))
   user  system elapsed 
  4.592   0.459   5.046 
> system.time(se <- eigs(A, k = 5, sigma = 0, opts = list(retvec = FALSE)))
   user  system elapsed 
  6.589   0.733   7.313 
> le$values[1]
[1] 138045534935
> se$values[se$nconv]
[1] 4.84661e-09

This has significant improvement over the previous version.

A lot thanks to you!

[swajnautcz]

Great, that's a significant improvement indeed. There is one more issue with my script, which is more on the R side I suppose. When trying to call eigs_sym() as:

library(Matrix)
f <- file("file://./lambda.mtx", "r")
str(m <- readMM(f))
close(f)
# read sparse matrix

str(A <- sparseMatrix(m@i, m@j, index1 = FALSE, x = m@x, giveCsparse = FALSE, symmetric = TRUE))
# or giveCsparse = TRUE

library(rARPACK)
str(le <- eigs_sym(A, k = 5, which = "LM", opts = list(retvec = FALSE))) # or as dsyMatrix

Then depending on the combination of giveCsparse and whether using A or as(A, "dsyMatrix"), I get either:

Error in eigs_sym(A, k = 5, which = "LM", opts = list(retvec = FALSE)) :
  unsupported matrix type
Calls: str -> eigs_sym
Execution halted

or:

Error in asMethod(object) :
  Cholmod error 'problem too large' at file ../Core/cholmod_dense.c, line 105
Calls: str -> eigs_sym -> as -> asMethod -> as -> asMethod
Execution halted

or the infamous Windows message that "process R stopped working". I know that this probably does not concern the package implementation so much, but still - do you have an idea how to call eigs_sym() in this case?

[swajnautcz]

Also, not sure if you have noticed, the Eigen::SimplicialLDLT has a setShift() method, so instead of doing:

SpMat I(n, n);
I.setIdentity();
solver.compute(mat - I * sigma);

it is possible to:

solver.setShift(-sigma, 1.0);
solver.compute(mat);

which will probably be even slightly faster and require less memory.

[yixuan]

For now eigs_sym() only works with the ordinary dense matrix type (of class matrix), so it does not accept other type of matrices. I'll extend it later.

And dsyMatrix is a dense matrix type, so the conversion will fail since the whole matrix is too large.

[yixuan]

Thanks for pointing out. Actually I have already adopted this usage. :-)
https://github.com/yixuan/rARPACK/blob/master/inst/include/RMatOp/RealShift_sym_sparseMatrix.h

[swajnautcz]

Right, you have just proved my poor knowledge of R. I didn't realize I'm explicitly asking for a dense matrix, thought that 'd' stands for 'double' there. Thanks a lot for your help and good luck developing the package further.

[yixuan]

Also thanks for your enormous kind help. I'll let you know once I finish the extension of eigs_sym().
Your experience on this is a treasure for me. :-)

[yixuan]

I've added the support for dgCMatrix in eigs_sym() and now we can do

library(Matrix)
f <- file("file://./lambda.mtx", "r")
str(m <- readMM(f))
close(f)
# read sparse matrix

A <- sparseMatrix(m@i, m@j, index1 = FALSE, x = m@x, giveCsparse = TRUE)

library(rARPACK)
system.time(le <- eigs_sym(A, k = 5, which = "LM", opts = list(retvec = FALSE)))
##   user  system elapsed 
##  1.942   0.009   1.949 
system.time(le <- eigs_sym(A, k = 5, sigma = 0, opts = list(retvec = FALSE)))
##   user  system elapsed 
##  3.921   0.347   4.264