Fast Laplacian Solvers, in depth: condition numbers, preconditioned CG, spanning-tree preconditioners measured, the near-linear-time lineage and what to run in practice

The system and its null space

For an undirected graph with non-negative edge weights w, the Laplacian is L = D - A: the weighted degree on the diagonal and -w(u, v) off it. Every row sums to zero, so the all-ones vector is in the null space and L is singular. On a connected graph that is the only null direction, so L x = b has a solution exactly when b sums to zero, and the solution is unique up to adding a constant. On a graph with k components, b must sum to zero on each component and there are k free constants. Two standard ways to remove the ambiguity are to ground one vertex (delete its row and column, which makes the matrix positive definite) or to keep L, start from zero and project the mean out of the answer.

Laplacians matter beyond graphs because any symmetric, diagonally dominant (SDD) matrix can be reduced to a Laplacian of a graph at most twice the size, so a fast Laplacian solver is a fast SDD solver. That class includes the matrices from finite-difference and finite-element discretisations of diffusion problems.

Why it is hard: condition numbers and fill

A dense solve costs O(n³), so the question is how fast iterative methods go. Conjugate gradient needs about √κ log(1/ε) iterations, each one sparse matrix-vector product costing O(m), where κ is the condition number, the ratio of the largest to the smallest relevant eigenvalue. For Laplacians κ can be enormous. The smallest nonzero eigenvalue of a path or grid with side s shrinks like 1/s², and weights spanning several orders of magnitude multiply the problem: the largest eigenvalue tracks the heaviest weighted degree while the smallest is limited by the weakest cut.

Direct methods have their own limit. Sparse Cholesky with a nested-dissection ordering costs about O(n^1.5) on planar graphs and fits comfortably for 2D meshes of a few million vertices, but on expanders and social graphs every elimination creates dense fill, and the factor no longer fits in memory. A fast Laplacian solver has to beat both: few iterations even when κ is large, with only near-linear setup and memory.

Preconditioned conjugate gradient and tree preconditioners

Preconditioning replaces L x = b with an equivalent system whose condition number is small. Pick a matrix M that approximates L and is cheap to solve with; preconditioned CG then needs about √κ(M^-1L) iterations, each costing one product with L and one solve with M. The whole art is in M. The Jacobi preconditioner uses only the diagonal, which corrects for vertices with very different weighted degrees and costs nothing. A spanning-tree preconditioner, proposed by Vaidya in the early 1990s, uses the Laplacian of a spanning tree: a tree Laplacian can be solved exactly in O(n) by eliminating leaves, and if the tree keeps the heavy edges it captures most of the graph's stiffness. The code below implements PCG, the Laplacian and a maximum-weight spanning-tree preconditioner; the Kruskal's algorithm article explains the tree construction.

import numpy as np, scipy.sparse as sp, scipy.sparse.linalg as spla
from scipy.sparse.csgraph import minimum_spanning_tree

def laplacian(n, E, w):
    A = sp.coo_matrix((np.r_[w, w], (np.r_[E[:, 0], E[:, 1]], np.r_[E[:, 1], E[:, 0]])),
                      shape=(n, n)).tocsr()
    return (sp.diags(np.asarray(A.sum(1)).ravel()) - A).tocsr()

def pcg(L, b, M=None, tol=1e-8, maxit=20000):
    x = np.zeros_like(b); r = b.copy(); z = M(r) if M else r.copy()
    p = z.copy(); rz = r @ z; bn = np.linalg.norm(b)
    for k in range(1, maxit + 1):
        Lp = L @ p; a = rz / (p @ Lp); x += a * p; r -= a * Lp
        if np.linalg.norm(r) <= tol * bn:
            return x, k
        z = M(r) if M else r; rz2 = r @ z; p = z + (rz2 / rz) * p; rz = rz2
    return x, maxit

def max_tree_preconditioner(n, E, w):
    """Vaidya-style: exact solver for the maximum-weight spanning tree's Laplacian."""
    T = minimum_spanning_tree(sp.coo_matrix((1 / w, (E[:, 0], E[:, 1])),
                                            shape=(n, n)).tocsr()).tocoo()
    Lt = laplacian(n, np.c_[T.row, T.col], 1 / T.data)
    lu = spla.splu(Lt[1:, 1:].tocsc())         # ground vertex 0; trees allow no-fill elimination
    return lu.solve

# usage: Lg = laplacian(n, E, w)[1:, 1:]; x, iters = pcg(Lg, b, max_tree_preconditioner(n, E, w))

Worked example: iteration counts on grids

Grids of side 32, 64 and 128 were built with either unit weights or weights e^{U(-3, 3)}, which spread over a factor of up to about 400, as a heterogeneous conductance field would. Vertex 0 was grounded, b was a random normal vector, and each solver ran to a relative residual of 10^-8. The condition numbers are those of the grounded matrix.

