A preconditioner that is a tree

A preconditioner is a matrix that is cheap to solve with and resembles the one that is not. On a Laplacian there is a family of candidates with an unusual property: they are subgraphs, so the question “how good is this preconditioner” becomes a question about a graph, and it can be answered before anything is factorised.

The extreme member of the family is a spanning tree. A tree’s Laplacian is solved in linear time with no fill at all — eliminating a leaf creates no clique, because a leaf has one neighbour — so it is as cheap as a preconditioner gets.

A condition number with no arithmetic in it #

The theorem is two inequalities, and both are about a purely combinatorial quantity.

The stretch of an edge over a tree is the resistance of the path the tree makes it take. For an unweighted graph that is the number of tree edges between its endpoints, an integer. The total stretch is the sum over the graph’s edges.

Then, for the generalised eigenvalues of L_G against L_T:

1  ≤  λ  ≤  total stretch of G over T

The lower bound is immediate: a tree is a subgraph, so its quadratic form is never larger than the whole graph’s, and the ratio is at least one. The upper bound is a path-counting argument — each of G’s edges is routed along its tree path, and the cost of that routing is exactly its stretch.

Measured on four graphs and two trees each, the lower bound is attained at 1.000000 to nine digits in seven of eight cases, and the upper bound is loose by:

graph          tree             stretch   λmin      λmax     slack
grid 5×5       breadth-first      120     1.000000  34.62     3.47
grid 5×5       random             126     1.000000  55.17     2.28
cycle 24       breadth-first       46     1.000000  24.00     1.92
two blocks 30  breadth-first      220     1.000000  62.25     3.53
hypercube 4    breadth-first       80     1.000000  21.64     3.70
hypercube 4    random             100     1.037772  40.93     2.54

A bound loose by a factor of two to four is a useful bound — it is in the same class as Cheeger’s upper half, which was loose by 1.4 to 13.7 — and it costs a breadth-first search per vertex rather than an eigenvalue computation.

What the bound is for #

It is not for predicting how long a solve will take. A factor of three in the condition number is a factor of 1.7 in the iteration count, which is a real error in a prediction.

It is for choosing between trees, and there it is exactly right. Two trees can be compared by their total stretch without factorising either, without running a solve, and without touching a floating-point number. The tree with the smaller stretch is the better preconditioner, and on the grid above the ordering by stretch (116, 120, 126) matches the ordering by measured condition number (37.1, 34.6, 55.2) in two of three pairs — which is a caveat worth having, since the bound is a bound and not a ranking.

That is the general shape of a combinatorial bound on a spectral quantity, and it is the same shape the sparsity field’s ordering heuristics have: a cheap combinatorial proxy for an expensive numerical property, useful for choosing and not for predicting.

How the spectrum is computed, and why it is not a solve #

The measurement above is a generalised eigenvalue problem, and the way it is set up is worth a paragraph because the obvious way does not work.

L_G and L_T are both singular — both are Laplacians, both annihilate the vector of ones — so the pencil (L_G, L_T) has a shared null vector and the generalised problem is singular. The eigenvalues that matter live on the complement, and the arithmetic has to be arranged so that the shared kernel is removed exactly rather than being left for the eigensolver to discover.

Grounding does it. Delete the same vertex’s row and column from both, and what is left are two positive definite matrices whose generalised eigenvalues are precisely the ones wanted, with the shared kernel gone by construction rather than by rounding. Then a Cholesky factorisation of the grounded tree, L_T = CCᵀ, turns the problem into an ordinary symmetric one: the eigenvalues of C⁻¹L_G C⁻ᵀ.

The Cholesky is the step where the tree’s cheapness shows. A tree’s grounded Laplacian factorises with no fill whatsoever — it is a tree, so an ordering by leaves eliminates every vertex without creating a clique — and the factor has exactly as many nonzeros as the tree has edges. That is the whole reason a tree is a candidate preconditioner, and it is visible in the measurement rather than merely assumed: the same construction on a denser subgraph would cost a dense factorisation.

