How to use `ARSymStdEig` in the shift-invert mode?

[dburov190]

The problem I'm solving is a truncated SVD, so following an example in the documentation, I'm using ARSymStdEig with ARdsNonSymMatrix in the following way:

    typedef ARdsNonSymMatrix<double, double> ArMatrix;
    ArMatrix ar_matrix(n_rows, n_cols, raw_matrix);

    ARSymStdEig<double, ArMatrix>
      right_eigen_solver(n_cols, n_sing, &ar_matrix, &ArMatrix::MultMtMv);

The issue with this is that my singular values are clustered around 1.0, the first few are actually exactly 1.0, and then they decay very slowly, so I was thinking of using the shift-invert mode to speed up the convergence. Documentation says in that case the method you provide in the constructor should return (A - sigma * I)^(-1) * v instead of A*v, but there is no ready-made method for that in ARdsNonSymMatrix.

Should I inherit from ARdsNonSymMatrix and write my own method that solves a linear system (A - sigma * I) * x = v within itself? What are the arguments, and how do I know the value of sigma?

Thanks!

[m-reuter]

Hi, I used shift invert to solve a generalized symmetric sparse problem, similar to https://github.com/m-reuter/arpackpp/blob/master/examples/superlu/sym/lsymgshf.cc I had the first eigenvalue be 0, which lead to instabilities, so I moved it out of the way slightly (using sigma to shift the problem away from zero). In your case, I don't know. You can play around with inverting (so that eigenvalues increase) and shift so that all are above zero.
About the missing method in ARdsNonSymMatrix, not everything is completely implemented it seems. I have never used non symmetric matrices, so not sure if this was omitted on purpose or if it is just missing. Sorry I can't be more helpful.

[dburov190]

No worries! And thanks for a quick reply.

I understand why mathematically it should work, I'm just wondering if there is a ready-made solution (it seems not). I think they left out those implementations on purpose (which is good), but didn't provide a convenient interface to implement it yourself (not good).

The manual says the (A - sigma * I)^(-1) * v solver is automatically provided by ARPACK++ for ARluSymStdEig, but it uses LU-decomposition. This is too memory-hungry and inefficient for large matrices, I believe? I'd probably use some sort of an iterative solver.

I wonder if there's any theory for how to combine IRAM and an iterative solver that is also based on Krylov subspaces. Because it seems to me that something like MINRES could benefit greatly from reusing the already computed Krylov subspaces.