CAST: Canonical Approximate Schur Tree for Approximate Cholesky on Graphs

Meher Chaitanya Email: mcpi@kth.se Affiliation: KTH Royal Institute of Technology Affiliation: Stockholm, Sweden Cameron Musco Email: cmusco@cs.umass.edu Affiliation: University of Massachusetts Amherst Affiliation: Amherst, Massachusetts, USA Aristides Gionis Email: argioni@kth.se Affiliation: KTH Royal Institute of Technology Affiliation: Stockholm, Sweden

Abstract

Graph-data workloads such as diffusion estimation, ranking, semi-supervised learning, and network optimization often solve many Laplacian or symmetric diagonally dominant M-matrix (SDDM) systems with the same coefficient matrix. Approximate Cholesky preconditioners eliminate vertices one at a time and store the resulting sparse approximate factorization, the factor, whose construction cost is amortized across these solves. But eliminating a vertex, the pivot, creates a dense Schur-complement clique among its $d$ active neighbors. We introduce CAST (Canonical Approximate Schur Tree), which replaces this clique with a weighted random spanning tree sampled directly from it. Every realization is connected and contains exactly $d-1$ edges, while reweighting each selected edge by the reciprocal of its tree-inclusion probability makes the update unbiased. The distribution is independent of the ordering of the pivot neighbors, and we prove that its leverage-score marginals minimize the largest normalized reweighted-edge contribution among unbiased inverse-marginal one-tree estimators.

We also introduce CAST- $\rho$ , which replaces each pivot neighbor with $\rho$ copies, each carrying a $1/\rho$ share of that neighbor’s incident weight, samples a weighted random spanning tree on the expanded clique, and contracts the copies back to the original neighborhood. The resulting update remains unbiased and connected, can be sampled exactly in $O(\rho d)$ time, and satisfies a $1/\rho$ bound on the second moment of the normalized local Schur error. Increasing $\rho$ therefore reduces certified local sampling variability, but may increase construction cost and downstream fill.

We evaluate CAST on $201$ matrix systems, including $168$ Newton-step Laplacians from maximum-flow interior-point methods, constructing one factor and reusing it for $250$ right-hand sides. CAST- $\rho$ is $1.11$ – $4.43\times$ faster than the better of the state-of-the-art AC and AC2 baselines Gao et al. (2026). Overall, we observe that CAST-1 is the faster default, whereas CAST-2 is preferable when its additional edge contributions remain inexpensive.

1 Introduction

Many graph-data and network-optimization pipelines repeatedly solve linear systems defined by a fixed graph operator. Examples include ranking and diffusion estimation Page et al. (1999); Andersen et al. (2006), label propagation Zhu et al. (2003); Zhou et al. (2003), shift-and-invert spectral graph methods Huang et al. (2019), and network optimization Spielman and Teng (2004); Spielman and Srivastava (2008). These applications lead to Laplacian or symmetric diagonally dominant M-matrix (SDDM) systems in which the operator is reused across many right-hand sides. The relevant computational objective is therefore not only to solve one system quickly, but to construct a preconditioner whose cost can be amortized across repeated solves.

Approximate Cholesky factorization provides a practical approach to this problem. As in sparse Gaussian elimination, vertices are eliminated sequentially, but the dense fill created by exact elimination is replaced by sparse randomized updates Kyng and Sachdeva (2016); Gao et al. (2026). The resulting approximate factor is used as a conjugate-gradient preconditioner and reused across multiple right-hand sides. Eliminating a vertex $v$ of degree $d$ , with incident edge weights $a_{1},\ldots,a_{d}$ and total weight $a=\sum_{i=1}^{d}a_{i}$ , creates the Schur-complement clique

K_{v}\;=\;\sum_{1\leq i<j\leq d}\frac{a_{i}a_{j}}{a}\,(\mathbf{e}_{i}-\mathbf{e}_{j})(\mathbf{e}_{i}-\mathbf{e}_{j})^{\top},

where $\mathbf{e}_{i}$ denotes the standard basis vector of the $i$ -th neighbor. This clique contains $\Theta(d^{2})$ edges. Approximate Cholesky methods avoid this quadratic fill by replacing $K_{v}$ with a random surrogate containing only $\mathcal{O}(d)$ edges.

Kyng and Sachdeva sample clique edges independently and control the accumulated error through a matrix-martingale analysis Kyng and Sachdeva (2016). The practical AC and AC2 solvers Gao et al. (2026) take a complementary approach: every sampled local update has connected support, the graph on the $d$ pivot neighbors formed by its nonzero-weight edges. A graph Laplacian has nullspace $\operatorname{span}\{\mathbf{1}\}$ exactly when its support is connected, matching the exact Schur clique. Since a spanning tree is the sparsest connected support, with exactly $d-1$ edges, we first study surrogates consisting of a single spanning tree $T$ of the clique.

For a clique edge $e$ with weight $w_{e}$ , let $p_{e}=\mathbb{P}\left(e\in T\right)>0$ be its inclusion probability, and assign $e$ the weight $w_{e}/p_{e}$ whenever it is selected. This inverse-marginal reweighting makes the surrogate unbiased, since the expected contribution of $e$ is $p_{e}(w_{e}/p_{e})=w_{e}$ . The remaining choice is the distribution over trees. If $p_{e}$ is small, edge $e$ receives a correspondingly large weight when selected. Let $\tau_{e}$ be the leverage score of $e$ in the clique, equal to its weight times its effective resistance within $K_{v}$ Spielman and Srivastava (2008). The normalized size of the reweighted contribution—its size measured relative to the quadratic form of $K_{v}$ itself—is exactly $\tau_{e}/p_{e}$ ; see Section 2. A large ratio lets one selected edge dominate the local update and inflates the error bounds obtained from matrix-concentration analyses Kyng and Sachdeva (2016). We therefore ask:

Question 1.1.

Among unbiased inverse-marginal one-tree replacements of a Schur clique, which edge-inclusion probabilities minimize the worst normalized contribution $\max_{e}\tau_{e}/p_{e}$ ?

The connected updates of AC and AC2 are generated through randomized sequential edge pairing Gao et al. (2026). Although these methods perform strongly across many SDDM systems, their induced sampling distributions depend on the order in which the incident edges are processed and are not designed to minimize $\max_{e}\tau_{e}/p_{e}$ . The distinguishing feature of our method is therefore not connectivity alone, but the distribution used to generate the local update.

We introduce CAST (Canonical Approximate Schur Tree): at each elimination step, it samples a weighted random spanning tree of the Schur clique and applies inverse-marginal reweighting. By the transfer-current theorem Lyons (2003), each clique edge is then included with probability equal to its leverage score, $p_{e}=\tau_{e}$ . This answers Question 1.1: the leverage-score marginals are the unique locally minimax-optimal choice among unbiased inverse-marginal one-tree estimators (Theorem 4.5). We call the sampling rule canonical because it depends only on the weighted Schur clique and not on the ordering of the pivot neighbors.

We further introduce CAST- $\rho$ , with integer parameter $\rho\geq 1$ , which splits each pivot neighbor into $\rho$ equal-share auxiliary copies, samples one weighted random spanning tree from the expanded Schur clique, and contracts the copies back to the original neighbors. For $\rho=1$ it coincides with CAST; for $\rho>1$ the contracted update need not be a tree, but it remains unbiased and connected and satisfies a local second-moment bound of order $1/\rho$ , reducing sampling variability at the cost of additional construction work and possible downstream fill. We write CAST-1 and CAST-2 for the $\rho=1$ and $\rho=2$ configurations used in our evaluation.

Our guarantees concern the estimator associated with a single Schur-clique replacement; we do not claim a worst-case spectral guarantee for the complete approximate factorization.

In summary, we make the following contributions:

1.

Canonical tree sampling and its splitting extension. We introduce CAST and CAST- $\rho$ . At a degree- $d$ pivot, a tree is sampled exactly from the target weighted spanning-tree distribution and contracted in $\mathcal{O}(\rho d)$ time using weighted Prüfer codes, without materializing either the original or the expanded dense Schur clique.
2.

Local minimax and second-moment guarantees. For the base $\rho=1$ update, the leverage-score vector is the unique marginal vector minimizing the largest normalized reweighted-edge contribution among unbiased inverse-marginal one-tree estimators (Theorem 4.5). For general $\rho$ , CAST- $\rho$ remains unbiased and connected and satisfies a normalized local second-moment bound of order $1/\rho$ , quantifying the reduction in sampling variability obtained through splitting (Section 4).
3.

Reuse-aware evaluation. We evaluate CAST by constructing one preconditioner and reusing it across many right-hand sides. Across heterogeneous SDDM systems and max flow IPM sequences, our experiments identify when CAST reduces repeated-solve cost relative to AC and AC2. We further characterize when increasing $\rho$ improves convergence and when its additional edge contributions create excessive downstream fill.

Independent contemporaneous work.

We recently became aware of independent work by Baumann, Kyng, and Zöcklein Baumann et al. (2026), who propose, under the name VAC, the same base $\rho=1$ local sampling rule as CAST-1: a weighted random spanning tree of the product-form Schur clique sampled via a weighted Pr"ufer code and equipped with inverse-marginal reweighting. Their work was developed independently of ours and establishes linear-time sampling rule. Our work additionally gives a local minimax characterization of the leverage-score marginals, introduces the CAST- $\rho$ splitting construction and its $1/\rho$ local second-moment guarantee, and evaluates these estimators as reusable preconditioners across a broad collection of SDDM systems.

2 Preliminaries and Problem Setting

We formalize SDDM linear systems, the factor-reuse setting, and the local Schur-clique replacement problem underlying CAST. The preliminaries on Prüfer codes are provided in Appendix A.

Graph Laplacians and SDDM systems.

Let $G=(V,E,w)$ be an undirected weighted graph with $n=|V|$ vertices. For an oriented edge $e=(u,v)$ write $\mathbf{b}_{e}=\mathbf{e}_{u}-\mathbf{e}_{v}$ , where $\mathbf{e}_{u}$ is the $u$ -th standard basis vector. The Laplacian of $G$ is $L_{G}=\sum_{e\in E}w_{e}\,\mathbf{b}_{e}\mathbf{b}_{e}^{\top}$ . If $G$ is connected, then $\ker(L_{G})=\operatorname{span}\{\mathbf{1}\}$ . Since $L_{G}$ is symmetric, $\operatorname{range}(L_{G})=\ker(L_{G})^{\perp}=\mathbf{1}^{\perp}$ , so the system $L_{G}\mathbf{x}=\mathbf{b}$ is solvable if and only if $\mathbf{b}\perp\mathbf{1}$ .

Throughout, $A$ denotes the coefficient matrix of the system to be solved: either a graph Laplacian $L_{G}$ or a symmetric diagonally dominant M-matrix (SDDM matrix), i.e. a symmetric positive-definite matrix with nonpositive off-diagonal entries satisfying $A_{ii}\geq\sum_{j\neq i}|A_{ij}|$ . Equivalently, an SDDM matrix can be written as $A=L_{G}+\Gamma$ , where $L_{G}$ is a graph Laplacian and $\Gamma\succeq 0$ is a diagonal matrix. The matrix $A$ is positive definite provided that $\Gamma$ has at least one positive entry on each connected component of $G$ . This class of matrices includes many shifted and regularized graph operators. For example, let $W\in\mathbb{R}^{n\times n}$ be a symmetric, nonnegative weighted adjacency matrix, let $D=\operatorname{diag}(W\mathbf{1})$ be its degree matrix, and let $L=D-W$ be a graph Laplacian. For a damping parameter $\alpha\in(0,1)$ and a seed distribution $\mathbf{s}\in\mathbb{R}^{n}$ , the personalized PageRank vector $\mathbf{x}$ satisfies $(D-\alpha W)\mathbf{x}=(1-\alpha)D\mathbf{s}$ . The coefficient matrix is SDDM because $D-\alpha W=\alpha(D-W)+(1-\alpha)D=\alpha L+(1-\alpha)D$ is symmetric, has nonpositive off-diagonal entries, and is diagonally dominant Gleich (2015). If every vertex has positive weighted degree, the matrix is also positive definite. In repeated-query settings, the operator $D-\alpha W$ remains fixed while the seed distribution $\mathbf{s}$ varies across queries.

The factor-reuse setting.

Many applications solve a sequence of linear systems $A\,\mathbf{x}^{(r)}=\mathbf{b}^{(r)}$ , for $r=1,\ldots,q$ , in which the coefficient matrix $A$ remains fixed while the right-hand side changes: diffusion from different seed sets, label propagation across classes, effective-resistance queries for many source–sink pairs, and inner systems in iterative optimization. An approximate Cholesky method performs the elimination once and stores the result as a sparse approximate factorization. Recall that the exact Cholesky decomposition of a positive-definite matrix is $A=\mathcal{C}\mathcal{C}^{\top}$ with $\mathcal{C}$ lower triangular; approximate elimination instead produces a sparse $\widehat{\mathcal{C}}$ with $\widehat{\mathcal{C}}\widehat{\mathcal{C}}^{\top}\approx A$ , which we call the factor. The factor is stored implicitly as the elimination data recorded at each pivot and is applied through forward and backward substitutions as a preconditioner for each right-hand side. The factor is constructed at cost $T_{\mathrm{build}}$ and reused for all $q$ solves. Writing $\overline{T}_{\mathrm{solve}}$ for the average time of one preconditioned-CG solve using the factor, the total and amortized workload costs are

T_{q}^{\mathrm{total}}=T_{\mathrm{build}}+q\,\overline{T}_{\mathrm{solve}},\text{ and }T_{q}^{\mathrm{amort}}=\frac{T_{q}^{\mathrm{total}}}{q}=\frac{T_{\mathrm{build}}}{q}+\overline{T}_{\mathrm{solve}}.

Thus $q=1$ measures one-shot performance, while increasing $q$ amortizes construction and captures the benefit of factor reuse.

Exact elimination and Schur cliques

For the local elimination analysis, we work with a Laplacian representation; general SDDM systems can be reduced to this form using the standard Gremban expansion Gremban (1996). Approximate Cholesky methods eliminate vertices one at a time. Suppose a pivot $v$ has active neighbors $u_{1},\ldots,u_{d}$ with $a_{i}=w_{vu_{i}}$ and $a=\sum_{i}a_{i}$ . Eliminating $v$ exactly removes the star at $v$ and adds the Schur clique

K_{v}=\sum_{1\leq i<j\leq d}\frac{a_{i}a_{j}}{a}\,(\mathbf{e}_{i}-\mathbf{e}_{j})(\mathbf{e}_{i}-\mathbf{e}_{j})^{\top}

(1)

