Score-Schur Topological Sort in Causal Discovery

• Score-Schur Topological Sort (SSTS) is an optimization‐free method that algebraically recovers causal order from the geometric signature of a learned score function.
• It uses the score-Jacobian information matrix (SJIM) and Schur complement updates to identify leaf nodes and eliminate them in a reverse topological order.
• Block-SSTS extends this approach to non-linear systems by grouping nodes to control expectation-gap errors while maintaining scalability.

Searching arXiv for the primary SSTS paper and closely related topological-sorting papers mentioned in the source material. Score-Schur Topological Sort (SSTS) is an optimization-free method for extracting a topological order in continuous causal discovery by using the geometry of a learned score function rather than solving a DAG-constrained optimization problem. In the formulation of "Optimization-Free Topological Sort for Causal Discovery via the Schur Complement of Score Jacobians" (Wu et al., 28 Apr 2026), causal discovery is split into two stages: learning an unconstrained score model for the observational distribution and then recovering causal order algebraically from the Score-Jacobian Information Matrix (SJIM). In linear additive noise models, the method is exact because iterative graph marginalization is mathematically equivalent to repeated Schur complements of the SJIM; in non-linear systems it becomes approximate, and the paper introduces Block-SSTS to control the resulting expectation-gap error while retaining scalability (Wu et al., 28 Apr 2026).

1. Conceptual formulation

SSTS is presented as a response to a specific bottleneck in continuous DAG learning: methods such as NOTEARS and DAGMA are described as coupling representation learning with structural optimization via a non-convex acyclicity penalty, making them sensitive to local optima, hyperparameters, and scaling limitations (Wu et al., 28 Apr 2026). The central shift is therefore from constrained structure optimization to statistical score estimation.

The score function is defined as

$s(x)=\nabla_x \log p(x),$

and the paper’s claim is that the causal hierarchy leaves a geometric signature in this object. Once a score model has been trained, topological sorting is recast as a matrix-elimination problem on the SJIM rather than a search over DAGs (Wu et al., 28 Apr 2026).

The paper characterizes this shift in two contrasting views. The older view is to learn a DAG by minimizing a loss plus an acyclicity penalty. The SSTS view is to learn a score density model and then recover order by Schur complements. This suggests that SSTS is not primarily a graph-search algorithm; it is an algebraic extraction procedure whose correctness depends on structural information encoded in the score geometry.

A common misconception is to equate any topological-sort-based ordering method with SSTS. That is explicitly incorrect for "Topological Sort for Sentence Ordering" (Prabhumoye et al., 2020), which learns pairwise sentence precedence constraints and then applies a standard DFS-based topological sort. That paper contains no Schur decomposition, no score-Jacobian construction, and no Score-Schur mechanism (Prabhumoye et al., 2020). By contrast, "Optimizing the Learning Order of Chinese Characters Using a Novel Topological Sort Algorithm" (Loach et al., 2016) is closer in spirit because it first constructs a score-based priority order and then repairs it into a hierarchy-respecting topological order, but the method is not named SSTS and is not Schur-based (Loach et al., 2016).

2. Mathematical basis: score geometry, SJIM, and leaf identifiability

The theoretical setting is a continuous additive noise model (ANM),

$$x_i = f_i(x_{pa(i)}) + \epsilon_i,\quad i=1,\dots,d,$$

with independent Gaussian noise

$\epsilon_i \sim \mathcal{N}(0,\sigma^2).$

Within this setting, the SJIM is defined as the expected negative Hessian of the log-density,

$$I = \mathbb{E}_{p(x)}\left[-\nabla_x^2 \log p(x)\right].$$

The paper treats $I$ as a structural information matrix that encodes causal topology (Wu et al., 28 Apr 2026).

The principal identifiability result is a leaf-node criterion based on diagonal energy. Under the ANM with homoscedastic Gaussian noise,

$$I_{ii} = \frac{1}{\sigma^2} + \frac{1}{\sigma^2}\sum_{j\in ch(i)} \mathbb{E}_{p(x)}\left[\left(\frac{\partial f_j}{\partial x_i}\right)^2\right].$$

