Preconditioner for finding the smallest eigenpairs of a large, but structured, matrix

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.

Assuming your overall matrix is positive definite (it is definitely symmetric), then I would suggest looking into algebraic multigrid (AMG) methods as preconditioners. They compute hierarchies of sparsified matrices themselves. If you're already using PETSc, take a look at the hypre preconditioner. Using this may force you to actually multiply out the elements of the matrix, but this may be worth doing. (It may also be possible that hypre has a way of accepting a fine level matrix that is unassembled.)

• My matrix is positive definite: the Laplacian is PSD and adding cI to it makes it positive definite. I'm not using PETSc currently but I could switch to using it and SLEPc. – Jeff Sep 28 '13 at 3:58
• $M^{-1}$ is a 4,000,000 x 4,000,000 matrix. If I approximate $M^{-1}$ using AMG / hyper, can I avoid storing a dense 4,000,000 x 4,000,000 matrix in memory? – Jeff Sep 28 '13 at 4:03
• Absolutely. You can only work with matrices of this size if they are sparse or (in the case that you only know the inverse) you solve linear systems whenever you need to apply it. – Wolfgang Bangerth Sep 28 '13 at 10:47

If you are only interested in the smallest eigenvalue, the conjugate gradient method applied to the matrix $L$ gives you a good approximation after a reasonably small number of steps, and you won't have to solve any linear systems. The details are in Y. Saad's book on iterative methods, but here a short summary:

1. From the coefficients that are computed in the conjugate gradient (cg) iteration you can construct a tridiagonal matrix $T$ of dimension number of cg steps. This matrix is the Galerkin projection of the matrix $L$ into the current Krylov space.

2. When the number of cg steps increases, the two extremal eigenvalues of $T$ approximate those of $L$.

3. Since $T$ is small and tridiagonal, it is fairly cheap to compute its eigenvalues, for instance with LAPACK.

4. Once you have the eigenvalue, computing the corresponding eigenvector is fairly cheap.

As a side node: the cg solver implemented in deal.II has an option that allows you to compute the spectrum of $T$ in every step and thus to observe convergence of the eigenvalues.

• This looks similar to Lanczos algorithm, which I'm currently using. With Lanczos algorithm I could directly find the smallest eigenvalues of $L$ (without solving linear systems), but that's slow because the smallest eigenvalues are all near $0$---the eigengap between the smallest two eigenvalues relates inversely to the convergence rate. Instead I find the largest eigenvalues of $(L + cI)^{-1}$---as I understand it, the largest eigenvalues of this matrix are farther apart. It's faster, but still slow. – Jeff Sep 28 '13 at 4:16

You should definitely try the Jacobi-Davidson method. One advantage of this method is that an exact $M^{-1} v$ is not necessary. Only an approximate $M^{-1} v$ is needed. The convergence of the Jacobi-Davidson method may be improved by using a preconditioner. If it is possible to compute $L=D−SBS^T$ explicitly, then you can try a level '0' incomplete cholesky factorization of this matrix. Possibly you need a small shift to avoid that the incomplete factorization breaks down on a zero or negative pivot. You can also use Jaobi-Davidson without preconditioner, but the convergence may be less well in that case. There exists a matlab code of Jacobi-Davidson that is optimized for symmetric matrices:JDCG

• Thanks for posting. I can't compute $L$ explicitly because it is so large and not sparse ($S$ and $B$ and $D$ are sparse, but $L$ isn't). I read here, however, that "The Jacobi–Davidson (JD) algorithm is particularly attractive if [the LU or Cholesky] factorization is not feasible". So it sounds like Jacobi-Davidson might be a great solution for me. I'm using Python, so I'll try the "pysparse" implementation. – Jeff Oct 14 '13 at 19:04