Why can't I compute the interior eigenvalues of a sparse matrix with "eigs" without inversion in MATLAB?

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?

Answer 1

[MathWorks Support Team]

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-locally-optimal-block-preconditioned-conjugate-gradient

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.

[Christine Tobler]

This method will only work for very particular matrices: With 'sm', EIGS uses the LU decomposition to solve linear systems with A. If A is symmetric, it's possible that the LDL decomposition uses less memory than LU.

If that is the case, you can use the function handle syntax of EIGS to make it use LDL instead of LU. From the doc:

 [L,U,p] = lu(A,'vector');
 Afun = @(x) U\(L\(x(p)));
 d = eigs(Afun,100,6,'sm')

Using LDL instead:

 [L,D,P] = ldl(A);
 Afun = @(x) P*(L'\(D\(L\(P'*b))))
 d = eigs(Afun,100,6,'sm')

For many problems, if LU goes out-of-memory, so will LDL. In these cases, another of the methods above may help instead.

Answer 2

[Andrew Knyazev]

check https://www.mathworks.com/matlabcentral/fileexchange/48-locally-optimal-block-preconditioned-conjugate-gradient