The symmetrised matrix is checked for symmetry before the eigensolver sees it, because C⁻¹L_G C⁻ᵀ is symmetric in exact arithmetic and is assembled from two triangular solves that are not — and an eigensolver handed a nearly-symmetric matrix returns nearly-real eigenvalues without saying so.

The measurement that does not go the way it should #

Here is the number the essay exists for.

graph      κ⁺(L)     tree κ    ratio (L ÷ tree)
path 48     933.11     1.00       933.1
cycle 48    233.78    48.00         4.87
cycle 36    131.65    36.00         3.66
cycle 24     58.70    24.00         2.45
grid 49      38.39    71.70         0.535
grid 36      27.86    51.52         0.541
grid 25      18.94    34.62         0.547
grid 16      11.66    21.02         0.555

On every grid the preconditioner makes the conditioning worse, by a factor of about 1.85, and the factor is nearly constant across sizes.

That is not a bug in the measurement and it is not a failure of the theorem. The theorem says the preconditioned eigenvalues lie between 1 and the stretch; it does not say the interval is narrower than the original spectrum’s. On a grid, κ⁺(L) is about 27 at thirty-six vertices while the breadth-first tree’s stretch is 210, and the resulting κ of 51 is what a bound of 210 permits.

The two families where it does pay are the ones that are nearly trees. A path is a tree, so its tree preconditioner is itself and κ is exactly 1 — a saving of the entire condition number, 933 at forty-eight vertices. A cycle is a tree plus one edge, its tree κ is exactly n while κ⁺(L) grows like n², so the win grows linearly with the size: 1.24, 2.45, 3.66, 4.87 at n = 12, 24, 36, 48.

Why the grid is the hard case #

A breadth-first tree of a grid has long paths in it. Two vertices adjacent in the grid can be separated by a tree path that goes all the way to the root and back, so their stretch is of order the diameter, and there are many such edges.

The total stretch of a breadth-first tree on a k × k grid grows like k³ — 60, 120, 210, 336 at k = 4, 5, 6, 7, which is close to k³ divided by a constant — while the graph’s own condition number grows like k². So the bound gets worse relative to the thing it bounds as the graph grows, and the measured condition number does the same.

This is exactly the problem the literature on low-stretch spanning trees exists to solve. There are constructions whose total stretch is O(m log n log log n) rather than O(n^{3/2}), built by a recursive decomposition into balls rather than by a search, and with one of those the whole approach becomes competitive. The trees here — breadth-first and uniformly random — are the two anybody would try first, and they are exactly the two that do not work on a grid.

So the honest summary of this essay’s measurement is: the theorem is correct, the bound is tight enough to choose with, and the obvious trees are the wrong trees. That last part is the reason the subject has a literature.

The lower bound is exact, and once it is not #

λmin = 1 to nine digits in seven of the eight rows above. That is a theorem being observed rather than approximated, and it has a clean reason: the eigenvector for λ = 1 is any vector whose quadratic form is the same over the tree and over the graph, which exists as soon as the graph has an edge the tree does not use.

The eighth row is the exception and is worth explaining. The hypercube’s random tree gives λmin = 1.037772. The theorem says λ ≥ 1, and 1.0378 satisfies it — the bound is not attained, which happens when no vector is supported entirely on the tree’s own edges in the required way. A four-dimensional hypercube is dense enough that a random spanning tree touches every vertex in a way that leaves no such direction.

The distinction between “the bound holds” and “the bound is attained” is the whole difference between those seven rows and the eighth, and it is the kind of thing a measurement finds and a proof does not mention.

The other reading of stretch #

Stretch has a second identity that ties this essay back to effective resistance, and it is not a coincidence.