to the remaining graph, where $\mathbf{e}_{i}$ denotes the $i$ -th standard basis vector in the coordinate system indexed by $u_{1},\ldots,u_{d}$ ; edges already present among the neighbors simply accumulate the corresponding weight. The clique contains $\binom{d}{2}=\Theta(d^{2})$ edges, even when the original graph is sparse. To avoid “filling up” the matrix, approximate Cholesky methods replace $K_{v}$ by a sparse randomly-sampled update $\widehat{K}_{v}$ satisfying $\mathbb{E}[\widehat{K}_{v}]=K_{v}$ Kyng and Sachdeva (2016); Gao et al. (2026). The AC and AC2 solvers construct unbiased connected local updates by randomized edge pairing, with AC2 doubling the per-edge sampling budget for robustness Gao et al. (2026). These sampling distributions, which have been empirically studied, depend on the order in which the neighbor edges are processed. In contrast, the sampling distribution in CAST is determined by the weighted Schur clique and admits the local minimax and second-moment guarantees established in Section 4.

The local Schur-clique replacement problem.

We write $K^{+}$ for the Moore–Penrose pseudoinverse and $K^{+/2}=(K^{+})^{1/2}$ . Let $K=\sum_{e\in E_{K}}w_{e}\,\mathbf{b}_{e}\mathbf{b}_{e}^{\top}$ be a connected local graph Laplacian: a $d\times d$ matrix indexed by the $d$ active neighbors of a pivot rather than by all $n$ vertices, with edge set $E_{K}$ and weights $w_{e}>0$ . Its support graph is $([d],E_{K})$ , the graph formed by its positive-weight edges; since it is connected, $\ker(K)=\operatorname{span}\{\mathbf{1}\}$ and $\operatorname{rank}(K)=d-1$ . We measure a positive-semidefinite local contribution $H$ relative to $K$ by $\lVert K^{+/2}HK^{+/2}\rVert=\max_{\mathbf{x}\perp\mathbf{1}}(\mathbf{x}^{\top}H\mathbf{x})/(\mathbf{x}^{\top}K\mathbf{x})$ , the largest ratio of $H$ -energy to $K$ -energy, and call this normalization the energy geometry of $K$ . The leverage score of edge $e$ is $\tau_{e}=w_{e}\,\mathbf{b}_{e}^{\top}K^{+}\mathbf{b}_{e}$ , which is exactly the size of that edge’s own term $w_{e}\mathbf{b}_{e}\mathbf{b}_{e}^{\top}$ in this geometry Spielman and Srivastava (2008); Kyng and Sachdeva (2016); leverage scores satisfy $\sum_{e}\tau_{e}=\operatorname{rank}(K)=d-1$ .

For a distribution $\mathcal{D}$ over spanning trees of the support graph of $K$ , let $p_{e}=\mathbb{P}_{T\sim\mathcal{D}}(e\in T)$ , and assume $p_{e}>0$ for every $e\in E_{K}$ . The inverse-marginal estimator

\widehat{K}_{\mathcal{D}}(T)=\sum_{e\in T}\frac{w_{e}}{p_{e}}\,\mathbf{b}_{e}\mathbf{b}_{e}^{\top}

is unbiased, since $\mathbb{E}[\widehat{K}_{\mathcal{D}}(T)]=\sum_{e}p_{e}(w_{e}/p_{e})\,\mathbf{b}_{e}\mathbf{b}_{e}^{\top}=K$ . We measure the local risk of $\mathcal{D}$ by the largest normalized contribution of any sampled edge,

R_{K}(\mathcal{D})=\max_{e\in E_{K}}\left\lVert K^{+/2}\!\left(\frac{w_{e}}{p_{e}}\mathbf{b}_{e}\mathbf{b}_{e}^{\top}\right)\!K^{+/2}\right\rVert=\max_{e\in E_{K}}\frac{\tau_{e}}{p_{e}},

(2)

where the second equality holds because each normalized edge contribution is rank one. Question 1.1 is then the following general design problem: among such $\mathcal{D}$ , which minimizes $R_{K}(\mathcal{D})$ ? Since $R_{K}$ depends on $\mathcal{D}$ only through its marginals $(p_{e})_{e\in E_{K}}$ , this is equivalently a question about the feasible marginal vectors of spanning-tree distributions. CAST uses the weighted random spanning-tree distribution of $K$ , whose marginals satisfy $p_{e}=\tau_{e}$ by the transfer-current theorem Lyons (2003); this equalizes the normalized size of every sampled edge and, as we show in Theorem 4.5, minimizes $R_{K}$ .

3 CAST: Canonical Approximate Schur Trees

In this section we define the local CAST- $\rho$ update and the approximate-Cholesky preconditioner obtained by applying it at each pivot. The construction replaces the Schur clique of Eq. (1) with a sparse random update derived from a weighted spanning tree on a temporarily expanded neighborhood.

3.1 The CAST update

Consider an elimination step with pivot $v$ , active neighbors $u_{1},\ldots,u_{d}$ , and incident edge weights $a_{i}=w_{vu_{i}}>0$ , where $a=\sum_{i=1}^{d}a_{i}$ . We interpret these edge weights as electrical conductances, and refer to the pivot neighbors as terminals. Exact elimination of $v$ adds the Schur clique $K_{v}$ of Eq. (1), which contains $\Theta(d^{2})$ edges.

For an integer splitting factor $\rho\geq 1$ , CAST- $\rho$ temporarily replaces each terminal $u_{i}$ by $P_{i}=\{(i,1),\ldots,(i,\rho)\}$ , a block of $\rho$ auxiliary copies, each carrying conductance $c_{p}=a_{i}/\rho$ . Let $P=\bigcup_{i=1}^{d}P_{i}$ , and let $\phi(p)=i$ denote the terminal associated with $p\in P_{i}$ , so that $\sum_{p\in P_{i}}c_{p}=a_{i}$ and $\sum_{p\in P}c_{p}=a$ . On $P$ , CAST- $\rho$ defines the complete weighted expanded Schur clique with edge conductances $w^{\mathrm{exp}}_{pq}=c_{p}c_{q}/a$ for $p\neq q$ , and samples a weighted random spanning tree $T$ according to

\mathbb{P}(T)\propto\prod_{(p,q)\in T}w_{pq}^{\mathrm{exp}}.

(3)

Sampling from such a distribution is not straightforward for general edge weights; Section 3.2 shows that the product form of $w^{\mathrm{exp}}_{pq}$ admits exact sampling in $\mathcal{O}(\rho d)$ time. The expanded clique is therefore used only to define the distribution and is never materialized explicitly.

Each sampled auxiliary edge $(p,q)\in T$ is assigned conductance $c_{p}c_{q}/(c_{p}+c_{q})$ , which is precisely its inverse-marginal reweighting on the expanded clique (see Lemma 4.1). The auxiliary-copy blocks are then contracted back to their terminals. If $\phi(p)=i\neq j=\phi(q)$ , the edge $(p,q)$ contributes this conductance to the terminal edge $(u_{i},u_{j})$ ; if $\phi(p)=\phi(q)$ , it becomes a self-loop after contraction and contributes nothing to the terminal Laplacian. Writing $\widehat{K}_{v}$ for the random Laplacian on the pivot neighborhood obtained after contracting the blocks and discarding self-loops,

\widehat{K}_{v}=\sum_{\begin{subarray}{c}(p,q)\in T\\ \phi(p)\neq\phi(q)\end{subarray}}\frac{c_{p}c_{q}}{c_{p}+c_{q}}\bigl(\mathbf{e}_{\phi(p)}-\mathbf{e}_{\phi(q)}\bigr)\bigl(\mathbf{e}_{\phi(p)}-\mathbf{e}_{\phi(q)}\bigr)^{\top},

(4)

where $\mathbf{e}_{i}$ denotes the basis vector associated with terminal $u_{i}$ . The auxiliary tree contains exactly $\rho d-1$ edges, so after discarding within-block edges and aggregating parallel contracted edges, the update contains at most $\rho d-1$ nonzero terminal-edge contributions. For $\rho>1$ it need not itself be a tree on the terminals, although its support is always connected (Theorem 4.3).

CAST-1: the base $\rho=1$ update.

When $\rho=1$ , each block holds a single copy carrying its terminal’s full conductance, so the expanded clique is the Schur clique itself and contraction is the identity: CAST- $\rho$ reduces to sampling one weighted random spanning tree of $K_{v}$ . Because Schur complementation preserves effective resistances among the retained vertices, the effective resistance between terminals $u_{i}$ and $u_{j}$ in $K_{v}$ equals the series resistance of the path $u_{i}$ – $v$ – $u_{j}$ before elimination, $R_{ij}=1/a_{i}+1/a_{j}$ . The leverage score of edge $(i,j)$ is, $\tau_{ij}=\frac{a_{i}a_{j}}{a}\left(\frac{1}{a_{i}}+\frac{1}{a_{j}}\right)=\frac{a_{i}+a_{j}}{a}.$ By the transfer-current theorem the weighted random spanning tree selects edge $(i,j)$ with probability $\tau_{ij}$ , and inverse-marginal reweighting assigns it conductance $\frac{a_{i}a_{j}/a}{\tau_{ij}}=\frac{a_{i}a_{j}}{a_{i}+a_{j}}$ ,

exactly the series conductance of the two-hop path $u_{i}$ – $v$ – $u_{j}$ that the edge replaces. Thus CAST-1 is an order-independent one-tree replacement whose edge marginals equal the clique leverage scores; Theorem 4.5 shows that these marginals uniquely minimize the local risk $R_{K_{v}}$ .

3.2 Exact sampling in $\mathcal{O}(\rho d)$ time via Prüfer codes

Sampling the update requires a weighted random spanning tree of the expanded Schur clique, a complete graph on $m=\rho d$ auxiliary copies. Its conductances have the product form $w^{\mathrm{exp}}_{pq}=c_{p}c_{q}/a$ , and this structure permits exact sampling from (3) using an i.i.d. weighted Prüfer code, without materializing the dense clique. Background on Prüfer codes is provided in Appendix A.

Fix a spanning tree $T$ on the auxiliary-copy set $P$ . Since $T$ has $m-1$ edges,

\prod_{(p,q)\in T}w_{pq}^{\mathrm{exp}}=a^{-(m-1)}\prod_{p\in P}c_{p}^{\deg_{T}(p)},

so, up to a factor independent of $T$ , the target probability is determined entirely by the vertex degrees. For $m\geq 2$ the Prüfer correspondence is a bijection between labeled trees on $P$ and sequences in $P^{m-2}$ , under which $p$ appears exactly $\deg_{T}(p)-1$ times in the code of $T$ Aigner and Ziegler (1999); West (2001). Drawing the $m-2$ symbols independently with $\mathbb{P}\left(\text{symbol}=p\right)=c_{p}/a$ therefore generates $T$ with probability

\prod_{p\in P}\left(\frac{c_{p}}{a}\right)^{\deg_{T}(p)-1}=a^{-(m-2)}\prod_{p\in P}c_{p}^{\deg_{T}(p)-1}.

The ratio of this probability to the unnormalized target weight is $a\prod_{p\in P}c_{p}^{-1}$ , independent of $T$ . Hence the i.i.d. Prüfer construction samples exactly from the weighted random spanning-tree distribution of the expanded Schur clique.

The sampler needs neither the expanded clique nor an explicit array of auxiliary-copy weights. All copies in $P_{i}$ carry conductance $a_{i}/\rho$ , so each symbol can be drawn hierarchically: sample terminal $i$ with probability $a_{i}/a$ , then one of its $\rho$ copies uniformly. An alias table for the terminal distribution is built once in $\mathcal{O}(d)$ time, after which each of the $m-2$ symbols costs $\mathcal{O}(1)$ . Decoding the code into a tree takes $\mathcal{O}(m)$ time by the standard leaf-pointer algorithm. Contraction is a single pass over the $m-1$ auxiliary edges: cross-block edges are placed on the corresponding terminal edge with conductance $c_{p}c_{q}/(c_{p}+c_{q})$ , within-block edges become self-loops and are discarded, and parallel terminal edges are aggregated by summing conductances. Sampling, decoding, and contraction therefore cost $\mathcal{O}(d+m)=\mathcal{O}(\rho d)$ in total.

Proposition 3.1 (Sampling cost).

For every pivot with $d\geq 2$ active neighbors and every integer $\rho\geq 1$ , CAST-SchurTree $(a_{1},\ldots,a_{d},\rho)$ samples the exact weighted random spanning-tree distribution of the expanded Schur clique and contracts the sampled tree onto the original pivot neighborhood in $\mathcal{O}(\rho d)$ time, without materializing either the $\Theta(d^{2})$ -edge Schur clique or the expanded clique.

Proof.

See Appendix B.1. ∎

Summing Proposition 3.1 over pivots, the elimination takes $\mathcal{O}\!\left(\rho\sum_{v}d_{v}\right)$ time, where $d_{v}$ is the degree of $v$ in the residual graph when it is eliminated. The residual degrees depend on the fill realized by earlier updates, which is why we study the effect of $\rho$ on downstream fill empirically (Section 5).

Two structural features of the Schur clique make this possible. For $\rho=1$ , its star origin gives the leverage scores in closed form (Section 3.1), identifying the target distribution without any resistance computation. Under $\rho$ -way splitting, the expanded clique retains product-form conductances, so the same Prüfer argument applies verbatim. Splitting thus refines the update while preserving both the algebraic structure needed for exact sampling and the order-independence that makes the rule canonical.

The CAST preconditioner.

Applying the local CAST- $\rho$ update at every pivot yields an approximate Cholesky factorization, CAST-Chol (Algorithms 1–2). The method maintains a residual graph, the weighted graph on the not-yet-eliminated vertices, initialized to $A$ . At each step it selects a pivot $v$ , records the pivot star—the pivot, its active neighbors, and their incident conductances, which together form the column of $\widehat{\mathcal{C}}$ associated with $v$ —then deletes $v$ from the residual graph and inserts the contracted terminal-edge contributions returned by Algorithm CAST-SchurTree. In the reuse setting, the resulting factor is built once for the fixed operator $A$ and applied within PCG for the right-hand sides $\mathbf{b}^{(1)},\ldots,\mathbf{b}^{(q)}$ .

Algorithm 1 CAST-SchurTree

(a_{1},\ldots,a_{d},\rho)

1: if

d\leq 1

then

2: return an empty terminal-edge multiset

3: end if

a\leftarrow\sum_{i=1}^{d}a_{i}

m\leftarrow\rho d

5: Build an alias table for

\pi_{i}=a_{i}/a

i\in[d]

\triangleright

\mathcal{O}(d)

6: Initialize an empty Prüfer sequence

S

7: for

t=1,\ldots,m-2

8: Sample

i\sim\pi

and

k\sim\operatorname{Unif}([\rho])

9: Append the auxiliary copy

p=(i,k)

S

\triangleright

\phi(p)=i

10: end for

11: Decode

S

into its spanning tree

