MATLAB: Can’t I compute the interior eigenvalues of a sparse matrix with “eigs” without inversion in MATLAB

arpackeigsMATLABsparse

I have a large sparse matrix and would like to compute its interior eigenvalues (i.e., the eigenvalues of smallest magnitude). The matrix is large enough that using the LU decomposition (x = A\b) to solve a linear system involving the matrix results in an out-of-memory error.   Since 'eigs' calls the ARPACK library, I would instead like to solve for the interior eigenvalues with the following approach:
  • Create a function to apply a "shift" of the matrix (A – s*I),
  • Call 'eigs', specifying 'sm' as the value of the input 'sigma' to compute the smallest eigenvalue of the shifted matrix (therefore getting the interior eigenvalues of A), but
  • Specify "regular mode" for ARPACK (that is, the inverse-free method that ARPACK provides) to avoid the out-of-memory error.
However, the 'sm' option for 'eigs' assumes that the user-supplied function returns A\x rather than A*x.  How can I force 'eigs' to use the inverse-free method in ARPACK?

Best Answer

While the ARPACK library that 'eigs' uses allows users to specify both "regular mode" and the 'sm' option, in the general case this combination will lead to poor (or no) convergence on the eigenvalues. Therefore, it would prove best to pursue an alternate method of computing the values.
 
If the matrix is symmetric, one possibility is to use the "folded spectrum" method.  With this technique, 'eigs' is used to compute the smallest algebraic eigenvalues (and corresponding eigenvectors) of the matrix (A-s*I)^2 rather than (A-s*I).  With this approach, the user would provide a function that applies the squared matrices and supply the 'sa' option for 'sigma' to 'eigs'. While the function requires a second application of the matrix in each iteration, it allows an inverse-free computation of the interior eigenvalues of (A - s*I). Please be advised, however, that the convergence of this approach will be slow (but will still occur) if the eigenvalues are clustered close to the shift s. Alternatively, the out-of-memory error encountered when solving a linear system with the matrix may be avoided by constructing an LDL factorization of the matrix with MATLAB's 'ldl' function and calling 'eigs' with the standard options for interior eigenvalues.
If the matrix is not symmetric, an alternate approach such as the "JDQR" method may be necessary.  While this algorithm is not built into MATLAB, there do exist external submissions on File Exchange, eg: 
https://www.mathworks.com/matlabcentral/fileexchange/48-lobpcg-m
Please be advised that these implementations are not authored by MathWorks, and we therefore cannot provide technical support for them. If you have any questions on these approaches or the code that implements them, please contact their respective authors.
Related Question