Mapping and Advancing the Scalability-Accuracy Frontier of Nonlinear Causal Discovery

• The comparison covers combinatorial score-based search, score matching [17] [24] [26], differentiable structure learning [10] [19] [33] [37], and amortized discovery [9] [11] [21].Source: §1 · pp. 1–2
• The paper empirically maps the scalability–accuracy frontier of four major nonlinear causal-discovery method families as variables and samples increase.Source: §1 · p. 2
• Repeated nonlinear local-score evaluation is identified as a key bottleneck limiting otherwise accurate combinatorial methods.Source: §1 · p. 2
• Spade combines a spline-accelerated local score with Topic search, reusing sufficient statistics to reduce repeated nonlinear scoring and improve the observed scalability–accuracy frontier.Source: §1 · p. 2

• The setup assumes an i.i.d. sample from a causally sufficient structural causal model whose DAG edges represent direct causes.Source: §2 · p. 2
• A score-based search can select acyclic parent sets, while a causal order constrains edges to run from earlier to later variables and limits each node's candidates to its predecessors.Source: §2 · p. 2
• In an additive-noise model, each variable is a function of its parents plus an exogenous noise term; under standard regularity conditions, nonlinear ANMs are generically identifiable from observational data.Source: §2 · p. 2
• The paper focuses on ANMs because scalability remains difficult even in this simpler, fully identifiable nonlinear setting, despite broader identifiable functional causal models [7].Source: §2 · p. 2

3 SCALABILITY–ACCURACY TRADE-OFFS

• The comparison covers combinatorial scoring, differentiable structure learning, score matching, and amortized inference; these replace graph search with continuous optimization, density-derivative estimation, or pretrained prediction in different ways.Combinatorial scores evaluate how proposed parents explain a child; differentiable methods optimize a relaxed adjacency matrix and nonlinear mechanisms subject to acyclicity.Source: §3 · pp. 2–3
• Amortized methods predict graph distributions from datasets after simulator-based pretraining, shifting test-time computation to a single forward pass.The simulator distribution partly determines their inductive biases; the evaluated methods are Avici and CauScale.Source: §3 · p. 3

3.1 Experimental Setup

• The controlled nonlinear additive-noise experiments vary dimension at n = 2 500 and sample size at d = 100, using five random DAGs per feature dimension.The feature experiment uses expected degree three; runs exceeding 24 hours time out, and GPU runs exceeding 40 GB are recorded as OOM.Source: §3.1 · p. 3

3.2 Results

• Topic achieves the highest F1 and degrades only mildly with dimension, while SDCD and CauScale reach the largest dimensions but remain well below the best F1 scores.Topic's cubic number of local-score evaluations, each requiring a nonlinear model fit, makes it increasingly expensive as graphs grow.Source: Figure 2 · p. 4, §3.2 · p. 4
• The compared representatives are SDCD, DAS, CauScale, and Topic for differentiable, score-matching, amortized, and combinatorial methods, respectively.The paper reports full metric tables in Appendix A.6.Source: §3.2 · p. 4
• DAS leads score-matching methods in feature scaling and is competitive at d = 25, but its accuracy deteriorates rapidly as dimension increases.The paper suggests high-dimensional estimation of score-derived quantities is the bottleneck.Source: §3.2 · p. 4
• As sample size grows, Topic remains the strongest structural performer, while DAS reaches its time limit at 10k samples and CauScale runs out of memory at 20K.DAS uses a Stein estimator with cubic runtime in sample size; CauScale's accuracy worsens as sample size increases.Source: §3.2 · p. 4
• Repeated nonlinear local-score evaluation, rather than a lack of signal, limits Topic; caching sample-dependent computations for additive spline scores enables faster search.The paper identifies this as an opportunity to improve scalability without sacrificing structural accuracy.Source: §3.2 · p. 5

4.1 Spline-Accelerated Score Computation

