I'm trying to find the eigenvector corresponding to the second smallest eigenvalue of a large $(4,000,000 \times 4,000,000)$ matrix $L$. $L$ is a graph Laplacian, with the following structure: $L = D - SBS^\intercal$, where $D$ is a diagonal matrix and $S$ and $B$ are a large sparse matrices. $S$ has dimension $4,000,000 \times 10,000,000$, but only about $40,000,000$ non-zero entries. $B$ has similar scale and sparsity. So I can rapidly perform matrix-vector multiplication: $Lv = Dv - S(B(S^\intercal v))$.

Currently I'm using Scipy (which calls an Arnoldi algorithm implemented in ARPACK) to find the smallest eigenvalues and corresponding eigenvectors of $L$. Rather than directly finding the smallest eigenpairs of $L$, I find the largest eigenpairs of $M^{-1}$, where $M = L + cI$. (Adding $cI$ and inverting doesn't change the eigenvectors.) To compute $M^{-1}v$, where $v$ is an arbitrary vector, I use the conjugate gradient algorithm.

But this approach is slow: the Arnoldi algorithm often takes hundreds of iterations to converge, and each Arnoldi iteration requires a few hundred conjugate gradient (CG) iterations. Since I'm solving hundreds of CG problems with the same design matrix $M$, I was thinking there could be much efficiency gained by preconditioning. But what preconditioner would be suitable for this problem? I've read that provably good sparse approximations of Laplacians exists ("ultrasparsifiers"), but can they be shifted and inverted, without destroying their sparsity, to form good approximations of $M^{-1}$? Or are there techniques for preconditioning that exploit the structure of $M$?

Thanks much.