Solving Lx = b for big sparse Laplacian matrices

Question

What algorithm is more practically suited in terms of performance for solving the $\mathbf{Lx=b}$ equation, where $\mathbf{L}$ is a generic Laplacian matrix (associated to a strongly connected graph, its first eigenvalue being null in most cases). $\mathbf{L}$ can have more than 1 million rows and less than 15 non-null entries per row.

In my specific case, I tried using C++ with the Eigen library for a small example. Seeing that most numerical solver libraries tend to offer both conjugate gradient and Cholesky methods, I still need a more educated advice regarding the pros and cons of such sparse solvers. As far as my experience goes, CG should be better, but I suspect the shear size of the problem can be overwhelming in the end.

Answer 1

Which one was faster when you tried using Eigen? You may also consider Trilinos if you're using C++.

By "generic Laplacian matrix", do you mean a finite difference / finite element discretization of the Laplace operator on some domain, or do you mean the Laplacian of some graph $G$?

Whether or not your solver is overwhelmed as you scale up the number of unknowns depends strongly on the underlying matrix. For example, if $L$ were a band matrix, then direct solution methods like Cholesky require no more space than the original matrix takes up to begin with. On the other hand, if the matrix has a less regular structure -- e.g. it comes from a 5-point stencil in 2D or a random graph -- then matrix factorizations will have a lot of fill-in and CG becomes preferable. In my own experience, the latter case is more common. Generally speaking, for a PDE in 2 dimensions or more, direct methods are faster for problems with less than ~25,000 unknowns, after which iterative methods are the clear winner.

In either case, there's a lot to gain by permuting the matrix so that its bandwidth is reduced. For direct solution methods, this reduces the amount of fill-in; for iterative solution methods, it can enhance locality of reference when accessing the elements of $x$.

Finally, for iterative methods, choosing a good preconditioner is as important, if not more so, than choosing a good solver. Both Trilinos and PETSc have implementations of multigrid preconditioners, for which your problem is the textbook application.

Answer 2

Daniel Shapero's answer is excellent, but I felt I should add the following:

Properly preconditioned iterative methods will almost certainly win here, unless your systems have a very, very special structure.

There's a considerable recent literature, unfortunately mostly of the "asymptotic a priori bounds on running time" variety, on solvers for graph Laplacian systems. See Dan Spielman's page for references to papers.

Spielman's page suggests using the algebraic multigrid method, in particular Oren Livne and Achi Brandt's LAMG to solve Laplacian systems in practice. Since it's already implemented and many of the recent Laplacian solvers that come with provable running time bounds aren't, I'd suggest giving it a shot.