Lanczos Restarting question

[peekxc]

Hi @yixuan,

I'm trying to understand the Lanczos method, mostly via Golub & Van Loan's Matrix Computations and various other articles.

Following this breakdown, my understanding is that the basic Lanczos (Algorithm 1) can compute the full spectrum of eigenvalues of a (n x n) sparse symmetric matrix using just O(n) workspace memory via the three-term recurrence, but only in exact arithmetic. That reference calls the 3-term-recurrence approach local orthogonalization.

Rounding errors that occur in the eigenvalue estimation procedure are due to loss of orthogonality in the Lanczos vectors, so we must either use more Lanczos vectors or some other means of retaining orthogonality using finite-precision arithmetic. Reorthogonalization via every Lanczos vector is possible, but this implies an O(n^2) memory overhead per-iteration. Other approaches seem to include, but are not limited to, block-Lanczos methods which effectively use O(nb) memory, polynomial filtering schemes, and the so-called implicit restart scheme used here.

I'm wondering whether it's possible to use Spectra (or ARPACK) to do either an explicit restart-like scheme, similar to Algorithm 4 from here, or possibly the O(nb) block Lanczos approach. Essentially, I want to have the ability to fix the number of Lanczos vectors used to acquire the Ritz values. Both Spectra and e.g. the SciPy port of ARPACK has a parameter option ncv, but it must be larger than the number of requested eigenvalues! I'm would like to be able to fix some value k, such that the implementation doesn't exceed O(nk) memory, at the possible sacrifice of needing more iterations to obtain all n eigenvalues up to tol precision.

Is this possible to do in Spectra, possibly via some loop that uses a shift-and-invert solver?

EDIT: To give more context, I know that shifting can be used obtain eigenvalues in the vicinity of the chosen shift. For example, if I wanted an iterative-solver-based approach to obtaining the full spectrum of some psd matrix that used at most k=2 lanczos vectors, I could use something like:

from scipy.sparse import random
from scipy.sparse.linalg import eigsh, aslinearoperator
np.set_printoptions(precision=2, suppress=True)

M = random(20,20, density=0.15, random_state=1234)
A = aslinearoperator(M @ M.T)

print(eigsh(A, k=19, return_eigenvectors=False))
# 0.   0.   0.01 0.03 0.06 0.07 0.12 0.2  0.3  0.52 0.62 0.67 0.89 1.13 1.3 1.59 2.47 2.97 4.63

print(eigsh(A, sigma=4.0, k=2, return_eigenvectors=False, which='LM'))
# [2.97 4.63]
print(eigsh(A, sigma=2.40, k=2, return_eigenvectors=False, which='LM'))
# [2.47 2.97]
print(eigsh(A, sigma=1.50, k=2, return_eigenvectors=False, which='LM'))
# [1.3  1.59]
print(eigsh(A, sigma=1.20, k=2, return_eigenvectors=False, which='LM'))
# [1.13 1.3 ]
# ...

But of course the shifts were chosen manually here, and I'm sure there exists better ways of doing this. It seems quite similar to the explicit restart method that partitions the lanczos vectors into "locked" and "active" sets, but it's not clear to me how everything connects together. Do you have any good references on this? Or is the ability to do this already embedded in spectra, just not at the interface level?

[yixuan]

Hi @peekxc, sorry for the late reply. I was not able to get the time looking at the issues here.

Let me try to understand your question. Basically you want to compute all eigenvalues of a sparse matrix with only $O(nk)$ memory use. Is that correct?

[peekxc]

Yes. I was confused at the time as the material I was reading suggested this was possible, but none of the Lanczos software I could find seemed to support this, which I thought seemed like a huge discrepancy.

I ended up finding PRIMME actually achieves this with variations of the Jacobi-Davidson method; you can set parameters like maxBasisSize

[yixuan]

Right. The reason is that implicitly restarted Arnoldi/Lanczos algorithm is a Krylov space method, and the restart needs to maintain the upper Hessenberg (or tridiagonal in the Lanczos case) structure of the $H$ matrix such that $AV\approx VH$. As a result, extra space is needed to compute the "unwanted eigenvalues", i.e., the shifts for the QR iteration.

Jacobi-Davidson no longer maintains the upper Hessenberg structure of the $H$ matrix, so its restart is easy. The cost is that $H$ is now a general matrix, thus it is more time-consuming to compute the Schur decomposition.

In my experience, the ncv parameter is typically at the same order of k, so it does not result in much larger memory footprint. The real problem is to compute all eigenvalues using iterative methods like Lanczos and Jacobi-Davidson. In my opinion, they are designed to compute just a few eigenvalues.

[peekxc]

The real problem is to compute all eigenvalues using iterative methods like Lanczos and Jacobi-Davidson. In my opinion, they are designed to compute just a few eigenvalueMs.

Right I realize most people just require a few eigenvalues but I encountered a situation where I needed to estimate all of the non-zero eigenvalues (but not the eigenvectors) of a nearly-full-rank LinearOperator which admitted a highly structured matvec operation (e.g. its a Laplacian type operation, $O(n)$ to evaluate).

In my opinion, they are designed to compute just a few eigenvalues.

Yes I've had quite a difficult time figuring out how to adjust the hyper-parameters to keep them competitive time-wise with direct methods or with the ARPACK implementation (e.g. havent had any luck with preconditioning, really). Nonetheless, it seems if the matvec is cheap enough and the restart parameters are configured well, I can get all of the eigenvalues in a time thats within a small constant factor of scipy.sparse.linalg.eigh