• The local score combines a conditional-mean fit with residual negative log-likelihood and a BIC complexity penalty.The Gaussian likelihood uses cached quadratic quantities; the histogram likelihood estimates a more flexible residual density but adds sample-dependent scoring.Source: §4.1 · pp. 5–6
• Fixed-basis spline regression lets Spade compile sample-dependent algebra once and reuse it across repeated local-score queries during combinatorial search.The cached formulation is equivalent to direct fitting up to numerical precision for fixed bases and penalties, while avoiding repeated passes over observations.Source: §4.1 · pp. 5–6, §4 · p. 5

4.2 Search

• Spade uses Topic’s topological-order search, which selects the next node, adds score-improving outgoing edges, and prunes redundant incoming edges.The paper replaces Topic’s separate cubic-spline fits and two-part MDL score with a shared fixed spline basis, smoothing penalty, and BIC-type criterion.Source: §4.2 · p. 6

4.3 Complexity

• With fixed basis size K, cache construction costs O(nd^2) time and O(d^2) storage; Spade-Hist additionally stores O(nd) spline-design data.For candidate parent-set size p, cached Gaussian scoring costs O(p^3), while histogram scoring adds O(np) residual scoring per query.Source: §4.3 · p. 6
• Under bounded indegree, Spade-Gauss reduces total score-evaluation cost from O(nd^3) to O(nd^2 + d^3), while Spade-Hist remains O(nd^3).Caching still removes repeated mean-fit computations for both variants.Source: §4.3 · p. 6

5.1 Scalability

• At d = 800, Spade-Gauss finishes in 99.3 seconds and Spade-Hist in 386.9 seconds, versus Topic’s 12.5 hours; at d = 1600, Topic exceeds 24 hours while both Spade variants finish within the limit.In feature scaling, both Spade variants achieve higher F1 than Topic; SDCD and CauScale are faster at the largest feature sizes but have substantially lower F1.Source: §5.1 · p. 7
• At the largest sample size where Topic terminates, Spade-Gauss and Spade-Hist are 2708× and 291× faster, respectively; at n = 160 000, both retain the strongest structural accuracy.Spade variants continue improving with additional samples while Topic plateaus.Source: §5.1 · p. 7

5.2 Robustness Across Data-Generating Regimes

• Across marginalized synthetic regimes, one of the Spade variants obtains the best F1; Spade-Hist is strongest overall and more robust to non-Gaussian noise, while Spade-Gauss performs best in Gaussian-noise settings.The evaluated grid varies graph density, mechanism class, and noise family.Source: §5.2 · pp. 7–8
• Across synthetic settings, Spade-Hist averages 0.748 F1 versus Topic’s 0.616, while reducing mean runtime from 75.35 seconds to 1.24 seconds, roughly 61× faster.It also exceeds CauScale’s 0.542 F1 while running about 5× faster than CauScale’s 6.23-second average.Source: §5.2 · p. 8

5.3 Real-World Data

• On the Causal Chamber benchmark, Spade-Hist reaches 0.50 F1 in 0.97 seconds, compared with SDCD’s 0.49 F1 in 520 seconds and CauScale’s 0.27 F1 in 10 seconds.The paper reports that Spade remains among the strongest and fastest methods on the Sachs protein-signaling dataset as well.Source: §5.3 · p. 8

6 CONCLUSION

• Across the considered benchmarks, combinatorial methods provide the strongest signal but are limited by repeated nonlinear local-score evaluations.Spade’s fixed-basis spline representation compiles and reuses sample-dependent computations during sparse order search.Source: §6 · p. 8
• The authors argue that nonlinear causal-discovery scalability should be assessed across the full scalability–accuracy frontier and all major algorithmic families.Such comparisons can localize where improvements occur and make trade-offs explicit.Source: §6 · p. 8

6.1 Limitations & Future Work

• The comparison is empirical, so the observed strengths and bottlenecks describe the studied benchmarks rather than universal properties of each method family.The authors identify further simulators and application domains as future replication settings.Source: §6.1 · p. 8
• The study focuses on causally sufficient, fully identifiable additive-noise models; extensions to partially identifiable CPDAGs, heteroscedastic mechanisms, latent confounding, and broader functional models remain open.Source: §6.1 · p. 8

AI use statement

