Matrix-Tree Theorem, in depth: counting spanning trees with an exact determinant, weighted and directed forms, edge marginals and dependency parsing

Statement: the Laplacian and its cofactors

Let G be an undirected graph on n vertices, possibly with parallel edges, and with a weight w(e) on each edge (use 1 for an unweighted graph). The Laplacian L is the n by n matrix with L[i][i] equal to the total weight of edges at vertex i and L[i][j] equal to minus the total weight of edges between i and j. Equivalently L = D - A, the degree matrix minus the adjacency matrix. Self-loops are ignored, because no spanning tree can use one.

Every row of L sums to zero, so L is singular. The theorem says that if you delete row i and column i for any i, the determinant of the remaining (n-1) by (n-1) matrix equals the sum, over all spanning trees T, of the product of the weights of T's edges. With unit weights that sum is simply the number of spanning trees, written t(G). Two equivalent forms are worth knowing: t(G) equals the product of the n-1 nonzero eigenvalues of L divided by n, and the graph is connected exactly when t(G) is positive, which is the same as the second-smallest Laplacian eigenvalue being nonzero.

A quick sanity check is Cayley's formula. The complete graph K_n has n^n-2 spanning trees: 16 for K₄ and 1,296 for K₆. The code below reproduces both.

Why it is true

The cleanest proof goes through the incidence matrix. Orient each edge arbitrarily and let B be the n by m matrix with +1 at the edge's head, -1 at its tail and 0 elsewhere. Then L = B W B^T, where W is the diagonal matrix of edge weights. Delete vertex i's row from B to get B_i; the reduced Laplacian is B_i W B_i^T.

The Cauchy-Binet formula expands the determinant of a product of an (n-1) by m and an m by (n-1) matrix as a sum over all choices S of n-1 columns: the product of the two corresponding (n-1) by (n-1) minors. Each choice S is a set of n-1 edges. A classical lemma says the minor of B_i on S is plus or minus 1 if S forms a spanning tree and 0 otherwise (if S contains a cycle, the columns are linearly dependent; if it is a tree, you can peel leaves to triangularise it). The two minors carry the same sign, so each tree contributes +1 times its weight product and everything else contributes 0. That is the theorem.

Computing the determinant exactly

The obvious way to evaluate the determinant, Gaussian elimination in floating point, is the wrong tool here. Tree counts grow exponentially, so a double runs out of precision quickly and returns a rounded number that looks plausible. You want exact integers, and there are two good ways to get them. Fraction-free Bareiss elimination keeps every intermediate value an integer by dividing each update by the previous pivot, a division that is always exact. It costs O(n³) arithmetic operations on integers whose size grows only linearly. Alternatively, compute the determinant modulo several large primes and reconstruct with the Chinese remainder theorem, or stop at one prime when you only need the count mod p, which is common in programming contests. The pivoting background is in Gaussian elimination.

def bareiss_det(M):
    """Exact determinant of an integer matrix, fraction-free (Bareiss)."""
    A = [row[:] for row in M]
    n = len(A)
    if n == 0:
        return 1
    sign, prev = 1, 1
    for k in range(n - 1):
        if A[k][k] == 0:                       # swap in a nonzero pivot row
            for i in range(k + 1, n):
                if A[i][k] != 0:
                    A[k], A[i] = A[i], A[k]
                    sign = -sign
                    break
            else:
                return 0
        for i in range(k + 1, n):
            for j in range(k + 1, n):
                A[i][j] = (A[i][j] * A[k][k] - A[i][k] * A[k][j]) // prev   # exact
        prev = A[k][k]
    return sign * A[n - 1][n - 1]

def laplacian(n, edges):                       # edges: (u, v, weight)
    L = [[0] * n for _ in range(n)]
    for u, v, w in edges:
        if u != v:
            L[u][u] += w; L[v][v] += w
            L[u][v] -= w; L[v][u] -= w
    return L

def count_spanning_trees(n, edges, drop=0):
    L = laplacian(n, edges)
    minor = [r[:drop] + r[drop + 1:] for i, r in enumerate(L) if i != drop]
    return bareiss_det(minor)

Always test a counter against brute force on small graphs: enumerate every subset of n-1 edges, keep those that are acyclic using a union-find, and sum their weight products. The same union-find drives Kruskal's algorithm. On every graph in this article the two methods agree.

Worked example: a square with a chord

Take a square on vertices 0, 1, 2, 3 with an extra chord from 0 to 2. Vertex 0 and vertex 2 have degree 3, the others degree 2, which gives the Laplacian in the figure. Deleting row and column 0 leaves [[2, -1, 0], [-1, 3, -1], [0, -1, 2]]. Its determinant is 2(6 - 1) - (-1)(-2 - 0) = 10 - 2 = 8.