T

on the auxiliary-copy set

\triangleright

\mathcal{O}(m)

12: Initialize an empty multiset

\mathcal{E}

13: for each auxiliary-tree edge

(p,q)\in T

with

i=\phi(p)\neq\phi(q)=j

14: Insert

\bigl(i,\,j,\;a_{i}a_{j}/(\rho(a_{i}+a_{j}))\bigr)

into

\mathcal{E}

15: end for

16: return

\mathcal{E}

Algorithm 2 CAST-Chol

(A,\rho)

1: Initialize the residual graph to

A

and the factor

F

to empty

2: while active vertices remain do

3: Select the next pivot

v

by the elimination-ordering rule

4: Read its active neighbors

u_{1},\ldots,u_{d}

and incident conductances

a_{1},\ldots,a_{d}

5: Record the pivot star of

v

F

6: Delete

v

and its incident edges from the residual graph

\mathcal{E}_{v}\leftarrow\textsc{CAST-SchurTree}(a_{1},\ldots,a_{d},\rho)

8: for each

(i,j,h)\in\mathcal{E}_{v}

9: Add conductance

h

to the residual edge

(u_{i},u_{j})

\triangleright

parallel contributions accumulate

10: end for

11: end while

12: return

F

4 Local Theory of the CAST Update

We analyze the CAST- $\rho$ update at a single elimination step. For the base $\rho=1$ update, we prove that the weighted spanning-tree distribution uniquely minimizes the local risk defined in Eq. (2) among unbiased inverse-marginal one-tree estimators. For general integer $\rho\geq 1$ , we establish unbiasedness, connected support, and a normalized local second-moment bound of order $1/\rho$ , quantifying the reduction in local sampling variability obtained through finer splitting.

Fix a pivot $v$ with active neighbors $u_{1},\ldots,u_{d}$ , where $d\geq 2$ , and let $a_{i}=w_{vu_{i}}>0$ , $a=\sum_{i=1}^{d}a_{i}$ . Exact elimination of $v$ creates the Schur clique $K_{v}$ of (1). We identify the pivot neighborhood with $[d]$ , and all Loewner-order comparisons between local Laplacians on the pivot neighborhood are understood on $\mathbf{1}_{d}^{\perp}=\{\mathbf{x}\in\mathbb{R}^{d}:\mathbf{x}^{\top}\mathbf{1}_{d}=0\}$ .

Auxiliary-copy notation.

We use the notation of Section 3.1: the blocks $P_{i}=\{(i,1),\ldots,(i,\rho)\}$ of auxiliary copies with conductances $c_{p}=a_{i}/\rho$ , the set $P=\bigcup_{i=1}^{d}P_{i}$ , the map $\phi(p)=i$ , the expanded Schur clique with conductances $w^{\mathrm{exp}}_{pq}=c_{p}c_{q}/a$ , and the tree $T$ drawn from (3). For an auxiliary edge $e=(p,q)$ we write $H_{e}=H_{pq}$ for its contracted contribution to the terminal Laplacian: setting $i=\phi(p)$ and $j=\phi(q)$ ,

H_{pq}=\frac{c_{p}c_{q}}{c_{p}+c_{q}}(\mathbf{e}_{i}-\mathbf{e}_{j})(\mathbf{e}_{i}-\mathbf{e}_{j})^{\top}\quad\text{if }i\neq j,\qquad H_{pq}=0\quad\text{if }i=j,

where $\mathbf{e}_{i}$ is the basis vector of terminal $u_{i}$ , so that $\widehat{K}_{v}=\sum_{(p,q)\in T}H_{pq}$ as in (4). Throughout this section $p_{e}=\mathbb{P}\left(e\in T\right)$ denotes the inclusion probability of an edge $e$ of the expanded clique; when $\rho=1$ the expanded clique is $K_{v}$ itself, and $p_{e}$ agrees with the clique-edge marginals of Section 2.

Unbiasedness.

We first show that contracting the sampled auxiliary tree preserves the exact Schur update in expectation.

Lemma 4.1 (Auxiliary-edge marginals).

For every auxiliary edge $(p,q)$ of the expanded Schur clique, $\mathbb{P}\left((p,q)\in T\right)=\frac{c_{p}+c_{q}}{a}$ .

Proof.

See Appendix B.2. ∎

In particular, the conductance assigned to a selected auxiliary edge is exactly its inverse-marginal reweighting:

\frac{w^{\mathrm{exp}}_{pq}}{\mathbb{P}\left((p,q)\in T\right)}=\frac{c_{p}c_{q}/a}{(c_{p}+c_{q})/a}=\frac{c_{p}c_{q}}{c_{p}+c_{q}}.

Theorem 4.2 (Unbiasedness).

For every pivot $v$ and every integer $\rho\geq 1$ , $\mathbb{E}[\widehat{K}_{v}]=K_{v}$ .

Proof.

See Appendix B.3. ∎

Connected support and correct local nullspace.

Unbiasedness controls the local update in expectation. CAST- $\rho$ also has a deterministic structural property: every realization has connected support on the pivot neighborhood.

Theorem 4.3 (Connected local support).

For every integer $\rho\geq 1$ , the terminal support graph of $\widehat{K}_{v}$ is connected on $\{u_{1},\ldots,u_{d}\}$ .

Proof.

See Appendix B.4. ∎

Corollary 4.4 (Correct local nullspace and rank).

For every integer $\rho\geq 1$ , $\ker(\widehat{K}_{v})=\operatorname{span}\{\mathbf{1}_{d}\}$ and $\operatorname{rank}(\widehat{K}_{v})=d-1$ .

Proof.

See Appendix B.5. ∎

Connectivity alone does not distinguish CAST from the practical AC and AC2 updates, which are also connected; the distinction lies in the distribution from which the local update is drawn. We turn to that distribution.

4.1 Local minimax optimality of the base $\rho=1$ update

When $\rho=1$ , CAST samples a weighted random spanning tree directly from the Schur clique, so the update is an inverse-marginal one-tree estimator (Section 2). We now answer Question 1.1: which edge marginals minimize the local risk $R_{K}(\mathcal{D})=\max_{e\in E_{K}}\tau_{e}/p_{e}$ , the largest normalized contribution of any sampled edge? We answer it for an arbitrary connected local Laplacian $K$ ; the Schur clique $K_{v}$ is the case relevant to CAST.

Theorem 4.5 (Local minimax optimality of leverage-score marginals).

For every spanning-tree distribution $\mathcal{D}$ satisfying $p_{e}>0$ for all $e\in E_{K}$ , $R_{K}(\mathcal{D})\geq 1$ , with equality if and only if $p_{e}=\tau_{e}$ for every $e\in E_{K}$ . In particular, the weighted random spanning-tree distribution of $K$ , whose edge marginals satisfy $p_{e}=\tau_{e}$ , is locally minimax-optimal among unbiased inverse-marginal one-tree estimators.

Proof.

See Appendix B.6. ∎

The proof is a short averaging argument: every spanning tree has $d-1$ edges, so $\sum_{e}p_{e}=d-1=\sum_{e}\tau_{e}$ , and a maximum is at least a weighted average.

Theorem 4.5 identifies the optimal edge-marginal vector uniquely, namely $p_{e}=\tau_{e}$ ; it does not assert uniqueness of the full distribution over spanning trees. CAST-1 realizes these marginals with the weighted random spanning-tree distribution of the clique, which for the Schur clique $K_{v}$ has the closed form derived in Section 3.1: edge $(i,j)$ is included with probability $\tau_{ij}=(a_{i}+a_{j})/a$ and reweighted to conductance $a_{i}a_{j}/(a_{i}+a_{j})$ . Hence CAST-1 attains $R_{K_{v}}=1$ , the smallest value achievable in this class.

4.2 Local second-moment certificate for CAST- $\rho$

Theorem 4.5 concerns estimators supported on a single spanning tree of the terminal clique, and so applies directly to the base $\rho=1$ update. For $\rho>1$ , CAST- $\rho$ samples one spanning tree on the expanded Schur clique and contracts the auxiliary-copy blocks back to the terminals; the resulting terminal update need not itself be a tree, and is therefore not a competing tree distribution on $K_{v}$ .

We analyze the general $\rho$ construction through a different local certificate: the normalized second moment of the one-pivot Schur error. The bound below scales as $1/\rho$ , quantifying how splitting decomposes the local estimator into smaller normalized contributions.

Let $E^{\mathrm{exp}}$ denote the edge set of the expanded Schur clique. For each auxiliary edge $e=(p,q)\in E^{\mathrm{exp}}$ , let $H_{e}$ be its contracted terminal contribution, with $H_{e}=0$ when $\phi(p)=\phi(q)$ , and set $A_{e}=K_{v}^{+/2}H_{e}K_{v}^{+/2}$ , its size in the energy geometry of $K_{v}$ . Let $\Pi_{K_{v}}=K_{v}^{+/2}K_{v}K_{v}^{+/2}$ be the orthogonal projection onto $\operatorname{range}(K_{v})=\mathbf{1}_{d}^{\perp}$ .

Lemma 4.6 (Normalized auxiliary-edge contribution).

For every auxiliary edge $e\in E^{\mathrm{exp}}$ , $0\preceq A_{e}\preceq\frac{1}{\rho}\Pi_{K_{v}}$ . Moreover, $A_{e}^{2}=\frac{1}{\rho}A_{e}$ .

Proof.

See Appendix B.7. ∎

Each auxiliary edge therefore contributes at most $1/\rho$ in the energy geometry of $K_{v}$ , independently of the pivot degree and of the incident weights: splitting shrinks the largest possible single contribution in direct proportion to $\rho$ . The second-moment analysis combines this with the determinantal negative dependence of weighted random spanning-tree edge indicators.

Lemma 4.7 (Covariance domination for spanning-tree indicators).

Let $T$ be a weighted random spanning tree of a connected weighted graph, and for each edge $e$ let $X_{e}=\mathbf{1}_{\{e\in T\}}$ and $p_{e}=\mathbb{E}[X_{e}]$ . If $C$ is the covariance matrix of $X$ , then $C\preceq\operatorname{diag}(p)$ .

Proof.

See Appendix B.8. ∎

Negative dependence is what makes the tree structure work in our favor: the sampled edges are not independent, but their covariance never exceeds what independent sampling with the same marginals would give, so the second moment can be bounded edgewise.

Theorem 4.8 (Local second-moment bound for CAST- $\rho$ ).

Let $Y_{v}=K_{v}^{+/2}(\widehat{K}_{v}-K_{v})K_{v}^{+/2}$ be the normalized local Schur error of the CAST- $\rho$ update. For every integer $\rho\geq 1$ ,

\mathbb{E}[Y_{v}]=0,\qquad\mathbb{E}[Y_{v}^{2}]\preceq\frac{1}{\rho}\,\Pi_{K_{v}}.

Proof.

See Appendix B.9. ∎

At $\rho=1$ the bound reads $\mathbb{E}[Y_{v}^{2}]\preceq\Pi_{K_{v}}$ , matching the unit local risk attained by CAST-1 (Theorem 4.5); each further doubling of $\rho$ halves the certified bound.

Scope of the guarantees.

The results above certify the CAST primitive at a single elimination step: for every $\rho$ the update is unbiased with connected support, at $\rho=1$ its marginals uniquely minimize the local risk, and for general $\rho$ its normalized second moment obeys a $1/\rho$ bound. The second moment is the quantity that matrix-concentration analyses of approximate elimination control alongside the largest single increment, so the bound identifies $\rho$ as a principled control on sampling variability, uniform over pivot degrees and incident weights. It remains a local certificate: it does not track how these errors accumulate across eliminations, and it says nothing about the magnitude of the resulting gain or about the construction cost and downstream fill that splitting introduces. Section 5 measures these effects across four benchmark families.

5 Experimental Evaluation

Methods.

We compare four preconditioners within a common elimination and PCG framework. AC Gao et al. (2026) is the practical approximate Cholesky factorization based on the elimination estimator of Kyng and Sachdeva Kyng and Sachdeva (2016), implemented in Laplacians.jl¹¹ 1 https://github.com/danspielman/Laplacians.jl; at each pivot it replaces the exact Schur-complement clique by a connected update generated through randomized sequential edge pairing. AC2 is its doubled-budget variant (split $=$ merge $=2$ ): each edge is represented by up to two half-weight multiedges, with fill multiplicity capped at two. These are the state of the art on the families our corpus is drawn from, where multigrid and incomplete-Cholesky solvers Koutis et al. (2011b); Falgout and Yang (2002); Balay et al. (2019); Livne and Brandt (2012); Meijerink and van der Vorst (1977) each fail on some instance while AC and AC2 converge throughout Gao et al. (2026), as our own runs on these collections confirmed; we therefore confine the comparison to the randomized approximate-elimination family.

CAST-1 applies the base $\rho=1$ update, replacing each elimination clique by a weighted random spanning tree whose edge marginals equal the clique leverage scores (Section 3.1). CAST-2 applies two-way splitting: each pivot neighbor is temporarily replaced by two equal-conductance auxiliary copies, one weighted random spanning tree is sampled from the resulting expanded Schur clique, and the copies are contracted back to the original neighborhood. We compare AC with CAST-1 as the base-granularity methods and AC2 with CAST-2 as their doubled-granularity variants; all four share the same elimination ordering (greedy on approximate minimum unweighted degree, following Gao et al. (2026)), factor representation, and PCG implementation, so the comparison isolates the local update.

AC2 versus CAST-2.

The two robustness variants increase sampling granularity in different ways. In AC2, edge replication persists throughout the factorization: the residual operator is represented as a multigraph with at most two multiedges per vertex pair, and a current neighbor represented by $t\in\{1,2\}$ multiedges contributes $t$ samples to the sequential clique update. Thus AC2 is not equivalent to averaging two independent AC factorizations.

In CAST-2, splitting is local to the current pivot. Eliminating a degree- $d$ pivot creates $2d$ temporary auxiliary copies, one weighted random spanning tree is sampled on the expanded Schur clique, and its copy blocks are contracted onto the $d$ original neighbors. The auxiliary copies are then discarded, although the contracted terminal edges remain in the residual graph and may increase downstream fill. Hence AC2 and CAST-2 are comparable doubled-granularity variants, but they differ in sample dependence, persistence, and the distribution of fill.

Metrics.

For $q=250$ right-hand sides, we report the total reuse workload

T_{250}=T_{\mathrm{build}}+\sum_{r=1}^{250}T_{\mathrm{solve},r}.

The two primary paired speedups are

S_{1}=\frac{T_{250}(\textit{AC})}{T_{250}(\textit{CAST-1})},\qquad S_{2}=\frac{T_{250}(\textit{AC2})}{T_{250}(\textit{CAST-2})},