Three things stand out. First, plain CG doubles its iterations each time the side doubles, consistent with the √κ growth expected from a condition number that rises a little faster than the square of the side. Second, Jacobi does nothing on unit weights, where every interior degree is already 4, and saves about 40 percent on spread weights. Third, the tree preconditioner loses on unit weights, where no edge is more important than another and a tree throws away half the edges, but wins by a factor of nine to ten over plain CG on spread weights, because the heavy edges it keeps are the ones that make the system stiff. Its iteration count still grows with size, which is the gap the theory set out to close. A threshold-based incomplete LU was also tried as a preconditioner and did not converge in 20,000 iterations; CG requires a symmetric positive definite M, and an ILU factorisation is not symmetric in general. Its symmetric counterpart, incomplete Cholesky, is the classic CG preconditioner for these matrices, but SciPy does not ship one.

From trees to near-linear time

The near-linear-time solvers make the tree idea quantitative. The stretch of an edge (u, v) with respect to a spanning tree is its weight times the resistance of the tree path between u and v; total stretch bounds how badly the tree approximates the graph. Low-stretch spanning trees with total stretch about m log n, up to log log factors, can be built in near-linear time, and preconditioning with one gives roughly √(m log n) iterations.

Spielman and Teng (STOC 2004) closed the rest of the gap with recursive preconditioning: augment the tree with a small number of sampled off-tree edges, eliminate the many vertices of degree one and two exactly, which shrinks the system, and solve the smaller system recursively with the same method. They also introduced spectral sparsifiers, sparse graphs whose Laplacians approximate the original within a factor of 1 ± ε on every vector. Spielman and Srivastava later showed that sampling edges with probability proportional to weight times effective resistance yields sparsifiers with O(n log n / ε²) edges. Koutis, Miller and Peng simplified and sped up the recursion around 2010 and 2011, and later work by Cohen and co-authors pushed the running time to roughly m times a square root of log n.

Two simpler algorithms came next. Kelner, Orecchia, Sidford and Zhu (STOC 2013) treat the solve as finding an electrical flow: route b along a low-stretch tree, then repeatedly pick a random off-tree edge and push flow around the cycle it closes until Kirchhoff's voltage law holds. Kyng and Sachdeva (FOCS 2016) gave approximate Gaussian elimination: eliminating a vertex of a Laplacian creates a weighted clique on its neighbours, and instead of adding the whole clique they add a few random edges whose expectation equals it. The result is a sparse approximate Cholesky factor used as a preconditioner, with no trees, sparsifiers or expanders in the algorithm. The connection to Gaussian elimination is direct: exact elimination of a Laplacian vertex is the Schur complement, and that is again a Laplacian.

A caution on practice: an implementation study by Hoske and colleagues (2015) confirmed that the Kelner et al. solver behaves near-linearly but found its constant factors made it much slower than conventional methods on realistic inputs. Asymptotics are a guide, not a benchmark.

What to use in practice

Expander-like graphs deserve emphasis: a random graph or a social network has a large spectral gap, so even plain CG converges in tens of iterations, and the elaborate solvers exist for the opposite case, graphs with long thin regions and wild weights. Measure κ or simply the plain-CG iteration count on a sample before choosing. For eigenvector work on the same graphs, such as the Fiedler vector in spectral clustering, the same preconditioners plug into LOBPCG.

Failure modes

• Right-hand side not orthogonal to the null space. If b does not sum to zero on each component, CG on the singular L drifts and never converges. Project b, or ground one vertex per component.
• Disconnected input. A solver handed two components returns constants of its own choosing on each; compare potentials only within a component.
• Residual is not error. With κ around 10⁶, a residual of 10^-8 can hide a much larger error in x. Set tolerances from the quantity you need, such as voltage differences or effective resistances, and check one.
• Non-SPD preconditioners. ILU factors, one-sided Gauss-Seidel and unsymmetric smoothers break CG's assumptions, as the ILU run showed. Use incomplete Cholesky or symmetric smoothers, or switch to GMRES knowingly.
• Zero or negative weights. Zero weights can disconnect the graph; negative ones make the matrix indefinite, and nothing in this article applies.
• Float32 accumulation. High weight contrast plus single precision stalls convergence well above the requested tolerance. Use float64 for the outer loop.

Trade-offs

The design space is setup cost against iteration count. Plain CG has no setup and is ideal when the graph is well conditioned. Jacobi is free and should always be on. Tree preconditioners are cheap to build and help exactly when weights are heterogeneous, but they lose when the graph is uniform. Multigrid and approximate Cholesky cost more to set up and give iteration counts that barely grow with size, which pays off for big graphs or many solves. Direct factorisation wins on moderate planar problems and loses badly on graphs with poor separators. The theoretical near-linear solvers are the reason we know Laplacian systems are fundamentally easy; the practical descendants are what you deploy.

What to do next

1. Write down how your Laplacian is built, whether the graph is connected, and how you remove the null space: grounding or projection.
2. Run plain CG and Jacobi PCG with the code above on a representative instance and record iteration counts; that one number tells you which row of the table you are in.
3. If weights vary by orders of magnitude, try the maximum-spanning-tree preconditioner and compare iteration counts as this article did.
4. For large or repeated solves, benchmark pyamg and an approximate-Cholesky solver against sparse Cholesky on your real data, including setup time.
5. Validate against a direct solve on a small instance, checking the error in the quantity you report, not just the residual.
6. Read the electrical networks article to connect potentials, currents and effective resistance to the solve.