Check by hand. There are 5 edges and C(5, 3) = 10 ways to choose 3 of them. Exactly two of those choices contain a cycle: the triangle 0-1-2 and the triangle 0-2-3. That leaves 8 trees, matching the determinant. Now give the edges weights 0-1: 2, 1-2: 1, 2-3: 3, 3-0: 1 and 0-2: 1. The code returns a weighted total of 29, and brute force agrees. Scaling up, the 3 by 3 grid graph has 192 spanning trees, and the 10 by 10 grid has 5,694,319,004,079,097,795,957,215,725,765,328,371,712,000, a 43-digit number that Bareiss computes exactly in a fraction of a second; its residue modulo the prime 2⁶¹ - 1 matches an independent mod-p elimination.

Directed graphs: Tutte's theorem

Tutte extended the theorem to directed graphs. An arborescence rooted at r is a spanning tree in which every vertex other than r has exactly one incoming arc and every vertex is reachable from r. Build the in-degree Laplacian: put the total weight of arcs entering v on the diagonal at v, and minus the weight of arc u to v at position (u, v). Delete r's row and column; the determinant counts arborescences rooted at r. Unlike the undirected case, the answer depends on which row you delete, because it depends on the root.

def count_arborescences(n, arcs, root):
    """arcs: (u, v, w) meaning u -> v. Counts spanning arborescences rooted at root."""
    L = [[0] * n for _ in range(n)]
    for u, v, w in arcs:
        if u != v:
            L[v][v] += w            # in-degree on the diagonal
            L[u][v] -= w
    minor = [r[:root] + r[root + 1:] for i, r in enumerate(L) if i != root]
    return bareiss_det(minor)

arcs = [(0, 1, 1), (0, 2, 1), (1, 2, 1), (2, 1, 1), (1, 3, 1), (2, 3, 1)]
count_arborescences(4, arcs, 0)     # 6, confirmed by enumerating parent choices

Finding the single best arborescence is a different problem, solved by Chu-Liu-Edmonds. The matrix-tree theorem gives you the sum over all of them, which is what probabilistic models need.

Edge marginals and dependency parsing

Suppose you pick a spanning tree with probability proportional to its weight product. What is the probability that edge e = (u, v) is in it? The answer is w(e) times the effective resistance between u and v when each edge is a resistor of conductance w(e). You can read effective resistance off the inverse of the reduced Laplacian: with vertex 0 grounded and G = (L with row and column 0 removed)^-1, padded with zeros for vertex 0, R(u, v) = G[u][u] + G[v][v] - 2G[u][v].

In the square with a chord, each side lies in 5 of the 8 trees and the chord in 4. The exact computation gives marginals 5/8 for the four sides and 1/2 for the chord, and they sum to 3 = n - 1, as every spanning tree has n - 1 edges. That sum is a useful assertion in any implementation. The connection to sampling, including Wilson's algorithm, is developed in random spanning trees.

This is exactly how non-projective dependency parsers are trained. A sentence's words are vertices, a neural scorer gives each head-to-dependent arc a score s, and the arc weight is exp(s). The log of Tutte's determinant is the log partition function of a distribution over all dependency trees, and its gradient with respect to each arc score is that arc's marginal probability. Koo et al. and Smith and Smith both showed this in 2007. In practice, compute log-det with an LU factorisation in floating point (not Bareiss), let automatic differentiation produce the marginals, and handle the single-root constraint by a modification of the root row.

Failure modes

• Floating-point determinants for counting. Exact counts overflow double precision quickly; use Bareiss or modular arithmetic. Floating point is fine for log-det in learning, where you want gradients, not exact counts.
• Forgetting parallel edges or loops. Parallel edges add their weights into the off-diagonal entry; loops must be skipped or the diagonal is wrong.
• Disconnected graphs. The determinant is 0; check connectivity first if a zero would be confusing to callers.
• Wrong Laplacian for directed graphs. Out-degree versus in-degree decides whether you count trees pointing away from or toward the root.
• Ill-conditioned inverses. For marginals on large graphs with very unequal weights, solve linear systems instead of forming an explicit inverse, and work in log space for the weights.

Trade-offs: choosing a method

Special structure beats all of these. Series-parallel graphs and cactus graphs have counts given by simple products over their cycles, and families such as complete, complete bipartite and grid graphs have closed forms from their eigenvalues.

What to do next

1. Implement laplacian and bareiss_det, and check them against Cayley's formula for K₃ to K₇.
2. Write a brute-force counter with union-find and compare on 50 random small graphs, including parallel edges and loops.
3. Add the directed version and verify it against enumeration of parent choices.
4. Compute edge marginals via the grounded inverse and assert that they sum to n - 1.
5. For an ML use case, rewrite the computation with torch.logdet and confirm that autograd gradients match the marginals from the inverse.