where values greater than one favor the corresponding CAST variant. Aggregate speedups are arithmetic means of per-system ratios, computed in the direction stated for each comparison; we also report the number of systems on which each CAST variant is faster, which is independent of this convention. Following Gao et al. (2026), we also report factor-construction and per-right-hand-side solve costs normalized by the number of input nonzeros,

C_{\mathrm{build}}=10^{6}\,\frac{T_{\mathrm{build}}}{\operatorname{nnz}(A)},\qquad C_{\mathrm{solve}}=10^{6}\,\frac{\overline{T}_{\mathrm{solve}}}{\operatorname{nnz}(A)},

with $\overline{T}_{\mathrm{solve}}$ the average per-solve time of Section 2. Both normalized costs are reported in $\mu\mathrm{s}/\mathrm{nnz}$ .

Setup.

Unless otherwise stated, every factorization is evaluated only as a preconditioner for the same PCG implementation, which uses a recurrence-based relative-residual stopping threshold of $10^{-8}$ and an iteration cap of $10^{3}$ . After each solve, we explicitly compute $\|A\mathbf{x}-\mathbf{b}\|_{2}/\|\mathbf{b}\|_{2}$ to verify the requested tolerance. We report any discrepancy between the recurrence-based stopping test and this explicit residual check; the treatment of failed solves in each aggregate is stated alongside the corresponding result.

For each matrix and method, one factor is constructed and reused for $q=250$ right-hand sides. The right-hand sides are independent standard Gaussian vectors projected onto $\mathbf{1}^{\perp}$ , as required for Laplacian compatibility. A fixed right-hand-side sequence, generated independently of the factor-sampling seeds, is presented to every method and every factor draw, so comparisons are paired at the level of individual right-hand sides. The number of independently seeded factor draws and the aggregation rule are stated for each collection in the corresponding appendix; all methods entering a paired comparison use the same right-hand-side sequence and the same number of factor draws.

Matrix collections.

The evaluation corpus contains $201$ SDDM and Laplacian systems from the SDDM2023 benchmark suite²² 2 https://rjkyng.github.io/SDDM2023/ of Gao, Kyng, and Spielman: $28$ SuiteSparse matrices, $128$ Chimera-IPM systems, $40$ Spielman-IPM systems, and five Sachdeva-star instances. Eleven SuiteSparse matrices are excluded from the comparative aggregates as described below, leaving $190$ systems with paired measurements. Collection-specific protocols and per-system results appear in the appendix. All experiments run in a single process on an Apple Silicon MacBook Pro with $24$ GB of RAM under Julia 1.12. Preconditioner construction and application are single-threaded, and dense operations in the shared PCG implementation use the same OpenBLAS configuration for all methods.

The reported Spielman corpus contains the four scales $k\in\{100,200,300,400\}$ . The larger $k\in\{500,600\}$ sequences, containing approximately $1.5\times 10^{8}$ and $2.2\times 10^{8}$ nonzeros, exceed the memory capacity of the benchmark machine under the multi-solve protocol and are not included.

Within the SuiteSparse collection, ten matrices are retained as correctness checks but excluded from comparative aggregates because all four methods solve them in one PCG iteration and their timings are dominated by fixed construction and application overhead. We also exclude bcsstm25: at the requested tolerance of $10^{-8}$ , all methods reach the iteration cap for at least some right-hand sides, whereas at tolerance $10^{-5}$ every method converges in one iteration, indicating that the failure is associated with the matrix–tolerance pair rather than a particular preconditioner. The SuiteSparse comparative aggregate therefore contains the remaining $17$ matrices.

5.1 Cross-collection results

Table 1 summarizes the paired comparative subsets at $q=250$ right-hand sides per factor. Across the $190$ systems with paired measurements, selecting the faster of CAST-1 and CAST-2 for each system in hindsight gives an arithmetic-mean speedup of $1.9\times$ over the correspondingly faster of AC and AC2, with CAST faster on $178$ systems. Per-system results are in Appendix C.

Base granularity.

CAST-1 improves on AC across every collection, and does so uniformly rather than on average: it is faster on all $128$ Chimera-IPM systems, with per-system ratios between $1.08\times$ and $1.31\times$ and a mean of $1.18\times$ , and on $14$ of $17$ SuiteSparse matrices with a mean of $1.17\times$ . Against a tuned implementation that is the most reliable available solver on these families, a consistent $1.17$ – $1.18\times$ speedup in total reuse workload—concentrated in the solve phase, and therefore compounding with the number of right-hand sides—is a substantive gain from changing only the local sampling distribution. On the Sachdeva-star ladder the margin instead grows with instance size, from $1.58\times$ at $k=50$ to $5.30\times$ at $k=600$ ( $3.09\times$ on average), as AC’s iteration count degrades faster than CAST-1’s.

Doubled granularity.

The $\rho=2$ comparison is bimodal, and the elimination-degree profile predicts which regime holds. When pivot degrees are uniformly small, splitting is nearly free and the gain is large: on the Spielman-IPM systems, where every pivot has degree at most five, CAST-2 reaches the verified $10^{-8}$ tolerance in one to two PCG iterations and is $4.43\times$ faster than AC2. The Sachdeva-star family benefits through a different mechanism. Its neighborhoods are clique-dominated rather than low-degree and CAST-2 needs multiple iterations for convergence, but it holds iteration growth to $25$ – $45$ across the ladder against $28$ – $65$ for AC2. When elimination produces a heavy tail of high-degree pivots the additional contributions instead dominate: on the Chimera collection AC2 is faster on $101$ of $128$ systems, CAST-2 is faster only on instance $i_{3}$ , whose maximum pivot degree of $48$ is the smallest in the family. On SuiteSparse both are close, CAST-2 leading by $1.07\times$ on $11$ of $17$ matrices.

Table 1: Cross-collection performance at

q=250

right-hand sides per factor. Entries are per-collection aggregate speedups; the aggregation rule is stated in the corresponding Appendix C. Parentheses give the number of systems on which the CAST variant is faster. The final column is a post hoc comparison that selects the faster configuration within each method family separately for each system:

\min(T_{\textit{AC}},T_{\textit{AC2}})/\min(T_{\textit{CAST-1}},T_{\textit{CAST-2}})

Collection	Elimination profile	$\frac{\textit{AC}}{\textit{CAST-1}}$	$\frac{\textit{AC2}}{\textit{CAST-2}}$	$\frac{\min(\textit{AC},\textit{AC2})}{\min(\textit{CAST-1},\textit{CAST-2})}$
SuiteSparse	mostly low-degree, $\bar{d}\lesssim 7$	$1.17\times$ (14/17)	$1.07\times$ (11/17)	$1.11\times$ (12/17)
Chimera-IPM	$\bar{d}=6$ – $27$ , $d_{\max}=48$ – $363$	$1.18\times$ (128/128)	$0.65\times$ (27/128)	$1.16\times$ (121/128)
Spielman-IPM	all pivots satisfy $d\leq 5$	$1.27\times$ (38/40)	$4.43\times$ (40/40)	$4.43\times$ (40/40)
Sachdeva-star	clique-dominated	$3.09\times$ (5/5)	$1.90\times$ (5/5)	$1.90\times$ (5/5)

Choosing $\rho$ .

Two conclusions follow. First, the sampling distribution matters on its own. All four methods share an elimination ordering, factor format, and PCG path; AC and CAST-1 differ only in how the local update is drawn. That difference alone cuts total reuse workload by $14$ – $15\%$ on the heterogeneous collections, almost entirely in the solve phase.

Second, $\rho$ is a structural choice, not a strictly better setting. Splitting pays when its extra terminal-edge contributions stay confined: under small pivot degrees, as on Spielman-IPM, or in clique-dominated neighborhoods where iterations would otherwise grow, as on Sachdeva-star. It costs when a heavy degree tail carries them into later eliminations, compounding as downstream fill. Degree alone does not decide this; the clique-dominated case shows why. What matters is whether the contributions propagate. A CAST-1 elimination profile predicts this on our corpus, but reading it requires a first factorization, and we do not evaluate it as an adaptive rule. The post hoc column of Table 1 therefore bounds what any per-system choice between $\rho=1$ and $\rho=2$ could achieve.

6 Related Work

Laplacian and SDDM solvers.

Nearly-linear-time Laplacian solvers originated from support-graph preconditioning, low-stretch graph constructions, and spectral sparsification Spielman and Teng (2004); Koutis et al. (2011a); Cohen et al. (2014); Jambulapati and Sidford (2025), with effective-resistance sampling an especially influential route to the last Spielman and Srivastava (2008). These works give strong global approximation and running-time guarantees, but their recursive preconditioning structures differ from the sparse approximate factorizations studied here. Practical solvers for these systems also include incomplete Cholesky and algebraic multigrid, among them graph-specialized variants such as LAMG Livne and Brandt (2012). CAST belongs instead to the randomized approximate-elimination line, which builds a sparse factor by eliminating vertices sequentially and replacing each dense Schur-complement clique by a sparse random update.

Randomized approximate elimination.

Kyng and Sachdeva introduced a nearly-linear-time approximate Gaussian elimination algorithm for graph Laplacians, replacing elimination fill by unbiased random samples and controlling the accumulated error through a matrix-martingale analysis Kyng and Sachdeva (2016). Related sparsified-Cholesky constructions extend the approach to connection Laplacians Kyng et al. (2016). On the practical side, RCHOL adds shared-memory parallelization Chen et al. (2021), and more recent work studies parallel frameworks and CPU/GPU construction of randomized approximate-Cholesky preconditioners Baumann and Kyng (2024); Liang et al. (2025). These address global guarantees, alternative sampling schemes, or parallel implementation. CAST instead studies the distribution used for a single local Schur-clique replacement, then applies that primitive within a sequential factorization.

Practical approximate-Cholesky solvers.

The closest predecessor is the AC( $k$ ) framework of Gao, Kyng, and Spielman Gao et al. (2026), which turns randomized approximate elimination into a practical SDDM solver in Laplacians.jl. Its local updates preserve connected support, and its configurations AC and AC2 use one and two samples per entry to trade construction cost against robustness. CAST adopts the same connected-update principle but differs in the distribution. AC( $k$ ) generates its trees by randomized sequential edge pairing, so the induced edge marginals depend on the order in which incident edges are processed. CAST samples from the weighted random spanning-tree distribution of the clique itself, which is order-independent and attains the minimax edge marginals of Theorem 4.5.

Leverage scores and random spanning trees.

Edge leverage scores—equivalently, edge conductance times effective resistance—are central to spectral sparsification Spielman and Srivastava (2008). For a weighted random spanning tree, the transfer-current theorem identifies the inclusion probability of each edge with its leverage score Lyons (2003), and random spanning trees have accordingly been studied as spectral sparsifiers, including from unions of a few independent trees Kyng and Song (2018). Sampling such trees on general graphs requires nontrivial graph-algorithmic machinery Durfee et al. (2017).

CAST uses this distribution locally rather than globally, which is what makes it cheap. Sampling one tree per elimination clique avoids the general machinery entirely: the star origin of the clique gives its leverage scores in closed form, and the product-form conductances of the expanded clique admit exact sampling in $\mathcal{O}(\rho d)$ time via weighted Prüfer codes Aigner and Ziegler (1999); West (2001). The union-of-trees results of Kyng and Song (2018) are also not the right analogy for $\rho>1$ : CAST- $\rho$ samples a single tree on an expanded vertex set, not $\rho$ independent trees on the original terminals.

Conclusion

We introduced CAST, a canonical sparse replacement for the Schur clique created during approximate Cholesky elimination. We show that among unbiased inverse-marginal one-tree estimators, leverage-score marginals uniquely minimize the largest normalized sampled-edge contribution. Because the clique is induced by a star, CAST samples the corresponding weighted random spanning tree exactly in $\mathcal{O}(\rho d)$ time. Empirically, CAST-1 reduces total reuse workload by $14$ – $15\%$ relative to AC on SuiteSparse and Chimera-IPM systems. When additional fill remains controlled, CAST-2 improves robustness: on Spielman-IPM systems it reaches the verified $10^{-8}$ residual tolerance in one to two PCG iterations and is, on average, $4.43\times$ faster than AC2.

References

