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.

  • 1
    $\begingroup$ I would be careful with the statement that $L$ is definitely SPD, since a graph Laplacian has a nullspace of size of the number of connected components. Similarly, a finite discretization of the Laplacian operator can also have a nullspace. $\endgroup$ – Victor Liu Aug 26 '14 at 19:34
  • $\begingroup$ I agree.. I'll edit the question to be more specific. $\endgroup$ – teodron Aug 26 '14 at 20:36
  • 1
    $\begingroup$ Algebraic multigrid might be something to look into $\endgroup$ – James Aug 27 '14 at 14:02

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.

  • $\begingroup$ Thanks for the nice pointers. I'm indeed considering graph Laplacians (planar but irregular strongly connected graphs). $\endgroup$ – teodron Aug 26 '14 at 20:43
  • $\begingroup$ I've tried to benchmark LDLT vs CG in Eigen (since I already had an application that used this library) and LDLT crashed at runtime even for a very simple 3x3 matrix that was positive definite. To make sure it was doable, I ran a Cholesky decomposition using SciLab and it worked like a charm. Maybe I'm using a version of Eigen (from 2013) that has a horrendous bug.. $\endgroup$ – teodron Aug 27 '14 at 5:39

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.

  • $\begingroup$ Very helpful information. Although I cannot accept two answers, your also deserves it. $\endgroup$ – teodron Aug 29 '14 at 6:49
  • $\begingroup$ @teodron: No, you accepted the right answer. Maybe after a few more years of development in the field something like this one will be the right one to accept. But we're not there yet. $\endgroup$ – tmyklebu Aug 29 '14 at 6:54
  • $\begingroup$ What about HSC - cs.huji.ac.il/~raananf/projects/hsc? $\endgroup$ – Royi Mar 23 '15 at 11:16
  • $\begingroup$ @Drazick: I haven't used it. Give it a shot on whichever problems you're interested in. $\endgroup$ – tmyklebu Mar 24 '15 at 2:11
  • $\begingroup$ I just thought it would be great of you added it to your answer. As you seemed to create a list. $\endgroup$ – Royi Mar 24 '15 at 7:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.