Hence, if $i$ is a leaf node, then $ch(i)=\emptyset$ and

$I_{ii} = \frac{1}{\sigma^2},$

whereas any node with children has strictly larger diagonal energy (Wu et al., 28 Apr 2026).

This yields the operational rule that causal sinks or leaves are exactly the nodes with minimum diagonal SJIM value. The paper’s intuition is that the diagonal of the Hessian or expected Hessian measures how much the score bends with respect to a variable. A leaf contributes no extra curvature through descendants because it appears only in its own equation and not in any child mechanism. A plausible implication is that SSTS converts a combinatorial ordering question into repeated identification of locally minimal curvature coordinates.

3. Exact linear-case mechanism: Schur complements as topological elimination

The exactness claim is established for linear ANMs. Writing

$x = Bx + \epsilon,$

with $$x_i = f_i(x_{pa(i)}) + \epsilon_i,\quad i=1,\dots,d,$$0 strictly upper triangular after topological ordering, one has

$$x_i = f_i(x_{pa(i)}) + \epsilon_i,\quad i=1,\dots,d,$$1

and

$$x_i = f_i(x_{pa(i)}) + \epsilon_i,\quad i=1,\dots,d,$$2

Because the distribution is Gaussian, the SJIM equals the precision matrix,

$$x_i = f_i(x_{pa(i)}) + \epsilon_i,\quad i=1,\dots,d,$$3

This is the matrix on which SSTS operates (Wu et al., 28 Apr 2026).

If $$x_i = f_i(x_{pa(i)}) + \epsilon_i,\quad i=1,\dots,d,$$4 is a leaf node, the paper proves that exact marginalization of that leaf corresponds to the Schur complement

$$x_i = f_i(x_{pa(i)}) + \epsilon_i,\quad i=1,\dots,d,$$5

Equivalently,

$$x_i = f_i(x_{pa(i)}) + \epsilon_i,\quad i=1,\dots,d,$$6

The significance is explicit: removing a leaf from the graph and updating the remaining structure is mathematically identical to applying a Schur complement to the SJIM (Wu et al., 28 Apr 2026).

The resulting elimination process is a topological sort in reverse leaf order. Repeatedly identify a minimum-diagonal leaf, eliminate it via a Schur complement, and continue on the reduced matrix. In the linear Gaussian case, the paper reports zero edge violations across all tested sizes for both empirical and population matrices, and treats this as validation of the exact Schur-complement theorem (Wu et al., 28 Apr 2026).

This exactness is what separates SSTS from earlier score-guided or priority-aware topological ordering heuristics. In the Chinese-character ordering method, for example, a ranking score is repaired into a valid hierarchy order by moving violated prerequisites leftward as little as possible (Loach et al., 2016). SSTS instead derives leaf structure from the SJIM and updates that structure algebraically under marginalization.

4. Algorithmic pipeline and Block-SSTS

The full procedure begins by fitting an unconstrained score network $$x_i = f_i(x_{pa(i)}) + \epsilon_i,\quad i=1,\dots,d,$$7 to observational data by score matching,

$$x_i = f_i(x_{pa(i)}) + \epsilon_i,\quad i=1,\dots,d,$$8

No DAG constraint is imposed at this stage (Wu et al., 28 Apr 2026).

The empirical Jacobian is then accumulated over the dataset,

$$x_i = f_i(x_{pa(i)}) + \epsilon_i,\quad i=1,\dots,d,$$9

and symmetrized to obtain

$\epsilon_i \sim \mathcal{N}(0,\sigma^2).$0

For high-dimensional data, the paper also allows group lasso on input-layer weights to suppress weak spurious connections and stabilize the empirical SJIM (Wu et al., 28 Apr 2026).

At each iteration, SSTS identifies the current leaf block by selecting nodes whose diagonal energy lies within a tolerance band of the minimum:

$\epsilon_i \sim \mathcal{N}(0,\sigma^2).$1

The nodes in $\epsilon_i \sim \mathcal{N}(0,\sigma^2).$2 are sorted inside the block by diagonal energy and prepended to the output order $\epsilon_i \sim \mathcal{N}(0,\sigma^2).$3. The matrix is then updated by a block Schur complement,