[1] M. Aigner and G. M. Ziegler (1999) Proofs from the book. Vol. 1, Springer. Cited by: Appendix A, §3.2, §6.
[2] R. Andersen, F. Chung, and K. Lang (2006) Local graph partitioning using pagerank vectors. In 2006 47th annual IEEE symposium on foundations of computer science (FOCS’06), pp. 475–486. Cited by: §1.
[3] S. Balay, S. Abhyankar, M. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Eijkhout, W. Gropp, et al. (2019) PETSc users manual. Technical report Argonne National Laboratory. Cited by: §5.
[4] Y. Baumann, R. Kyng, and G. Zöcklein (2026) VAC: a volume-sampling-based elimination rule for approximate cholesky factorization. Note: Manuscript, July 28, 2026 Cited by: §1.
[5] Y. Baumann and R. Kyng (2024) A framework for parallelizing approximate gaussian elimination. In Proceedings of the 36th ACM Symposium on Parallelism in Algorithms and Architectures, pp. 195–206. Cited by: §6.
[6] C. Chen, T. Liang, and G. Biros (2021) RCHOL: randomized cholesky factorization for solving sdd linear systems. SIAM Journal on Scientific Computing 43 (6), pp. C411–C438. Cited by: §6.
[7] M. B. Cohen, R. Kyng, G. L. Miller, J. W. Pachocki, R. Peng, A. B. Rao, and S. C. Xu (2014) Solving sdd linear systems in nearly m log1/2 n time. In Proceedings of the forty-sixth annual ACM symposium on Theory of computing, pp. 343–352. Cited by: §6.
[8] D. Durfee, R. Kyng, J. Peebles, A. B. Rao, and S. Sachdeva (2017) Sampling random spanning trees faster than matrix multiplication. In Proceedings of the 49th Annual ACM SIGACT Symposium on Theory of Computing, pp. 730–742. Cited by: §6.
[9] R. D. Falgout and U. M. Yang (2002) Hypre: a library of high performance preconditioners. In International Conference on computational science, pp. 632–641. Cited by: §5.
[10] Y. Gao, R. Kyng, and D. A. Spielman (2026) AC(k): robust solution of laplacian equations by randomized approximate cholesky factorization. SIAM Journal on Scientific Computing. Cited by: §C.1, §C.2, §C.3, §C.4, §C.4, §1, §1, §1, §2, §5, §5, §5, §6, Abstract.
[11] D. F. Gleich (2015) PageRank beyond the web. siam REVIEW 57 (3), pp. 321–363. Cited by: §2.
[12] K. D. Gremban (1996) Combinatorial preconditioners for sparse, symmetric, diagonally dominant linear systems. Ph.D. Thesis, Carnegie Mellon University Pittsburgh. Cited by: §C.1, §2.
[13] W. Huang, W. Lin, H. H. Lu, and S. Yau (2019) ISIRA: integrated shift–invert residual arnoldi method for graph laplacian matrices from big data. Journal of Computational and Applied Mathematics 346, pp. 518–531. Cited by: §1.
[14] A. Jambulapati and A. Sidford (2025) Ultrasparse ultrasparsifiers and faster laplacian system solvers. ACM Transactions on Algorithms 21 (3), pp. 1–49. Cited by: §6.
[15] I. Koutis, G. L. Miller, and R. Peng (2011) A nearly-m log n time solver for sdd linear systems. In 2011 IEEE 52nd Annual Symposium on Foundations of Computer Science, pp. 590–598. Cited by: §6.
[16] I. Koutis, G. L. Miller, and D. Tolliver (2011) Combinatorial preconditioners and multilevel solvers for problems in computer vision and image processing. Computer Vision and Image Understanding 115 (12), pp. 1638–1646. Cited by: §5.
[17] R. Kyng, Y. T. Lee, R. Peng, S. Sachdeva, and D. A. Spielman (2016) Sparsified cholesky and multigrid solvers for connection laplacians. In Proceedings of the forty-eighth annual ACM symposium on Theory of Computing, pp. 842–850. Cited by: §6.
[18] R. Kyng and S. Sachdeva (2016) Approximate gaussian elimination for laplacians-fast, sparse, and simple. In 2016 IEEE 57th Annual Symposium on Foundations of Computer Science (FOCS), pp. 573–582. Cited by: §1, §1, §1, §2, §2, §5, §6.
[19] R. Kyng and Z. Song (2018) A matrix chernoff bound for strongly rayleigh distributions and spectral sparsifiers from a few random spanning trees. In 2018 IEEE 59th Annual Symposium on Foundations of Computer Science (FOCS), pp. 373–384. Cited by: §6, §6.
[20] T. Liang, C. Chen, Y. Yaniv, H. Luo, D. Tench, X. S. Li, A. Buluc, and J. Demmel (2025) Parallel gpu-accelerated randomized construction of approximate cholesky preconditioners. arXiv preprint arXiv:2505.02977. Cited by: §6.
[21] O. E. Livne and A. Brandt (2012) Lean algebraic multigrid (lamg): fast graph laplacian linear solver. SIAM Journal on Scientific Computing 34 (4), pp. B499–B522. Cited by: §5, §6.
[22] R. Lyons (2003) Determinantal probability measures. Publications Mathématiques de l’IHÉS 98, pp. 167–212. Cited by: §B.2, §B.6, §B.8, §1, §2, §6.
[23] J. A. Meijerink and H. A. van der Vorst (1977) An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix. Mathematics of Computation 31 (137), pp. 148–162. Cited by: §5.
[24] L. Page, S. Brin, R. Motwani, and T. Winograd (1999) The pagerank citation ranking: bring order to the web. In Proc. of the 7th International World Wide Web Conf.–1998, Cited by: §1.
[25] D. A. Spielman and N. Srivastava (2008) Graph sparsification by effective resistances. In Proceedings of the fortieth annual ACM symposium on Theory of computing, pp. 563–568. Cited by: §1, §1, §2, §6, §6.
[26] D. A. Spielman and S. Teng (2004) Nearly-linear time algorithms for graph partitioning, graph sparsification, and solving linear systems. In Proceedings of the thirty-sixth annual ACM symposium on Theory of computing, pp. 81–90. Cited by: §1, §6.
[27] D. B. West (2001) Introduction to graph theory. 2nd edition, Prentice Hall. Cited by: Appendix A, §3.2, §6.
[28] D. Zhou, O. Bousquet, T. Lal, J. Weston, and B. Schölkopf (2003) Learning with local and global consistency. Advances in neural information processing systems 16. Cited by: §1.
[29] X. Zhu, Z. Ghahramani, and J. D. Lafferty (2003) Semi-supervised learning using gaussian fields and harmonic functions. In Proceedings of the 20th International conference on Machine learning (ICML-03), pp. 912–919. Cited by: §1.

Appendix A Background on Prüfer Codes

This appendix collects the facts about Prüfer codes underlying the exact sampler of Section 3.2.

The Prüfer correspondence.

Let $P$ be a set of $m\geq 2$ labeled vertices with an arbitrary total ordering. A Prüfer code is a sequence $S=(S_{1},\ldots,S_{m-2})\in P^{m-2}$ . The Prüfer correspondence is a bijection between such sequences and labeled spanning trees on $P$ [1, 27]: every labeled tree has a unique code, and every sequence of length $m-2$ decodes to a unique tree. To encode a tree, repeatedly remove the smallest-labeled leaf and append its unique neighbor to the sequence, stopping when two vertices remain; one vertex is removed per step, so the code has exactly $m-2$ symbols.

For example, the tree on $\{1,2,3,4,5\}$ with edges $\{(1,3),(2,3)\\ ,(3,4),(4,5)\}$ has smallest leaf $1$ , whose neighbor is $3$ , so the first symbol is $3$ . After removing $1$ the smallest leaf is $2$ , again with neighbor $3$ ; after removing $2$ it is $3$ , whose remaining neighbor is $4$ . The code is $(3,3,4)$ .

The inverse procedure reconstructs the tree. For each label $p$ initialize $r_{p}=1+\lvert\{t:S_{t}=p\}\rvert$ . At each step choose the smallest label $\ell$ with $r_{\ell}=1$ , connect $\ell$ to the next code symbol $S_{t}$ , remove $\ell$ , and decrease $r_{S_{t}}$ by one. Two labels remain once all symbols are processed; connecting them completes the tree. When $m=2$ the code is empty and the tree is the single edge between the two labels.

Degrees from code multiplicities.

The property CAST relies on is that $p$ appears in the code exactly $\deg_{T}(p)-1$ times, i.e.

\deg_{T}(p)=1+\bigl\lvert\{t:S_{t}=p\}\bigr\rvert.

Intuitively, $p$ is recorded whenever a neighboring leaf is removed while $p$ remains, which happens once per incident edge except the edge through which $p$ is itself removed, or the final edge if $p$ is one of the last two vertices. Consistently, $\sum_{p\in P}(\deg_{T}(p)-1)=2(m-1)-m=m-2$ , the length of the code.

Nonuniform Prüfer sampling.

The correspondence is purely combinatorial. A weighted tree distribution can nonetheless be induced by drawing the code symbols nonuniformly: if $S_{1},\ldots,S_{m-2}$ are independent with $\mathbb{P}\left(S_{t}=p\right)=\pi_{p}$ , then, because the correspondence is a bijection, the decoded tree $T$ has probability $\mathbb{P}\left(T\right)=\prod_{p\in P}\pi_{p}^{\deg_{T}(p)-1}$ , which depends on $T$ only through its degree sequence. This matches the structure of complete graphs whose edge weights factor over endpoints.

Lemma A.1 (Weighted Prüfer sampling).

Let the complete graph on $P$ have edge weights $w_{pq}=\beta\,\theta_{p}\theta_{q}$ for $p\neq q$ , with $\beta>0$ and $\theta_{p}>0$ . Draw $S_{1},\ldots,S_{m-2}$ independently with $\mathbb{P}\left(S_{t}=p\right)=\theta_{p}/\sum_{r\in P}\theta_{r}$ and let $T$ be the decoded tree. Then $T$ follows the weighted random spanning-tree distribution, $\mathbb{P}\left(T\right)\propto\prod_{(p,q)\in T}w_{pq}$ .

Proof.

A spanning tree has $m-1$ edges and each $p$ appears in $\deg_{T}(p)$ of them, so $\prod_{(p,q)\in T}w_{pq}=\beta^{m-1}\prod_{p\in P}\theta_{p}^{\deg_{T}(p)}$ . By the degree-multiplicity property,

\mathbb{P}\left(T\right)=\prod_{p\in P}\left(\frac{\theta_{p}}{\sum_{r}\theta_{r}}\right)^{\deg_{T}(p)-1}=\Bigl(\sum_{r}\theta_{r}\Bigr)^{-(m-2)}\prod_{p\in P}\theta_{p}^{\deg_{T}(p)-1}.

The ratio of the two expressions is $\beta^{m-1}\bigl(\sum_{r}\theta_{r}\bigr)^{m-2}\prod_{p}\theta_{p}$ , which does not depend on $T$ . Both are distributions over the same finite set, so they coincide. ∎

The endpoint-product form is essential: independent Prüfer symbols do not sample the correct weighted spanning-tree distribution for general edge weights.

Application to CAST- $\rho$ .

At a pivot with $d$ terminals, the expanded Schur clique is the complete graph on the $m=\rho d$ auxiliary copies with conductances $w^{\mathrm{exp}}_{pq}=c_{p}c_{q}/a$ . This is Lemma A.1 with $\theta_{p}=c_{p}$ and $\beta=1/a$ , and since $\sum_{p\in P}c_{p}=a$ , the symbol distribution is $\mathbb{P}\left(S_{t}=p\right)=c_{p}/a$ . Because every copy of terminal $i$ carries the same conductance $a_{i}/\rho$ , a symbol can be drawn hierarchically: sample $i$ with probability $a_{i}/a$ , then one of its $\rho$ copies uniformly. Section 3.2 gives the resulting sampler and its $\mathcal{O}(\rho d)$ cost.

Appendix B Proofs

B.1 Proof of Proposition 3.1

Proof.

Let $m=\rho d$ . A spanning tree $T$ of the expanded Schur clique has $m-1$ edges, and each copy $p$ appears in $\deg_{T}(p)$ of them, so

\prod_{(p,q)\in T}w_{pq}^{\mathrm{exp}}=a^{-(m-1)}\prod_{p\in P}c_{p}^{\deg_{T}(p)}.

The Prüfer correspondence is a bijection between spanning trees of the complete graph on $P$ and sequences in $P^{m-2}$ , under which $p$ appears exactly $\deg_{T}(p)-1$ times in the code of $T$ . Hence drawing the $m-2$ symbols independently with $\mathbb{P}\left(\text{symbol}=p\right)=c_{p}/a$ generates $T$ with probability

\prod_{p\in P}\left(\frac{c_{p}}{a}\right)^{\deg_{T}(p)-1}=a^{-(m-2)}\prod_{p\in P}c_{p}^{\deg_{T}(p)-1},

where we used $\sum_{p\in P}\bigl(\deg_{T}(p)-1\bigr)=m-2$ . The ratio of this probability to the unnormalized tree weight above is $a\prod_{p\in P}c_{p}^{-1}$ , independent of $T$ . Two probability distributions on the same finite set whose ratio is constant are equal, so the sampler is exact.

For the running time, all copies of terminal $i$ carry the same conductance $a_{i}/\rho$ , so a symbol is drawn by sampling a terminal from an alias table for $(a_{1}/a,\ldots,a_{d}/a)$ and then one of its $\rho$ copies uniformly. The table is built once in $\mathcal{O}(d)$ time, after which each of the $m-2$ symbols costs $\mathcal{O}(1)$ . Decoding takes $\mathcal{O}(m)$ time by the standard leaf-pointer algorithm. Contraction is one pass over the $m-1$ auxiliary edges with $\mathcal{O}(1)$ work each: a cross-block edge $(p,q)$ is emitted as the terminal-edge contribution $(\phi(p),\phi(q),c_{p}c_{q}/(c_{p}+c_{q}))$ , and a within-block edge is discarded. The total is $\mathcal{O}(d+m)=\mathcal{O}(\rho d)$ .

Finally, the only structures allocated are the alias table over the $d$ terminals, the code of length $m-2$ , and the decoded tree on $m$ copies. Neither the $\Theta(d^{2})$ -edge Schur clique nor the expanded clique is ever formed. ∎

B.2 Proof of Lemma 4.1

Proof.

Let $\mathrm{S}$ be the star on $P\cup\{z\}$ in which the center $z$ is joined to each auxiliary copy $r\in P$ by an edge of conductance $c_{r}$ , so that the total conductance at $z$ is $\sum_{r\in P}c_{r}=a$ . Eliminating $z$ from $\mathrm{S}$ produces, by Eq. (1), the clique on $P$ with edge conductances $c_{p}c_{q}/a=w^{\mathrm{exp}}_{pq}$ . The expanded Schur clique is therefore the Schur complement of $\mathrm{S}$ onto $P$ .

Schur complementation preserves effective resistances among the retained vertices, so the effective resistance between $p$ and $q$ in the expanded clique equals their effective resistance in $\mathrm{S}$ . Since $\mathrm{S}$ is a tree, the unique $p$ – $q$ path $p$ – $z$ – $q$ determines this resistance, giving

R_{pq}^{\mathrm{exp}}=\frac{1}{c_{p}}+\frac{1}{c_{q}}.

By the transfer-current theorem [22], the inclusion probability of an edge in a weighted random spanning tree is its conductance times its effective resistance, whence

\mathbb{P}\left((p,q)\in T\right)=w_{pq}^{\mathrm{exp}}R_{pq}^{\mathrm{exp}}=\frac{c_{p}c_{q}}{a}\left(\frac{1}{c_{p}}+\frac{1}{c_{q}}\right)=\frac{c_{p}+c_{q}}{a}.\qed

B.3 Proof of Theorem 4.2

Proof.

Write $\widehat{K}_{v}=\sum_{e\in E^{\mathrm{exp}}}X_{e}H_{e}$ with $X_{e}=\mathbf{1}_{\{e\in T\}}$ , so that by linearity $\mathbb{E}[\widehat{K}_{v}]=\sum_{e}\mathbb{P}\left(e\in T\right)\,H_{e}$ . Within-block edges contract to self-loops and contribute $H_{e}=0$ , so only pairs with endpoints in distinct blocks remain.

Fix distinct terminals $i\neq j$ . Collecting the terms carrying $(\mathbf{e}_{i}-\mathbf{e}_{j})(\mathbf{e}_{i}-\mathbf{e}_{j})^{\top}$ and applying Lemma 4.1,

	$\displaystyle\sum_{p\in P_{i}}\sum_{q\in P_{j}}\mathbb{P}\left((p,q)\in T\right)\,\frac{c_{p}c_{q}}{c_{p}+c_{q}}$	$\displaystyle=\sum_{p\in P_{i}}\sum_{q\in P_{j}}\frac{c_{p}+c_{q}}{a}\cdot\frac{c_{p}c_{q}}{c_{p}+c_{q}}$
		$\displaystyle=\frac{1}{a}\Bigl(\sum_{p\in P_{i}}c_{p}\Bigr)\Bigl(\sum_{q\in P_{j}}c_{q}\Bigr)=\frac{a_{i}a_{j}}{a},$