• The authors report using generative AI tools for discussion, formulation, interpretation, writing, readability, and coding, with review and revision of AI-assisted material.They state that they take responsibility for the final work.Source: AI use statement · p. 9

CHECKLIST

• The paper states that its formal claims are computational-complexity derivations rather than novel theorem or recovery-guarantee claims.The checklist says assumptions for the complexity analysis are stated in Section 4.3.Source: CHECKLIST · p. 11
• The supplementary reproduction package includes experimental code, a containerized environment, paper presets, seeds, data-loading instructions, and reproduction commands.Source: CHECKLIST · p. 11
• Appendix materials report data-generation procedures, metric definitions, seed counts, standard errors, hyperparameter selection, and compute infrastructure.Source: CHECKLIST · p. 11
• The implementation is included in the supplementary package, but the manuscript does not list license information for external datasets, code, and pretrained models.Source: CHECKLIST · p. 11

Combinatorial additive-score methods

• CAM, Topic, and Globe use nonlinear local scores, but differ in whether they separate order search from edge selection, interleave those operations, or search directly over directed graphs.CAM separates order search and edge selection; Topic interleaves ordering, insertion, and pruning; Globe greedily modifies parent sets.Source: §A.1 · p. 13

Differentiable structure learning

• NoTears [38], Dagma [1], and SDCD [19] replace discrete graph search with continuous adjacency optimization and enforce acyclicity through differentiable formulations.NoTears uses a smooth acyclicity equality, Dagma a log-determinant characterization, and SDCD a quadratic-structure constraint with staged training.Source: §A.1 · p. 13

Score-matching and ordering methods

• Score [24] and NoGAM [17] infer causal orders from density-score derivatives, with NoGAM relaxing the Gaussian-noise restriction of Score.Both use score-derived quantities to infer an order; NoGAM accommodates arbitrary noise distributions under additive-model assumptions.Source: §A.1 · p. 13
• DiffAN [26] and SciNO estimate score-derived Hessian information for ordering, while Caps [32] uses parent scores to refine orders for mixed linear and nonlinear relations.DiffAN uses a diffusion-model score estimator; SciNO approximates the Hessian diagonal with a neural operator; Caps prunes and adds candidate parents after order estimation.Source: §A.1 · p. 14

Amortized causal discovery

• Avici [11] and CauScale [21] predict graph structure through neural inference after pretraining on simulated structural causal models.Avici maps a dataset to edge probabilities in a forward pass; CauScale compresses embeddings and ties attention weights to improve memory and runtime.Source: §A.1 · p. 14

Constraint-based references

• PC and FCI use conditional-independence tests to recover equivalence-class or partial-ancestral-graph structures rather than generally identifying a single fully oriented DAG.PC assumes causal sufficiency; FCI allows latent confounding and selection effects.Source: §A.1 · p. 14

Synthetic data

• Synthetic datasets use causally sufficient SCMs with independent additive noise and order-compatible Erdős–Rényi DAGs, with each eligible edge sampled at a probability set by target expected degree.Variables are generated in topological order, and observed variables are standardized to zero empirical mean and unit empirical variance.Source: §A.2 · p. 14

Controlled scalability

• The controlled scalability setting fixes expected degree at 3 and uses Gaussian-process mechanisms with Gaussian noise, varying features or samples across the specified grids.The feature-scaling experiment fixes n = 2,500; sample sizes range from 2,500 to 160,000.Source: §A.2 · p. 15

Heterogeneous scalability

• The heterogeneous scalability setting samples node-wise mechanism classes and noise families, including linear, additive, non-additive, and interaction-based mechanisms.Its graph-density and problem-size grids match the controlled experiment, and each setting uses five independent seeds.Source: §A.2 · p. 15

Qualitative synthetic comparison

• The qualitative comparison covers 96 synthetic data-generating settings and 480 datasets, varying graph, mechanism, and noise factors.Each setting has five independently generated datasets, with homoscedastic noise and column-standardized variables.Source: §A.2 · p. 15

Differentiable-mechanism diagnostic

