Lanczos algorithm for finding top eigenvalues of a matrix sum

Question

I am trying to find the top k leading eigenvalues of a NumPy matrix (using python dot product notation) L@L + a*[email protected], where $L$ and $X$ are a symmetric $n \times n$ and an $n \times d$ matrix, respectively, and $a$ is a scalar.

According to the below text from this paper, I should be able to calculate these leading eigenvalues with L@(L@v) + a*X@(X.T@v), where I guess v is an arbitrary vector. The Lanczos paper they cite is here.

I'm not quite sure where to start. I know that scipy has scipy.sparse.linalg.eigsh here, and from the notes it looks like it uses the Lanczos algorithm - but I am at a loss as to whether it's possible to use sparse.linalg.eigsh for my specific use case. I googled around and didn't find a Python implementation for this very quickly -- does anybody know if I can use sparse.linalg.eigsh to calculate this somehow? I definitely don't want to write this algorithm out myself if I can avoid it.

---

Answer 1

You could define a linear opearator and pass it to the function eigsh. Ideally, your matrices $L$ and $X$ are sparse so you can take advantage of the matrix-vector product.

In your case, you would have something like the following.

import numpy as np
from scipy.sparse.linalg import LinearOperator, eigsh

def mv(v):
    a = 2.3
    return L@(L@v) + a*X@(X.T@v)

n = L.shape[0]
k = 10
A = LinearOperator((n, n), matvec=mv)
eigenvalues, eigenvectors = eigsh(A, k=k, which="LM")
```