using $\sum_{p\in P_{i}}c_{p}=a_{i}$ . This is exactly the coefficient of the terminal edge $(i,j)$ in $K_{v}$ . Both $\mathbb{E}[\widehat{K}_{v}]$ and $K_{v}$ are graph Laplacians on $\{u_{1},\ldots,u_{d}\}$ , and a graph Laplacian is determined by its edge coefficients, so

\mathbb{E}[\widehat{K}_{v}]=K_{v}.

∎

B.4 Proof of Theorem 4.3

Proof.

The sampled tree $T$ is connected on the auxiliary-copy set $P$ . Contract each block $P_{i}$ to its terminal $u_{i}$ . Contraction preserves connectedness: given terminals $u_{i}$ and $u_{j}$ , pick $p\in P_{i}$ and $q\in P_{j}$ ; the $p$ – $q$ path in $T$ maps to a walk from $u_{i}$ to $u_{j}$ in the quotient graph.

It remains to check that the quotient edges are exactly the edges of the support graph of $\widehat{K}_{v}$ . An auxiliary edge with both endpoints in one block becomes a self-loop and is discarded, which does not affect connectivity. A cross-block edge $(p,q)$ with $\phi(p)=i\neq j=\phi(q)$ contributes conductance $c_{p}c_{q}/(c_{p}+c_{q})>0$ to the terminal edge $(u_{i},u_{j})$ , and parallel contributions are summed. Since all $c_{p}>0$ , every such contribution is strictly positive and no cancellation occurs, so each quotient edge carries positive weight in $\widehat{K}_{v}$ . The support graph of $\widehat{K}_{v}$ is therefore the quotient graph with self-loops removed, and is connected.

When $\rho=1$ each block is a single copy, so the contraction is the identity and $\widehat{K}_{v}$ is itself a spanning tree of the terminal neighborhood; for $\rho>1$ it need not be a tree. ∎

B.5 Proof of Corollary 4.4

Proof.

The nullspace of a weighted graph Laplacian has dimension equal to the number of connected components of its support graph, and is spanned by the indicator vectors of those components. By Theorem 4.3 the support graph of $\widehat{K}_{v}$ is connected on $\{u_{1},\ldots,u_{d}\}$ , so $\ker(\widehat{K}_{v})=\operatorname{span}\{\mathbf{1}_{d}\}$ and, by rank–nullity, $\operatorname{rank}(\widehat{K}_{v})=d-1$ . ∎

B.6 Proof of Theorem 4.5

Proof.

Every spanning tree of the support graph of $K$ has exactly $d-1$ edges, so $\sum_{e\in E_{K}}\mathbf{1}_{\{e\in T\}}=d-1$ for every realization; taking expectations gives $\sum_{e\in E_{K}}p_{e}=d-1$ . The leverage scores satisfy $\sum_{e\in E_{K}}\tau_{e}=\operatorname{rank}(K)=d-1$ as well, so

1=\frac{\sum_{e\in E_{K}}\tau_{e}}{\sum_{e\in E_{K}}p_{e}}=\frac{\sum_{e\in E_{K}}p_{e}\,(\tau_{e}/p_{e})}{\sum_{e\in E_{K}}p_{e}}\leq\max_{e\in E_{K}}\frac{\tau_{e}}{p_{e}}=R_{K}(\mathcal{D}),

since the middle expression is a weighted average of the ratios $\tau_{e}/p_{e}$ with positive weights $p_{e}$ . Hence $R_{K}(\mathcal{D})\geq 1$ .

A weighted average with positive weights attains its maximum only when every term equals that maximum. Equality therefore forces $\tau_{e}/p_{e}=1$ , i.e. $p_{e}=\tau_{e}$ , for every $e\in E_{K}$ ; conversely, if $p_{e}=\tau_{e}$ throughout then $R_{K}(\mathcal{D})=1$ . Finally, for the weighted random spanning-tree distribution of $K$ the transfer-current theorem [22] gives $p_{e}=\mathbb{P}\left(e\in T\right)=w_{e}\,\mathbf{b}_{e}^{\top}K^{+}\mathbf{b}_{e}=\tau_{e}$ , so this distribution attains the bound. ∎

B.7 Proof of Lemma 4.6

Proof.

If $e$ is a within-block edge then $H_{e}=0$ , hence $A_{e}=0$ and both claims hold trivially. Otherwise let $e=(p,q)$ with $\phi(p)=i\neq j=\phi(q)$ , and write $\mathbf{b}=\mathbf{e}_{i}-\mathbf{e}_{j}$ , so that $H_{e}=\frac{c_{p}c_{q}}{c_{p}+c_{q}}\,\mathbf{b}\mathbf{b}^{\top}$ and $A_{e}=\frac{c_{p}c_{q}}{c_{p}+c_{q}}\,(K_{v}^{+/2}\mathbf{b})(K_{v}^{+/2}\mathbf{b})^{\top}$ . Thus $A_{e}$ is positive semidefinite and rank one, and its only nonzero eigenvalue is

\|A_{e}\|=\frac{c_{p}c_{q}}{c_{p}+c_{q}}\,\bigl\|K_{v}^{+/2}\mathbf{b}\bigr\|^{2}=\frac{c_{p}c_{q}}{c_{p}+c_{q}}\,\mathbf{b}^{\top}K_{v}^{+}\mathbf{b}.

By the argument of Lemma 4.1 applied to the pivot star, the effective resistance between $u_{i}$ and $u_{j}$ in $K_{v}$ is $\mathbf{b}^{\top}K_{v}^{+}\mathbf{b}=1/a_{i}+1/a_{j}$ . Substituting $c_{p}=a_{i}/\rho$ and $c_{q}=a_{j}/\rho$ ,

\|A_{e}\|=\frac{a_{i}a_{j}}{\rho\,(a_{i}+a_{j})}\left(\frac{1}{a_{i}}+\frac{1}{a_{j}}\right)=\frac{1}{\rho},

independently of $d$ and of the incident weights.

Write $A_{e}=\frac{1}{\rho}\mathbf{v}\mathbf{v}^{\top}$ with $\mathbf{v}=K_{v}^{+/2}\mathbf{b}/\|K_{v}^{+/2}\mathbf{b}\|$ . Since $\mathbf{v}\in\operatorname{range}(K_{v}^{+/2})=\operatorname{range}(K_{v})$ and $\Pi_{K_{v}}$ is the orthogonal projection onto that subspace, $\mathbf{v}\mathbf{v}^{\top}\preceq\Pi_{K_{v}}$ , giving $0\preceq A_{e}\preceq\frac{1}{\rho}\Pi_{K_{v}}$ . Finally, a rank-one matrix $\lambda\mathbf{v}\mathbf{v}^{\top}$ with $\|\mathbf{v}\|=1$ satisfies $(\lambda\mathbf{v}\mathbf{v}^{\top})^{2}=\lambda^{2}\mathbf{v}\mathbf{v}^{\top}=\lambda\cdot\lambda\mathbf{v}\mathbf{v}^{\top}$ , so $A_{e}^{2}=\frac{1}{\rho}A_{e}$ . ∎

B.8 Proof of Lemma 4.7

Proof.

Fix an arbitrary orientation of the edges. The edge indicators of a weighted random spanning tree then form a determinantal point process whose kernel is the transfer-current matrix $M$ , which is symmetric and positive semidefinite, with $M_{ee}=p_{e}$ [22]. Determinantal pair correlations give, for $e\neq f$ ,

\mathbb{P}\left(e\in T,\,f\in T\right)=\det\begin{pmatrix}M_{ee}&M_{ef}\\ M_{fe}&M_{ff}\end{pmatrix}=p_{e}p_{f}-M_{ef}^{2},

so $C_{ef}=\mathbb{P}\left(e\in T,f\in T\right)-p_{e}p_{f}=-M_{ef}^{2}$ . On the diagonal, $X_{e}$ is a Bernoulli indicator, so $C_{ee}=\mathbb{E}[X_{e}^{2}]-p_{e}^{2}=p_{e}-p_{e}^{2}=p_{e}-M_{ee}^{2}$ . Hence

C=\operatorname{diag}(p)-M\circ M,

where $\circ$ is the entrywise product; note that $C$ does not depend on the chosen orientation, since only the squares $M_{ef}^{2}$ appear. Since $M\succeq 0$ , the Schur product theorem gives $M\circ M\succeq 0$ , and therefore $C\preceq\operatorname{diag}(p)$ . ∎

B.9 Proof of Theorem 4.8

Proof.

The mean-zero statement follows from Theorem 4.2:

\mathbb{E}[Y_{v}]=K_{v}^{+/2}\bigl(\mathbb{E}[\widehat{K}_{v}]-K_{v}\bigr)K_{v}^{+/2}=0.

For $e\in E^{\mathrm{exp}}$ let $X_{e}=\mathbf{1}_{\{e\in T\}}$ and $p_{e}=\mathbb{E}[X_{e}]$ . Since $\widehat{K}_{v}=\sum_{e}X_{e}H_{e}$ and, by unbiasedness, $K_{v}=\sum_{e}p_{e}H_{e}$ , conjugating by $K_{v}^{+/2}$ gives

Y_{v}=\sum_{e\in E^{\mathrm{exp}}}(X_{e}-p_{e})A_{e}.

Fix $\mathbf{z}\in\mathbb{R}^{d}$ . As $Y_{v}$ is symmetric, $\mathbf{z}^{\top}\mathbb{E}[Y_{v}^{2}]\mathbf{z}=\mathbb{E}\bigl[\lVert Y_{v}\mathbf{z}\rVert_{2}^{2}\bigr]$ , and expanding the square,

\mathbb{E}\bigl[\lVert Y_{v}\mathbf{z}\rVert_{2}^{2}\bigr]=\sum_{e,f}\mathbb{E}\bigl[(X_{e}-p_{e})(X_{f}-p_{f})\bigr]\langle A_{e}\mathbf{z},A_{f}\mathbf{z}\rangle=\sum_{e,f}C_{ef}\,G_{ef},

where $G_{ef}=\langle A_{e}\mathbf{z},A_{f}\mathbf{z}\rangle$ is the Gram matrix of the vectors $\{A_{e}\mathbf{z}\}$ and is therefore positive semidefinite. The right-hand side is the Frobenius inner product $\langle C,G\rangle_{F}$ .

By Lemma 4.7, $\operatorname{diag}(p)-C\succeq 0$ . For positive-semidefinite $A$ and $B$ one has $\langle A,B\rangle_{F}=\operatorname{tr}(AB)=\operatorname{tr}(A^{1/2}BA^{1/2})\geq 0$ , so $\langle\operatorname{diag}(p)-C,\,G\rangle_{F}\geq 0$ and hence

\mathbf{z}^{\top}\mathbb{E}[Y_{v}^{2}]\mathbf{z}\;\leq\;\langle\operatorname{diag}(p),G\rangle_{F}=\sum_{e}p_{e}\lVert A_{e}\mathbf{z}\rVert_{2}^{2}=\mathbf{z}^{\top}\Bigl(\sum_{e}p_{e}A_{e}^{2}\Bigr)\mathbf{z},

the last equality because each $A_{e}$ is symmetric. As $\mathbf{z}$ was arbitrary, $\mathbb{E}[Y_{v}^{2}]\preceq\sum_{e}p_{e}A_{e}^{2}$ .

Lemma 4.6 gives $A_{e}^{2}=\frac{1}{\rho}A_{e}$ , so $\sum_{e}p_{e}A_{e}^{2}=\frac{1}{\rho}\sum_{e}p_{e}A_{e}$ , and unbiasedness gives

\sum_{e\in E^{\mathrm{exp}}}p_{e}A_{e}=K_{v}^{+/2}\Bigl(\sum_{e}p_{e}H_{e}\Bigr)K_{v}^{+/2}=K_{v}^{+/2}K_{v}K_{v}^{+/2}=\Pi_{K_{v}}.

Therefore $\mathbb{E}[Y_{v}^{2}]\preceq\frac{1}{\rho}\Pi_{K_{v}}$ . ∎

Appendix C Empirical Evaluation

C.1 SuiteSparse benchmark

Benchmark.

We evaluate on 28 symmetric diagonally dominant M-matrices from the SuiteSparse Matrix Collection, following the benchmark selection of Gao, Kyng, and Spielman [10]. The collection includes grid and mesh Laplacians, finite-element and finite-volume discretizations, shallow-water models, ill-conditioned structural systems, and irregular graph problems, with the largest matrices containing up to $4.8$ million nonzeros.

Ten matrices are diagonal or effectively diagonal, and every method converges in one PCG iteration for every right-hand side. These instances do not exercise the randomized clique estimator; their timings primarily reflect fixed construction and application overhead. We therefore retain them as correctness checks but exclude them from aggregate performance comparisons.

Of the remaining 18 matrices, bcsstm25 is reported as did not finish (DNF) and excluded from the aggregates. Its conditioning places the requested tolerance of $10^{-8}$ below the attainable double-precision accuracy for some right-hand sides, causing every method to reach the iteration limit on those right-hand sides. At tolerance $10^{-5}$ , all methods converge in one iteration, indicating that the failure is attributable to the matrix–tolerance pair rather than to a particular solver. The aggregate comparison therefore contains $17$ matrices.

The SDDM inputs are reduced to Laplacian form using the standard Gremban expansion [12]. This preprocessing is applied once per matrix and shared identically by all methods.

Results.

Tables 2 and 3 report the complete per-matrix results for the base- and doubled-granularity comparisons, AC versus CAST-1 and AC2 versus CAST-2, respectively. Matrices are ordered by total workload time. Table 4 reports the ten trivial instances excluded from the comparative aggregates.

Base-granularity comparison.

CAST-1 is faster than AC on 14 of the 17 matrices, with an average speedup of $1.17\times$ . The largest gains occur on nos7 ( $1.55\times$ ), shallow_water2 ( $1.54\times$ ), and nos6 ( $1.48\times$ ). On the two largest matrices, ecology1 and ecology2, CAST-1 is $1.25\times$ and $1.26\times$ faster, respectively.

Build costs are similar for the two methods on most matrices, whereas solve costs generally favor CAST-1. The aggregate improvement therefore arises primarily from the solve phase rather than from cheaper factor construction. The exceptions are apache1 ( $0.73\times$ ), jnlbrng1 ( $0.77\times$ ), and bcsstm24 ( $0.85\times$ ). On apache1 and jnlbrng1, CAST-1 requires more PCG iterations; CAST-1 requires 54 iterations per solve on apache1, compared with 28 for AC.

Doubled-granularity comparison.

CAST-2 is faster than AC2 on 11 of the 17 matrices, with an average speedup of $1.069\times$ . It performs particularly well on shallow_water2 ( $1.42\times$ ), the two ecology matrices ( $1.27$ – $1.31\times$ ), and apache1 ( $1.22\times$ ). Thus, the $\rho=2$ construction recovers the loss observed for CAST-1 on apache1, consistent with improved robustness under finer splitting.