$\epsilon_i \sim \mathcal{N}(0,\sigma^2).$4

This repeats until the active set is empty (Wu et al., 28 Apr 2026).

After the topological order is obtained, edge recovery is performed by regressing each node only on earlier nodes in the order,

$\epsilon_i \sim \mathcal{N}(0,\sigma^2).$5

where $\epsilon_i \sim \mathcal{N}(0,\sigma^2).$6 is a sparsity penalty such as Lasso (Wu et al., 28 Apr 2026).

Block-SSTS is introduced because the non-linear case does not preserve exact Schur marginalization under expectation. The grouping of approximately equal minimal diagonals serves two functions stated in the paper: it avoids forcing a spurious sequential order among nodes in the same topological stratum, and it reduces the number of Schur steps from potentially $\epsilon_i \sim \mathcal{N}(0,\sigma^2).$7 to fewer block iterations (Wu et al., 28 Apr 2026). This suggests that Block-SSTS is both a statistical tolerance mechanism and a depth-compression mechanism.

5. Non-linear approximation, expectation gap, and failure modes

For non-linear ANMs, the Hessian

$\epsilon_i \sim \mathcal{N}(0,\sigma^2).$8

depends on the sample, so the key equality from the linear case fails:

$\epsilon_i \sim \mathcal{N}(0,\sigma^2).$9

The paper defines the discrepancy as

$$I = \mathbb{E}_{p(x)}\left[-\nabla_x^2 \log p(x)\right].$$0

For a leaf node $$I = \mathbb{E}_{p(x)}\left[-\nabla_x^2 \log p(x)\right].$$1, Proposition 1 gives

$$I = \mathbb{E}_{p(x)}\left[-\nabla_x^2 \log p(x)\right].$$2

The paper emphasizes that this non-linear mismatch is a covariance term of the leaf’s gradient and that it is localized to the parent block $$I = \mathbb{E}_{p(x)}\left[-\nabla_x^2 \log p(x)\right].$$3, rather than being globally distributed across the graph (Wu et al., 28 Apr 2026).

The practical consequence is that SSTS in non-linear settings is only approximately correct, and repeated marginalization can accumulate error. The paper explicitly describes a tradeoff in the block tolerance $$I = \mathbb{E}_{p(x)}\left[-\nabla_x^2 \log p(x)\right].$$4: if the tolerance is too small, the method performs too many sequential Schur updates and accumulates more error; if it is too large, nodes may be incorrectly grouped as parallel, causing structural errors (Wu et al., 28 Apr 2026).

The limitations identified by the paper are fourfold. First, performance depends on score estimation quality: if the score network is inaccurate, the empirical SJIM is noisy and the diagonal-energy rule can misidentify leaves. Second, finite-sample variance in the empirical Jacobian can generate spurious cross-derivatives; the paper states that this variance, rather than optimization failure, ultimately limits structural fidelity at large $$I = \mathbb{E}_{p(x)}\left[-\nabla_x^2 \log p(x)\right].$$5. Third, the non-linear expectation gap accumulates through elimination. Fourth, the exact theory is sensitive to the ANM assumption with Gaussian homoscedastic noise. The failure-mode analysis reported in the paper states that multiplicative noise models can still work well, whereas post-nonlinear models break the method badly, with very high SHD and very low TPR (Wu et al., 28 Apr 2026).

A further caveat is computational exactness. The appendix reportedly includes a covariance-patching mechanism that can theoretically restore exactness via covariance terms, but the paper states that it is too expensive in memory and compute and is therefore not used in practice (Wu et al., 28 Apr 2026).

6. Empirical behavior, scaling, and relation to adjacent methods

The paper’s dominant-cost claim is that the extraction phase is reduced to matrix operations with worst-case cost $$I = \mathbb{E}_{p(x)}\left[-\nabla_x^2 \log p(x)\right].$$6 for dense matrices, while memory is roughly $$I = \mathbb{E}_{p(x)}\left[-\nabla_x^2 \log p(x)\right].$$7 for storing the SJIM because Jacobians are accumulated in streaming form rather than stored per sample (Wu et al., 28 Apr 2026). This is the computational basis for the method’s scaling claims.