The stretch of an edge over a tree is the resistance of its tree path. Summed over the graph’s edges and divided by n − 1, it is the average resistance a random edge is routed through. And by Foster’s identity, the sum of effective resistances over a graph’s own edges is exactly n − 1 — so total stretch is the same quantity computed against the tree instead of against the graph, and their ratio measures how much worse the tree is at conducting than the graph is.

That reading makes the bound less mysterious. A tree with total stretch S routes the graph’s currents along paths that are, on average, S/(n − 1) times as resistive as the graph’s own — and the preconditioned condition number is bounded by exactly that inflation.

It also explains why a low-stretch tree is possible at all. Foster’s identity says the graph’s resistances sum to n − 1 whatever the graph is, so a tree that gets close to that sum is one that routes most edges along short paths, and the existence of such trees is a statement about how much of a graph’s connectivity a spanning subgraph can capture.

What a tree cannot do at all #

One limitation is structural rather than quantitative, and it is why combinatorial preconditioning is a graph subject rather than a general one.

The construction needs the matrix to be a Laplacian, or close enough to one that a subgraph makes sense. Symmetric diagonally dominant matrices qualify — they can be written as a Laplacian plus a nonnegative diagonal, so the whole apparatus transfers — and that class covers a useful amount: resistive networks, some finite-element discretisations, the normal equations of certain least-squares problems, image-processing operators.

It does not cover a general symmetric positive definite matrix, nor anything the sparsifier of an earlier anchor would be applied to outside this class. There is no subgraph of a dense covariance matrix, no path along which to route an entry, and no stretch to compute. The moment the off-diagonal signs are mixed, the quadratic form stops being a sum over edges of squared differences and every argument in this essay fails at its first line.

That is a sharper boundary than most preconditioning families have. An incomplete Cholesky applies to anything positive definite; a multigrid hierarchy applies to anything with a notion of smoothness. A tree preconditioner applies to a class defined by a sign condition, and inside that class it has a bound nothing else offers — a condition number known before any arithmetic, from a walk on a graph.

What to do with all this #

Compute the stretch, always. It is a breadth-first search per vertex, it needs no arithmetic, and it is the only number available before a solve that says anything about how the solve will go — in the way a count of spanning trees is available from the same object and says something else.

Compare the tree’s κ against the graph’s own. The measurement above is two eigenvalue computations on a small instance and it is the check that the preconditioner is not making things worse — which, on four of the eight rows here, it was.

Do not use a breadth-first or a random tree on a mesh. Both have stretch growing faster than the graph’s condition number, and the result is a preconditioner that is cheap, correct, and slower than no preconditioner at all.

And remember what the tree buys when it does work. A tree solve is linear and exact, with no fill and no dropped entries, so a tree-preconditioned iteration has a cost per step that no incomplete factorisation can match. The trade is entirely in the iteration count, and the iteration count is bounded by a number that can be computed by walking the graph.

Where the whole idea came from #

A closing note on why a subject exists at all, because the shape of it is instructive.

Solving a Laplacian system is the inner loop of a great deal — every finite-element solve on a diffusion problem, every step of an interior-point method on a network flow, every maximum-flow algorithm of the last decade. So the question of how fast a Laplacian system can be solved has a sharp answer worth chasing, and the answer turned out to be combinatorial rather than numerical: the fastest known methods are built from spanning trees, sparsifiers and recursive graph decompositions, and the numerical linear algebra in them is a preconditioned conjugate gradient iteration that has not changed since 1952.

That is unusual. Almost everywhere else on this site, a faster method is a better factorisation, a better ordering of the same arithmetic, or a lower precision. Here the improvements are theorems about graphs, and the arithmetic is a constant.

It also explains why the two obvious trees fail. The subject’s whole difficulty is concentrated in constructing a tree whose stretch is small, and that construction is a graph algorithm with no numerical content — a recursive decomposition into balls of bounded radius, chosen so that few edges have their endpoints in different balls. Nothing about it would occur to somebody thinking about matrices, and the reason it is needed is entirely visible in the grid rows of the table above.

At other settings #