On jnlbrng1, AC2 requires 14 iterations per solve, compared with 27 for CAST-2, yielding a speedup ratio of $0.64\times$ . The largest loss occurs on Andrews ( $0.50\times$ ), whose highly skewed degree distribution leads to a CAST-2 build cost of $2.09\,\mu\mathrm{s}/\mathrm{nnz}$ , compared with $0.63\,\mu\mathrm{s}/\mathrm{nnz}$ for AC2. This construction overhead makes CAST-1, rather than CAST-2, the preferable CAST configuration on this matrix.

Cost decomposition and factor reuse.

The performance differences reflect both iteration count and preconditioner-application cost. On ecology1, for example, CAST-2 requires 26.3 iterations per solve, compared with 29.5 for AC2, while also having a lower per-iteration application cost. Together, these effects produce the observed $1.31\times$ total-time speedup.

Factor reuse is particularly important for the doubled-granularity variants. In a separate solve-count sweep, at $q=1$ , AC2 is $51\%$ more expensive than AC and is slower on all 17 matrices, while each CAST variant remains within a few percent of its corresponding baseline. The additional construction cost is amortized as the number of right-hand sides increases; at $q=250$ , the aggregate comparison modestly favors CAST-2.

Table 2: Base-granularity comparison (AC versus CAST-1) on the 17 nontrivial SuiteSparse matrices.

T_{250}

is the total time for factor construction and 250 PCG solves to relative residual

10^{-8}

. Each entry is reported in seconds as the median over five independently seeded factor draws, and the speedup is the ratio of these medians,

T_{250}(\textit{AC})/T_{250}(\textit{CAST-1})

. Build and per-right-hand-side solve costs are normalized by the number of input nonzeros.

Excluded (DNF): bcsstm25 — its conditioning places the $10^{-8}$ target below attainable double-precision accuracy,
		$T_{250}$ (s)			solve $\mu$ s/nnz		build $\mu$ s/nnz
Matrix	nnz	AC	CAST-1	speedup	AC	CAST-1	AC	CAST-1
ecology2	4,995,991	227.96	181.51	1.26 $\times$	0.174	0.139	0.12	0.12
ecology1	4,996,000	226.49	181.90	1.25 $\times$	0.171	0.138	0.12	0.12
apache1	542,184	14.26	19.59	0.73 $\times$	0.097	0.137	0.09	0.09
Andrews	760,154	13.84	11.29	1.23 $\times$	0.070	0.059	0.19	0.23
shallow_water2	327,680	5.46	3.55	1.54 $\times$	0.059	0.044	0.07	0.10
torsion1	197,608	3.80	3.29	1.16 $\times$	0.074	0.066	0.07	0.08
obstclae	197,608	3.79	3.24	1.17 $\times$	0.073	0.065	0.07	0.08
shallow_water1	327,680	3.46	2.43	1.42 $\times$	0.039	0.029	0.06	0.10
jnlbrng1	199,200	3.35	4.36	0.77 $\times$	0.067	0.087	0.07	0.07
nopoly	70,842	1.48	1.29	1.14 $\times$	0.081	0.073	0.05	0.06
fv3	87,025	1.26	1.23	1.03 $\times$	0.057	0.056	0.05	0.06
fv2	87,025	0.72	0.64	1.12 $\times$	0.033	0.029	0.05	0.05
fv1	85,264	0.71	0.62	1.14 $\times$	0.033	0.029	0.05	0.05
nos7	4,617	0.102	0.066	1.55 $\times$	0.086	0.055	0.07	0.08
gr_30_30	7,744	0.085	0.078	1.09 $\times$	0.043	0.039	0.05	0.06
bcsstm24	3,562	0.084	0.100	0.85 $\times$	0.087	0.103	0.04	0.07
nos6	3,255	0.076	0.051	1.48 $\times$	0.092	0.062	0.06	0.07
Arithmetic mean of speedups				1.172 $\times$	wins 14/17
so a subset of right-hand sides never converges for any method.

Table 3: Doubled-granularity comparison (AC2 versus CAST-2) on the same 17 SuiteSparse matrices and under the same protocol as Table 2. Speedup is the ratio of these medians defined as

T_{250}(\textit{AC2})/T_{250}(\textit{CAST-2})

Excluded (DNF): bcsstm25 — its conditioning places the $10^{-8}$ target below attainable double-precision accuracy,
		$T_{250}$ (s)			solve $\mu$ s/nnz		build $\mu$ s/nnz
Matrix	nnz	AC2	CAST-2	speedup	AC2	CAST-2	AC2	CAST-2
ecology1	4,996,000	164.76	125.98	1.31 $\times$	0.124	0.095	0.25	0.28
ecology2	4,995,991	163.20	128.32	1.27 $\times$	0.124	0.097	0.27	0.28
Andrews	760,154	14.77	29.58	0.50 $\times$	0.075	0.144	0.63	2.09
apache1	542,184	12.01	9.85	1.22 $\times$	0.082	0.069	0.24	0.23
shallow_water2	327,680	4.71	3.31	1.42 $\times$	0.051	0.040	0.14	0.12
obstclae	197,608	3.24	2.92	1.11 $\times$	0.065	0.058	0.18	0.15
torsion1	197,608	3.23	2.92	1.11 $\times$	0.065	0.058	0.18	0.15
shallow_water1	327,680	2.95	2.45	1.21 $\times$	0.035	0.030	0.10	0.12
jnlbrng1	199,200	2.70	4.21	0.64 $\times$	0.053	0.083	0.19	0.14
nopoly	70,842	1.06	1.03	1.02 $\times$	0.059	0.058	0.10	0.11
fv3	87,025	0.94	0.98	0.96 $\times$	0.043	0.044	0.12	0.15
fv2	87,025	0.61	0.63	0.97 $\times$	0.027	0.029	0.11	0.11
fv1	85,264	0.60	0.61	0.99 $\times$	0.027	0.028	0.12	0.11
bcsstm24	3,562	0.085	0.099	0.85 $\times$	0.086	0.103	0.07	0.07
nos7	4,617	0.077	0.061	1.26 $\times$	0.065	0.051	0.14	0.16
gr_30_30	7,744	0.068	0.063	1.09 $\times$	0.034	0.032	0.12	0.13
nos6	3,255	0.061	0.049	1.25 $\times$	0.074	0.059	0.08	0.10
Arithmetic mean of speedups				1.069 $\times$	wins 11/17
so a subset of right-hand sides never converges for any method.

Table 4: The ten trivial SuiteSparse matrices, retained as correctness checks and excluded from aggregate comparisons. Every method converges in exactly one PCG iteration on every right-hand side. The reported totals, in seconds, include factor construction and 250 solves to relative residual

10^{-8}

and therefore primarily measure fixed construction and application overhead rather than preconditioner quality.

Matrix	nnz	AC	AC2	CAST-1	CAST-2
bcsstm39	46,772	0.141	0.142	0.153	0.153
t3dl_e	20,360	0.058	0.059	0.063	0.063
t2dal_e	4,257	0.012	0.012	0.013	0.013
bcsstm21	3,600	0.012	0.012	0.012	0.012
bibd_81_2	3,240	0.010	0.010	0.011	0.011
bcsstm23	3,134	0.009	0.008	0.009	0.010
bcsstm26	1,922	0.006	0.006	0.006	0.006
bcsstm11	1,473	0.005	0.004	0.005	0.005
bcsstm08	1,074	0.003	0.003	0.003	0.003
bcsstm09	1,083	0.003	0.003	0.004	0.004

C.2 IPM Sequences on Chimera Graphs

Collection.

We evaluate on the maximum-flow interior-point-method sequences from the SDDM2023 benchmark suite of Gao, Kyng, and Spielman [10]. Each system is a weighted graph Laplacian with $n=100{,}000$ vertices and arises from a Newton step of an interior-point method for undirected maximum flow. The underlying graphs are five independent draws, denoted $i_{1},\ldots,i_{5}$ , from the Chimera generator in Laplacians.jl ³³ 3 https://github.com/danspielman/Laplacians.jl. Chimera graphs combine heterogeneous structures, including grid-like components, star joins, and graph products, and are designed to stress Laplacian solvers.

The five graph instances differ in their elimination geometry. Under CAST-1 elimination, instance $i_{1}$ has mean pivot degree $\bar{d}=26.7$ and maximum pivot degree $d_{\max}=363$ , whereas $i_{3}$ has $\bar{d}=6.4$ and $d_{\max}=48$ (Table 5). For each instance, the benchmark provides IPM runs at five duality-gap targets,

\varepsilon\in\{10^{-1},\ldots,10^{-5}\},

with each run contributing between three and six Newton-step systems. The resulting collection contains $23$ – $28$ matrix systems per instance and $128$ systems in total. Relative performance is stable across the duality-gap targets. We therefore aggregate over Newton steps and targets within each underlying graph instance.

Configuration.

We use the same protocol, metrics, and hardware described in Section 5. Every timed solve satisfies the explicitly verified residual criterion

\frac{\lVert A\mathbf{x}-\mathbf{b}\rVert_{2}}{\lVert\mathbf{b}\rVert_{2}}\leq 10^{-8},

and no method reaches the iteration cap. Within each instance we report the mean total workload over its systems and take speedups as ratios of these means, which measures the aggregate cost of processing an entire IPM sequence. Win counts give the complementary per-system view (Table 5).

Results.

CAST-1 is faster than AC on all $128$ systems (Tables 5 and 6.). The per-instance speedups range from $1.13\times$ to $1.24\times$ , and the per-system ratios range from $1.08\times$ to $1.31\times$ . The normalized solve cost is lower for CAST-1 on all five instances. Build costs are equal or nearly equal on $i_{1}$ and $i_{2}$ and higher for CAST-1 on $i_{3}$ – $i_{5}$ ; at $q=250$ , however, the reduction in solve cost dominates these construction differences.

The behavior of CAST-2 depends more strongly on the elimination-degree profile. On $i_{3}$ , which has the smallest mean and maximum pivot degrees, CAST-2 is the fastest method: it is $1.16\times$ faster than AC2 and wins on all $27$ systems. On the other four instances, the AC2/CAST-2 ratios are $0.35\times$ , $0.76\times$ , $0.45\times$ , and $0.52\times$ , respectively.

This difference is associated with construction fill. At a degree- $d$ pivot, CAST-2 samples one spanning tree on $2d$ auxiliary copies. The tree contains $2d-1$ auxiliary edges before within-block edges are discarded and parallel terminal-edge contributions are aggregated, compared with $d-1$ edges for CAST-1. On instances with high-degree pivots, the additional terminal-edge contributions propagate through subsequent eliminations. Accordingly, the median CAST-2 build cost relative to CAST-1 increases from $3.0\times$ on $i_{3}$ and $4.6\times$ on $i_{2}$ to approximately $10\times$ on $i_{4}$ and $i_{5}$ , and $13.1\times$ on $i_{1}$ . The doubled-granularity baseline is also structure dependent: AC2 has lower mean total cost than AC on $i_{2}$ and $i_{3}$ , but higher cost on $i_{1}$ , $i_{4}$ , and $i_{5}$ .

Choosing the splitting factor.

At $q=250$ , CAST-2 outperforms CAST-1 only on $i_{3}$ . Its factor is $3.0\times$ more expensive to construct on this instance, so CAST-1 remains $1.38\times$ faster at $q=4$ . The two variants cross between $q=10$ and $q=50$ , and CAST-2 is $1.17\times$ faster by $q=250$ . On the four instances with $d_{\max}\geq 92$ , CAST-2 is $1.3$ – $4.6\times$ slower than CAST-1 at $q=250$ . Within this collection, the results support a degree-aware empirical choice of $\rho$ . A light upper tail in the pivot-degree distribution favors $\rho=2$ when the factor is reused sufficiently to amortize its higher construction cost. In contrast, $\rho=1$ is preferable when high-degree pivots make downstream fill and factor construction dominant. The pivot-degree profile observed during a CAST-1 construction may inform this choice, although we treat this as an empirical heuristic rather than as an evaluated adaptive selection algorithm.

Table 5: Chimera-IPM total workload at

q=250

right-hand sides per factor.

T_{250}

includes factor construction and 250 PCG solves to verified relative residual

10^{-8}

. For each system–method pair we report the median over five independently seeded factor draws; displayed times are arithmetic means of these medians over the

23

–

28

systems associated with each underlying graph instance. Parentheses give the number of systems on which the corresponding CAST variant is faster.

				$T_{250}$ (s)				Speedup (mean)
Instance	nnz	$\bar{d}$	$d_{\max}$	AC	AC2	CAST-1	CAST-2	$\textit{AC}/\textit{CAST-1}$	$\textit{AC2}/\textit{CAST-2}$
$i_{1}$	1,100,592	26.7	363	30.6	40.1	24.6	113.8	$1.24\times$ (28/28)	$0.35\times$ (0/28)
$i_{2}$	797,974	9.5	92	19.8	17.6	17.5	23.1	$1.13\times$ (24/24)	$0.76\times$ (0/24)
$i_{3}$	499,696	6.4	48	17.9	15.5	15.7	13.4	$1.14\times$ (27/27)	$1.16\times$ (27/27)
$i_{4}$	814,436	15.2	229	23.1	25.8	19.5	57.4	$1.19\times$ (26/26)	$0.45\times$ (0/26)
$i_{5}$	499,982	10.7	246	13.9	15.6	11.9	30.1	$1.17\times$ (23/23)	$0.52\times$ (0/23)
Across instances: AC/CAST-1 arithmetic mean $\mathbf{1.18\times}$ ; CAST-1 is faster on $128/128$ systems.

Table 6: Chimera-IPM costs normalized by the number of input nonzeros. The relative construction cost of CAST-2 broadly follows the upper tail of the elimination-degree distribution: its build cost is

3.0\times

that of CAST-1 on

i_{3}

, approximately

10\times

i_{4}

and

i_{5}

, and

13.1\times

i_{1}

	Solve ( $\mu\mathrm{s}/\mathrm{nnz}$ )				Build ( $\mu\mathrm{s}/\mathrm{nnz}$ )
Instance	AC	AC2	CAST-1	CAST-2	AC	AC2	CAST-1	CAST-2
$i_{1}$	0.105	0.139	0.087	0.392	0.40	1.15	0.41	5.37
$i_{2}$	0.093	0.083	0.083	0.110	0.16	0.37	0.16	0.73
$i_{3}$	0.133	0.115	0.116	0.099	0.10	0.32	0.15	0.45
$i_{4}$	0.109	0.119	0.091	0.269	0.25	0.68	0.28	2.85
$i_{5}$	0.101	0.114	0.089	0.230	0.18	0.59	0.28	2.80

Table 7: Spielman-IPM total cost at

q=250

right-hand sides per factor.

T_{250}

includes factor construction and

250

PCG solves to verified relative residual

