- Approximate Cholesky
- Graph Laplacian
- Random Spanning Trees
- Schur Complement
- Leverage Scores
- Sparse Solvers
Every time you eliminate a vertex from a graph, its neighbours collapse into a dense little cluster where none may have existed before. Do that thousands of times, as any direct solver for a large sparse system must, and the fill explodes. A new preprint from KTH Royal Institute of Technology and the University of Massachusetts Amherst proposes a stubbornly simple answer. Instead of storing that dense cluster, keep one tree drawn at random from it, and reweight the surviving edges so the tree carries the same behaviour on average. The method is called CAST, short for Canonical Approximate Schur Tree, and on a corpus of two hundred and one real linear systems it runs between roughly one and four times faster than the strongest published competitor.
Key points
- CAST replaces the dense Schur complement clique formed by eliminating a pivot vertex with a single weighted random spanning tree sampled from that clique.
- Each surviving edge is reweighted by the reciprocal of its probability of landing in the tree, which makes the sparse update statistically unbiased.
- Because the clique is induced by a star, the correct spanning tree can be drawn exactly in time proportional to the pivot degree using weighted Prufer codes, so the dense clique is never actually built.
- Leverage score edge marginals are the unique choice that minimises the largest reweighted single edge contribution, giving a clean minimax guarantee.
- A robustness variant called CAST rho splits each neighbour into rho copies, halving a certified error bound for every doubling of rho.
- Across two hundred and one Laplacian and SDDM systems, CAST cut total reuse cost by fourteen to fifteen percent over the state of the art and reached tolerance far faster on the hardest families.
Why eliminating a vertex is so expensive
Start with the object at the centre of this work, a graph Laplacian. Take any weighted undirected graph. Its Laplacian is a matrix that places the sum of a vertex incident weights on the diagonal and the negated edge weight off the diagonal for every pair of connected vertices. Laplacians and their close relatives, the symmetric diagonally dominant M matrices known as SDDM matrices, sit underneath an enormous range of computation. They show up in the Newton steps of interior point methods for maximum flow, in finite element discretisations, in spectral clustering, in semi supervised learning on graphs, and in physical simulation. Solving a system built on one of these matrices quickly is a workhorse operation, not a niche one.
The classical way to solve such a system exactly is Cholesky factorisation, which is Gaussian elimination specialised to symmetric positive matrices. You pick a pivot vertex, eliminate it, and update the rest of the matrix to account for its removal. Here is where the graph view makes the pain visible. When you eliminate a vertex from a graph, every pair of its former neighbours becomes connected, because the paths that used to run through the eliminated vertex now have to be represented directly. A vertex with degree d turns its d neighbours into a clique with on the order of d squared new edges. On a sparse graph that starts with only a handful of edges per vertex, repeated elimination can fill the matrix until it is effectively dense, and the memory and arithmetic both blow up.
The escape route that has driven twenty years of nearly linear time Laplacian solvers is to stop insisting on an exact factor. Approximate elimination keeps a sparse stand in for each dense clique, accepting a small controlled error in exchange for a factor that stays sparse and can be applied cheaply as a preconditioner. The preconditioner does not need to be exact. It only needs to be close enough that an iterative method such as preconditioned conjugate gradient converges in a small number of steps.
The one move that decides everything
Approximate Cholesky methods differ almost entirely in one decision. When a pivot vertex is eliminated and its neighbours would form a dense clique, what sparse object do you put in its place? The exact clique has an edge for every pair of neighbours, and each of those edges has a specific weight determined by the pivot. Concretely, if the pivot is joined to neighbour p with conductance c sub p, then the exact Schur complement places an edge between neighbours p and q whose weight is the product of their conductances divided by the total conductance at the pivot.
The leading practical family, the AC method of Gao, Kyng and Spielman that ships inside the widely used Laplacians.jl library, replaces that clique with a random connected structure grown by pairing edges together as they are processed. It works well and comes with strong guarantees, but it has a subtle property. Because the tree is built by sequentially pairing incident edges, the probability that any given edge ends up in the sparse update depends on the order in which the pivot neighbours happen to be processed. The distribution the method samples from is a side effect of the traversal, not a quantity chosen on purpose.
CAST asks a sharper question. Forget how the tree is grown. What is the ideal distribution over sparse replacements, chosen deliberately, and can we sample from it directly without ever forming the clique? The answer the authors give is that the replacement should be a single weighted random spanning tree of the clique, and the marginals, meaning the inclusion probability of each edge, should equal the edge leverage scores.
The dense clique is never built, never stored, and never touched. CAST samples one tree straight from the star that would have created it.
Reweighting is what keeps it honest
A spanning tree of the neighbour clique has only as many edges as there are neighbours minus one, which is a dramatic reduction from the roughly d squared edges of the full clique. But a tree is a biased picture of the clique unless you correct for it. Some edges are far more likely to be chosen than others, and if you simply keep the chosen edges at their natural weight, the update systematically misrepresents the matrix. CAST fixes this with inverse marginal reweighting. Every edge that lands in the sampled tree gets its weight divided by the probability that it would have been selected.
This is the same idea that makes importance sampling unbiased. A rare event, when it does occur, has to count for more, and the reciprocal of its probability is exactly the right amount more. The consequence is that the expected value of the sparse update is the true dense clique. Each individual elimination introduces zero bias, and the only thing left to control is variance, which is where the choice of marginals and the robustness variant come in.
There is a beautiful reason the marginals have a closed form here. The clique among the neighbours is not an arbitrary weighted graph. It is the Schur complement of a star, the pivot at the centre joined to each neighbour. For a star the effective resistance between any two leaves is just the sum of the two reciprocal conductances along the unique path, and the transfer current theorem then says the probability that an edge appears in a weighted random spanning tree equals its conductance times that effective resistance. Everything collapses into quantities you already have.
Sampling a tree in linear time with Prufer codes
Sampling a weighted random spanning tree on a general graph is not cheap. The best known methods lean on nontrivial graph algorithmic machinery and cost far more than a single elimination step can afford. CAST sidesteps the general problem entirely by exploiting the fact that the clique it needs to sample from is complete, meaning every neighbour is joined to every other, and its edge weights factor into a product over endpoints.
For complete graphs with product form edge weights there is a classical encoding called the Prufer code. A Prufer code is a short sequence that is in exact one to one correspondence with the labelled spanning trees of a complete graph. Encode a tree by repeatedly removing its smallest leaf and recording that leaf neighbour, and you get a sequence of length two less than the number of vertices. Decode by reversing the process. The key fact the authors use is that a vertex appears in the code exactly one fewer time than its degree in the tree, which ties the code statistics directly to tree structure.
Because of that correspondence, drawing each symbol of the code independently from the right distribution over neighbours produces a tree drawn from exactly the weighted spanning tree distribution you want. There is no rejection, no iterative construction over the clique, and no clique. The cost is proportional to the pivot degree, which is the best you could hope for, since you have to write down a tree with that many edges regardless.
The robustness dial, CAST rho
A single spanning tree is a lean estimator, and on some matrices leanness costs you. The variance of a one tree update can be large when a pivot has a heavy tailed neighbourhood, and large variance means the preconditioner is a rougher approximation, which means more iterations to converge. CAST adds a dial for exactly this situation. Rather than sampling one tree on the original neighbours, CAST rho splits each neighbour into rho equal conductance copies, samples a single spanning tree on this expanded vertex set, and then contracts the copies back onto the original neighbours.
The effect is to spread each neighbour contribution across rho smaller pieces, which shrinks the largest single edge contribution in direct proportion to rho. The authors prove a local second moment bound that scales as one over rho, so every doubling of rho halves the certified variance bound. The base case rho equal to one recovers the plain single tree method and its unit local risk.
Crucially the expanded sampling stays cheap. The Prufer sampler on rho copies of d neighbours still runs in time proportional to rho times d, because a symbol can be drawn hierarchically, first picking which original neighbour and then which of its rho copies uniformly. The robustness never breaks the linear time budget that made the base method attractive.
A useful way to read CAST rho is as a knob between two familiar regimes. At rho equal to one you get the sparsest possible unbiased update and the lowest construction cost. As rho grows you spend more work per pivot to buy a smoother preconditioner, which pays off precisely when extra iterations would otherwise pile up. The paper is careful to note this is a structural choice rather than a strictly better setting, since splitting can increase downstream fill when a pivot has a heavy degree tail.
Why leverage scores are the right marginals
One of the more satisfying results in the paper is that the leverage score marginals are not just convenient, they are optimal in a precise sense. Among all unbiased single tree estimators built by inverse marginal reweighting, the leverage score choice uniquely minimises the largest normalised reweighted edge contribution. That is a minimax statement. It picks the distribution that makes the worst case single edge as small as possible, which is exactly the quantity that matrix concentration analyses of approximate elimination care about.
The intuition is a balancing act. If you make an edge very unlikely to be selected, then on the rare occasions it is selected its reweighted contribution is huge, because you divide by a tiny probability. If you make it very likely, it barely gets scaled up but crowds out other edges. The leverage scores are the sweet spot where no single edge can blow up the update, and the star structure hands them to you in closed form. The combination of an optimality guarantee and a closed form is unusual, and it is the theoretical heart of why CAST is more than an engineering tweak.
A reference implementation of the CAST update
The clearest way to see how little machinery CAST needs is to write the single pivot update in plain code. The routine below takes the conductances from a pivot to its neighbours, draws a weighted random spanning tree of the induced clique using a Prufer code, and reweights each selected edge by the reciprocal of its inclusion probability. Notice that the dense clique is never formed. The whole thing runs in time proportional to the number of neighbours.
import numpy as np
def cast_schur_tree(neighbors, cond, rng):
"""One CAST elimination step for a star induced clique.
neighbors : array of the pivot neighbour ids, length d
cond : conductances from the pivot to each neighbour, length d
returns : list of (p, q, reweighted_weight) sparse update edges
"""
d = len(neighbors)
if d < 2:
return [] # nothing to connect
a = cond.sum() # total conductance at the pivot
p_select = cond / a # Prufer symbol distribution
# A star induced complete graph on d leaves needs a code of length d - 2.
code = rng.choice(d, size=max(d - 2, 0), p=p_select)
# Decode the Prufer code into tree edges by the leaf pointer algorithm.
degree = np.ones(d, dtype=int)
for s in code:
degree[s] += 1
edges = []
import heapq
leaves = [i for i in range(d) if degree[i] == 1]
heapq.heapify(leaves)
for s in code:
leaf = heapq.heappop(leaves)
edges.append((leaf, s))
degree[leaf] -= 1
degree[s] -= 1
if degree[s] == 1:
heapq.heappush(leaves, s)
u = heapq.heappop(leaves)
v = heapq.heappop(leaves)
edges.append((u, v)) # last two labels close the tree
# Reweight every selected edge by the reciprocal of its inclusion prob.
update = []
for i, j in edges:
ci, cj = cond[i], cond[j]
w_exact = ci * cj / a # true Schur weight for this pair
p_incl = (ci + cj) / a # closed form inclusion probability
w_hat = w_exact / p_incl # unbiased inverse marginal reweight
update.append((neighbors[i], neighbors[j], float(w_hat)))
return update
def demo():
rng = np.random.default_rng(0)
neighbors = np.array([11, 12, 13, 14, 15])
cond = np.array([0.9, 0.4, 1.3, 0.6, 0.8])
# Average many CAST updates and confirm the expected clique is recovered.
d = len(neighbors)
a = cond.sum()
accum = np.zeros((d, d))
trials = 200000
idx = {v: i for i, v in enumerate(neighbors)}
for _ in range(trials):
for p, q, w in cast_schur_tree(neighbors, cond, rng):
accum[idx[p], idx[q]] += w
accum[idx[q], idx[p]] += w
accum /= trials
exact = np.zeros((d, d))
for i in range(d):
for j in range(i + 1, d):
w = cond[i] * cond[j] / a
exact[i, j] = exact[j, i] = w
print("max abs error vs exact clique:", np.abs(accum - exact).max())
if __name__ == "__main__":
demo()
Run the demo and the averaged sparse updates converge to the exact dense clique, which is the unbiasedness property made concrete. In a real solver this routine is called once per pivot inside the elimination loop, the reweighted edges are added to the residual graph, and the eliminated pivot rows and columns become part of the sparse Cholesky factor. The factor is then reused as a preconditioner across many right hand sides.
How it performs on real systems
The evaluation is the part that turns a clean idea into a claim worth repeating. The authors assembled a corpus of two hundred and one SDDM and Laplacian systems drawn from a public benchmark suite, spanning four families. There are SuiteSparse matrices from real applications, a large batch of Newton step Laplacians produced inside maximum flow interior point solvers, another set of Spielman interior point matrices, and a handful of adversarial star structured instances designed to stress the method. Every method in the comparison shares the same elimination ordering, factor format, and conjugate gradient implementation, so the comparison isolates exactly one thing, the local update.
The headline is that CAST at base granularity improves on AC across every collection, and does so uniformly rather than on average. It reaches a fourteen to fifteen percent reduction in total reuse cost on the heterogeneous collections, and that saving concentrates in the solve phase, which means it compounds as the number of right hand sides grows. When you factor a matrix once and then solve against it hundreds of times, a cheaper solve is the gift that keeps giving.
| Collection | Elimination profile | CAST base vs AC | CAST rho vs AC2 |
|---|---|---|---|
| SuiteSparse | Mostly low degree neighbourhoods | 1.17 times faster | 1.07 times faster |
| Chimera interior point | Moderate degree, some heavy tails | 1.18 times faster | 0.65 times, AC2 wins |
| Spielman interior point | All pivots small degree | 1.27 times faster | 4.43 times faster |
| Sachdeva star | Clique dominated neighbourhoods | 3.09 times faster | 1.90 times faster |
The doubled granularity story is more nuanced, and the paper is admirably honest about it. On the Spielman family, where every pivot has small degree, splitting is nearly free and the payoff is enormous, with CAST rho reaching tolerance in one or two iterations and running four times faster than the corresponding robust baseline. On the star structured Sachdeva family, where neighbourhoods are clique dominated, splitting again helps by holding iteration growth flat as the instances scale. But on the Chimera family, where elimination produces a heavy tail of high degree pivots, the extra terminal contributions dominate and the robust baseline wins on most instances. The lesson the authors draw is that rho is a structural choice tied to the degree profile, not a universally better setting.
One honest limitation deserves a flag. The guarantees in the paper are local certificates. They bound the error introduced by a single elimination step, and they say nothing directly about how those errors accumulate across the full sequence of eliminations, nor about the downstream fill that splitting can create. The strong end to end speedups are empirical. The theory explains why each step is well behaved and why leverage marginals are optimal locally, but the global behaviour is established by measurement on the corpus rather than proved.
Where CAST sits in the landscape
It helps to place this work against its neighbours. The nearly linear time Laplacian solver line began with support graph preconditioners and low stretch spanning trees, and later spectral sparsification through effective resistance sampling gave strong global approximation guarantees. Those methods build recursive global structures. CAST is not in that family. It lives inside the randomised approximate elimination line, which builds a sparse factor by eliminating vertices one at a time and replacing each dense clique with a sparse random update.
Within that line, the closest predecessor is the AC framework that CAST is measured against. Both preserve connected support and both trade construction cost against robustness with a granularity setting. The difference is the distribution. AC generates its trees by randomised sequential edge pairing, so the induced edge marginals depend on processing order. CAST samples from the weighted random spanning tree distribution of the clique itself, which is order independent and attains the minimax optimal leverage score marginals. That single conceptual shift, from a distribution that emerges from an algorithm to one chosen on purpose and sampled exactly, is the contribution.
If you have followed related threads on scaling graph computation, the flavour will feel familiar. Work on universal graph coarsening with hashing chases the same goal of shrinking a graph while preserving its essential shape, and studies of the graph Laplacian in spectral clustering lean on the very same spectral objects that CAST manipulates. The through line across all of these is that the Laplacian spectrum, effective resistances, and leverage scores are the currency of fast graph computation, and whoever samples from them most cleanly wins.
What to take away
CAST is a reminder that in numerical computation the distribution you sample from is a design decision worth agonising over, not an accident of implementation. The authors identified that approximate Cholesky methods had been sampling their sparse updates from whatever distribution their tree building procedure happened to induce, asked what the ideal distribution would be, and discovered that for the star induced cliques these methods actually face, the ideal is a leverage score weighted random spanning tree that can be sampled exactly in linear time. The result is a method that is unbiased at every step, provably optimal in its local risk, and consistently faster on a large corpus of the linear systems that solvers meet in practice.
For a practitioner the practical read is straightforward. If you are solving many right hand sides against a fixed Laplacian or SDDM matrix, and your neighbourhoods are not dominated by a heavy tail of very high degree pivots, CAST at base granularity is a low risk drop in that should shave a meaningful slice off your reuse cost, mostly in the solve phase where it compounds. If your pivots are uniformly small degree or your cliques are dense, the robustness variant can turn a modest gain into a large one. And if your elimination produces a heavy degree tail, measure before committing, because that is the regime where the existing robust baseline can still come out ahead.
Follow the research that actually ships
We turn dense preprints like this one into clear, faithful explainers you can act on. Subscribe for weekly breakdowns of the graph and numerical methods quietly making large scale computation cheaper.
Read more on aitrendblendFrequently asked questions
What does CAST actually replace inside approximate Cholesky?
When Cholesky elimination removes a pivot vertex, its neighbours would form a dense clique with a weight on every pair. CAST replaces that entire dense clique with a single weighted random spanning tree drawn from it, and reweights each surviving edge by the reciprocal of its inclusion probability so the sparse update matches the dense clique in expectation.
Why is CAST able to sample a spanning tree so cheaply?
The clique among a pivot neighbours is induced by a star, so it is a complete graph whose edge weights factor into a product over endpoints. Complete graphs with product form weights can be encoded by Prufer codes, and drawing each code symbol independently from the endpoint weight distribution decodes to exactly the weighted random spanning tree you want. The cost is proportional to the pivot degree.
What is the point of the CAST rho variant?
CAST rho controls variance. It splits each neighbour into rho equal conductance copies, samples one spanning tree on the expanded set, then contracts the copies back. This shrinks the largest single edge contribution in proportion to rho, and the paper proves a local second moment bound that scales as one over rho, so each doubling of rho halves the certified error bound while keeping the sampler linear in cost.
How much faster is CAST than existing methods?
On a corpus of two hundred and one Laplacian and SDDM systems, CAST at base granularity beat the AC baseline on every collection, cutting total reuse cost by fourteen to fifteen percent on the heterogeneous families. The robustness variant reached tolerance up to about four times faster on small degree families, though the older robust baseline still won on collections with a heavy tail of high degree pivots.
Are the speed claims backed by proofs or by experiments?
Both, at different levels. The proofs are local certificates. They show each elimination step is unbiased with connected support and that leverage score marginals uniquely minimise the local risk, plus a one over rho second moment bound. The end to end speedups on full systems are empirical, measured on the corpus, because the theory does not directly bound how per step errors accumulate or how much fill splitting adds.
Should I be cautious that this is a preprint?
Yes, treat the results as promising rather than settled. This work is posted as an arXiv preprint and has not yet completed peer review, so the numbers and proofs have not been independently vetted by referees. The method builds on well established ideas in approximate elimination and the evaluation is on a public benchmark, which is reassuring, but independent replication is the right standard before relying on it in production.
Source paper. Meher Chaitanya, Cameron Musco, and Aristides Gionis. CAST, Canonical Approximate Schur Tree for Approximate Cholesky on Graphs. arXiv preprint arXiv:2609.09255v1, 8 September 2026. Affiliations, KTH Royal Institute of Technology and University of Massachusetts Amherst. This article explains a preprint that has not yet been peer reviewed and is provided for general technical education, not as engineering or financial advice.