On the Number of Conditional Independence Tests in Constraint-based Causal Discovery
Authors: Marc Franquesa Monés, Jiaqi Zhang, Caroline Uhler
Organizations: Eric and Wendy Schmidt Center, Broad Institute of MIT and Harvard · Centre de Formació Interdisciplinària Superior, Universitat Politècnica de Catalunya · Laboratory for Information and Decision Systems, Massachusetts Institute of Technology
Learning causal relations from observational data is a fundamental problem with wide-ranging applications across many fields. Constraint-based methods infer the underlying causal structure by performing conditional independence tests. However, existing algorithms such as the prominent PC algorithm need to perform a large number of independence tests, which in the worst case is exponential in the maximum degree of the causal graph. Despite extensive research, it remains unclear if there exist algorithms with better complexity without additional assumptions. Here, we establish an algorithm that achieves a better complexity of pO(s) tests, where p is the number of nodes in the graph and s denotes the maximum undirected clique size of the underlying essential graph. Complementing this result, we prove that any constraint-based algorithm must perform at least 2Ω(s) conditional independence tests, establishing that our proposed algorithm achieves exponent-optimality up to a logarithmic factor in terms of the number of conditional independence tests needed. Finally, we validate our theoretical findings through simulations, on semi-synthetic gene-expression data, and real-world data, demonstrating the efficiency of our algorithm compared to existing methods in terms of number of conditional independence tests needed.
Figures & tables
Figure 1: A graph with p nodes, where d=p−2 , and s=2 .
Figure 2: Meek rules ( Meek,, 1995 ) . The dashed lines partition nodes into sets based on ancestral relationships. Only Meek rule 1 establishes an ancestor for a node that previously had none, forcing a component to split.
Figure 3: Example graph illustrating the classification of nodes into sets DS1 and FS1 given the prefix node set S={0,1} .
Algorithm 1 Greedy Ancestral Search ( GAS )
Figure 4: Comparison across Erdős-Rényi graphs on 50 nodes and increasing expected neighborhood size. Results are averaged over 3 runs. (a) Execution time (seconds) presented on a logarithmic scale. (b) Accuracy comparison measured by Structural Hamming Distance (SHD) between predicted and ground-truth graphs normalized by the number of possible edges. Shaded regions show the standard deviation.
Figure 5: Comparison on semi-synthetic data generated using SERGIO. Results are averaged over 5 runs. (a) Execution time (seconds) presented on a logarithmic scale. (b) Accuracy comparison measured by Structural Hamming Distance (SHD) between predicted and ground-truth graphs normalized by the number of possible edges.
Figure 6: Partially directed graphs learned in the Airfoil dataset. GAS , GAS+ , and PC performed 46 , 50 , and 104 distinct CI tests, respectively.
Figure 12: An example demonstrating how Algorithm 1 can recover the correct causal graph when the generating probability distribution is Markov but not faithful to that graph.
Algorithm 2 Algorithm used to remove edges from E in Algorithm 1
Algorithm 3 Algorithm used to calculate DSm
Algorithm 4 Algorithm used to calculate FSm
Algorithm 5 Modification of the final loop in Algorithm 1 (Lines 24 – 25 )
Figure 13: Comparison across Erdős-Rényi graphs on 50 nodes and increasing expected neighborhood size. Results are averaged over 3 runs. (a) Execution time (seconds) presented on a logarithmic scale. (b) Accuracy comparison measured by Structural Hamming Distance (SHD) between predicted and ground-truth graphs normalized by the number of possible edges. Shaded regions show the standard deviation.
Figure 14: Number of distinct CI tests performed across Erdős-Rényi graphs of 50 nodes and increasing expected neighborhood size. Results are averaged over 10 runs.
Figure 15: Comparison across Erdős-Rényi graphs on 50 nodes and increasing expected neighborhood size. Results are averaged over 10 runs. (a) Number of missing edges between predicted and ground-truth skeletons. (b) Number of extra edges between predicted and ground-truth skeletons. Shaded regions show the standard deviation.
Exp. Neigh. Size
Precision
Recall
F1-Score
GAS+
PC
GAS+
PC
GAS+
PC
15
0.42
0.83
0.92
0.29
0.58
0.43
25
0.57
0.81
0.81
0.16
0.67
0.27
35
0.73
0.87
0.77
0.12
0.75
0.21
45
0.93
0.96
0.75
0.09
0.83
0.17
Appendix
Table 1: Additional skeleton metrics for GAS+ and PC across Erdős-Rényi graphs of 50 nodes and increasing expected neighborhood size.
Figure 16: Comparison across Erdős-Rényi graphs of increasing number of nodes and edge probability p=0.5 . Results are averaged over 3 runs. (a), (c) Execution time (seconds) presented on a logarithmic scale. (b), (d) Accuracy comparison measured by the Structural Hamming Distance between predicted and ground-truth graphs normalized by the number of possible edges. Shaded regions show the standard deviation.
Figure 17: Comparison of the number of distinct conditional independence tests across Erdős-Rényi graphs of increasing size. Results are averaged over 10 runs; shaded regions indicate the standard deviation.
m
Exp. Neigh. Size
Structural Hamming Distance
Execution Time (s)
PC
GAS
GAS+
PC
GAS
GAS+
5
8.57
105.00
380.67
283.33
1.7576
0.0377
0.0432
10
14.29
225.00
426.00
289.67
1.4815
0.0085
0.0138
15
17.14
297.00
477.33
296.33
1.7844
0.0074
0.0192
20
17.14
305.00
450.00
294.67
3.3219
0.0071
0.0150
Appendix
Table 2: Comparison across Barabási–Albert graphs of 40 nodes and increasing parameter m . Results are averaged over 3 runs.
Number of Samples
PC
GAS
GAS+
50
568.33
585.00
598.00
100
567.00
579.67
573.33
500
540.00
597.67
514.33
1000
554.33
596.33
500.33
5000
523.00
607.00
492.33
10000
538.00
604.00
476.33
Appendix
Table 3: Skeleton Structural Hamming Distance between the predicted and ground-truth graphs for increasing sample sizes. Data was generated from 50 -node Erdős-Rényi graphs with an edge probability of p=0.5 . Results are averaged over 3 runs.
Conditional Independence (CI) tests are the statistical engine of constraint-based causal discovery: in algorithms such as PC (Peter-Clark) and FCI (Fast Causal Inference), skeleton pruning and key orientations follow directly from CI decisions. This survey reviews CI testing with emphasis on assumptions, robustness, and scalability in high-dimensional and mixed-type settings common in biomedical domains. The survey organizes widely used CI methods into six families: partial-correlation, contingency-table, regression, nearest-neighbor, kernel, and machine-learning-based. Special emphasis is provided on the robustness layers that address the limitations of these families. For each family, the survey examines when CI decisions reflect the data-generating distribution and when they fail. By this, we link test-level properties, including power decay with conditioning set size and asymmetric type I/II error consequences, to graph-level errors in skeleton recovery and v-structure orientation. The survey also compares adoption across major R and Python libraries and summarizes open challenges, including mixed-type CI testing without discretization, small-sample error control, and strategies for improving scalability of CI-testing.
Pavel Averin, Theodoros Moysiadis, Ioannis Katakis
Department of Computer Science School of Sciences and Engineering University of Nicosia 2417, Nicosia, Cyprus
Constraint-based causal discovery relies on repeated conditional independence tests, but fast nonparametric tests often sacrifice calibration, especially when variables depend on the conditioning set through nonlinear relationships. We introduce BLITZ (Broad-to-Local Independence Testing via residualiZation), a nonparametric conditional independence test designed to run well under a second while maintaining the accuracy needed for the thousands of queries performed by constraint-based causal discovery algorithms. BLITZ first removes broad smooth dependence on the conditioning set using low-order polynomial regression, then applies a small nonlinear feature map and residualizes those features with shallow tree regressions. The resulting statistic tests residual cross-covariance, with a moment-matched chi-square approximation to the null distribution. We show theoretically that the two-stage design reduces the effective complexity faced by the tree residualizers, allowing shallow trees to control residual conditional-mean bias while avoiding excessive overfitting. In simulations, BLITZ provides better null calibration than fast kernel, random-feature, and regression-based competitors while remaining among the fastest methods tested. In causal discovery experiments on synthetic graphs and flow-cytometry data, BLITZ yields more reliable endpoint orientations among retained adjacencies and competitive structural recovery. These results suggest that broad-to-local residualization is a practical route to calibrated, scalable nonparametric conditional independence testing for causal discovery.
Eric V. Strobl
Department of Biomedical Informatics, University of Pittsburgh, Pittsburgh, PA 15213, USA
Causal-discovery algorithms return a directed graph, yet provide no principled means of distinguishing edge directions identified by the data from those assigned without an identifying assumption. Under the standard Markov and faithfulness conditions, the observational distribution identifies only a Markov equivalence class; orientations within that class are not determined by the joint distribution and cannot be recovered from additional samples alone, but require either a functional restriction or an intervention. We introduce a protocol for observational causal discovery on continuous data that attaches to each candidate edge a discrete impossibility certificate: a RESOLVED code records the identifiability theorem under which the direction was committed, while an IMPOSSIBLE code records the failure mode together with the specific question a domain expert must answer to resolve it. The bivariate cascade is extended with five gated identifiability tiers LSNM, IGCI, Stein, MDL, and PEIT that abstain when their precondition test rejects. Two oracle primitives, the meta-hub query and the node-children query, jointly establish an upper bound of 1+K expert interactions sufficient to recover any DAG, where K denotes the number of non-leaf vertices. Under an ideal-oracle assumption, the bound is met exactly on the asia, sachs, child, and alarm benchmarks.