10^{-8}

. Each scale has

10

systems. For each system–method pair we report the median over five independently seeded factor draws; displayed times are arithmetic means of these medians over the systems in each instance, and speedups are ratios of the displayed means, so they measure the aggregate cost of solving an entire IPM sequence. Parentheses report the number of systems on which the CAST variant is faster. AC and AC2 use the buffered application kernel.

Across scales: $\textit{AC}/\textit{CAST-1}$ mean $1.27\times$ ; $\textit{AC2}/\textit{CAST-2}$ mean $\mathbf{4.43\times}$ , with CAST-2 faster on $40/40$ systems.
		$T_{250}$ (s)				Speedup
Scale	nnz	AC	AC2	CAST-1	CAST-2	$\textit{AC}/\textit{CAST-1}$	$\textit{AC2}/\textit{CAST-2}$
$k=100$	1,025,404	7.35	6.84	8.23	1.45	$0.89\times$ (9/10)	$4.72\times$ (10/10)
$k=200$	8,080,604	59.68	54.65	48.19	12.30	$1.24\times$ (9/10)	$4.44\times$ (10/10)
$k=300$	27,226,204	204.60	182.78	132.75	40.42	$1.54\times$ (10/10)	$4.52\times$ (10/10)
$k=400$	64,321,604	409.21	403.72	290.87	99.81	$1.41\times$ (10/10)	$4.04\times$ (10/10)

Table 8: Spielman IPM ladder normalized by the number of input nonzeros, under the same protocol as Table 7.

	solve $\mu$ s/nnz				build $\mu$ s/nnz
Scale	AC	AC2	CAST-1	CAST-2	AC	AC2	CAST-1	CAST-2
$k=100$	0.0248	0.0247	0.0204	0.0054	0.032	0.061	0.049	0.049
$k=200$	0.0248	0.0247	0.0157	0.0055	0.043	0.072	0.070	0.069
$k=300$	0.0247	0.0246	0.0161	0.0055	0.042	0.103	0.073	0.077
$k=400$	0.0247	0.0248	0.0168	0.0055	0.044	0.121	0.076	0.081

Table 9: Sachdeva-star results at

q=250

right-hand sides per factor. Iteration counts are means over the

250

solves, and

T_{250}

includes factor construction and

250

PCG solves to relative residual

10^{-8}

. For every instance–method pair, each entry is the median over three independently seeded factor draws. Speedups are ratios of the displayed total times. The PCG iteration cap is

5\times 10^{3}

and is never reached.

		Iterations per solve				$T_{250}$ (s)				Speedup
Scale	nnz	AC	AC2	CAST-1	CAST-2	AC	AC2	CAST-1	CAST-2	$\frac{AC}{\textit{CAST-1}}$	$\frac{ACtwo}{\textit{CAST-2}}$
$k=50$	62,551	50.8	27.6	37.9	25.2	0.88	0.46	0.55	0.37	$1.58\times$	$1.23\times$
$k=100$	500,101	108.9	37.9	68.3	28.9	15.8	5.49	7.77	3.32	$2.04\times$	$1.65\times$
$k=200$	4,000,201	279.1	51.9	139.0	36.3	349	66.6	127	33.2	$2.75\times$	$2.00\times$
$k=400$	32,000,401	693.4	58.6	287.9	42.8	7,519	638	1,986	298	$3.79\times$	$2.14\times$
$k=600$	108,000,601	1,120.8	65.0	352.9	44.7	42,804	2,539	8,072	1,033	$5.30\times$	$2.46\times$
CAST-2 is the fastest method on every instance and every solve count measured.

C.3 IPM Sequences on Spielman Graphs

Collection.

The second IPM family from the SDDM2023 benchmark suite [10] forms a scaling ladder with one Spielman graph for each size parameter $k\in\{100,200,300,400\}$ . Each graph contributes the Laplacians from the ten Newton steps of a maximum-flow interior-point-method run, yielding $40$ systems in total. The systems range from $n=3.4\times 10^{5}$ to $2.1\times 10^{7}$ vertices and contain up to $6.4\times 10^{7}$ nonzeros.

These graphs are structurally close to trees but have highly heterogeneous edge weights. Most of their vertices have degree at most two, while the number of edges beyond a spanning tree ranges from $5{,}100$ to $80{,}400$ . Each graph also contains $k$ hub vertices of input degree $2k$ . Across the IPM trajectory, edge weights span between $9.1$ and $14.3$ orders of magnitude, increasing primarily during the first approximately four Newton steps and then stabilizing.

Despite the high input degrees of the hubs, minimum-degree elimination produces no large pivots. Path vertices are eliminated first, and every hub has degree at most five when it is eventually eliminated. The mean pivot degree is $2.00$ at every scale. Both CAST variants produce factors of the same measured size.

Configuration.

We use the protocol and metrics of Section 5. Unless otherwise stated, each factor is evaluated on $q=250$ right-hand sides. For every system–method pair, the reported value is the median over five independently seeded factor draws. For each scale ( $10$ systems) we report the mean total workload over its systems and take speedups as ratios of these means, which measures the aggregate cost of processing an entire IPM sequence. Win counts give the complementary per-system view. The larger $k\in\{500,600\}$ sequences, containing approximately $1.5\times 10^{8}$ and $2.2\times 10^{8}$ nonzeros, exceed the memory capacity of the benchmark machine for every method under the multi-solve protocol and are therefore excluded.

$\rho=2$ : fast iteration convergence.

CAST-2 converges in one to two PCG iterations on every recorded solve across all four problem scales. The worst verified relative residual is $5.5\times 10^{-9}$ , which is below the requested tolerance of $10^{-8}$ . At $q=250$ , the available comparisons give mean speedups of $4.04$ – $4.72\times$ over AC2.

The normalized solve cost of CAST-2 remains between $0.0054$ and $0.0055\,\mu\mathrm{s}/\mathrm{nnz}$ on the reported scales. The improved convergence does not increase the measured factor size: CAST-1 and CAST-2 both produce factors containing approximately $0.66$ times as many nonzeros as the input, and their construction costs are similar at the reported scales.

This behavior is associated with the uniformly small pivot degrees. At a degree- $d$ pivot, where $d\leq 5$ throughout this family, CAST-1 samples a spanning tree with $d-1$ edges. In contrast, CAST-2 samples one spanning tree on $2d$ auxiliary copies, containing $2d-1$ auxiliary edges before within-block edges are discarded and parallel terminal-edge contributions are aggregated. Two-way splitting therefore increases local sampling granularity while remaining inexpensive when $d$ is small.

The low-iteration convergence is an empirical property of this graph family. It does not imply that CAST-2 reproduces every Schur clique exactly, nor is CAST-2 the union of two independent trees on the original terminal neighborhood. Rather, the resulting factor is sufficiently accurate for PCG to satisfy the requested tolerance after one to two iterations on every tested Spielman system.

$\rho=1$ : reliability on extreme-weight systems.

The base $\rho=1$ estimator is less consistent on these extreme-weight, near-tree systems, where an unlucky draw can produce a noticeably worse factor. At $k=100$ this is visible in the aggregate, which favors AC despite CAST-1 being faster on nine of ten systems: one system dominates the scale total (Table 7).

Comparison with Chimera-IPM.

The Spielman and Chimera families exhibit complementary behavior. On the Chimera instances with heavy elimination-degree tails (Section C.2), the additional local contributions introduced by $\rho=2$ propagate through subsequent eliminations and substantially increase construction cost. On the Spielman family, every pivot has degree at most five, so the same increase in local sampling granularity remains inexpensive and yields convergence in one to two iterations.

Together, these results support a degree-aware empirical choice of the splitting factor. A uniformly light elimination-degree distribution favors $\rho=2$ when the factor is reused, whereas $\rho=1$ is preferable when high-degree pivots make downstream fill and construction cost dominant. We treat this as an empirical selection heuristic rather than as an evaluated adaptive algorithm.

C.4 Sachdeva-star stress tests

Collection.

The Sachdeva-star family is a synthetic construction on which the AC estimator is known to require rapidly increasing iteration counts [10]. For an even parameter $k$ , each instance consists of a hub connected by unit-weight gateway edges to $k/2$ disjoint cliques of size $k$ . All edges have unit weight. The resulting Laplacian has

n=1+\frac{k^{2}}{2}\qquad\text{and}\qquad\mathrm{nnz}=\frac{k^{3}}{2}+k+1.

We evaluate

k\in\{50,100,200,400,600\},

corresponding to $n=1{,}251$ through $180{,}001$ vertices and up to $108{,}000{,}601$ stored nonzeros.

Unlike the Spielman family, whose difficulty arises from extreme weights on near-tree graphs, the Sachdeva-star instances are unweighted and structurally challenging. Eliminations within the clique blocks create dense local neighborhoods, and the quality of the base-granularity estimators deteriorates as $k$ increases.

Configuration.

We use the protocol and metrics of Section 5. Each factor is evaluated on $q=250$ right-hand sides, and for every instance–method pair the reported value is the median over three independently seeded factor draws. Because iteration growth is a primary quantity of interest on this family, we raise the PCG iteration cap from $10^{3}$ to $5\times 10^{3}$ . No method reaches this cap. Every solve satisfies the verified residual criterion

\frac{\lVert A\,\mathbf{x}-\mathbf{b}\rVert_{2}}{\lVert\mathbf{b}\rVert_{2}}\leq 10^{-8}.

Results.

Table 9 shows a clear separation between the base- and doubled-granularity variants. The mean iteration count of AC increases from $50.8$ at $k=50$ to $1{,}120.8$ at $k=600$ . CAST-1 also becomes less effective as the instances grow, but its iteration count increases more slowly, from $37.9$ to $352.9$ . The doubled-granularity variants are substantially more stable. Across the same range, AC2 increases from $27.6$ to $65.0$ iterations per solve, while CAST-2 increases from $25.2$ to $44.7$ . At $q=250$ , CAST-2 is faster than AC2 on every instance, with speedup increasing monotonically from $1.23\times$ at $k=50$ to $2.46\times$ at $k=600$ . Relative to AC, its total-time advantage at the largest instance is

\frac{42{,}804}{1{,}033}\approx 41.4\times.

A separate solve-count sweep shows that CAST-2 is already the fastest method at $q=1$ , so its advantage on this family does not depend on amortizing construction cost.

Interpretation.

Unlike on the Spielman family, CAST-2 does not reduce these systems to one or two PCG iterations. The clique-dominated elimination structure remains nontrivial after contraction, and its mean iteration count grows from $25.2$ to $44.7$ along the ladder. Nevertheless, the growth is mild compared with the base-granularity methods. Although clique blocks produce high-degree pivots, the additional contributions of $\rho=2$ fall within blocks that are already dense, so they do not propagate as new fill — the opposite of the Chimera instances with heavy degree tails, where the extra contributions land on sparse neighborhoods and compound through later eliminations. This behavior is consistent with the local second-moment improvement of the $\rho=2$ construction: finer splitting reduces local sampling variability, although the theory does not by itself imply a global iteration bound. The Spielman and Sachdeva families therefore illustrate two distinct benefits of $\rho=2$ . On the low-degree Spielman eliminations, CAST-2 empirically produces an almost exact preconditioner and converges in one to two iterations. On the clique-dominated Sachdeva instances, it does not eliminate iteration growth, but it limits that growth sufficiently to provide the best total time throughout the tested ladder. These results agree qualitatively with those reported in [10]: the method ordering and the growth of the base-granularity iteration counts are unchanged, although absolute iteration counts here are approximately $15\%$ – $30\%$ higher, since our right-hand sides are Gaussians projected onto $\mathbf{1}^{\perp}$ rather than the range-restricted vectors used there.

CAST: Canonical Approximate Schur Tree for Approximate Cholesky on Graphs

Abstract

1 Introduction

Question 1.1.

Independent contemporaneous work.

2 Preliminaries and Problem Setting

Graph Laplacians and SDDM systems.

The factor-reuse setting.

Exact elimination and Schur cliques

The local Schur-clique replacement problem.

3 CAST: Canonical Approximate Schur Trees

3.1 The CAST update

CAST-1: the base ρ=1\rho=1 update.

3.2 Exact sampling in 𝒪⁡(ρ​d)\mathcal{O}(\rho d) time via Prüfer codes

Proposition 3.1 (Sampling cost).

Proof.

The CAST preconditioner.

4 Local Theory of the CAST Update

Auxiliary-copy notation.

Unbiasedness.

Lemma 4.1 (Auxiliary-edge marginals).

Proof.

Theorem 4.2 (Unbiasedness).

Proof.

Connected support and correct local nullspace.

Theorem 4.3 (Connected local support).

Proof.

Corollary 4.4 (Correct local nullspace and rank).

Proof.

4.1 Local minimax optimality of the base ρ=1\rho=1 update

Theorem 4.5 (Local minimax optimality of leverage-score marginals).

Proof.

4.2 Local second-moment certificate for CAST-ρ\rho

Lemma 4.6 (Normalized auxiliary-edge contribution).

Proof.

Lemma 4.7 (Covariance domination for spanning-tree indicators).

Proof.

Theorem 4.8 (Local second-moment bound for CAST-ρ\rho).

Proof.

Scope of the guarantees.

5 Experimental Evaluation

Methods.

AC2 versus CAST-2.

Metrics.

Setup.

Matrix collections.

5.1 Cross-collection results

Base granularity.

Doubled granularity.

Choosing ρ\rho.

6 Related Work

Laplacian and SDDM solvers.

Randomized approximate elimination.

Practical approximate-Cholesky solvers.

Leverage scores and random spanning trees.

Conclusion

References

Appendix A Background on Prüfer Codes

The Prüfer correspondence.

Degrees from code multiplicities.

Nonuniform Prüfer sampling.

Lemma A.1 (Weighted Prüfer sampling).

Proof.

Application to CAST-ρ\rho.

Appendix B Proofs

B.1 Proof of Proposition 3.1

Proof.

B.2 Proof of Lemma 4.1

Proof.

B.3 Proof of Theorem 4.2

Proof.

B.4 Proof of Theorem 4.3

Proof.

B.5 Proof of Corollary 4.4

Proof.

B.6 Proof of Theorem 4.5

Proof.

B.7 Proof of Lemma 4.6

Proof.

B.8 Proof of Lemma 4.7

CAST-1: the base $\rho=1$ update.

3.2 Exact sampling in $\mathcal{O}(\rho d)$ time via Prüfer codes

4.1 Local minimax optimality of the base $\rho=1$ update

4.2 Local second-moment certificate for CAST- $\rho$

Theorem 4.8 (Local second-moment bound for CAST- $\rho$ ).

Choosing $\rho$ .

Application to CAST- $\rho$ .

$\rho=2$ : fast iteration convergence.

$\rho=1$ : reliability on extreme-weight systems.