• The differentiable-mechanism diagnostic uses neural SCMs at d ∈ {10, 25, 50, 100}, with n = 1000 and non-additive parent interactions when nodes have multiple parents.The neural mechanisms use one hidden layer; each parent pair contributes an additional random neural interaction.Source: §A.2 · pp. 15–16

Real-world data

• The real-world benchmarks are Causal Chamber with 10,000 observations of 23 variables and Sachs with 853 observations of 11 variables.Causal Chamber has 39 reference edges; Sachs has 17, and both settings use five bootstrap samples.Source: Real-world data · p. 16

A.3 Hardware

• GPU methods used one NVIDIA A100 40 GB GPU and four CPU threads; CPU methods used four threads on AMD EPYC 7773X nodes, with no GPU.Reported runtimes are wall-clock runtimes of the benchmark invocation.Source: §A.3 · p. 16

Heterogeneous benchmark

• At d = 1600 in the heterogeneous feature-scaling experiment, Topic does not finish within 24 hours, while Spade-Gauss and Spade-Hist finish in 918 and 2416 seconds.Spade-Hist has the strongest accuracy; in sample scaling it reaches F1 0.829 versus 0.724 for Spade-Gauss and 0.622 for Topic.Source: §A.4 · pp. 16–17

Differentiable Structure Learning Experiment

• No differentiable method reliably ranks the true DAG highest among 1000 random-order DAGs, whereas the spline score ranks it highest consistently despite non-additive generating mechanisms.The diagnostic uses frozen learned mechanisms to score candidate DAGs; the authors interpret the failure as suggesting poor local optima in the joint graph–mechanism space.Source: §A.4 · p. 17

Runtime Ablation

• Caching compiled Gram matrices and feature–response cross-products significantly speeds up SPADE, especially Spade-Gauss.The runtime ablation compares cached variants with versions that refit the additive-spline score from scratch.Source: §A.4 · p. 17, Figure 7 · p. 18

Qualitative results at n = 2000

• At n = 2000, the three strongest competitors differ from those at n = 1000, while SPADE variants show similar patterns across the two sample sizes.The n = 2000 comparison is presented in Fig. 8.Source: §A.4 · p. 17

Sachs

• On Sachs, SPADE variants rank among the three highest F1 scores in all three data settings while achieving the lowest SHD and runtimes.The settings use raw data and log and arcsinh transformations.Source: §A.4 · pp. 17–18

A.6 Additional Metrics

• Additional metrics include directed and skeleton precision, recall, and F1, structural Hamming distance, normalized topological divergence, and runtime in seconds.SHD counts missing and extra directed edges while treating a reversed edge as one operation; normalized topological divergence is the fraction of true edges with incorrectly ordered endpoints.Source: §A.6 · p. 18

A.5 Cached Sufficient Statistics

• Algorithm 1 caches centered spline sufficient statistics, then evaluates parent-set scores using Gaussian residuals or a histogram residual likelihood.The cache stores G = B^T B, H = B^T Y, and response sums of squares; local scoring extracts parent-set blocks and solves a penalized regression.Source: §A.6 · p. 19

A.6 Additional Metrics

• Appendix tables report means over successful rows, excluding failed, timed-out, or out-of-memory runs.The tables cover controlled feature-size and sample-size experiments, heterogeneous scalability, and the Causal Chamber benchmark.Source: §A.6 · p. 19

A.7 Implementation Details

• Methods use default or author-provided hyperparameters without tuning against benchmark ground truth; computational settings change only when needed for GPU memory.Stochastic methods use one predetermined seed per method–dataset pair, without random restarts or best-run selection.Source: §A.7 · pp. 19–20
• Spade uses centered cubic B-spline bases with ridge regularization, a second-difference smoothing penalty, and deviance-scale scores with a BIC penalty.Spade-Hist replaces Spade-Gauss’s Gaussian residual term with a histogram residual likelihood; both use the same spline basis and search settings.Source: §A.7 · p. 20
• The appendix specifies implementation settings for the compared baselines, including model variants, optimization parameters, thresholds, and conditional-independence tests.For example, PC and FCI use kernel conditional-independence tests, while NoTears and Dagma use nonlinear MLP variants.Source: §A.7 · pp. 20–21