The reported empirical results span linear ANMs, non-linear synthetic benchmarks, and the Sachs protein-signaling dataset. The paper explicitly uses Edge Violations (EV), Structural Hamming Distance (SHD), and True Positive Rate (TPR), and argues that rank-correlation metrics such as Kendall’s $$I = \mathbb{E}_{p(x)}\left[-\nabla_x^2 \log p(x)\right].$$8 are inappropriate for topologically ambiguous DAGs because many valid orders exist (Wu et al., 28 Apr 2026).

Setting	Reported result	Notes
Linear ANMs	Zero edge violations	Up to $$I = \mathbb{E}_{p(x)}\left[-\nabla_x^2 \log p(x)\right].$$9
Non-linear Block-SSTS	$I$0, extraction time $I$1 s	About $I$2 block inversions
Memory-throttled variant	$I$3, extraction time about $I$4 s	Peak VRAM around $I$5 GB
Sachs dataset	SHD $I$6, TPR $I$7, total time about $I$8 s	$I$9, $$I_{ii} = \frac{1}{\sigma^2} + \frac{1}{\sigma^2}\sum_{j\in ch(i)} \mathbb{E}_{p(x)}\left[\left(\frac{\partial f_j}{\partial x_i}\right)^2\right].$$0

In the non-linear synthetic benchmarks with $$I_{ii} = \frac{1}{\sigma^2} + \frac{1}{\sigma^2}\sum_{j\in ch(i)} \mathbb{E}_{p(x)}\left[\left(\frac{\partial f_j}{\partial x_i}\right)^2\right].$$1-based ANMs and $$I_{ii} = \frac{1}{\sigma^2} + \frac{1}{\sigma^2}\sum_{j\in ch(i)} \mathbb{E}_{p(x)}\left[\left(\frac{\partial f_j}{\partial x_i}\right)^2\right].$$2, the paper reports that SSTS is strong and much faster than DAGMA and NOTEARS at $$I_{ii} = \frac{1}{\sigma^2} + \frac{1}{\sigma^2}\sum_{j\in ch(i)} \mathbb{E}_{p(x)}\left[\left(\frac{\partial f_j}{\partial x_i}\right)^2\right].$$3, remains very fast in extraction at $$I_{ii} = \frac{1}{\sigma^2} + \frac{1}{\sigma^2}\sum_{j\in ch(i)} \mathbb{E}_{p(x)}\left[\left(\frac{\partial f_j}{\partial x_i}\right)^2\right].$$4, but begins to lose structural accuracy relative to some masked score methods as the non-linear expectation gap and finite-sample score error become more prominent (Wu et al., 28 Apr 2026). At $$I_{ii} = \frac{1}{\sigma^2} + \frac{1}{\sigma^2}\sum_{j\in ch(i)} \mathbb{E}_{p(x)}\left[\left(\frac{\partial f_j}{\partial x_i}\right)^2\right].$$5, the paper states that EV approaches a near-random baseline, interpreting this as evidence that score-model estimation no longer cleanly separates the hierarchy. The article’s broader conclusion is therefore not that causal discovery becomes trivial once optimization is removed, but that the principal bottleneck changes from constrained optimization to statistical estimation.

Relative to neighboring literature, SSTS occupies a distinct position. The sentence-ordering method of (Prabhumoye et al., 2020) is a constraint-solving use of topological sort over pairwise sentence relations, but it is not SSTS because it contains no score matrix and no Schur-complement elimination. The Chinese-character ordering method of (Loach et al., 2016) is a weighted priority order repaired into a valid topological sequence with minimal disturbance to the original ranking; this is a score-guided topological ordering heuristic, but not an optimization-free causal-order extraction method. SSTS is therefore best understood as a specialized causal-discovery procedure in which topological order is read from score geometry through SJIM diagonal energies and Schur complement updates (Wu et al., 28 Apr 2026).