Paper deep dive
Test-Time Learning of Causal Structure from Interventional Data
Wei Chen, Rui Ding, Bojun Huang, Yang Zhang, Qiang Fu, Yuxuan Liang, Han Shi, Dongmei Zhang
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 91%
Last extracted: 7/20/2026, 8:33:44 PM
Summary
The paper introduces TICL (Test-time Interventional Causal Learning), a novel method for causal discovery from interventional data. TICL combines Test-Time Training (TTT) with Joint Causal Inference (JCI) to address challenges in generalization and versatility when intervention targets are unknown. It employs a self-augmentation strategy using an IS-MCMC algorithm to generate instance-specific training data at test time, avoiding distribution shifts. A PC-inspired two-phase supervised learning scheme is used to predict identifiable causal structures (I-CPDAG) and detect intervention targets, demonstrating superior performance on bnlearn benchmarks.
Entities (8)
Relation Signals (7)
TICL → uses → Test-time Training
confidence 95% · TICL ... synergizes Test-Time Training with Joint Causal Inference.
TICL → uses → Joint Causal Inference
confidence 95% · TICL ... synergizes Test-Time Training with Joint Causal Inference.
TICL → employs → IS-MCMC
confidence 92% · We introduce the IS-MCMC algorithm ... to fuel the self-augmentation process.
TICL → predicts → I-CPDAG
confidence 90% · we focus on identifying the identifiable components in the I-CPDAG
IS-MCMC → generates → Self-Augmentation
confidence 88% · We introduce the IS-MCMC algorithm ... to fuel the self-augmentation process.
PC Algorithm → inspires → Two-Phase SCL
confidence 87% · Inspired by the PC algorithm ... we implement the supervised causal learning as a two-phase process
TICL → outperforms → bnlearn
confidence 85% · Extensive experiments on bnlearn benchmarks demonstrate TICL's superiority
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Supervised causal learning has shown promise in causal discovery, yet it often struggles with generalization across diverse interventional settings, particularly when intervention targets are unknown. To address this, we propose TICL (Test-time Interventional Causal Learning), a novel method that synergizes Test-Time Training with Joint Causal Inference. Specifically, we design a self-augmentation strategy to generate instance-specific training data at test time, effectively avoiding distribution shifts. Furthermore, by integrating joint causal inference, we developed a PC-inspired two-phase supervised learning scheme, which effectively leverages self-augmented training data while ensuring theoretical identifiability. Extensive experiments on bnlearn benchmarks demonstrate TICL's superiority in multiple aspects of causal discovery and intervention target detection.
Tags
Links
- Source: https://arxiv.org/abs/2602.19131v1
- Canonical: https://arxiv.org/abs/2602.19131v1
Trouble viewing inline? Open PDF directly →
Full Text
148,822 characters extracted from source content.
Expand or collapse full text
Test-Time Learning of Causal Structure from Interventional Data Wei Chen ∗ 1,2,3 , Rui Ding 3 , Bojun Huang 4 , Yang Zhang 5 , Qiang Fu 3 , Yuxuan Liang 2 , Han Shi 3 , Dongmei Zhang 3 1 HKUST, 2 HKUST(GZ), 3 Microsoft Research, 4 Sony Research, 5 NUS ∗ Work done during an internship at Microsoft Research. Supervised causal learning has shown promise in causal discovery, yet it often struggles with generalization across diverse interventional settings, particularly when intervention targets are unknown. To address this, we proposeTICL(Test-time Interventional Causal Learning), a novel method that synergizes Test-Time Training with Joint Causal Inference. Specifically, we design a self-augmentation strategy to generate instance-specific training data at test time, effectively avoiding distribution shifts. Furthermore, by integrating joint causal inference, we developed a PC-inspired two-phase supervised learning scheme, which effectively leverages self-augmented training data while ensuring theoretical identifiability. Extensive experiments on bnlearn benchmarks demonstrateTICL’s superiority in multiple aspects of causal discovery and intervention target detection. Correspondence: onedeanxxx@gmail.com, juding@microsoft.com Date: Oct 01, 2024 1 Introduction Causal discovery, the identification of causal relations from data, underpins modern scientific progress (Re- ichenbach, 1956). Recently, Supervised Causal Learning (SCL) has emerged as a promising modeling paradigm for causal discovery (Lorch et al., 2022; Ke et al., 2023a,b; Wu et al., 2025). By training on synthetic datasets spanning diverse causal structures, SCL models aim to operationalize causal discovery as a structured prediction task where empirical accuracy under realistic conditions is the central objective. Despite these successes, existing SCL works has primarily focused on observational data. However, experimen- tation remains the gold standard for causality (Hill, 1952; Pearl, 2009). Interventional data, in particular, allows for the identification of more causal relations with fewer assumptions (Hauser and Bühlmann, 2012). Applying SCL under interventional settings is thus of considerable importance. Unfortunately, the SCL paradigm faces two main challenges when interventions are involved, particularly if the specific intervention actions are unknown, a common scenario in real world third-party experiments (Eaton and Murphy, 2007a): Challenge I: Versatility. Different from the observational setting (which enjoys a simple and uniform problem formulation), interventional settings may vary widely, including factors such as known / unknown targets, hard / soft interventions, and single / multiple-variable interventions. This diversity complicates learning versatile models, and existing SCL methods (Lorch et al., 2022; Ke et al., 2023b) are typically limited to specific settings (e.g., hard + known), limiting their applicability to diverse practical scenarios. Challenge I: Generalizability. As a supervised learning approach, SCL methods fundamentally face the core challenge of generalization – models trained on data following a specific distribution (e.g., one generated by a configured simulator) may perform poorly on real-world test data. The formulation diversity as mentioned above may further exacerbate the generalization issue in the interventional setting. To circumvent the generalization bottleneck, we investigate Test-Time Training (T) paradigm (Sun, 2023), an emerging that trains instance-specific models at test/inference time in an on-the-fly manner rather than seeking a single universally generalizable model. Thus, we explore three key questions: (1) When should 1 arXiv:2602.19131v1 [cs.LG] 22 Feb 2026 +58.0% +13.7% ℐ-CPDAG ℐ: Int. Targets ℐ: Int. Targets +20.0% +43.0% +60.0%+63.0% +10.5% +6.3% +8.0% ℐ-CPDAG Figure 1 (Empirical dominance across interventional SCL tasks).TICLconsistently outperforms all SoTA methods on bothI-CPDAG discovery and intervention targets detection, across diverse intervention families and limited-sample regimes. ! Test Data ! ! ,# ! ,...,! " ,# " Training set ! ! ,! # ,...,! " Candidates ℳ SCL Model IS-MCMC Self-augmentation 3 Te s t-Time Tr a i n i n g 4 Inference 5 2 1 " # ∗ ⭐ output Obtain candidate graphs set to represent $(#|!) Forward Sampling Each # " fits parameter from ! and forward-samples the ! " 1 3 2 2,3 1 ! ℐ " 1 3 2 2,3 1 ! ℐ # 1 3 2 2,3 1 ! ℐ $ Space of DAGsMEC ℐ-MEC " ∗ Identifiability Figure 2 Left: The workflow of test-time learning of causal structure. Right: An identifiability example. training data be acquired? (2) What kind of training data is effective for SCL? (3) How can such data be acquired in our setting? Specifically, When: We identify and exploit the window between accessing test data and performing the actual inference. In this window of timing, we generate free training data via test data and train the model specifically for the final testing. What: We observe that the posterior estimation of causal graphs,P(G|D), provides highly effective training data. By sampling graphs (G 1 ,...,G n ) from the posterior and generating compatible datasets via forward sampling, we create paired instances ⟨G i ,D i ⟩ for self-augmentation (Fig. 2 Left). How: We introduce the IS-MCMC algorithm, which constructs a Markov chain over the augmented graph structure space with intervention constraints and performs multiple optimizations. This enables efficient sampling from the posterior P (G|D) to fuel the self-augmentation process. To address the versatility challenge, we advocate for the adoption of Joint Causal Inference (JCI) frame- work (Mooij et al., 2020). JCI unifies diverse interventional settings by pooling datasets into an augmented representation, allowing algorithms to treat them as observational. While JCI simplifies the input formulation, adapting SCL to this framework still requires to answer two critical questions: (1) What should be the appropriate learning target? (2) How should the learning process be designed? Specifically, What: We prioritize theoretical identifiability, emphasizing that models should predict identifiable causal structures, i.e.,I-CPDAG. This ensures the model only predicts relations that are structurally identifiable from the training data. How: We introduce the two-phase SCL algorithm. Inspired by the PC algorithm (Spirtes and Glymour, 1991), we focus on identifying the identifiable components of theI-CPDAG, namely the skeleton and v- structures. This approach ensures asymptotic correctness while promoting systematic feature characterization and improved classification mechanisms. 2 By incorporating these ideas together, we proposeTICL, a novel method for interventional causal discovery, focusing on discrete data to illustrate its effectiveness. Specifically,TICLestablishes a T + JCI paradigm, employing a self-augmentation algorithm for training data acquisition and a two-phase SCL algorithm for designing training targets. Evaluation Highlights: Figure 1 summarizes our main empirical finding:TICLestablishes a new state-of- the-art for interventional SCL. These improvements are hold consistently across diverse benchmark causal graphs and multiple evaluation criteria. In summary, our key contributions are: •We proposeTICL, a novel method for SCL under interventions. OurTICLconsistently outperforms existing state-of-the-art methods in experiments on two tasks: causal discovery and intervention target detection. • We introduce a specific test-time training technique to the SCL domain, which (self-)augments the training data by efficient sampling causal graphs from the posterior distribution via an optimized IS-MCMC process. •Our systematic study highlight JCI as a promising direction for unifying interventional settings in SCL. Within this framework, our two-phase learning algorithm provides a concrete solution for what and how to learn. 2 Preliminaries 2.1 Interventional Causal Discovery A Causal Graphical Model (CGM)M=<G,P >overdrandom variables X :=X 1 ,...,X d comprises: (i) a directed acyclic graph (DAG)Gwith nodes corresponding to the variables X and edges encoding direct causal relations between them, and (i) a joint probability distributionP X that is Markov compatible withG, i.e., P X = Q d i=1 P (X i |pa(X i )), where pa(X i ) are the parents of X i . Given an unknown causal model <G,P >, the Observational Causal Discovery problem asks to infer about the causal graphGbased on an i.i.d. sample of the joint distributionP X . It is however well known that under the observational setting, in the general case, we can only identify a causal graph up to its Markov Equivalent Class (MEC), even if with the faithfulness assumption (thatX⊥ P Y|Z ⇒ X⊥ G Y|Z) and with the Markov compatibility condition , where the MEC of theGis the set of graphs with the same skeleton and v-structures with G (Verma and Pearl, 1990). The identifiability limit in the observational setting can be effectively mitigated by conditioning the observations upon interventions, i.e., actions that purposefully perturb the causal model. In general, an intervention can apply to a subset of variablesI ⊂X, called the intervention targets, which can either be a single variable or contain multiple variables. For each target variableX i inI, the intervention replaces the conditional distributionP(X i |pa(X i )) with a new one:P (I) (X i |pa(X i )). The joint distribution after the intervention thus becomesP (I) X = Q i/∈I P (X i |pa(X i )) Q j∈I P (I) (X j |pa(X j )). The intervention is called a hard (a.k.a. perfect, structural) intervention (Eaton and Murphy, 2007b) if it eliminates the intervention targets’ causal dependence on their parents entirely, i.e., ifP (I) (X i |pa(X i )) =P (I) (X i ),∀X i ∈ I; otherwise it is called a soft (a.k.a. imperfect, parametric) intervention (Tian and Pearl, 2001) as it maintains at least part of the original causal dependence. It is an unknown intervention if the target set I is unknown. The Interventional Causal Discovery problem asks to infer about the causal graphGbased on a collection of data samplesD=D 0 ,D 1 ,...,D K , each obtained under a different intervention (Cooper and Yoo, 1999; Hauser and Bühlmann, 2012; Brouillard et al., 2020). LetI:=I 0 ,I 1 ,...,I K be the intervention targets for each of the interventions here.Iis called the intervention family of the datasetD. It is often useful to obtain D 0 as a sample of the observational distributionP X without any actual intervention (Ke et al., 2023a), and in this case we denoteI 0 =∅. In this paper, we consider the situation that the interventions are unknown, and we need to infer both about the causal graphGand about the intervention familyIfrom the given data collection D. 3 2.2 Joint Causal Inference Framework Our method is based on the JCI framework (Mooij et al., 2020), which reduces an interventional causal discovery problem into an observational causal discovery problem over an augmented causal model. The basic idea is to treat a data sample under intervention as an observation under an imposed condition. The data-sample collectionDin the intervention setting can then be seen as a single data sample under a variety of observation conditions. More specifically, Augmented Data. Given data-sample collectionD=D 0 ,D 1 ,...,D K obtained under intervention family I:=I 0 ,I 1 ,...,I K , an augmented dataD I can be constructed by stacking theK+ 1 data samples inD together by rows, then appendingKcolumns, each corresponding to a newly added environment variable (or intervention variable)X I k which takes binary value andX I k = 1 in and only in data points (= rows) originally from D k . Please refer to the data table at the top of Figure. 3 for an example. Augmented Graph. Accordingly, an augmented causal graphG I can be constructed by adding nodesX I k , and adding edgesX I k → X i for allX i ∈ I k if the intervention targetsI k are known. In theG I , nodes corresponding to the original variablesX i ∈X are called system nodes, and those corresponding toX I k are called environment nodes. In the JCI framework, we first convert data collectionDto augmented dataD I , optionally incorporating edges from known interventions regimeI, then infer the augmented graphG I , which encodes both the information of the original causal graphGover the system variables X, and reveals unknown intervention targets as edges from environment nodes to system nodes are identified. Moreover, the causal inference over the augmented graph in the JCI framework may be facilitated by some a priori constraints when prior knowledge / assumptions about the interventions are available. For example, exogeneity assumption requires no causal edges from any system variableX i ∈X to environment nodes. Similarly, complete randomized context and generic context assumption assume no confounding between system and environment nodes, and no causation between environment nodes, respectively. For more details, please refer to the paper (Mooij et al., 2020). 2.3 Identifiability with Intervention Data Two causal DAGsG 1 andG 2 are indistinguishable under an intervention familyIwithI 0 =∅if and only if their interventional graphsG 1 I k andG 2 I k have the same skeleton and v-structures for allI k ∈I(Yang et al., 2018). It is thus clear that with interventional data, we may identify the true causal graph up to a smaller equivalence class – called theI-MEC. Fig. 2 (Right) gives a simple example that illustrates how interventions can help with causal identifiability: The three DAGs on the left belong to the same MEC. Considering the intervention familyI=∅,1,2,3,G 1 I is not in the sameI-MEC asG 2 I andG 3 I due to the absence of the v-structureX 2 → X 1 ←1. The schematic on the right further elucidates the relationships among the true DAG G ⋆ , the I-MEC, the MEC, and the space of all possible DAGs. 3 Test-Time Learning of Causal Structure In this section, we present the details ofTICLmethod (see Figure 3). We first formally summarize the problem and outline its solution, then delve into the two essential components of ourTICL: self-augmentation strategy and two-phase supervised causal structure learning. Problem Summary. Given a data collectionDgenerated by an unknown causal model (G,P) under an intervention familyI(potentially multi-variable, soft, and unknown), we adopt standard causal discovery assumptions: i)Pis Markovian and faithful w.r.t the causal graph and causal sufficiency, and i) exogeneity, complete-randomized-context, and generic-context property for the intervention family. We aim to predict all causal relations entailed by the given dataD(subject to the above assumptions), which correspond to the invariant causal structures in theI-MEC set of the causal graphGbehindD, as explained in Section 2.3. Such invariant causal structures can be computationally encapsulated as a partial DAG, called the Interventional-Complete Partial Directed Acyclic Graph (I-CPDAG), in which each directed edge indicates an 4 Intervention Data Observation Data Hard Soft Single Multiple Know Unknow A C B 퐼 ! 퐼 " ABC퐼 ! 퐼 " 퐷 # 00 퐷 ! 10 퐷 " 01 Augmented Graph Augmented Data JCI-Style Step 1: JCIPooling Step 2: Self-Augmentation Step 3: Two-PhaseSCL forward sampling Good Initial State ⭐ Random Seed Proxy Seed Efficient Parallel Chains ⭐ 풢 (%&!) 풢 (%) 훩 (%&!) 훩 (%) IS-MCMC Intervention Constraints CPT Reuse Training Pair 풟 (%) 풢 # All Edges Feature train Skeleton Model Orientation Model All UTs Feature ✅ ❌ Intervention Target Detection JCI Prior train Meek Rules 퓘 - CPDAG Discovery 6⃣ 1⃣ 2⃣ 3⃣ 4⃣ 5⃣ ... Figure 3 The overall workflow of TICL for test-time learning of causal structure from interventional data. invariant causal relation in theI-MEC set. Besides theI-CPDAG Discovery task thus discussed, we also want to identify the unknown intervention targets in the intervention family I (Intervention Target Detection). Solution Outline.TICLproceeds in three steps: (1) JCI Pooling: First, we convert the given data collection Dinto an augmented dataD I following the JCI protocol. (2) Self-Augmentation: Then, we use the self- augmentation strategy to construct a Markov chain over the space of augmented graph structures constrained by interventions, where eachG i fits the parameters inD I , and forward sampling to getD i as the training instance. (3) Two-Phase SCL: Last, we extract edge and triplet features of training data to enable two-phase SCL for the skeleton and orientation models. The former leverages JCI priors (i.e., prior knowledge of environment-system variable relations) to identify unknown targets, while the latter employs Meek Rules to enhance causal discovery. 3.1 Training Data Acquisition via Self-Augmentation We propose leveraging the posterior estimation of causal graphs,P(G I |D I ) to generate training data. Specif- ically, we sample causal graphs (G 1 ,...,G n ) fromP(G I |D I ) using a tailored Markov Chain Monte Carlo (MCMC) designed for interventional data. For eachG i , the parameters governing the conditional distributions of variables, given their parents inG i , are re-estimated using the datasetD. Forward sampling is then applied with these parameters to generate a datasetD i compatible withG i These paired instances⟨G i ,D i ⟩are used for subsequent two-phase SCL. This self-augmented data offers two key advantages: First, since the true causal graphG ∗ is unknown, the posterior estimation (G 1 ,...,G n ) provides a diverse range of plausible causal structures, capturing epistemic uncertainty. Second, by re-estimating parameters and forward sampling compatible datasetsD i from eachG i , we generate accurately labeled instances⟨G i ,D i ⟩, whereG i likely retains similarity to the augmented graph G I , which entails the similarity properties of theD i . More importantly, this generation mechanism is free, and we can generate as much training data as we want. Implementation Details. We summarize the key steps here: I. Initialization: The initial graph G (0) is randomly initialized with the augmented graph’s node number. I. Mutation: In each iteration, a candidate graphG cand is generated fromG (t−1) by mutating the graph structure (adding, deleting, or reversing an edge). I. Parameters Re-Estimation: For eachG cand , the conditional probability table (CPT) is re-estimated via maximum likelihood estimation (MLE), ensuring alignment with the augmented dataD I . If a node’s parents remain unchanged, the corresponding CPT is inherited from the previous graph. IV. Evaluation: The candidate graph G cand is evaluated using a goodness-of-fit score, such as log-likelihood, computed based on its alignment with the observed data. V. Acceptance: The acceptance probability ofG cand is computed as:α(G cand |G (t−1) ) =min 1, score(G cand ) score(G (t−1) ) , 5 G cand is accepted as G (t) with probability α, else it retains G (t−1) . VI. Forward Sampling: OnceG cand is accepted, the correspondingD (t) is generated by forward sampling from the graph, and the pair ⟨G (t) ,D (t) ⟩ is stored for future use. There are, however, several critical considerations to address. First, since the IS-MCMC process operates on augmented graphs rather than standard causal graph, additional constraints must be imposed to ensure validity. Furthermore, managing time complexity is crucial for efficient sampling. Accordingly, we implement the following optimizations: •Good Initial State: Using a proxy algorithm to generate an initial graph significantly reduces convergence time compared to starting from a purely random graph. •Intervention Constraints: The mutation process must respect system-environment (sys – env node) con- straints that are specific to the interventional setup. Specifically, edges from system to environment nodes (sys→env) are not permitted, and interactions between environment nodes are excluded. This ensures that the augmented graphs remain consistent under the JCI framework. • Parameters Reuse: Only the nodes affected by structural changes need to have their CPTs re-estimated. The rest of the graph can retain its parameters from the previous. • Efficient Parallel Chains: Running multiple parallel IS-MCMC chains allows for faster exploration of the augmented graph space, improving the sampling efficiency. Remark. By addressing these factors, we ensure that the generation of training data is both effective and efficient, thereby enhancing the task performance. We detail optimization strategies in Appendix B.1 and summarize the procedure in Algorithms 1 and 2. Theoretically, our IS-MCMC follows a standard structured MCMC framework (Madigan et al., 1995; Su and Borsuk, 2016; Kuipers and Moffa, 2017) and converges to the posterior distributionP(G I |D I ), thus avoiding convergence to incorrect graph structures. (See formal discussion and convergence visualization in Appendix C). Empirically, we also demonstrate its task effectiveness in Section 4.2 with experimental evidence. 3.2 Two-Phase Supervised Causal Learning With the generated training data, a straightforward idea is to train a model that directly predicts the causal graphs. However, we argue that this may not be the best choice, and that the model should instead predict the identifiable causal structures, i.e., theI-CPDAG. Recall that a directed causal edge is identifiable if and only if it appears consistently in every causal graph compatible with the given data. In other words, there does not exist stable association at all between the given data and the unidentifiable causal edges, despite their presence in the true causal graph. For this reason, we choose to focus on identifying the identifiable components in the I-CPDAG, namely, skeletons and v-structures. Inspired by the PC algorithm (Spirtes and Glymour, 1991), we implement the supervised causal learning as a two-phase process, with phase 1 focusing on identifying the skeleton, followed by orientation predictions in phase 2. This shift implies that for each phase, we must define the learning target, the feature set, and the classification mechanism. Implementation Details. We first review the PC algorithm, which consists of two main stages: Phase 1 identifies the skeleton and the separating setsSS, determining the existence of edges. Starting from an undirected complete graph, edges are removed iteratively through conditional independence (CI) tests. Specifically, an edgeX i − X j is removed ifX i is conditionally independent ofX j given a subsetSof other variables in the currentk-order graph. Phase 2 orients the unshielded triplesUin the skeleton based on the S, assigning edge directionality. Here, the triple⟨X a ,X c ,X b ⟩is oriented into a v-structureX a → X c ← X b if X c is not in the separating set of X a and X b . Then, we establish a formal connection between PC and SCL. Due to limited space, we focus on phase 2 as an example. Phase 1 and more details are provided in B.2. Task: For all unshielded triples U, classify whether each triple ⟨X a ,X c ,X b ⟩ is a v-structure. Featurization:QuerySSsatisfyingX a ⊥ X b |S, and calculate existence feature:F ⟨X a ,X c ,X b ⟩ = 1 X c ∈S 0 X c /∈S 6 Classifier:Train a binary classifier using features: C o (F ⟨X a ,X c ,X b ⟩ ) := v-structure F ⟨X a ,X c ,X b ⟩ = 0 non-v-structure F ⟨X a ,X c ,X b ⟩ ̸= 0 In summary, detecting v-structures can be framed as a binary classification task, with the PC algorithm viewed as feature engineering combined with a static classifier. It’s clear that this approach can be extended by enriching the feature set, such as considering all possible separating sets or more conditional dependency information. Crucially, it enables us to replace heuristic searches and potentially erroneous CI tests with robust classification mechanisms trained on data. Remark. Building on this insight, we introduce two classifiers: one to detect the existence of edges between nodes, and the other to determine whether an unshielded triple is a v-structure. We detail feature engineering in Appendix B.2 and summarize the learnable components in Algorithms 4 and 5. Moreover, the interventional setting (augmented graph) offers further benefits. During inference, prior knowledge from the augmented graph allow pre-identifying certain edges (e.g., env→sys node), facilitating Meek rules (Meek, 1995) to orient additional edges. Theoretically, our method preserves the PC algorithm’s asymptotic properties: as sample size grows, the classifiers converges to a theoretically plausible solution, i.e., correct CPDAG (orI-CPDAG in intervention), thus ensuring identifiability (Proof in Appendix D). Empirically, we also demonstrate its effectiveness and improvements in Sections 4.1 and 4.4. 4 Experiments In this section, we conduct extensive experiments to investigate the following core research questions (RQs): RQ1: How does theTICLperform in causal structure identification from interventional data? (Effectiveness). RQ2: How the quality and quantity of training data from self-augmentation affectTICL? (Generalizability). RQ3: How does the sampling and complete running time ofTICLcompare to other methods? (Efficiency). RQ4: How does theTICLperform compared to others across different intervention settings? (Versatility). Benchmarks & Baselines. We use 14 (semi-)real causal graph datasets from the bnlearn (Scutari, 2010) repository as benchmarks. We compareTICLwith methods designed for intervention data, including score-based methods: GIES (Hauser and Bühlmann, 2012), IGSP (Wang et al., 2017; Yang et al., 2018), UT-IGSP (Squires et al., 2020), and continuous optimization-based methods: ENCO (Lippe et al., 2022), AVICI (Lorch et al., 2022). We also extend observational data methods using the JCI framework to intervention settings, including constraint-based methods like PC (Spirtes and Glymour, 1991), score-based methods like HC (Tsamardinos et al., 2006), BLIP (Scanagatta et al., 2015), and gradient-based methods like GOLEM (Ng et al., 2020). Note that some methods require known intervention targets, marked with *; DCDI (Brouillard et al., 2020) and BaCaDI (Hägele et al., 2023) failed due to their focus on continuous data, and results from CSIvA (Ke et al., 2023b) and SDI (Ke et al., 2023a) are omitted due to inaccessible code. More details are provided in Appendix E.1 and E.2. Training Datasets & Metrics. For the training data, by default we use 400 synthetic causal graph instances sampled via the IS-MCMC and forward sample 10k instances from re-parameterized conditional probability tables. XGBoost (Chen and Guestrin, 2016) is used as the classifier. The F1-Score evaluates intervention target detection andI-CPDAGdiscovery. For causal assessment, we also use Structural Hamming Distance (SHD) and Structural Intervention Distance (SID) (Peters and Bühlmann, 2015). More details in Appendix E.3 and E.4. Moreover, to enable unified and fair comparison, we implemented extensive engineering optimizations (e.g., specialized data structures, intervention mechanisms). The resulting codebase (>10k lines of Python) will be open-sourced to support future research. 4.1 Causal Structure Identification Performance (RQ1) Following DCDI and GIES, we conduct multiple intervention experiments for each graph, with the number of experiments being 20% of the nodes. To enhance baseline diversity, we use single-node interventions (e.g., ENCO limited to single-node interventions). Since hard interventions are a special case of soft ones, we choose soft interventions. Finally, we forward-sample 10k samples for both observational and intervention cases to generate test data. 7 Table 1 Performance comparison ofI-CPDAG(SHD↓/ SID↓/ F1↑).: Best,: Next. /: Crashes, –: Timeout. Methods TICLJCI-PCJCI-BLIPJCI-HCJCI-GOLEMAVICI ∗ ENCO ∗ IGSP ∗ UT-IGSPGIES SHD SID F1 SHD SID F1 SHD SID F1 SHD SID F1 SHD SID F1 SHD SID F1 SHD SID F1 SHD SID F1 SHD SID F1 SHD SID F1 Earthquake001.00001.00001.00 5 14 0.40 6 19 0.0036 0.403 12 0.67001.00001.00 100.89 Survey001.00 5 18 0.40140.91001.00 7 23 0.40 8 18 0.22 5 21 0.29140.91140.91140.91 Asia2130.86 4 19 0.753180.77 6 28 0.46 15 36 0.00 7 28 0.31 613 0.74 4 17 0.62 5 24 0.73 4 140.77 Sachs360.90 19 56 0.52 7 41 0.675 290.74 18 61 0.11 37 46 0.22 41 42 0.26 1016 0.70 13 21 0.61 15 38 0.48 Child7670.81 40 299 0.6210 1880.78 16 211 0.65 41 330 0.21 / / / 110 251 0.06 29124 0.63 33 214 0.45 46 197 0.29 Insurance172950.78 66 537 0.5324 3600.68 32 342 0.60 63 687 0.23 / / / 172 505 0.14 80 455 0.34 82 442 0.31 87261 0.52 Water454700.53 – – –46 5380.49 48503 0.46 81 564 0.22 / / / 58 527 0.46 / / / / / / / / / Mildew292390.76 58 7660.7635 5230.51 41 5000.51 – – – / / / 262 820 0.10 130 482 0.22 123 430 0.25 11681 0.41 Alarm9960.91 67 4810.7618152 0.75 38 416 0.61 – – – / / / 236 612 0.11 65 247 0.43 70 274 0.37 88 219 0.39 Barley559480.63 127 1789 0.57399190.70 78 1328 0.39 – – – / / / 98 1059 0.55 / / / / / / – – – Hailfinder425090.66 157 1152 0.33435580.70 81 935 0.38 – – – / / / 98 730 0.61 232 870 0.14 235 910 0.14 66 973 0.00 Hepar242 12160.79 – – – 50 1228 0.72438550.81 – – – / / / 75 1693 0.57 70 1663 0.62 67 1606 0.64 811085 0.67 Win95pts1016110.58 – – –109 900 0.53 – – – – – – / / / 1358930.54 / / / / / / 112 992 0.00 Pathfinder187 9473 0.27 – – – 207 /0.30 – – – – – – / / /15655280.47 13457261 0.10 1348 9226 0.08 – – – Rank (SHD)1.14 ± 0.125.36 ± 3.232.14 ± 0.413.86 ± 4.41 7.43 ± 3.67 8.07 ± 2.78 6.36 ± 5.804.64 ± 2.665.21 ± 2.745.50 ± 3.11 Rank (SID)1.57 ± 0.826.43 ± 3.823.57 ± 2.534.64 ± 5.37 8.14 ± 3.69 7.64 ± 2.66 5.29 ± 7.923.93 ± 2.784.57 ± 2.673.43 ± 2.96 Rank (F1)1.36 ± 0.374.43 ± 4.242.29 ± 0.783.93 ± 4.49 7.71 ± 3.35 8.00 ± 3.14 5.93 ± 7.494.86 ± 2.845.21 ± 2.455.36 ± 2.52 Table 2 Performance comparison of Intervention Targets Detection (F1↑).: Best,: Next. /: Crashes, –: Timeout. Datasets Earthquake Survey Asia Sachs Child Insurance Water Mildew Alarm Barley Hailfinder Hepar2 Win95pts Pathfinder Rank (F1) UT-IGSP0.501.00 0.44 0.33 0.350.17/0.220.15/0.270.23/0.045.14±3.27 CITE1.001.00 0.67 0.360.670.32/0.380.56/0.670.61/0.603.07±2.21 PreDITEr1.001.00 0.67 0.50–/–/–//5.21±5.17 JCI-GOLEM0.670.40 0.40 0.22 0.180.110.23–6.29±2.20 JCI-HC1.000.50 0.67 0.57 0.400.290.330.340.290.260.230.23–4.14±1.98 JCI-BLIP1.001.000.800.67 0.570.450.440.350.380.340.340.820.290.232.43±0.67 JCI-PC1.001.000.80 0.57 0.280.50–0.260.640.700.25–3.43±3.10 TICL1.001.001.001.001.000.830.860.580.820.830.760.900.500.741.00±0.00 As shown in Tables 1 and 2,TICLachieves competitive results across 14 causal graph datasets. Its advantages are highlighted in two aspects:❶Diverse and Challenging Benchmarking: Unlike methods tested only on synthetic datasets with 10-30 nodes (Brouillard et al., 2020; Lorch et al., 2022; Hägele et al., 2023; Ke et al., 2023b), we evaluate on the bnlearn benchmark, which includes real-world-inspired causal graphs ranging from small to over 100 nodes. Most methods fail as DAGs grow super exponentially, whileTICL demonstrates superior scalability and consistent performance.❷Highly Competitive Results: In intervention target detection,TICLoutperforms all methods, with an average F1-score improvement of 50.21% over the second-best method. ForI-CPDAGdiscovery,TICLimproves the F1 score by 13.62% over the best baseline. The exception is ENCO on Pathfinder, where prior knowledge of known intervention targets yields comparable performance. Overall,TICLleads with an average F1 rank of 1.36 ± 0.37 for causal discovery and 1.00 ± 0.00 for intervention target detection.❸Beyond this, we find methods based on the JCI framework also show competitive performance, suggesting it is a promising direction for intervention causal discovery. 4.2 Training Data Study: Quantity and Quality (RQ2) We maintain consistent setups and assess the impact of self-augmented data quality and quantity for SCL. Data Quality. To evaluate training data quality, we consider three generation strategies: Purely Random, generating random graphs (Erdős-Rényi, Scale-Free) with 10-50 nodes; IS-MCMC w/ Random Seed, using 8 F1@-CPDAG F1 @Target 0.50 0.75 1.00 Asia Purely Random w/ Random Seed w/ Proxy Seed F1@-CPDAG F1 @Target 0.50 0.75 1.00 Sachs 02468 SHD 10 14 18 22 26 SID Asia Purely Random w/ Random Seed w/ Proxy Seed 036912 SHD 0 5 10 15 20 SID Sachs Figure 4 The impact of different synthetic data strategies on the performance of causal structure learning. 15101520 0.50 0.75 1.00 F1 @ -CPDAG Random Proxy 15101520 0.50 0.75 1.00 F1 @ Target Random Proxy 15101520 0 3 6 SHD Random Proxy 15101520 0 5 10 SID Random Proxy 15101520 Step 0.50 0.75 1.00 F1 @ -CPDAG Random Proxy 15101520 Step 0.50 0.75 1.00 F1 @ Target Random Proxy 15101520 Step 0 7 14 SHD Random Proxy 15101520 Step 0 12 24 SID Random Proxy Earthquake Sachs Figure 5 The impact of different initialization strategies on the performance of IS-MCMC convergence. random graphs as initial seeds for IS-MCMC; and IS-MCMC w/ Proxy Seed, using proxy algorithms (e.g., JCI- BLIP) to generate seed graphs for IS-MCMC. We compare performance in Figure 4 and IS-MCMC convergence in Figure 5. We observe that:❶High-quality in-domain data is crucial for supervised causal models. The three generation strategies can be viewed as progressively approach the test data domain. As shown in Figure 4, both IS-MCMC variants outperform the purely random baseline across all metrics. Among them, the proxy seed method yields the best results, highlighting the value of self-augmentation strategy for test-time training. ❷Optimizing posterior-induced data generation process is key to self-augmentation. Although both IS-MCMC strategies target the same posterior, empirical results (Figure 5) show that proxy seeds converge faster and excel in both causal discovery and target detection. This confirms that our optimization strategy, via effective initialization and intervention constraints, ensures superior efficiency and convergence. 1002004008001600 # Training Instance 0 1 2 Distance SHD SID 0.8 0.9 1.0 Score Earthquake F1@-CPDAG F1 @Target 1002004008001600 # Training Instance 0 1 2 Distance SHD SID 0.8 0.9 1.0 Score Survey F1@-CPDAG F1 @Target 1002004008001600 # Training Instance 0 15 30 Distance SHD SID 0.50 0.75 1.00 Score Asia F1@-CPDAG F1 @Target 1002004008001600 # Training Instance 0 5 10 Distance SHD SID 0.6 0.8 1.0 Score Sachs F1@-CPDAG F1 @Target Figure 6 Effect of training data instances on the performance of causal structure learning. Data Quantity. We evaluateTICL’s scaling property with training data size: instance count and sample size per instance. Figure 6 shows the impact of instance counts (100-1600), while Figure 7 plots performance with sample sizes (1k–20k) under default intervention settings. We observe that: ❶ Our TICL benefits from more training instances. As shown in Figure 6, performance improves with more instances across all datasets. Smaller graphs require fewer instances, while larger ones gain more due to increased complexity.❷ TICLbenefits from 9 1k2k5k10k20k # Sample Size 0 1 2 Distance SHD SID 0.8 0.9 1.0 Score Earthquake F1@-CPDAG F1 @Target 1k2k5k10k20k # Sample Size 0 7 14 Distance SHD SID 0.6 0.8 1.0 Score Survey F1@-CPDAG F1 @Target 1k2k5k10k20k # Sample Size 0 15 30 Distance SHD SID 0.50 0.75 1.00 Score Asia F1@-CPDAG F1 @Target 1k2k5k10k20k # Sample Size 0 15 30 Distance SHD SID 0.4 0.7 1.0 Score Sachs F1@-CPDAG F1 @Target Figure 7 Effect of training sample size on the performance of causal structure learning. larger sample sizes. As shown in Figure 7, increasing sample size stabilizes independence tests, improving various metrics across graphs. These properties further highlight the promise, as our self-augmentation enables the generation of abundant free data. 4.3 Sampling and Running Efficiency Study (RQ3) Figure 8 Efficiency study of our method. We maintain the basic setup and assess the efficiency of our method from two perspectives. Sampling Time. Since our initialization en- sures rapid convergence, we evaluate the multi- chain MCMC strategy. As shown in Figure 8 (a), multi-chain MCMC shows minimal time increase as graphs scale, offering significant acceleration. This stems from IS-MCMC’s parallel adaptability, confirming the strategy’s efficiency. Running Time. We compare the wall-clock times (including proxy time) of all baseline methods. As shown in Figure 8 (b):❶Other learning-based methods are inefficient, typically requiring 100-500 seconds due to slow gradient optimization, while non-SCL methods finish within 50 seconds.❷ TICLmatches non-SCL efficiency. While execution time involves proxy algorithm, this is optional; in high-risk causal discovery, time cost becomes secondary. 4.4 Different Intervention Settings Study (RQ4) We further conducted experiments and analyses on RQ4 regarding intervention types, intervention ratios, and test sample sizes. Due to the complex diversity of different situations, this has been rarely considered comprehensively before. However, as shown in Figure 1,TICLalways maintains its superiority under various intervention combinations and a limited number of samples. Due to page limitation, for more results and analysis, please refer to Appendix F. 5 Conclusion and Future work In this paper, we presentTICLfor causal structural learning from interventional data, particularly when faced with generalization challenges in the real-world. We introduced a specific test-time training technique and demonstrated the prospects of JCI for SCL. In the future, we further explore the performance ofTICLin continuous setting and under fewer assumptions to pursue a causal foundation model. 10 References Alan Agresti.Categoricaldataanalysis, volume 792. John Wiley & Sons, 2012. Yashas Annadani, Nick Pawlowski, Joel Jennings, Stefan Bauer, Cheng Zhang, and Wenbo Gong. Bayesdag: Gradient- based posterior inference for causal discovery.AdvancesinNeuralInformationProcessingSystems, 36:1738–1763, 2023. Philippe Brouillard, Sébastien Lachapelle, Alexandre Lacoste, Simon Lacoste-Julien, and Alexandre Drouin. Dif- ferentiable causal discovery from interventional data.AdvancesinNeuralInformationProcessingSystems, 33: 21865–21877, 2020. Tianqi Chen and Carlos Guestrin. Xgboost: A scalable tree boosting system. InProceedingsofthe22ndacmsigkdd internationalconferenceonknowledgediscoveryanddatamining, pages 785–794, 2016. Wei Chen and Yuxuan Liang. Learning with calibration: Exploring test-time computing of spatio-temporal forecasting. InTheThirty-ninthAnnualConferenceonNeuralInformationProcessingSystems, 2025. J Cheng, R Greiner, J Kelly, D Bell, and W Liu. Learning bayesian networks from data: An information-theory based approach department of computing sciences.Artificialintelligence, 137:43–90, 2002. David Maxwell Chickering. Learning equivalence classes of bayesian-network structures.TheJournalofMachine LearningResearch, 2:445–498, 2002a. David Maxwell Chickering. Optimal structure identification with greedy search.Journalofmachinelearningresearch, 3(Nov):507–554, 2002b. Junsouk Choi, Robert Chapkin, and Yang Ni. Bayesian causal structural learning with zero-inflated poisson bayesian networks.Advancesinneuralinformationprocessingsystems, 33:5887–5897, 2020. Diego Colombo, Marloes H Maathuis, et al. Order-independent constraint-based causal structure learning.J.Mach. Learn.Res., 15(1):3741–3782, 2014. Pierre Comon. Independent component analysis, a new concept? Signalprocessing, 36(3):287–314, 1994. Gregory F Cooper and Changwon Yoo. Causal discovery from a mixture of experimental and observational data. In ProceedingsoftheFifteenthconferenceonUncertaintyinartificialintelligence, pages 116–125, 1999. Thomas M Cover.Elementsofinformationtheory. John Wiley & Sons, 1999. James Cussens. Bayesian network learning with cutting planes.arXivpreprintarXiv:1202.3713, 2012. Haoyue Dai, Rui Ding, Yuanyuan Jiang, Shi Han, and Dongmei Zhang. Ml4c: Seeing causality through latent vicinity. InProceedingsofthe2023SIAMInternationalConferenceonDataMining(SDM), pages 226–234. SIAM, 2023. Karan Dalal, Daniel Koceja, Jiarui Xu, Yue Zhao, Shihao Han, Ka Chun Cheung, Jan Kautz, Yejin Choi, Yu Sun, and Xiaolong Wang. One-minute video generation with test-time training. InProceedingsoftheComputerVisionand PatternRecognitionConference, pages 17702–17711, 2025. Tristan Deleu, António Góis, Chris Emezue, Mansi Rankawat, Simon Lacoste-Julien, Stefan Bauer, and Yoshua Bengio. Bayesian structure learning with generative flow networks. InUncertaintyinArtificialIntelligence, pages 518–528. PMLR, 2022. Payam Dibaeinia and Saurabh Sinha. Sergio: a single-cell expression simulator guided by gene regulatory networks. Cellsystems, 11(3):252–271, 2020. Rui Ding, Yanzhi Liu, Jingjing Tian, Zhouyu Fu, Shi Han, and Dongmei Zhang. Reliable and efficient anytime skeleton learning. InProceedingsoftheAAAIConferenceonArtificialIntelligence, volume 34, pages 10101–10109, 2020. Daniel Eaton and Kevin Murphy. Belief net structure learning from uncertain interventions.JMachLearnRes, 1: 1–48, 2007a. Daniel Eaton and Kevin Murphy. Exact bayesian structure learning from uncertain interventions. InArtificial intelligenceandstatistics, pages 107–114. PMLR, 2007b. Byron Ellis and Wing Hung Wong. Learning causal bayesian network structures from experimental data.Journalof theAmericanStatisticalAssociation, 103(482):778–789, 2008. 11 Arthur Gretton, Kenji Fukumizu, Choon Teo, Le Song, Bernhard Schölkopf, and Alex Smola. A kernel statistical test of independence.Advancesinneuralinformationprocessingsystems, 20, 2007. Wiebke Günther, Urmi Ninad, Jonas Wahl, and Jakob Runge. Conditional independence testing with heteroskedastic data and applications to causal discovery.AdvancesinNeuralInformationProcessingSystems, 35:16191–16202, 2022. Alexander Hägele, Jonas Rothfuss, Lars Lorch, Vignesh Ram Somnath, Bernhard Schölkopf, and Andreas Krause. Bacadi: Bayesian causal discovery with unknown interventions. InInternationalConferenceonArtificialIntelligence andStatistics, pages 1411–1436. PMLR, 2023. Moritz Hardt and Yu Sun. Test-time training on nearest neighbors for large language models. InTheTwelfth InternationalConferenceonLearningRepresentations, 2024. Alain Hauser and Peter Bühlmann. Characterization and greedy learning of interventional markov equivalence classes of directed acyclic graphs.TheJournalofMachineLearningResearch, 13(1):2409–2464, 2012. A Bradford Hill. The clinical trial. NewEnglandJournalofMedicine, 247(4):113–119, 1952. Steven M Hill, Chris J Oates, Duncan A Blythe, and Sach Mukherjee. Causal learning via manifold regularization. 2019. Patrik Hoyer, Dominik Janzing, Joris M Mooij, Jonas Peters, and Bernhard Schölkopf. Nonlinear causal discovery with additive noise models.Advancesinneuralinformationprocessingsystems, 21, 2008. Markus Kalisch and Peter Bühlman. Estimating high-dimensional directed acyclic graphs with the pc-algorithm. JournalofMachineLearningResearch, 8(3), 2007. Nan Rosemary Ke, Olexa Bilaniuk, Anirudh Goyal, Stefan Bauer, Hugo Larochelle, Bernhard Schölkopf, Michael C Mozer, Chris Pal, and Yoshua Bengio. Learning neural causal models from unknown interventions.arXivpreprint arXiv:1910.01075, 2019. Nan Rosemary Ke, Olexa Bilaniuk, Anirudh Goyal, Stefan Bauer, Hugo Larochelle, Bernhard Schölkopf, Michael Curtis Mozer, Christopher Pal, and Yoshua Bengio. Neural causal structure discovery from interventions.Transactionson MachineLearningResearch, 2023a. Nan Rosemary Ke, Silvia Chiappa, Jane Wang, Anirudh Goyal, Jorg Bornschein, Melanie Rey, Theophane Weber, Matthew Botvinic, Michael Mozer, and Danilo Jimenez Rezende. Learning to induce causal structure.TheEleventh InternationalConferenceonLearningRepresentations, 2023b. Daphne Koller and Nir Friedman.Probabilisticgraphicalmodels:principlesandtechniques. MIT press, 2009. Jack Kuipers and Giusi Moffa. Partition mcmc for inference on acyclic digraphs.JournaloftheAmericanStatistical Association, 112(517):282–299, 2017. Adam Li, Amin Jaber, and Elias Bareinboim. Causal discovery from observational and interventional data across multiple environments. InThirty-seventhConferenceonNeuralInformationProcessingSystems, 2023. Hebi Li, Qi Xiao, and Jin Tian. Supervised whole dag causal discovery.arXivpreprintarXiv:2006.04697, 2020. Wenqian Li, Yinchuan Li, Shengyu Zhu, Yunfeng Shao, Jianye Hao, and Yan Pang. Gflowcausal: Generative flow networks for causal discovery.arXivpreprintarXiv:2210.08185, 2022. Phillip Lippe, Taco Cohen, and Efstratios Gavves. Efficient neural causal discovery without acyclicity constraints.The TenthInternationalConferenceonLearningRepresentations, 2022. Yuejiang Liu, Parth Kothari, Bastien Van Delft, Baptiste Bellot-Gurlet, Taylor Mordan, and Alexandre Alahi. Ttt++: When does self-supervised test-time training fail or thrive?AdvancesinNeuralInformationProcessingSystems, 34:21808–21820, 2021. David Lopez-Paz, Krikamol Muandet, and Benjamin Recht. The randomized causation coefficient.J.Mach.Learn. Res., 16:2901–2907, 2015. Lars Lorch, Scott Sussex, Jonas Rothfuss, Andreas Krause, and Bernhard Schölkopf. Amortized inference for causal structure learning.AdvancesinNeuralInformationProcessingSystems, 35, 2022. Pingchuan Ma, Rui Ding, Haoyue Dai, Yuanyuan Jiang, Shuai Wang, Shi Han, and Dongmei Zhang. Ml4s: Learning causal skeleton from vicinal graphs. InProceedingsofthe28thACMSIGKDDConferenceonKnowledgeDiscovery andDataMining, pages 1213–1223, 2022. 12 David Madigan, Jeremy York, and Denis Allard. Bayesian graphical models for discrete data.InternationalStatistical Review/RevueInternationaledeStatistique, pages 215–232, 1995. Hakim Manghwar, Bo Li, Xiao Ding, Amjad Hussain, Keith Lindsey, Xianlong Zhang, and Shuangxia Jin. Crispr/cas systems in genome editing: methodologies and tools for sgrna design, off-target evaluation, and strategies to mitigate off-target effects.Advancedscience, 7(6):1902312, 2020. Alessandro Mascaro and Federico Castelletti. Bayesian causal discovery from unknown general interventions.arXiv preprintarXiv:2312.00509, 2023. Christopher Meek. Causal inference and causal explanation with background knowledge. InProceedingsoftheEleventh conferenceonUncertaintyinartificialintelligence, pages 403–410, 1995. Thomas Minka. Estimating a dirichlet distribution, 2000. Joris M Mooij, Sara Magliacane, and Tom Claassen. Joint causal inference from multiple contexts.TheJournalof MachineLearningResearch, 21(1):3919–4026, 2020. In Jae Myung. Tutorial on maximum likelihood estimation.JournalofmathematicalPsychology, 47(1):90–100, 2003. Ignavier Ng, Shengyu Zhu, Zhitang Chen, and Zhuangyan Fang. A graph autoencoder approach to causal structure learning.arXivpreprintarXiv:1911.07420, 2019. Ignavier Ng, AmirEmad Ghassami, and Kun Zhang. On the role of sparsity and dag constraints for learning linear dags.AdvancesinNeuralInformationProcessingSystems, 33:17943–17954, 2020. Judea Pearl.Causality. Cambridge university press, 2009. Jonas Peters and Peter Bühlmann. Structural intervention distance for evaluating causal graphs.Neuralcomputation, 27(3):771–799, 2015. Anne Helby Petersen, Joseph Ramsey, Claus Thorn Ekstrøm, and Peter Spirtes. Causal discovery for observational sciences using supervised machine learning.JournalofDataScience, 21(2):255–280, 2023. Joseph Ramsey, Jiji Zhang, and Peter L Spirtes. Adjacency-faithfulness and conservative causal inference.arXiv preprintarXiv:1206.6843, 2012. Hans Reichenbach. The direction of time.Philosophy, 34(128):65–66, 1956. Andreas Sauter, Nicolò Botteghi, Erman Acar, and Aske Plaat. Core: Towards scalable and efficient causal discovery with reinforcement learning. InProceedingsofthe23rdInternationalConferenceonAutonomousAgentsand MultiagentSystems, page 1664–1672. International Foundation for Autonomous Agents and Multiagent Systems, 2024. Mauro Scanagatta, Cassio P de Campos, Giorgio Corani, and Marco Zaffalon. Learning bayesian networks with thousands of variables.Advancesinneuralinformationprocessingsystems, 28, 2015. Marco Scutari. Learning bayesian networks with the bnlearn r package.JournalofStatisticalSoftware, 2010. Alex Smola, Arthur Gretton, Le Song, and Bernhard Schölkopf. A hilbert space embedding for distributions. In Internationalconferenceonalgorithmiclearningtheory, pages 13–31. Springer, 2007. Peter Spirtes and Clark Glymour. An algorithm for fast recovery of sparse causal graphs.Socialsciencecomputer review, 9(1):62–72, 1991. Peter Spirtes, Clark N Glymour, and Richard Scheines. Causation,prediction,andsearch. MIT press, 2000. Chandler Squires, Yuhao Wang, and Caroline Uhler. Permutation-based causal structure learning with unknown intervention targets. InConferenceonUncertaintyinArtificialIntelligence, pages 1039–1048. PMLR, 2020. Chengwei Su and Mark E Borsuk. Improving structure mcmc for bayesian networks through markov blanket resampling. TheJournalofMachineLearningResearch, 17(1):4042–4061, 2016. Yu Sun.Test-TimeTraining. PhD thesis, University of California, Berkeley, 2023. Yu Sun, Xiaolong Wang, Zhuang Liu, John Miller, Alexei Efros, and Moritz Hardt. Test-time training with self- supervision for generalization under distribution shifts. InInternationalconferenceonmachinelearning, pages 9229–9248. PMLR, 2020. 13 Yu Sun, Xinhao Li, Karan Dalal, Jiarui Xu, Arjun Vikram, Genghan Zhang, Yann Dubois, Xinlei Chen, Xiaolong Wang, Sanmi Koyejo, et al. Learning to (learn at test time): Rnns with expressive hidden states. InForty-second InternationalConferenceonMachineLearning, 2025. Arnuv Tandon, Karan Dalal, Xinhao Li, Daniel Koceja, Marcel Rød, Sam Buchanan, Xiaolong Wang, Jure Leskovec, Sanmi Koyejo, Tatsunori Hashimoto, et al. End-to-end test-time training for long context.arXivpreprint arXiv:2512.23675, 2025. Jin Tian and Judea Pearl. Causal discovery from changes. InProceedingsoftheSeventeenthconferenceonUncertainty inartificialintelligence, pages 512–521, 2001. Jean-François Ton, Dino Sejdinovic, and Kenji Fukumizu. Meta learning for causal direction. InProceedingsofthe AAAIconferenceonartificialintelligence, volume 35, pages 9897–9905, 2021. Ioannis Tsamardinos, Laura E Brown, and Constantin F Aliferis. The max-min hill-climbing bayesian network structure learning algorithm.Machinelearning, 65:31–78, 2006. Burak Varici, Karthikeyan Shanmugam, Prasanna Sattigeri, and Ali Tajer. Scalable intervention target estimation in linear models.AdvancesinNeuralInformationProcessingSystems, 34:1494–1505, 2021. Burak Varici, Karthikeyan Shanmugam, Prasanna Sattigeri, and Ali Tajer. Intervention target estimation in the presence of latent variables. InUncertaintyinArtificialIntelligence, pages 2013–2023. PMLR, 2022. Thomas Verma and Judea Pearl. Equivalence and synthesis of causal models. InProceedingsoftheSixthAnnual ConferenceonUncertaintyinArtificialIntelligence, pages 255–270, 1990. MK Vijaymeena and K Kavitha. A survey on similarity measures in text mining.MachineLearningandApplications: AnInternationalJournal, 3(2):19–28, 2016. Julius von Kügelgen, Michel Besserve, Liang Wendong, Luigi Gresele, Armin Kekić, Elias Bareinboim, David Blei, and Bernhard Schölkopf. Nonparametric identifiability of causal representations from unknown interventions.Advances inNeuralInformationProcessingSystems, 36, 2024. Dequan Wang, Evan Shelhamer, Shaoteng Liu, Bruno Olshausen, and Trevor Darrell. Tent: Fully test-time adaptation by entropy minimization.arXivpreprintarXiv:2006.10726, 2020. Qin Wang, Olga Fink, Luc Van Gool, and Dengxin Dai. Continual test-time domain adaptation. InProceedingsofthe IEEE/CVFConferenceonComputerVisionandPatternRecognition, pages 7201–7211, 2022. Xiaoqiang Wang, Yali Du, Shengyu Zhu, Liangjun Ke, Zhitang Chen, Jianye Hao, and Jun Wang. Ordering-based causal discovery with reinforcement learning.arXivpreprintarXiv:2105.06631, 2021. Yuhao Wang, Liam Solus, Karren Yang, and Caroline Uhler. Permutation-based causal inference algorithms with interventions.AdvancesinNeuralInformationProcessingSystems, 30, 2017. Menghua Wu, Yujia Bao, Regina Barzilay, and Tommi Jaakkola. Sample, estimate, aggregate: A recipe for causal discovery foundation models.TransactionsonMachineLearningResearch, 2025. Jing Xiang and Seyoung Kim. A* lasso for learning a sparse bayesian network structure for continuous variables. Advancesinneuralinformationprocessingsystems, 26, 2013. Karren Yang, Abigail Katcoff, and Caroline Uhler. Characterizing and learning equivalence classes of causal dags under interventions. InInternationalConferenceonMachineLearning, pages 5541–5550. PMLR, 2018. Yuqin Yang, Saber Salehkaleybar, and Negar Kiyavash. Learning unknown intervention targets in structural causal models from heterogeneous data. InInternationalConferenceonArtificialIntelligenceandStatistics, pages 3187– 3195. PMLR, 2024. Tjalling J Ypma. Historical development of the newton–raphson method. SIAMreview, 37(4):531–551, 1995. Kui Yu, Jiuyong Li, and Lin Liu. A review on algorithms for constraint-based causal discovery.arXivpreprint arXiv:1611.03977, 2016. Yue Yu, Jie Chen, Tian Gao, and Mo Yu. Dag-gnn: Dag structure learning with graph neural networks. InInternational ConferenceonMachineLearning, pages 7154–7163. PMLR, 2019. Alessio Zanga, Elif Ozkirimli, and Fabio Stella. A survey on causal discovery: Theory and practice.International JournalofApproximateReasoning, 151:101–129, 2022. 14 Jiaqi Zhang, Kristjan Greenewald, Chandler Squires, Akash Srivastava, Karthikeyan Shanmugam, and Caroline Uhler. Identifiability guarantees for causal disentanglement from soft interventions.AdvancesinNeuralInformation ProcessingSystems, 36, 2024. Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric P Xing. Dags with no tears: Continuous optimization for structure learning.Advancesinneuralinformationprocessingsystems, 31, 2018. 15 Appendix Table of Contents A Background Knowledge and Related Work18 A.1 Background Knowledge . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 A.1.1 Causal Graph-related Concept . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 A.1.2 JCI Assumption-related Concept . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 A.2 Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 A.2.1 Interventional Causal Discovery . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 A.2.2 Supervised Causal Learning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 A.2.3 Test-Time Training . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 A.2.4 Joint Causal Inference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 B Inplementation Details21 B.1 Training Data Acquisition via Self-Augmentation . . . . . . . . . . . . . . . . . . . . . . . . . 21 B.1.1 Parameter Re-Estimation for CPT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 B.1.2 Interventional Structural-Markov Chain Monte Carlo . . . . . . . . . . . . . . . . . . . 22 B.2 Two-phase Supervised Causal Learning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 B.2.1 Revisiting the PC Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 B.2.2 Connecting PC and Supervised Machine Learning . . . . . . . . . . . . . . . . . . . . 26 B.2.3 Skeleton and Orientation Feature Engineering . . . . . . . . . . . . . . . . . . . . . . . 26 C Convergence of Inteventional Structural-MCMC29 C.1 Guarantee of Posterior Distribution Convergence . . . . . . . . . . . . . . . . . . . . . . . . . 29 C.2 Optimization Considerations for IS-MCMC . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 C.3 Convergence Visualization Analysis of IS-MCMC . . . . . . . . . . . . . . . . . . . . . . . . . 30 D Identifiability Theory of Two-Phase SCL30 D.1 Intuitive Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 D.1.1 Observation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 D.1.2 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 D.2 Advantage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 D.3 Theoretical Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 D.3.1 Restating the PC-like SCL Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 D.3.2 Asymptotic Correctness of the PC-like SCL Approach . . . . . . . . . . . . . . . . . . 32 Go to Appendix Index16 E Experimental Details33 E.1 Benchmark . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 E.2 Baselines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 E.3 Metrics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 E.4 Experiments Setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 F Additional Experiments (RQ4)39 F.1 Effect of Intervention Type . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 F.1.1 Prefect Interventions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 F.1.2 Multiple Interventions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 F.2 Effect of Sample Size . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 F.2.1 Limited interventional data sample sizes . . . . . . . . . . . . . . . . . . . . . . . . . . 40 F.2.2 Limited observational data sample sizes . . . . . . . . . . . . . . . . . . . . . . . . . . 40 F.3 Effect of Intervention Experiment Size . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 G More Related Work41 H More Discussion / Limitation / Future Work45 H.1 Supervised Causal Discovery . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 H.2 Self-augmentation vs. Pre-training . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 H.3 Discrete vs. Continuous Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 I Broader Impacts45 Go to Appendix Index17 A Background Knowledge and Related Work A.1 Background Knowledge In this section, we give basic concepts related to causal graphs and joint causal inference. A.1.1 Causal Graph-related Concept Definition A.1 (Directed Acyclic Graph). A directed acyclic graph (DAG) is a directed graphGthat has no cycles, i.e. no directed paths starting and ending at the same vertex. Definition A.2 (Skeleton). A undirected graphK= (V K ,E K ) represents the skeleton of a causal DAG G = (V G ,E G ) if (X → Y )∈ E G ∪ (Y → X)∈ E G ⇐⇒ (X − Y )∈ E K . Definition A.3 (UT and V-structures). A triple of variables⟨X,T,Y⟩in a skeleton is an unshielded triple, or short for UT, ifXandYare adjacent toTbut are not adjacent to each other.⟨X,T,Y⟩can be further oriented to become a v-structure X → T ← Y , in which T is called the collider. Definition A.4 (PC). Denote the set of parents and children ofXin a skeleton asPC X , in other words, PC X are the neighbors ofXin the skeleton. For convenience, if we discussPC X in the context of a UT ⟨X,T,Y⟩, we intentionally mean the set of parents and children ofXbut excludeT. Similarly,PC T excludes X,Y . Definition A.5 (Vicinity). We define the vicinity of a UT⟨X,T,Y⟩asV ⟨X,T,Y⟩ :=X,T,Y∪ PC X ∪ PC Y ∪ PC T . Vicinity is a generalized version of PC, i.e., the neighbors of X,T,Y in the skeleton. Definition A.6 (Sepsets). Sepsets S can be define: S : X⊥Y|S,S ⊂ PC X ∪ T, or S ⊂ PC Y ∪ T. Under faithfulness assumption, sepsetsSis an ensemble where each item is a subset of variables within the vicinity that d-separates X and Y . Definition A.7 (Causal Graph). A causal graphGis a graphical description of a system in terms of cause-effect relationships, i.e. the causal mechanism. Specifically, for each edge (X,Y)∈E,Xis the direct cause ofY, andYis the direct effect ofX, satisfying the causal edge assumption, i.e., the value assigned to each variableXis completely determined by the functionFgiven its parent. Formally, this can be expressed as: X i : = f (Pa(X i )), ∀X i ∈V. As natural consequence of such definitions, we can define models that entail both the structural representation and the set of functions that regulate the underlying causal mechanism. Definition A.8 (Structural Causal Model). A structural causal model (SCM) is defined by the tuple M = (V,U,F,P ), where: • V is a set of endogenous variables, i.e. observable variables, • U is a set of exogenous variables, i.e. unobservable variables, where U ∪V =∅ • F is a set of functions, where each functionf i ∈ Fis defined asfi : = (V ∪U) p →V, withpthe ariety of f i , so that f i determines completely the value of V i , • P is a joint probability distribution over the exogenous variables P (U ) = Q i P (U i ). Assumption 1 (Causally Sufficiency). The set of variables V is said to be causally sufficient if and only if every cause of any subset of V is contained in V itself. Note that, in our setup, due to the causality sufficiency, exogenous variables are ignored. A.1.2 JCI Assumption-related Concept Assumption 0 ("Joint SCM"). The data-generating mechanism is described by a simple SCMMof the form: M : C k = f k (X pa H (k)∩I ,C pa H (k)∩K ,E pa H (k)∩J ),k ∈K, X i = f i (X pa H (i)∩I ,C pa H (i)∩K ,E pa H (i)∩J ),i∈I, P(E) = Q j∈J P(E j ), (1) Go to Appendix Index18 that jointly models the system and the context. Its graphG(M) has nodesI ∪K(corresponding to system variables X i i∈I and context variables C k k∈K ). It will always make this assumption in order to facilitate the formulation of JCI, the following three assumptions that we discuss are optional, and their applicability has to be decided based on a case-by-case basis. Typically, when a modeler decides to distinguish asystemfrom itscontext, the modeler possesses background knowledge that expresses that the context isexogenous to the system: Assumption 1. ("Exogeneity") No system variable causes any context variable, i.e., ∀k ∈K,∀i∈I : i→ k /∈G(M). Assumption 2. ("Complete randomized context") No context variable is confounded with a system variable, i.e., ∀k ∈K,∀i∈I : i↔ k /∈G(M). JCI Assumption 1 is often easily justifiable, but the applicability of JCI Assumption 2 may be less obvious in practice. Then, Assumption 3 is further stated, which can be useful whenever both JCI Assumptions 1 and 2 have been made as well. Assumption 3. ("Generic context model") The context graph 1 G(M) K is of the following special form: ∀k ̸= k ′ ∈K : k ↔ k ′ ∈G(M) ∧ k → k ′ /∈G(M). The following key result essentially states that when one is only interested in modeling the causal relations involving the system variables (under JCI Assumptions 1 and 2), one does not need to care about the causal relations between the context variables, as long as one correctly models the context distribution. Theorem 1. Assume that JCI Assumptions 0, 1 and 2 hold for SCM M: M : C k = f k (C pa H (k)∩K ,E pa H (k)∩J ),k ∈K, X i = f i (X pa H (i)∩I ,C pa H (i)∩K ,E pa H (i)∩J ),i∈I, P(E) = Q j∈J P(E j ), For any other SCM ̃ Msatisfying JCI Assumptions 0, 1 and 2 that is the same asMexcept that it models the context differently, i.e., of the form ̃ M : C k = ̃ f k (C pa ̃ H (k)∩K ,E pa ̃ H (k)∩ ̃ J ),k ∈K, X i = f i (X pa H (i)∩I ,C pa H (i)∩K ,E pa H (i)∩J ),i∈I, P(E) = Q j∈ ̃ J P(E j ), with J ⊆ ̃ J and pa H (i) = pa ̃ H (i) for all i∈I, we have that (i) the conditional system graphs coincide: G(M) do(K) =G( ̃ M) do(K) ; (i)if ̃ MandMinduce the same context distribution, i.e.,P M (C) =P ̃ M (C), then for any perfect intervention on the system variablesdo(I,ξ I ) withI ⊆ I(including the non-interventionI=∅), ̃ M do(I,ξ I ) is observationally equivalent to M do(I,ξ I ) . (i) if the context graphsG( ̃ M ) K andG(M) K induce the same separations, then alsoG( ̃ M ) andG(M) induce the same separations (where “separations” can refer to either d-separations or σ-separations). The following corollary of Theorem 1 states that JCI Assumption 3 can be made without loss of generality for the purposes of constraint-based causal discovery if the context distribution contains no conditional independences: 1 Remember that G(M) K denotes the subgraph on the context variables K induced by the causal graph G(M). Go to Appendix Index19 Corollary 2. Assume that JCI Assumptions 0, 1 and 2 hold for SCMM. Then there exists an SCM ̃ M that satisfies JCI Assumptions 0, 1 and 2 and 3, such that (i) the conditional system graphs coincide: G(M) do(K) =G( ̃ M) do(K) ; (i) for any perfect intervention on the system variablesdo(I,ξ I ) withI ⊆I(including the non-intervention I =∅), ̃ M do(I,ξ I ) is observationally equivalent to M do(I,ξ I ) ; (i)if the context distributionP M (C) contains no conditional or marginal independences, then the same σ-separations hold inG( ̃ M ) as inG(M); if in addition, the Directed Global Markov Property holds for M, then also the same d-separations hold in G( ̃ M) as in G(M). A.2 Related Work In this section, we present some of the most relevant work to the key points of this paper. A.2.1 Interventional Causal Discovery Recovering the underlying causal structure from observational and interventional data is a fundamental research problem (Spirtes et al., 2000; Pearl, 2009). When only observational data are available, constraint- based methods identify a directed acyclic graph (DAG) consistent with conditional independence constraints. Notable examples include the PC algorithm (Spirtes and Glymour, 1991) and its variants, such as Conservative- PC (Ramsey et al., 2012) and PC-stable (Colombo et al., 2014). Score-based methods, including GES (Chick- ering, 2002a) and hill-climbing (Koller and Friedman, 2009), search for the optimal DAG under a predefined scoring function combined with constraints. Gradient-based methods further extend score-based approaches by transforming discrete searches into continuous equality constraints, such as NOTEARS (Zheng et al., 2018), GAE (Ng et al., 2019), and GraN-DAG (Yu et al., 2019). However, without specific assumptions, the identifiability of structures solely derived from observational data is theoretically limited. Some works have extended causal discovery algorithms to intervention settings (Hauser and Bühlmann, 2012; Wang et al., 2017; Yang et al., 2018; Squires et al., 2020; Zhang et al., 2024; von Kügelgen et al., 2024). For example, GIES (Hauser and Bühlmann, 2012) pioneered the case with known hard interventions, while IGSP (Wang et al., 2017; Yang et al., 2018) introduced a greedy sparse ordering method to handle known general interventions. Building on this, UT-IGSP (Squires et al., 2020) partially addresses the scenario of unknown target interventions. However, given the diversity of intervention data types and experimental strategies, these methods still fail to universally address various intervention scenarios, limiting effective causal identification. In contrast, our approach offers a unified solution for causal discovery from intervention data by leveraging a joint causal inference framework combined with supervised causal learning methods, achieving practical empirical performance. A.2.2 Supervised Causal Learning Supervised Causal Learning (SCL) is an emerging paradigm of causal discovery, and has demonstrated strong empirical performance. The strength of SCL lies in its ability to learn complex classification mechanisms, contrasting with traditional rule-based logics for detecting causal relations. Early SCL work focused on pairwise causal discovery, such as RCC (Lopez-Paz et al., 2015) and MRCL (Hill et al., 2019). Further efforts have shifted toward multivariate causal discovery, with models like DAG-EQ (Li et al., 2020) and SLdisco (Petersen et al., 2023), which apply to linear causal models. ML4C (Dai et al., 2023) and ML4S (Ma et al., 2022) focus on v-structure detection and skeleton learning, respectively. For interventional data, methods such as CSIvA (Ke et al., 2023b) and AVICI (Lorch et al., 2022) extended SCL to handle known, hard interventions. However, these models are limited in generalization. In contrast, our approach adopts a test-time training paradigm, which introduces a process that acquires training data during test time. This allows the model to ’overfit’ to the specific biases of the test data, addressing the generalization issue and significantly improving performance in real-world settings. Go to Appendix Index20 A.2.3 Test-Time Training A body of work (Sun et al., 2020; Wang et al., 2020, 2022; Liu et al., 2021; Hardt and Sun, 2024; Sun et al., 2025; Dalal et al., 2025; Chen and Liang, 2025; Tandon et al., 2025) has explored test-time training paradigm to address the challenge of distribution shift in test data. A common approach is to identify an auxiliary task that aids the model in better adapting to the test data. For instance, T (Sun et al., 2020) jointly trains a model for rotation prediction and image classification. T++ (Liu et al., 2021) extends this by employing a contrastive learning approach as an auxiliary task for adaptation. In the context of causal discovery, test-time training exhibits distinctive characteristics. Unlike prior works, causal discovery benefits from the ability to generate test-time training data in a self-supervised manner. We refer to this approach as self-augmentation. This self-augmentation allows test-time training data to capture nuances and biases inherent to the test data, thereby making it well-suited for supservised causal learning. Currently, there is limited work in this area, with the most relevant being ML4S (Ma et al., 2022), which proposes a heuristic for generating vicinal graphs at test time for skeleton learning. To the best of our knowledge, we are the first to present a systematic approach to test-time adaptation in SCL, and we extend it to a broader context of causal discovery in general interventional settings. A.2.4 Joint Causal Inference Joint Causal Inference (JCI) (Mooij et al., 2020) presents a joint causal inference framework aiming to integrate multiple observed outcomes collected during different experiments (i.e., contexts). In this framework, the observed variable set is divided into two disjoint sets: system variables and context variables. After that, S-FCI (Li et al., 2023) builds upon JCI by introducing a new constraint-based algorithm, enabling learning from observational and intervention data across multiple domains. (Mascaro and Castelletti, 2023) also provides a graphical representation theory of I-MEC under general interventions and designs compatible priors for Bayesian inference to ensure score equivalence of indistinguishable structures. Although JCI provides a novel framework for addressing intervention problems and simplifies the unification of different intervention settings for supervised causal learning, it does not directly address the critical questions of what the appropriate learning objectives are when applying supervised learning, and how the learning process should be designed. To this end, we highlight the overlooked connection between identifiability in ICD and SCL. Through the modified PC-like SCL algorithm, a two-phase learning method is developed for predicting identifiable causal structures (I-CPDAG) while ensuring theoretical identifiability. B Inplementation Details B.1 Training Data Acquisition via Self-Augmentation B.1.1 Parameter Re-Estimation for CPT The joint distribution represented by a randomly generated conditional probability table (CPT) may sig- nificantly differ from the true joint distribution corresponding to the augmented data, thus undermining the forward-sampling process. Therefore, we propose approximating the joint posterior not only on the structure of the causal graph but also on the parameters of its conditional probability distributions. This approach ensures that the distribution of the augmented graph maintains in-domain "similarity" with the true underlying distribution. Specifically, we first use maximum likelihood estimation (Myung, 2003) to estimate the conditional probability tables (CPTs) of the initial seed graph. Secondly, during the IS-MCMC iteration process, we propose that the CPTs of the current step’s graph are determined by the CPTs from the previous step. For each node in the current step’s graph, if its parent nodes remain the same as those in the corresponding node from the previous step’s graph, its CPT is directly inherited from the previous step. Otherwise, nodes may lose previous parent nodes due to edge deletion or gain new parent nodes due to edge addition as a result of graph structure transformation operations. For edge deletions that result in the loss of previous parent nodes, we can adjust the CPT by marginalization. For instance, in stept−1, if nodeXhas parent nodesYandZ, the corresponding CPT encodes the Go to Appendix Index21 distributionP(X|Y,Z). In stept, if the edgeY → Xis deleted, the new conditional distribution should encodeP(X|Z). This can be naturally achieved by marginalizing out the nodeYusing the law of total probability, i.e., P (X|Z) = P y P (X|Z,Y = y). For edge additions that result in new parent nodes, the situation becomes more complicated. For instance, in stept, if a new edgeU → Xis added, introducing a new parent nodeU, the corresponding conditional probability distribution should encodeP(X|Y,Z,U). Considering we primarily deal with discrete data, and the Dirichlet-multinomial distribution is a natural choice for modeling categorical distributions, we use the Dirichlet distribution to sample different CPTs. Specifically, for each conditional probability distribution P(X|Y=y,Z=z,U=u), we estimate its parametersα i by maximizing the log-likelihood function of the data, which is given by: F (α) = logp(D|α) =log Y i p(p i |α) =log Y i Γ ( P k α k ) Q k Γ(α k ) Y k p α k −1 ik = N log Γ X k α k ! − X k log Γ (α k ) + X k (α k − 1) log ˆp k ! , wherelog ˆp k = 1 N P i logp ik represents the observed sufficient statistics, Γ(x) denotes the Gamma function and is defined to be: R ∞ 0 t x−1 e −t dt. Since this function does not have a closed-form solution, we employ a fixed-point iteration technique (Minka, 2000) to estimate the parameters. The core idea is to find an initial guess forαand iteratively refine it to converge to the maximum likelihood estimate. It seek a functionF(·) that provides a lower bound for the log-likelihood function: F (α)≥ N X k α k ! Ψ X k α old k ! − X k log Γ(α k ) + X k α k log ˆp k + C ! , whereCis a constant dependent onα, ensuring that the function is optimized atα. By setting the gradient to zero, a new iterative value for α is derived: α i k = Ψ −1 Ψ X k α i−1 k ! + log ˆp k ! , where inverse digamma function Ψ −1 can be efficiently solved using the Newton-Raphson method (Ypma, 1995). Subsequently, for each possible valuexof the variableX, we sample fromDirichlet(βα i ) to obtain P(X|Y=y,Z=z,U=u), thereby forming our target distributionP(X|Y,Z,U). Here,βis a hyperparameter that adjusts the variance, set to 0.25. B.1.2 Interventional Structural-Markov Chain Monte Carlo We first present the standard version of the IS-MCMC algorithm, as outlined in Algorithm 1. The first step involves initializing the seed graph, which is typically sampled from the prior distribution of the augmented graph. The second step estimates the parameters of the conditional probability tables using maximum likelihood estimation based on the initial graph structure and the augmented data. The main loop of the algorithm consists of three parts: (1) generating a proposed (or candidate) sample graph from the proposal distribution with intervention constraint; (2) calculating the acceptance probability using the acceptance functionα(·), based on the proposal distribution and the posterior probability; (3) accepting the candidate sample with probabilityα(the acceptance probability), or rejecting the candidate sample with probability 1− α. However, as the number of nodes increases, the computational efficiency of the standard IS-MCMC algorithm becomes untenable. Therefore, for the IS-MCMC algorithm applied to large graphs, we enhance it from two aspects. For the initial random seed graph in the first step, we recommend using proxy algorithms such as JCI+BLIP, which can provide a good starting point that is closer to the target graph structure, thereby accelerating the convergence rate. For the iterative MCMC process in the second step, we suggest Go to Appendix Index22 Algorithm 1 Standard IS-MCMC Algorithm Workflow Input: Chain Length: T, Expected Training Sample Size: N Output: Synthetic training data set: (G (1) ,D (1) ), (G (2) ,D (2) ),..., (G (N) ,D (N) ) Pipeline: 1. Initialize Seed Graph G (0) ∼G I 2. Estimating CPT parameters by maximum likelihood Θ (0) = MLE(G (0) ,D I ) 3. for iteration t = 1, 2,... do: (a) Propose: G cand ∼ q(G (t) |G (t−1) ) with Intervention Constraint. (b) Parameter estimate for current graph: Θ cand = Estimation(G cand ,D I ) (c) Calculate Acceptance Probability:α(G cand |G (t−1) ) = min n 1, q(G (t−1) |G cand )P(G cand |D I ) q(G cand |G (t−1) )P(G (t−1) |D I ) o (d) Randomly Sample: u∼ Uniform[0, 1] (e) Accept the proposal: G (t) ←G cand , if u < α then Reject the proposal: G (t) ←G t−1 , Θ (t) ← Θ (t−1) (f) Forward Sampling data: D (t) ← Sampling(G (t) , Θ (t) ) Algorithm 2 Optimized IS-MCMC Algorithm Workflow ❶ Initialize Seed Graph with Proxy Algorithm: G (0) ← proxy(D I ) ▷ Good Initial State Estimating CPT parameters by maximum likelihood: Θ (0) = MLE(G (0) ,D I ) Parallel i∈ 1, 2,...,N do: ▷ Efficient Parallel Chains for iteration t = 1, 2,...,L do: ❷ Propose: G cand ∼ q(G (t) |G (t−1) ) with Intervention Constraint. ▷ Intervention Constraints ❸ Re-parameterize according to the previous graph ( B.1.1): Θ cand =Estimation(G cand ,Θ (t−1) )▷Parameters Reuse ❹ Calculate Acceptance Probability:α(G cand |G (t−1) ) = min 1, q(G (t−1) |G cand )P(G cand |D I ) q(G cand |G (t−1) )P(G (t−1) |D I ) ❺ Randomly Sample: u∼ Uniform[0, 1] Accept proposal: G (t) ←G cand , Θ (t) ← Θ cand if u < α else Reject proposal: G (t) ←G t−1 , Θ (t) ← Θ (t−1) ❻ Forward Sampling data: D (t) ← Sampling(G (t) , Θ (t) ) parallelizing the computation by simultaneously executing multiple chains. Each chain performs consecutive edge perturbations, followed by sequential accept-reject sampling, to conduct approximate computations. The optimized IS-MCMC workflow is shown in Algorithm 2. If not specified otherwise, we use the optimized IS-MCMC method as the default setting. train infer Marginal Independence Test Intervention Observation Hard Soft Single Multiple Know Unknow A C B 퐼 1 퐼 2 ABC퐼 1 퐼 2 퐷 0 00 퐷 1 10 퐷 2 01 Augmented Graph Augmented Data Data Generation for Test-Time Training JCI-Style Good Initial State Random Seed Proxy Seed Efficient Parallel Chains 풢 (푡−1) 풢 (푡) 훩 (푡−1) 훩 (푡) IS-MCMC Intervention Constraints CPT Reuse Training Pair 풟 (푡) 풢 0 Skeleton Learning Orientation Learning Undirected Complete Graph 0-order Graph All 0-order features of edge train infer 1-order Graph ... ...... 1-order CI Test 2-order CI Test k-order CI Test All (k-1)-order features of edge train infer k-order Graph 풌→퐦퐚퐱(풅) Skeleton Extra Triple Vicinity Query c b a 푃퐶 푋 푎 푃퐶 푋 푏 푃퐶 푋 푐 All vicinity features of UT Unshielded Triples (UT) Identifiable Triples ... ... k-order Learner v-structure Learner Conflict filtering Meek Rules JCI Prior 퓘-푪푷푫푨푮 Figure B.1 The detailed framework of TICL. Go to Appendix Index23 B.2 Two-phase Supervised Causal Learning The two-phase supervised learning approach ofTICLis inspired by the PC algorithm (Spirtes and Glymour, 1991). We first describe the PC algorithm in detail and then reinterpret it from a machine learning perspective, providing a formalization. We also provide detailed feature engineering and techniques. As shown in Figure B.1, we further present a more detailed framework diagram, which includes the specific process of two-phase learning. B.2.1 Revisiting the PC Algorithm The PC algorithm consists of three main phases, as shown in Algorithm 3. Phase 1 identifies the skeleton and the separating sets, determining the existence of edges. Phase 2 orients the unshielded triples in the skeleton based on the separating sets, establishing edge directionality. Finally, phase 3 further refines the edge directionality using heuristic Meek Rules. Algorithm 3 Pipeline of PC Algorithm (Spirtes and Glymour, 1991) Require: Conditional Independence Information CI among Variables. Undirected Complete Graph C Obtained by Variables. Output: Separation Sets: S (Temporary Product) Skeleton: K (Intermediate Product) Unshielded Triple Sets of Skeleton: U (Temporary Product) Partially Directed Acyclic Graph, i.e., PDAG: P (Intermediate Product) Completed Partially Directed Acyclic Graph, i.e., CPDAG: G Pipeline: 1. Adjacency Determination: Find the skeletonKand separation setsSSbased on undirected complete graph C and conditional independence information CI using Algorithm 4. 2.Orientation Determination: Orient unshielded triplesUin the skeletonKto derive partially directed acyclic graph P based on the separation sets S using Algorithm 5. 3. InPorient as many of the remaining undirected edges as possible by repeated application of rules R1-R3 to derive completed partially directed acyclic graph G. In phase 1, we start with a complete undirected graphC. This graph is then sparsified through iterative conditional tests informationCI, where an edgeX i − X j is removed ifX i is conditionally independent ofX j given some subsetSof the remaining variables of the currentk-order graph. These conditional independence queries proceed in a cascading manner, making the algorithm computationally efficient for high-dimensional sparse graphs since we only need to query conditional independencies up to orderd−1, wheredis the maximal in degree of the underlying DAG. We summarize this process in Algorithm 4. In phase 2, it aims to identify V-structures. Specifically, it considers all unshielded triplesUin the skeletonK and orients an unshielded triple (X a ,X c ,X b ) into a V-structure if and only ifX c is not in the separating set of X a and X b . We also summarize this process in Algorithm 5. In phase 3, heuristic meek rules (Meek, 1995) are further applied iteratively to orient as many of the remaining undirected edges as possible. It contains the following three rules: • R1: If X a → X b − X c exists, change X b − X c to X b ← X c (to avoid creating a new V-structure). • R2: If X a → X b → X c exists, change X a − X c to X a → X c (otherwise a directed cycle is created). • R3: If there areX a − X c1 → X b ,X a − X c2 → X b , andX c1 ,X c2 are not adjacent, changeX a − X b to X a → X b (otherwise a new v-structure or a directed cycle is created). Go to Appendix Index24 Algorithm 4 Pipeline of Adjacency Determination / Phrase 1 of thePCAlgorithm (Kalisch and Bühlman, 2007) Require: Undirected Complete Graph C Induced by Variables. Sequential Order k. Output: Skeleton: K Separation Sets: S Pipeline: k = 1 repeat for each adjacent pair (X i ,X j ) in C do for each subset S ⊆ adj(C,X i )\X j or adj(C,X j )\X i with |S|= k do if X i and X j are conditionally independent given S then Delete edge X i − X j from C Let S (X i ,X j ) =S (X j ,X i ) =S end if end for end for▷ This process can be supervised for learning k = k + 1 unitl all adjacent pairs (X i ,X j ) in C satisfy |adj(C,X i )\X j |< k Algorithm 5 Pipeline of Orientation Determination / Phrase 2 of thePCAlgorithm (Kalisch and Bühlman, 2007) Require: Skeleton: K Separation Sets: S Unshielded Triple Sets of Skeleton: U Output: Partially Directed Acyclic Graph, i.e., PDAG: P Pipeline: for each non-adjacent pair (X a ,X b ) with common neighbour X c in U do if X c /∈S (X a ,X b ) then Replace X a − X c − X b in K by X a → X c ← X b end if end for▷ This process can be supervised for learning P =S Go to Appendix Index25 B.2.2 Connecting PC and Supervised Machine Learning In addition to phase 3, phase 1 and 2 can be considered as classification tasks regarding the determination of the presence of edges and the orientation of unshielded triples based on extracted conditionally independent features. We continue below with a detailed explanation. ■ From a machine learning perspective, we formalize phase 1 as follows: Task: For the current k-oreder graph, classify whether there is an edge between vertices X i and X j . Featurization: Query a subsetSof the current remaining variables ofk-oreder graph, and calculate the conditional dependence between X i and X j given S: F (k) (X i ,X j ) =min S⊆X\X i ,X j X i ∼ X j |S Classifier: Train a binary classifier for (k + 1)-order graph using current k-order features: C (k+1) skeleton (F (k) (X i ,X j ) ) := ( adjacent F (k) (X i ,X j ) ̸= 0 non− adjacent F (k) (X i ,X j ) = 0 ■ From a machine learning perspective, we formalize phase 2 as follows: Task: For all unshielded triples U, classify whether each triple < X a ,X c ,X b > is a v-structure. Featurization: Query all separating setsSSsatisfyingX a ⊥ X b |S, and calculate the existence Boolean feature: F <X a ,X c ,X b > = 1 X c ∈S 0 X c /∈S Classifier: Train a binary classifier using features: C orientation (F <X a ,X c ,X b > ) := v-structure F <X a ,X c ,X b > = 0 non− v-structure F <X a ,X c ,X b > ̸= 0 Theorem 3. (Spirtes and Glymour (Spirtes et al., 2000)) Let the distribution of X be faithful to a DAGG, and assume that we are given perfect conditional independence information about all pairs of variables (X i ,X j ) in X given subsetsS ⊆X\X i ,X j . The output of PC algorithm is the CPDAG that representsG. However, in practical applications, conditional independence tests relying on data estimates can encounter issues like limited sample size, data quality issues, and intricate dependency structures (Günther et al., 2022). Thus, it is advisable to employ more systematic featurization procedures to enhance classification performance in a more effective and robust manner. B.2.3 Skeleton and Orientation Feature Engineering In the realm of skeleton learning, drawing upon the experience of (Cheng et al., 2002; Ding et al., 2020; Xiang and Kim, 2013; Ma et al., 2022), we have extracted primary features of the two categories: quantitative k-order conditional dependencies and local structural information. For orientation learning, again relying on the experience of (Vijaymeena and Kavitha, 2016; Zanga et al., 2022; Dai et al., 2023), we have extracted main features of the two categories: quantitative unshielded triplet conditional dependencies and local structural information. The intuition behind it is detailed below. In skeleton inference, we first look for useful features of the type of conditional test information. Intuitively, for the numerical vectors resulting from the conditional tests performed on the current node pairs, which indicate the conditional dependency between nodes, we consider them as primary features. Additionally, when conducting higher-order conditional tests, if the dependency relationships decrease significantly, it suggests that they may be blocked by new nodes, so this residual conditional dependency can be seen as another Go to Appendix Index26 푋 1 푋 2 푋 3 푋 4 푋 5 푋 6 0.8 0.70.6 0.70.8 0.9 0.7 0.5 0.9 1-order Graph 푺 1 =1: 푋 1 ~푋 2 |푋 3 =0.7 푋 1 ~푋 2 |푋 4 =0.6 푋 1 ~푋 2 |푋 5 =0.5 푋 1 ~푋 2 |푋 6 =0.8 푺 2 =2: 푋 1 ~푋 2 |푋 3 ,푋 4 =0.6, 푋 1 ~푋 2 |푋 3 ,푋 5 =0.3, 푋 1 ~푋 2 |푋 3 ,푋 6 =0.7, 푋 1 ~푋 2 |푋 4 ,푋 5 =0.4, 푋 1 ~푋 2 |푋 4 ,푋 6 =0.6, 푋 1 ~푋 2 |푋 5 ,푋 6 =0.1, 푺 0 =0: 푋 1 ~푋 2 =0.8 푋 1 ~푋 3 =0.9,푋 1 ~푋 4 =0.7 푋 1 ~푋 5 =0.6,푋 1 ~푋 6 =0.8 푋 2 ~푋 3 =0.7,푋 2 ~푋 4 =0.7 푋 2 ~푋 5 =0.5,푋 1 ~푋 6 =0.9 Figure B.2 Feature extraction example of edge X 1 − X 2 in 1-order skeleton graph. Table B.1 The features extracted for skeletal learning. As shown in Figure B.2,k= 2, is used to illustrate the design intuition and computational examples. Skeleton Feature TypeNameCalculateExample Raw ResultsDimension CI Test Information k-order Conditional Dependence X i ∼ X j |S 2 , |S 2 |= 2 [0.6, 0.3, 0.7, 0.4, 0.6, 0.1] Not Fixed → 15 Residual Conditional Dependence X i ∼ X j |S 1 − min X q X i ∼ X j |S 2 ∪X q , |S 1 |= 1, |S 2 |= 2 [0.4, 0.2, 0.4 ,0.7]Not Fixed → 15 Structural Information Competitiveness |X q |X q ∈nbd(X i ),C k−1 (X i ,X j )>C k−1 (X i ,X q )| | nbd G k−1 (X i ) | −1 , |X q |X q ∈nbd(X j ),C k−1 (X j ,X i )>C k−1 (X j ,X q )| | nbd G k−1 (X j ) | −1 , C k−1 (X i ,X j ) [0.5, 0.75, 0.8]3 Degreedeg(X i ),deg(X j ) [4, 4]2 Density |nbd(X i )∩nbd(X j )| min (|nbd(X i )|,|nbd(X j )|) [1]1 complementary feature. Considering the structural aspect, conditional independence may be influenced by the information of the local graph structure in which the current node pairs are located. Generally, the sparser the local graph is, the less likely it is to be disconnected by conditional independence tests. Therefore, we select three main features for measurement, including the relative competitiveness of adjacent edges between node pairs, the degrees of node pairs, and the overlapping density of adjacent edges of node pairs. In Figure B.2, we provide an example of a 1-order skeleton graph, based on which we further summarize the types of feature variables and calculation methods in Table B.1. Go to Appendix Index27 푃퐶 푋 푎 푃퐶 푋 푏 푃퐶 푋 푐 푋 푎 푋 1 푋 푐 푋 푏 푋 2 푋 4 푋 3 푋 5 Variable Domain: 푋 푎 ~푋 푏 ,푋 푎 ~푃퐶 푋 푏 ,푃퐶 푋 푎 ~푋 푏 ,푃퐶 푋 푎 ~푃퐶 푋 푏 Condition Domain: ∅,푺,푋 푐 ,푋 푐 ∨푺,푃퐶 푋 푐 ,푃퐶 푋 푐 ∨푺 푋 푎 ~푋 푏 |∅=0.9 푋 푎 ~푋 푏 |푺=0.8 푋 푎 ~푋 푏 |푋 푐 ∨푺=0.6 푋 푎 ~푋 푏 |푃퐶 푋 푐 =0.8 푋 푎 ~푋 푏 |푃퐶 푋 푐 ∨푺=0.7 푋 푎 ~푃퐶 푋 푏 |∅=0.7 푋 푎 ~푃퐶 푋 푏 |푺=0.8 푋 푎 ~푃퐶 푋 푏 |푋 푐 ∨푺=0.7 푋 푎 ~푃퐶 푋 푏 |푃퐶 푋 푐 =0.5 푋 푎 ~푃퐶 푋 푏 |푃퐶 푋 푐 ∨푺=0.4 푃퐶 푋 푎 ~푋 푏 |∅=0.8 푃퐶 푋 푎 ~푋 푏 |푺=0.4 푃퐶 푋 푎 ~푋 푏 |푋 푐 ∨푺=0.3 푃퐶 푋 푎 ~푋 푏 |푃퐶 푋 푐 =0.6 푃퐶 푋 푎 ~푋 푏 |푃퐶 푋 푐 ∨푺=0.2 푃퐶 푋 푎 ~푃퐶 푋 푏 |∅=0.9 푃퐶 푋 푎 ~푃퐶 푋 푏 |푺=0.8 푃퐶 푋 푎 ~푃퐶 푋 푏 |푋 푐 ∨푺=0.7 푃퐶 푋 푎 ~푃퐶 푋 푏 |푃퐶 푋 푐 =0.6 푃퐶 푋 푎 ~푃퐶 푋 푏 |푃퐶 푋 푐 ∨푺=0.4 푺=푋 1 ,푋 5 ,푋 1 ,푋 5 푋 푎 ~푋 푏 |푋 푐 =0.8 푃퐶 푋 푎 ~푋 푏 |푋 푐 =0.4 푋 푎 ~푃퐶 푋 푏 |푋 푐 =0.8 푃퐶 푋 푎 ~푃퐶 푋 푏 |푋 푐 =0.9 Figure B.3 Orientation feature extraction example of unshielded triple X a − X c − X b . In the unshielded triplet orientation, we also seek useful features of the condition test information type. Unlike relying solely on specific triplet condition test information, we extend it to the neighborhood. Specifically, in addition to testing variablesX a andX b themselves, we also extend the set variables of their parents and children,PC X a andPC X b , as testing variable domains. For the condition variables, in addition to the empty set∅,S,X c , we further extend the set variables of the parents and children ofX c to the condition domain, including their element-wise unionSS∨ X c ,S∨ PC X c . Therefore, by selecting a pair of variables from the variable domain and a condition from the condition domain, we obtain a total of 4×6 = 24 neighborhood conditional dependency features. Furthermore, considering structural information, we overlap numerical and scale numerical measures of the size and overlap features of the sets of the variable domain and the condition domain to reflect the potential sparsity of the local graph structure. As shown in Figure B.3, we provide an example of directional feature extraction, based on which we further summarize the types of feature variables and calculation methods in Table B.2. Table B.2 The features extracted for orientation learning, as shown in Figure B.3, are used to illustrate the design intuition and computational examples. Orientation Feature TypeNameCalculateExample Raw Results Dimension CI Test Information Vicinity Conditional Dependence Variable| Condition, Variable: X a ∼X b ,X a ∼PC X b , PC X a ∼X b ,PC X a ∼PC X b Condition: ∅,S,X c ,X c ∨ S, PC X c ,PC X c ∨ S [0.9, 0.8, 0.8, 0.6, 0.8, 0.7, 0.7, 0.8, 0.8, 0.7, 0.5, 0.4, 0.8, 0.4, 0.4, 0.3, 0.6, 0.2, 0.9, 0.8, 0.9, 0.7, 0.6, 0.4] 24 Structural Information Overlap |set 1 ∩set 2 | min(|set 1 |,|set 2 |) , where set∈PC X a ,PC X b ,PC X c , S, and set 1 = X c , set 2 = S [0.5, 0, 0.5, 0.5, 0, 0.5, 1]7 Scaling #PC X a , #PC X b , #PC X c , #S, 1 #S P S i ∈S #S i [2, 3, 2, 3, 4 3 ]5 It is noted that the conditional dependence features mentioned above may be dynamically changing, which does not conform to the fixed-size features required by most traditional machine learning models. However, this can be easily addressed through classical kernel mean embedding techniques (Smola et al., 2007) to obtain Go to Appendix Index28 fixed-length embedding features: 1 |F| X z∈F (cos(< w j ,z > +b j )) m j=1 ∈R m . Here,m= 15, which signifies that each extended featureFnow possesses a fresh embedding dimension of 15. Furthermore, we incorporate an extra set of five statistics comprising maximum, minimum, mean, standard deviation, and set size. Consequently, we merge all these derived features to construct the input feature for the model. Additionally, various methods exist for assessing conditional dependence, such as determining the p-value through testing for conditional independence or using conditional mutual information (Cover, 1999). When dealing with categorical variables, the G2 test is a suitable option (Agresti, 2012). In our approach, we employ an approximate version of the G2 statistic and rely on p-values to evaluate conditional dependence. It is worth noting that the p-value can become insignificant due to the double-precision bit limit in computers. To overcome this issue, transforming p-values can be employed to mitigate such limitations and offer a means to evaluate conditional dependency. To begin with, we define the complementary error function as: g(z) = 1− 2 √ π Z z 0 e −t 2 dt, and utilize the quantity z as the inverse of g: z = g −1 (x). By applying the non-linear transformation ofg −1 to a given p-valuex, we derive a re-scaled quantity that enhances the assessment of conditional dependency. In essence,zcan be interpreted as a measure akin to z-scores in a standard normal distribution; for example, if the p-value is 0.01, thenz= 3, as a value of 3-sigma signifies that the probability of data falling within a 3-sigma range in a normal distribution is 0.99. C Convergence of Inteventional Structural-MCMC In this section, we first confirm the convergence of the posterior distribution in the IS-MCMC algorithm based on the Structure MCMC framework. Furthermore, we provide a detailed discussion of the key factors considered in optimizing the IS-MCMC process. Besides, we also provide visualization experiments on convergence analysis. C.1 Guarantee of Posterior Distribution Convergence Our IS-MCMC algorithm follows a standardized Structure MCMC process (Madigan et al., 1995; Su and Borsuk, 2016; Kuipers and Moffa, 2017) in the intervention-augmented graph space. Please refer to Algorithm 2, which corresponds to Steps ❶ - ❻ in the main text as a general introduction to the MCMC process. By employing the Metropolis-Hastings sampler within the standardized Structure MCMC framework, the Markov chain is guaranteed to have a stationary distribution equal to the posterior distributionP(G |D) (Su and Borsuk, 2016). It is important to note that, under the premise of this theoretical guarantee, the novelty of our work lies in proposing the use of the posterior distribution as the target for acquiring training samples and optimizing feasible solutions in intervention scenarios. C.2 Optimization Considerations for IS-MCMC However, two critical factors need to be addressed. First, managing the time complexity of this process is crucial for maintaining efficient inference. Second, ensuring the effectiveness of the IS-MCMC process on the augmented graph. Go to Appendix Index29 Asia Survey Figure C.1 Convergence Visualization Analysis on different datasets. •From the T-driven perspective on efficiency: Since we leverage the "test-time" phase, the convergence speed of MCMC is particularly important. We observed that a good initial state enables the Markov chain to converge faster, as evidenced by Figure 5 in the experimental section of our paper. Additionally, we modified the algorithm to use parallel multi-chain sampling, allowing for more efficient exploration of the graph space, as shown in Figure 8 of our experimental results. •From the JCI-driven perspective on adaptability: The MCMC process operates on augmented graphs rather than standard causal graphs. Thus, additional constraints are required to ensure the validity of the augmented graph. To address this, we introduced intervention constraints to ensure consistency of the augmented graph within the JCI framework. In summary, our IS-MCMC algorithm effectively addresses multiple challenges under the T + JCI framework and demonstrates that using posterior estimates as training data is highly beneficial for test-time SCL. C.3 Convergence Visualization Analysis of IS-MCMC As shown in Figure C.1, we plotted the goodness-of-fit scores using random seed and proxy seed strategies on the Asia and Survey datasets, respectively. By treating the initial state as low-quality, the IS-MCMC strategy gradually improves the score until it converges to a stationary state. The proxy seed accelerates this process, so in practice, our method can significantly benefit from this process, especially in test-time training. D Identifiability Theory of Two-Phase SCL D.1 Intuitive Analysis D.1.1 Observation We use the observations from phase 2 as an illustration. As mentioned in Section B.2, we formalize the v-structure orientation of PC as a specific classification problem, that is, determining whether the non-shielded triple ⟨X a ,X c ,X b ⟩ forms a v-structure by classifier C PC . The logic of PC’s orientation is clearly asymptotically correct (i.e., the CI test becomes fully accurate when Go to Appendix Index30 the number of records goes to infinity). However, in practice, when the number of records is finite, its empirical performance is not satisfactory. Thus, further enhancements, such as Majority-PC (MPC) (Colombo et al., 2014), have been developed to address this limitation. MPC is a sample-based enhancement of PC’s orientation, achieving better performance with finite samples. Instead of identifying only one separating setSS, MPC finds all possible separating setsSSand counts how many of them contain X c . The logic can be recast as follows: Featurization: Finds all separating sets S of X a ,X b , and defines a real-valued feature: F <X a ,X c ,X b > = |S i |X c ∈ S i ∈S| |S| Classifier: Train a binary classifier using features: C MPC (F <X a ,X c ,X b > ) := v-structure F <X a ,X c ,X b > ≤ 0.5 non− v-structure F <X a ,X c ,X b > > 0.5 Theoretically, both PC and MPC, as "hand-crafted" classifiers, are asymptotically correct. However, in practice, MPC exhibits greater complexity in its classification mechanism compared to PC, resulting in improved empirical performance. Nonetheless, from a machine learning perspective, both PC and MPC remain "simple" in terms of their feature representation and classification strategies. D.1.2 Motivation The primary motivation for modifying PC into a PC-like SCL method is to leverage both theoretical identifiability and enhanced empirical performance: •Theoretical Guarantees. Our PC-like SCL approach retains the asymptotic properties of the original PC algorithm. Specifically, the method detects the correct CPDAG (orI-CPDAG in the interventional setting) when the sample size approaches infinity. This ensures theoretical identifiability. •Empirical Performance. In finite-sample scenarios, where PC’s reliance on conditional independence (CI) tests can lead to errors, the SCL approach outperforms traditional methods. By combining feature- rich representations with a robust classification mechanism, SCL achieves superior empirical results, as demonstrated in our paper, by comparingTICLwith PC-JCI or other non-SCL methods. Similar empirical evidences are also presented in prior works (Li et al., 2020; Ma et al., 2022; Dai et al., 2023). D.2 Advantage Building on these observations and motivation analysis, our goal is to design a PC-like SCL algorithm that maintains theoretical asymptotic correctness while encouraging a more systematic feature representation and enabling the learning of more sophisticated classification mechanisms. This approach aims to make our method more "robust" than PC, achieving superior empirical performance. Specifically, Our PC-Like SCL method generalizes these ideas: •Featurization: SCL extracts a richer set of features that capture conditional dependencies and structural patterns around ⟨X a ,X c ,X b ⟩. •Classifier: Instead of relying on heuristic rules likeC PC orC MPC , we train a classifier on synthetic data. The model learns complex interactions among features, avoiding the limitations of error-prone conditional independence (CI) tests. In the asymptotic regime, the learned classifier converges to an equivalent, theoretically correct solution like C PC orC MPC . However, in practical settings, SCL’s broader feature set and data-driven optimization enable significantly better empirical performance. For example, Table 1 and 2 compare JCI+SCL with JCI-PC and other non-SCL methods, clearly showing SCL’s superiority. Go to Appendix Index31 D.3 Theoretical Analysis D.3.1 Restating the PC-like SCL Process The PC-like SCL approach proposed in our paper consists of two phases: •Phase 1: Skeleton Learning. A binary classifierC 1 takes the feature set (as detailed in Table B.1 and Figure B.2) as input and predicts the existence of edges between all pairs of nodes. •Phase 2: Orientation Learning. A binary classifierC 2 takes the feature set (as detailed in Table B.2 and Figure B.3) as input and predicts the v-structure for each unshielded triple (UT), using the skeleton learned from Phase 1. Finally, Meek’s rules are applied to orient as many causal directions as possible. We aim to prove that bothC 1 andC 2 are asymptotically correct (i.e., they output results equivalent to those of the PC algorithm when the sample size approaches infinity). For simplicity, we demonstrate the proof for C 2 , as the proof for C 1 follows the same principle. D.3.2 Asymptotic Correctness of the PC-like SCL Approach Definition D.1 (Overlap Coefficient).OLP (A, B):=|A∩B|/min (|A|,|B|), where A and B are two sets of variables. Definition D.2 (Discriminative Predicate). A discriminative predicate is a binary predicate function defined over the domain ofC 2 ’s feature vector. It can be regarded as a special classifier with a predefined mechanism. Definition D.3 (Strong Predicate). A strong discriminative predicate satisfies the following two criteria when applied to a UT’s feature vector: (a) It evaluates to ’true’ if the UT forms a v-structure. (b) It evaluates to ’false’ if the UT does not form a v-structure. A strong predicate is asymptotically correct because its output aligns with that of the PC algorithm when the sample size approaches infinity. The key statement here is that there exists at least one strong predicate, meaning we can construct a static classifier using the proposed feature set to achieve asymptotic correctness. Lemma 4 (Existence of a Strong Predicate). For a canonical dataset with infinite samples, the predicate OLP (X c ,S)= 0 is a strong discriminative predicate. It corresponds precisely to the v-structure detection logic of the PC algorithm. Proof. According to the PC algorithm’s asymptotic correctness theorem: For a canonical dataset with infinite samples (assuming faithfulness and perfect conditional independence information): •If an unshielded triple⟨X a ,X c ,X b ⟩forms a v-structure, thenX c does not belong to any separation set of (X a ,X b ). • If ⟨X a ,X c ,X b ⟩ does not form a v-structure, then X c belongs to all separation sets of (X a ,X b ). In this setting (infinite samples), there is no ambiguity:X c is either in all separation sets or in none of them. Now consider the predicate OLP (X c ,S) = 0, which evaluates to true if and only if X c /∈S for all S ∈S. • This implies that ⟨X a ,X c ,X b ⟩ is a v-structure. •Therefore,OLP (X c ,S)= 0 is a strong predicate because it perfectly discriminates v-structures from non-v-structures. Theorem 5 (Asymptotic Correctness ofC 2 ). By employing a learning model with the universal approximation property, the classifierC 2 is asymptotically correct when classifying a canonical dataset with infinite samples. Go to Appendix Index32 Proof.From Lemma 4, there exists a strong discriminative predicatePthat achieves zero loss on a canonical dataset with infinite samples. •GivenPas the ground truth for v-structure detection, a machine learning model with universal approximation capability can approximatePto arbitrary precision, achieving performance no worse than P. •Additionally, there may exist multiple strong predicates that satisfy the criteria of Definition D.3. Therefore, C 2 can converge to any one of these strong predicates in the asymptotic regime. Summary: The theoretical guarantees demonstrate that the PC-like SCL approach, by leveraging the universal approximation capabilities of learning-based classifiers, retains the asymptotic correctness of the original PC algorithm while enabling superior empirical performance on finite data. E Experimental Details E.1 Benchmark Our experiments are carried out on 14 different causal graph datasets inspired by real-world applications from bayesian network repository 2 . All discrete networks undergo thorough quality checks and necessary repairs, ensuring that the sum of all conditional probability distributions is 1, there are no single-level virtual nodes, and there are no dependencies. The statistics of these bayesian networks are shown in Table E.1. We also provide visualizations of some causal graphs, as shown in Figure E.1. Table E.1 Statistics and description of bayesian networks we used. Network #Nodes #Edges #Parameters Max in-degree Avg. degree Earthquake541021.60 Survey662122.00 Asia881822.00 Sachs111717833.09 Child202523022.50 Insurance27521,00833.85 Water326610,08354.12 Mildew3546540,15032.63 Alarm374650942.49 Barley4884114,00543.50 Hailfinder56662,65642.36 Hepar2701231,45363.51 Win95pts7611257472.95 Pathfinder10919572,07953.58 E.2 Baselines The development history of causal discovery from interventional data has been a gradual and incremental process. Here is a detailed introduction to some key works in this field: • GIES (Hauser and Bühlmann, 2012): This paper first extends the concept of Markov equivalence of Directed Acyclic Graphs for the first time to the case of interventional distributions arising from multiple 2 https://w.bnlearn.com/bnrepository/ Go to Appendix Index33 Burglary Earthquake Alarm JohnCalls MaryCalls (a) Earthquake AS E OR T (b) Survey asia tub smoke lung bronceither xray dysp (c) Asia PKC PKA JnkP38 Plcg PIP3 PIP2Raf Mek Erk Akt (d) Sachs BirthAsphyxia Disease LVH LVHreport DuctFlow CardiacMixing HypDistrib LungParench HypoxiaInO2 CO2 LowerBodyO2 RUQO2CO2Report LungFlow ChestXray XrayReport Sick Grunting GruntingReport Age (e) Child Age SocioEcon GoodStudentRiskAversion VehicleYearMakeModel RuggedAuto Mileage Antilock SeniorTrain DrivingSkill DrivQuality Accident ThisCarDam CarValue HomeBaseAntiTheft Theft ThisCarCost OtherCarCost PropCost OtherCar Airbag Cushioning MedCostILiCost DrivHist (f) Insurance Figure E.1 Graph Visualizations. interventional experiments. It further demonstrates that under reasonable assumptions of interventional experiments, the intervened Markov equivalence defines a more refined DAG partition than the observed Markov equivalence, thereby enhancing the identifiability of causal models. Moreover, they generalize the greedy equivalence search algorithm, proposing a greedy interventional equivalence search algorithm for regularized maximum likelihood estimation under such intervened conditions. • IGSP (Wang et al., 2017; Yang et al., 2018): In this paper, the authors propose two algorithms that utilize both observational and interventional data with consistency guarantees, and prove their consistency under the faithfulness assumption. These algorithms are intervention-adapted versions of the Greedy SP algorithm, and they are non-parametric, which makes them applicable to the analysis of non-Gaussian data as well. Subsequently, they first extend these identifiability results to general interventions that can modify the dependencies between the target variable and its causes without eliminating them, and propose the first consistent algorithm for learning a DAG in such an environment. •UT-IGSP (Squires et al., 2020): In this paper, the authors further extend interventions to unknown scenarios, that is, the problem of estimating causal Directed Acyclic Graph models from a mixture of observational and interventional data when the intervention targets are partially or completely unknown. They describe the intervened Markov equivalence classes of DAGs that can be identified from interventional data with unknown intervention targets. Additionally, they propose a provably consistent algorithm for learning the intervened Markov equivalence class from such data. The algorithm greedily searches the permutation space to minimize a novel scoring function. This algorithm is also non-parametric. • DCDI (Brouillard et al., 2020): In this paper, the authors first introduce a theoretically grounded method for learning causal structures based on neural networks that can utilize interventional data. They present two instances of this approach: one that relies on normalizing flows as a universal density approximator. They also demonstrate that the precise maximization of the proposed score will identify theI-Markov equivalence classes of causal graphs for both known and unknown target settings. • ENCO (Lippe et al., 2022): In this paper, the authors explore score-based methods for causal discovery from observational and interventional data. They argue that such methods often require constrained optimization to enforce acyclicity or lack convergence guarantees. Consequently, they formulate graph search as an optimization of the likelihood of independent edges, where the direction of edges is modeled as separate parameters. Thus, they provide convergence guarantees when interventions on all variables are Go to Appendix Index34 available, without the need to constrain the acyclicity of the scoring function. •BaCaDI (Hägele et al., 2023): In this paper, the authors aim to discuss the issue of causal discovery in the realistic scenario where interventional data is scarce. To address this shortcoming, they provide a principled Bayesian approach that operates within the continuous space of causal Bayesian networks (CBNs) and the latent probability representations of interventions. This enables them to approximate complex joint posteriors through efficient, gradient-based particle variational inference techniques, making it applicable to causal systems with many variables. Recently, Supervised Causal Discovery, which aims to compress data and map between causal relations using pre-synthesized causal graph datasets, has shown impressive performance in test data. We introduce some key works on supervised causal discovery from intervention data: •CSIvA (Ke et al., 2023b): In this paper, the authors believe that meta-learning enables models to generalize well to data from natural causal Bayesian networks, even with relatively few assumptions made during training on synthetic data. Therefore, they introduce for the first time a supervised approach to tackle the problem of causal structure induction. This method maps datasets composed of both observational and interventional samples to a structure. By introducing a novel transformer architecture, they aim to discover relations between variables across samples. •SDI (Ke et al., 2019, 2023a): In this paper, the authors introduce a novel neural network-based method for causal discovery from interventional data, capable of handling unknown interventions. By utilizing two sets of distinct parameters to model the causal mechanisms and the structure of the causal graph, experimental evidence indicates that this method can generalize to unseen interventions and can effectively perform partial graph discovery. •AVICI (Lorch et al., 2022): In this work, the authors posit that designing appropriate scores or tests that capture prior knowledge is challenging. Therefore, they propose amortizing the learning of causal structures, which involves training a variational inference model to directly predict causal structures from observational or interventional data. This approach leverages permutation invariance to exhibit robust generalization capabilities, particularly in the challenging field of genomics. Due to the limitations inherent in the type of methods and the constraints of the assumptions made, they are generally unable to effectively process interventional data across various experimental scenarios. Here, we further summarize, based on the experiments of their paper, the basic setup of their causal discovery from interventional data, as depicted in Table E.2: Table E.2 Summary of basic settings for various methods of causal discovery from intervention data Methods Groun-Truth DAGTesting DataIntervention Experiment Type#Variable Type#Obs. #Int. / per exp.Type#Int. Family#Int. Target GIES Random Graph (2007) (Linear Gaussian SEM) 10 ∼ 50 Continue 50 ∼ 10k50 ∼ 10k Single / Multiple Know Hard #Var. * .0 / .2 / .4 / .6 / .8 / 1.01 ∼ 4 IGSP ER (Linear Gaussian SEM) 10 ∼ 20 Continue 1k ∼ 100k1k ∼ 100k Single / Multiple Know Hard / Soft 1 ∼ 2#Var. * 0.4 / 0.5 UT-IGSP ER (Linear Gaussian SEM) 20Continue 1k ∼ 5k1 ∼ 5k Single / Multiple Know / Unknow Hard / Soft Know: (5) + Unknow: (5) Know: (1) + Unknow: (0 ∼ 3) DCDI ER (LGM / ANM / N) 10-20 Continue10k10k Single / Multiple Know / Unknow Hard / Soft #Variable#Var. * 0.1 SDI Random Graph (N) bnlearn 3 ∼ 13 Discrete 2560 (from code) 1111 (from code) Single Know Soft #Variable1 CSIvA ER (LGM / ANM / N) (MLP / Dirichlet) 30 ∼ 80 Continue 100 ∼ 1.5k100 ∼ 1.5k Single Know / Unknow Hard / Soft #Variable1 ENCO Random Graph (N) bnlearn 25Discrete 1k ∼ 100k20 ∼ 200 Single Know Soft #Variable1 AVICI ER / SF / WS / SBM / GRG (Linear / Random Fourier) 30Continue 200 (from code) 20 (from code) Single Know Hard / Soft #Var. * .51 BaCaDI ER / SF / SERGIO (2020) (Linear / Nolinear Gaussian) 20Continue10010 Single Know / Unknow Hard / Soft #Variable1 TICLbnlearn5-109 Discrete 2k ∼ 10k1k ∼ 10k Single / Multiple Know / Unknow Hard / Soft #Var. * .0 / .2 / .4 / .8 1∼ 3 Most existing intervention dataset learning methods have typically required a known set of intervention Go to Appendix Index35 targets, which is often a strong assumption. Therefore, when intervention targets are unknown, this naturally poses a challenge: how to identify intervention nodes individually becomes a key issue in structural learning. This task is of practical significance. For instance, in the realm of gene editing, techniques can cleave off-target genomic sites. Evaluating and detecting off-target effects accurately and devising corresponding strategies constitute significant research directions in current gene editing studies (Manghwar et al., 2020). Limited attention has been given to exploring methods for identifying intervention targets. We provide an introduction as follows: •UT-IGSP (Squires et al., 2020): As mentioned before, the UT-IGSP algorithm learns the intervention targets while learning the causal structure. Although the intervention refines the search space, greedy search of the sparsest permutations is very slow in the case of high-dimensional data, especially when using non-Gaussian conditional independence test. •CITE (Varici et al., 2021): In this paper, the authors primarily address the problem of estimating unknown intervention targets in causal directed acyclic graphs based on observational and interventional data. The focus lies on soft interventions within the framework of linear structural equation models. They propose a scalable and efficient algorithm capable of consistently identifying all intervention targets. The key idea is to estimate the intervention sites based on the difference between precision matrices associated with the observed data set and the intervened data set. • PreDITEr (Varici et al., 2022): After that, the authors of this paper further extend the previous method by eliminating the need to learn the entire causal model, focusing solely on learning intervention targets. The key point is to leverage the sparse changes imposed on the precision matrix of the linear model by intervention measures, composed of a series of precision difference estimation steps. Additionally, they infer the knowledge required to refine observational Markov equivalence classes into interventional MECs. •LIT (Yang et al., 2024): In this paper, the authors tackle for the first time the problem of identifying unknown intervention targets in a multi-environment setting. They further consider cases within the intervention target set that allow for potential confounding factors. They propose a two-stage algorithm to recover exogenous noise and match it with the corresponding endogenous variables. Under the assumption of causal sufficiency, intervention targets can be uniquely identified. In cases where potential confounding factors exist, a candidate intervention target set is provided, offering more information compared to previous works. As shown in Table E.3, we provide a comparison to briefly illustrate the types of unknown intervention experiments supported by the above intervention target detection algorithms. Table E.3 Comparison of situations (more suitable) supported by intervention target detection methods. Methods Causal Insufficiency Interv. Frequency Interv. MechanismData TypeDomain Single Multiple HardSoftContinue Discrete Single Multiple UT-IGSP✗✓∼✓ CITE✗✓∼✓∼✓✗ PreDITEr✗✓∼✓∼✓✗ LIT✓∼✓∼✓ To ensure fair experimentation, algorithms and hyperparameter tuning are involved in each of our experiments. For all baseline algorithms, we explore critical hyperparameters that govern their interactions, and report the best results. Moreover, we use open source code for all algorithms for evaluation, including code from the authors as well as various popular toolkits. The public code for the baseline is also included at the URL below. Firstly, we introduce a method specifically designed for causal discovery from observational and interventional data: • GIES: - Score Function: Gaussian BIC, https://github.com/juangamella/gies • IGSP: - Significanceα ∈ 1e3, 1e4, CI test∈ Gaussian, kci, hsic,https://github.com/uhlerlab/ graphical_model_learning Go to Appendix Index36 •UT-IGSP: - Significanceα∈1e3, 1e4, CI test∈Gaussian, kci, hsic,https://uhlerlab.github.io/ causaldag/utigsp.html •ENCO: - Sparsity Regularizerλ sparse ∈0.002,0.02, Epochs∈50,100,https://github.com/phlippe/ ENCO • AVICI: - Model ∈scm-v0, linear, rff, grn, https://github.com/larslorch/avici As mentioned in the main text, joint causal inference offers a robust theoretical framework that effectively combines interventional and observational data, simplifying causal discovery algorithms to operate solely on observational setting. We select four distinct algorithmic combinations, with parameter settings as follows: •PC - variant∈ original, stable, Significanceα ∈ 0.05, 0.01, CI test∈ fisherz, g2, chi2,https: //github.com/huawei-noah/trustworthyAI/tree/master/gcastle • HC - Score Function: Bdeu Score, https://github.com/pgmpy/pgmpy • BLIP - Time ∈ 60, 300, #Variable, https://cran.r-project.org/web/packages/r.blip/ • GOLEM - λ 1 ∈2e-2, 2e-3, ω ∈0.2, 0.3, https://github.com/ignavierng/golem Finally, for the specific task of intervention target detection, the configuration parameters for the comparative methods are as follows: •CITE -λ l1 ∈ 1e-1, 5e-1,Parent l1 ∈ 5e-3, 1e-2, 2e-2, 3e-2, 4e-2, 5e-2, 6e-2, 8e-2, 9e-2, 1e-1,https:// github.com/bvarici/intervention-estimation • PreDITEr - λ l1 ∈1e-1, 3e-1, λ pasp ∈1e-1, 2e-1, https://github.com/bvarici/uai2022-intervention-estimation-latents In addition to our selected baselines, we also enumerate other classic baselines relevant to causal discovery from intervention and explain why they were not chosen for our study: •DCDI: This method is a well-known and widely used approach for causal discovery using neural networks, relying on normalizing flows as a generic density estimator. Therefore, it is primarily applicable to continuous data types. We found in experimental testing that it almost never converges and operates at an unbearably slow pace in the case of discrete data. Hence, we disregard its results. •SDI: This method is the first to use neural networks and adapt to discrete data types for intervention causal discovery. However, due to its high complexity, it is often challenging to scale to large-scale node graphs. Additionally, according to the original repository code, its hybrid C-language requirement for compiler adaptation leads to abnormal errors. Hence, we disregard its results. •CSIvA: The lack of accessible code and complex parameter settings make it impossible to replicate the findings of the paper. • BaCaDI: Similar reasons to DCDI. • LIT: Similar reasons to CSIvA. E.3 Metrics We design different reasonable metrics for two different tasks, i.e.,I-CPDAG discovery and intervention target detection. The detailed introduction is as follows: I-CPDAG discovery. We first calculate Structural Hamming Distance (SHD) atI-CPDAG level. Specifically, SHD is computed between the learnedI-CPDAG( ˆ G) and ground truthI-CPDAG(G), i.e., the smallest number of edge additions, deletions, direction reversals and type changes (directed vs. undirected) to convert the outputI-CPDAG to ground truthI-CPDAG. As is shown in Table E.4, SHD is equal to the sum of the number of ✗s in the table. Different graphs may lead to different causal inference statements and different intervention distributions. To quantify this discrepancy, we use the (pre-) distance betweenI-CPDAG( ˆ G)s, known as Structural Intervention Distance (SID). The SID is solely based on graphical criteria and quantifies the closeness between two causal graphs according to their corresponding causal inference statements. Thus, it is highly suitable for assessing graphs used for the computation of interventions. Formally: SID : ( ˆ G,G) 7→ #(i,j),i̸= j| the intervention distribution from i to j is falsely estimated by ˆ G with respect to G Go to Appendix Index37 Table E.4 SHD calculation details. in Predict→ in Ground-Truth↓ identifiable (directed) unidentifiable (undirected)missing in skeleton rightwrong identifiable➀ ✓➁ ✗➂ ✗➃ ✗ unidentifiable➄ ✗➅ ✓➆ ✗ nonexist➇ ✗➈ ✗➉ ✓ For more detailed calculation algorithms, we recommend readers to refer to (Peters and Bühlmann, 2015) carefully. F1-score is then calculated based on the identifiable edges ofI-CPDAG( ˆ G) andI-CPDAG(G), where the accuracy (precision) is equal to True Positive Rate (TPR) and the recall (recall) is equal to 1 - False Discovery Rate (FDR). Details about the specific calculation can also refer to Table E.4: Precision = TPR = ➀ ➀ + ➁ + ➂ + ➃ , Recall = 1− FDR = ➀ ➀ + ➁ + ➄ + ➇ , F1-Score = 2× Precision× Recall Precision + Recall , Intervention Target Detection. We then use the F1-Score metric to measure the detection of intervention targets. Let the set of all intervention family edges be denoted as ˆ I, and all predicted intervention target edges be denoted as I. Then we can formally define: Precision = #(I∩ ˆ I) #I , Recall = #(I∩ ˆ I) # ˆ I , F1-Score = 2× Precision× Recall Precision + Recall . E.4 Experiments Setting Computational Resources. For conducting experiments in this work, we employed the IS-MCMC algorithm to pre-generate training data for each test dataset. Subsequently, we trained two types of models, namely the skeleton model and the orientation model, based on the data. The skeleton model comprises multiple cascaded models, with a default maximum of 4 order. Our benchmark environment consists of a Linux server equipped with a 1×AMD EPYC 7763 128-Core Processor CPU (512GB memory) and 4×NVIDIA RTX A6000 (48GB memory) GPUs. To carry out benchmark testing experiments, all baselines are set to run for a duration of 12 hours by default, with specific timings contingent upon the method. It is noteworthy that our method does not necessitate GPU computation, although optimization may further expedite processes. Other methods, depending on their implementation choices, may opt for GPU acceleration. Basic Configurations. As shown in Table E.5, we provide the basic parameter configurations for ourTICL, along with the fundamental details of the synthetic data and test data used in various types of experiments in this paper. Furthermore, we provide some basic conceptual explanations. The training data refers to the pre-sampled synthetic data, which includes pairs of instances used for training. These pairs consist of re-parameterized synthetic graphs and sample data tables obtained through forward sampling based on conditional probability tables. As for the test data, we directly retrieve the causal graph structure Go to Appendix Index38 Table E.5 Basic configuration of our experiments and methods CategoryDetailHyperparameters or Settings Training Data #Training Instances 400 (default) 100 ∼ 1600 (others, see 4.2) #Simples / per Instance 10k (default) 1k ∼ 20k (others, see 4.2) Testing Data #Observation Sample 10k (default) 2k ∼ 10k (others, see F.2.2) #Intervention Sample 10k (default) 1k ∼ 10k (others, see F.2.1) Intervention Experiment Intervention Type Single + Soft + Unknow (default) Single + Hard + Unknow (others, see F.1.1) Multiple + Soft + Unknow (others, see F.1.2) #Intervention Exp. #Var. × 0.2 (default) #Var. × 0.0, 0.4, 0.8 (others, see F.3) Initial Graph SeedProxy Algorithm JCI-BLIP (default) Random Graph (others, see 4.2) Model threshold Skeleton Threshold0.5 ∼ 0.75 Orientation Threshold0.1 and parameterized conditional probability tables from bnlearn. Then, the data obtained directly through forward sampling is referred to as observed samples. Otherwise, we can conduct various types of intervention experiments, including single-node intervention or multi-node intervention (i.e., intervening on one or multiple nodes at a time), soft intervention or hard intervention (i.e., whether to eliminate dependencies of parent nodes; for soft intervention, we replace probability distributions sampled from the Dirichlet distribution with parametersα∼ U[0.2,1.0]), and known intervention or unknown intervention targets (i.e., whether we have prior knowledge of intervention target nodes). For the starting point of IS-MCMC sampling, we default to selecting JCI+BLIP as the proxy algorithm to obtain the initial graph seed. We use thexgb.XGBClassifier() API provided byScikit-learn 3 as the classifier for both skeleton and orientation models. Different threshold parameters are set for skeleton identification and orientation identification, while all other hyper-parameters are set to default values. F Additional Experiments (RQ4) F.1 Effect of Intervention Type F.1.1 Prefect Interventions We follow the default parameter settings, replacing soft interventions with hard interventions, i.e., perfect interventions, to eliminate the dependency of the parent node. The results ofI-CPDAG structure discovery are shown in Table F.1, and the results of intervention target identification are shown in Table F.2. The results indicate that our method consistently outperforms on almost all datasets and metrics. F.1.2 Multiple Interventions Next, we consider conducting multi-node interventions, where for each intervention experiment, we randomly select 1 to 3 nodes for intervention, with the selection of nodes being without replacement, while the rest follow the default parameter settings. The results ofI-CPDAG structure discovery are shown in Table F.3, and the results of intervention target identification are shown in Table F.4. The results indicate that our method still maintains a consistent lead in almost all datasets and metrics. Additionally, it is worth noting that in the case of interventions on multiple nodes simultaneously, we observed that the performance does not necessarily improve compared to single-node interventions, which may be due to the complexity introduced by intervening on multiple nodes at the same time. 3 https://scikit-learn.org/stable/ Go to Appendix Index39 Table F.1 Performance comparison of I-CPDAG under perfect intervention (F1 Score ↑ / SHD ↓ / SID ↓). Datasets #nodes / #edges MetricTICL JCI-GOLEM JCI-PC JCI-BLIP JCI-HC AVICI ∗ ENCO ∗ IGSP ∗ GIES UT-IGSP Earthquake 5 / 4 F1 Score1.000.001.001.000.000.290.671.001.001.00 SHD0600442000 SID 0130020102000 Survey 6 / 6 F1 Score1.000.400.400.730.730.000.440.910.910.91 SHD0742264111 SID 02318991716444 Asia 8 / 8 F1 Score0.860.000.400.860.120.000.220.77 0.710.43 SHD2148213810346 SID 12363412473428101824 Sachs 11 / 17 F1 Score0.790.000.440.790.71/0.250.56 0.510.56 SHD6201768/36131813 SID 2761302729/35283432 Child 20 / 25 F1 Score0.760.240.250.750.65/0.070.49 0.380.53 SHD1034451116/105295329 SID 160301233188226/210125 140112 Insurance 27 / 52 F1 Score0.720.150.470.710.47/0.110.37 0.390.37 SHD2273642443/171779678 SID325671414390421/536442316452 Water 32 / 66 F1 Score0.580.14–0.500.30/0.24/// SHD 4277–4471/80/// SID419566–457642/508/// Table F.2 Performance comparison of Intervention Targets Detection under perfect intervention (F1 Score ↑). Datasets Earthquake Survey Asia Sachs Child Insurance Water 1122456 UT-IGSP0.400.500.440.400.360.22/ CITE1.001.001.000.570.400.55/ PreDITEr0.67/1.000.67–/ JCI-GOLEM0.671.000.000.180.220.280.21 JCI-HC0.501.000.360.670.470.340.32 JCI-BLIP1.001.000.800.800.530.440.42 JCI-PC1.000.670.670.470.470.45– TICL1.001.001.001.001.001.001.00 F.2 Effect of Sample Size F.2.1 Limited interventional data sample sizes Intervention experiments in the real world are often unrealistic or costly, such as gene knockout experiments, where only a small amount of experimental data may be available. Therefore, we further consider the case of limited intervention data sample size. Specifically, we reduce the sample size of each intervention experiment from 10,000 to 1,000 to simulate this scenario, while keeping other parameters at their default settings. The results ofI-CPDAG structure discovery are shown in Table F.5, and the results of intervention target identification are shown in Table F.8. The results indicate that most methods show a certain degree of performance decline. Nevertheless, our method still maintains a significant lead on almost all datasets and metrics, demonstrating the robustness and superiority of our TICL approach. F.2.2 Limited observational data sample sizes Similarly, we also consider the case of limited sample size of observational data. Specifically, we reduce the sample size of the observational data from 10,000 to 2,000 to simulate this scenario, while keeping the rest of the parameters at their default settings. The results ofI-CPDAG structure discovery are shown in Table F.7, and the results of intervention target identification are shown in Table F.8. The results are similar, indicating Go to Appendix Index40 Table F.3 Performance comparison of I-CPDAG under multiple interventions (F1 Score ↑ / SHD ↓ / SID ↓). Datasets #nodes / #edges MetricTICL JCI-GOLEM JCI-PC JCI-BLIP JCI-HC IGSP ∗ GIES UT-IGSP Earthquake 5 / 4 F1 Score0.600.000.551.000.251.000.401.00 SHD45406070 SID9183015060 Survey 6 / 6 F1 Score0.830.000.400.910.830.910.000.91 SHD29411161 SID103018464304 Asia 8 / 8 F1 Score0.860.080.560.950.560.460.440.77 SHD29717573 SID83317521202510 Sachs 11 / 17 F1 Score0.940.160.560.900.540.510.320.72 SHD21812312141910 SID861382148354315 Child 20 / 25 F1 Score 0.610.180.260.600.440.600.330.51 SHD1444451524254729 SID24330822925130472170117 Insurance 27 / 52 F1 Score0.730.160.360.690.520.340.450.33 SHD 2263722738818581 SID369624472398425470276453 Water 32 / 66 F1 Score0.620.20–0.600.24/// SHD3981–4271/// SID436605–478655/// Table F.4 Performance comparison of Intervention Targets Detection under multiple interventions (F1 Score ↑). Datasets Earthquake Survey Asia Sachs Child Insurance Water 33664512 UT-IGSP0.860.860.670.570.800.24/ CITE0.801.000.670.430.330.71/ PreDITEr0.800.50/0.29–/ JCI-GOLEM0.330.800.620.470.210.140.23 JCI-HC0.571.000.670.670.500.330.38 JCI-BLIP1.001.000.860.860.440.420.57 JCI-PC0.801.000.920.750.500.33– TICL0.801.001.000.861.000.861.00 the importance of observational data for causal structure identification. We also see that our method is still in a good leading position. F.3 Effect of Intervention Experiment Size Finally, we investigated the impact of different intervention trial frequencies on causal structure discovery. Specifically, we set intervention frequencies at 0%, 20%, 40%, and 80% of the target graph nodes, where 0% means inferring causal structure solely from observational data. We conducted experiments on the Child dataset (20 nodes) using seven methods, including ours, and the results of I-CPDAG structure discovery at different intervention frequencies are shown in Figure F.1. The experiments align with intuition, showing that with an increasing number of experiments, all methods exhibit improvements in various metrics, with our method standing out more prominently. G More Related Work Causal Discovery From Observational Data. Traditional causal discovery algorithms infer causal structures from static observational data and fall into four types: constraint-based, score-based, gradient- based, and function causal model methods. Constraint-based algorithms, such as PC (Spirtes and Glymour, 1991), FCI (Spirtes et al., 2000), and PC-Stable (Colombo et al., 2014), utilize conditional independence Go to Appendix Index41 Table F.5 Performance comparison ofI-CPDAGunder limited interventional data sample sizes (F1↑/ SHD↓/ SID ↓). Datasets #nodes / #edges MetricTICL JCI-GOLEM JCI-PC JCI-BLIP JCI-HC AVICI ∗ ENCO ∗ IGSP ∗ GIES UT-IGSP Earthquake 5 / 4 F1 Score1.000.000.891.001.000.000.671.001.001.00 SHD0710142000 SID 012000102000 Survey 6 / 6 F1 Score0.730.000.330.730.000.000.600.910.910.91 SHD21162753111 SID930229301714444 Asia 8 / 8 F1 Score0.800.200.330.800.330.000.700.77 0.710.77 SHD311731186343 SID838381934379101810 Sachs 11 / 17 F1 Score0.740.000.440.210.380.080.170.67 0.320.62 SHD7231715131741112112 SID31624255486047244030 Child 20 / 25 F1 Score0.800.160.250.750.540.560.140.51 0.470.49 SHD10413910201696283628 SID168315277185249188192128125135 Insurance 27 / 52 F1 Score0.640.060.370.600.570.170.120.30 0.530.38 SHD256566253751162825875 SID395677544404423625539451246438 Water 32 / 66 F1 Score0.470.19–0.400.45/0.47/// SHD4771–5152/53/// SID496545–540398/473/// Table F.6 Performance comparison of Intervention Targets Detection under limited interventional data sample sizes (F1 ↑). Datasets Earthquake Survey Asia Sachs Child Insurance Water 1122456 UT-IGSP0.331.000.571.000.360.27/ CITE1.000.670.671.000.400.22/ PreDITEr 0.671.00/0.400.33–/ JCI-GOLEM0.001.000.000.270.150.110.00 JCI-HC1.000.670.500.800.400.370.36 JCI-BLIP1.001.000.670.670.400.440.40 JCI-PC1.001.000.670.800.440.32– TICL1.001.001.001.001.000.911.00 relations inferred from the faithfulness assumption to establish graphical separations. Coupled with reliable conditional independence testing, these methods handle diverse data distributions and causal relations but are limited to identifying partial DAG. Score-based algorithms search for the optimal DAG under combination constraints using predefined scoring functions. Examples include GES (Chickering, 2002b), hill climbing (Koller and Friedman, 2009), and integer programming (Cussens, 2012). Gradient-based methods extend score-based approaches by transforming discrete search into continuous equality constraints. For example, NOTEARS (Zheng et al., 2018) utilizes algebraic features of DAGs for a smooth global search, but assumes linear relations. GAE (Ng et al., 2019), GraN-DAG (Yu et al., 2019) and SG-MCMC (Annadani et al., 2023) extend this to non-linear functions. Other models relax constraints (Ng et al., 2020) and use reinforcement learning (Wang et al., 2021) as a complement. Function causal model methods like LiNGAM (Comon, 1994) and ANM (Hoyer et al., 2008) describe causal relations in specific functional forms, differentiating between DAGs within the same equivalence class with additional assumptions on data Go to Appendix Index42 Table F.7 Performance comparison ofI-CPDAGunder limited observational data sample sizes (F1↑/ SHD↓/ SID↓). Datasets #nodes / #edges MetricTICL JCI-GOLEM JCI-PC JCI-BLIP JCI-HC AVICI ∗ ENCO ∗ IGSP ∗ GIES UT-IGSP Earthquake 5 / 4 F1 Score0.860.000.861.001.000.170.571.00 0.001.00 SHD1610083040 SID1161001320160 Survey 6 / 6 F1 Score0.730.000.400.730.290.000.250.50 0.360.50 SHD 21242675444 SID1130189292218112011 Asia 8 / 8 F1 Score0.770.000.290.670.530.380.710.77 0.400.77 SHD 41074584363 SID18463918313414104210 Sachs 11 / 17 F1 Score0.450.090.630.210.360.270.250.65 0.350.47 SHD 1217131514142992212 SID 49542455535745183739 Child 20 / 25 F1 Score 0.750.200.170.620.75/0.110.69 0.380.65 SHD 1044501711/110165317 SID 181290258266186/2206213999 Insurance 27 / 52 F1 Score 0.610.200.380.620.52/0.150.54 0.410.56 SHD 2969643038/127408939 SID 424630538451453/540364347428 Water 32 / 66 F1 Score0.500.21–0.490.35/0.40/// SHD4481–4661/57/// SID507620–510532/527/// Table F.8 Performance comparison of Intervention Targets Detection under limited observational data sample sizes (F1 ↑). Datasets EarthquakeSurveyAsiaSachsChildInsuranceWater 1122456 UT-IGSP0.330.670.400.670.810.24/ CITE/1.00/1.000.330.38/ PreDITEr/1.00//0.40–/ JCI-GOLEM0.000.000.000.270.150.270.20 JCI-HC1.000.670.440.500.570.370.30 JCI-BLIP1.001.000.670.670.440.440.37 JCI-PC1.000.670.400.570.290.30– TICL1.000.671.001.000.860.591.00 #Var. * 0.0#Var. * 0.2#Var. * 0.4#Var. * 0.8 Number of Intervention Experiments 0 10 20 30 40 50 SHD Score UT-IGSP GIES IGSP JCI-HC JCI-BLIP TICL #Var. * 0.0#Var. * 0.2#Var. * 0.4#Var. * 0.8 Number of Intervention Experiments 0 50 100 150 200 250 SID Score #Var. * 0.0#Var. * 0.2#Var. * 0.4#Var. * 0.8 Number of Intervention Experiments 0.3 0.4 0.5 0.6 0.7 0.8 0.9 F1 Score Figure F.1 Exploration of the impact of different size of intervention experiments on performance Go to Appendix Index43 distribution or function classes. Nevertheless, the identifiability of structures obtained solely from observational data without specific assumptions is theoretically limited. Table G.1 Comparison of different methods on causal discovery. MethodsLearning ParadigmProblem Space Intervention Targets Test-time Training Unknown Known Hard Soft Single Multiple DAG-GFlowNet (Deleu et al., 2022) Unsupervised Generative Learning Interventional (Passive)✗✓ ✓✗ ✓✗ GFlowCausal (Li et al., 2022)Unsupervised Generative LearningObservational✗ CORE (Sauter et al., 2024)Reinforcement LearningInterventional (Active)✗✓ ✓✗ ✓✗ meta-CGNN (Ton et al., 2021)Meta-LearningObservational (Bivariate)✗ DAG-GNN (Yu et al., 2019)Unsupervised LearningObservational✗ TICL (Our)Supervised LearningInterventional (Passive) ✓ ✓ ✓ ✓ Causal Discovery From Interventional Data. Despite the potential identifiability of Directed Acyclic Graphs (DAGs) under sufficient i.i.d. interventions and certain assumptions, limited literature explores this problems due to the diversity in types and strategies of intervention data. GIES (Hauser and Bühlmann, 2012) extends the GES algorithm to intervention settings by leveraging the similarity between observed causal graphs and their corresponding intervened graphs. It follows the same two-step approach as the original process, utilizing forward and backward stage traversal through the search space until achieving local maximum scores. Focusing on perfect interventions, IGSP (Wang et al., 2017; Yang et al., 2018) introduces a greedy sparse permutation method, offering an extension to general interventions. This involves optimizing the scoring function through a greedy approach and guiding permutation-based strategies for traversing the I-MEC space. Subsequently, UT-IGSP (Squires et al., 2020) partially addresses unknown target intervention scenarios. Recently, ENCO (Lippe et al., 2022), based on the concept of invariance, proposes a continuous optimization method to learn causal structures from intervention data, modeling the existence of edges and edge direction as separate parameters. Nonetheless, all of these methods typically focus on a limited set of intervention scenarios and fail to provide a unified approach for effectively addressing the problem of causal discovery from interventions. Table G.2 Comparison of IS-MCMC and Other MCMC Methods FeatureOther MCMC sampling methods (Ellis and Wong, 2008; Choi et al., 2020)IS-MCMC (Our) Sampling SpaceStandard graph space (for posterior estimation)Intervention-augmented graph space (JCI framework) Intervention HandlingCannot handle unknown intervention targetsHandles unknown intervention targets UsagePosterior inferenceGenerates diverse ⟨G i ,D i ⟩ training pairs for T ConstraintsStandard graph constraintsEnforces intervention constraints ObjectivePosterior estimationData generation for adaptation Supervised Causal Discovery. Supervised causal discovery aims to pre-access synthetic datasets of causal relations and learn causal directions in a supervised manner. Early research focused on pairwise relations, including RCC (Lopez-Paz et al., 2015) and MRCL (Hill et al., 2019). For multivariate causal learning, DAG-EQ (Li et al., 2020), based on permutation-equivariant edge models, investigated supervised causal discovery using synthetic data from linear causal models, showing promising results. Subsequently, SLdisco (Petersen et al., 2023) address shortcomings in sample size and graph density settings. It trained on synthetic linear Gaussian data to learn equivalence classes of causal graphs from observational datasets. Recently, CSIvA (Ke et al., 2023b) devised a transformer architecture with permutation invariance, extending the supervised learning paradigm to incorporate intervention data for increased flexibility. Concurrently, AVICI (Lorch et al., 2022), leveraging amortized variational inference optimization, circumvented the challenges of structural search to predict causal structures directly from the provided dataset. Considering that the aforementioned methods only evaluate a limited set of real causal model data and may not effectively capture fundamental causal information, such as persistent dependencies (Yu et al., 2016) and directional asymmetry. ML4S (Ma et al., 2022) and ML4C (Dai et al., 2023), operating under the supervised paradigm for skeleton and direction learning, achieved progressively significant identifiability. However, causal discovery for various types of intervention data scenarios remains inadequately addressed in the supervised paradigm. Our approach represents the first unified solution for supervised intervention causal discovery. In addition, we have also added two comparison tables, G.1 and G.2, to illustrate the conceptual and formal differences between some existing methods and our TICL method and IS-MCMC strategy with examples. Go to Appendix Index44 H More Discussion / Limitation / Future Work H.1 Supervised Causal Discovery We advocate for causal discovery from intervention data in supervised paradigm. Identifiable causal relation- ships between causal graphs can be treated as ground truth labels, and corresponding data can be obtained through forward sampling. Importantly, this sampling strategy is cost-free, allowing us to obtain an infinite number of samples for training a model to compress the mapping between data and causal relationships. As the main study on data quantity, typically, the more data, the better performance. Furthermore, although we can freely generate parameterized causal graphs, for specific systems of interest, such as the aforementioned study on data quality, we consider training data derived using posterior distributions to be useful for causal discovery of the data of interest. H.2 Self-augmentation vs. Pre-training We noticed that all the SCL methods in the literature, during the training phase, integrate or generate datasets containing various causal mechanisms and their related data samples through simulators (if available), however, such generation mechanisms are entirely dependent on manual specification. Ideally, it is preferable to predefine rich synthetic graph structures (possible examples include ER / SF / Low Rank / Random graphs, etc.) and synthetic function types (including linear, quadratic, nonlinear, etc.) as much as possible. We refer to the process of synthesizing data during the training phase and then training them together as "pre-training", similar to today’s language model pre-training, using massive training data to achieve generalization in different domains. However, causal structure learning is a high-risk field, where people often pay more attention to performance. The potential danger of this pre-training method is that if the causal mechanisms of the system of interest have never been seen before, there may be out-of-distribution generalization issues. Therefore, in this paper, we advocate acquiring training data more relevant to the test data after accessing the test data, providing high-quality training data by observing the posterior estimates of causal graphs to avoid the "domain shift" problem. We call this self-augmentation method, similar to in-context learning, adjusting on specific testing instance of interest in order to achieve high performance that is usable. H.3 Discrete vs. Continuous Data In this paper, we mainly focus on causal structure identification from discrete data. The Markov completeness theorem states that for discrete or linear Gaussian data, we can only identify causal graphs up to their CPDAG. For continuous data satisfying linear non-Gaussian mechanisms or additive noise assumptions, we can identify more causal directions. However, our goal is to properly handle intervention data by using a supervised learning paradigm to learn identifiable causal structures. Therefore, more considerations may be needed to design corresponding learning tasks for continuous data. Nevertheless, in practice, our method can be naturally extended to continuous data by utilizing standard algorithms from existing literature, such as the Hilbert-Schmidt Independence Criterion (Gretton et al., 2007) to compute conditional dependencies of numerical variables. Despite this, more attention may need to be paid, which we consider as future work. I Broader Impacts We proposeTICL, a new method for causal discovery from data with unknown intervention targets. Specifically, we introduce the concept of self-augmentation, using test data to obtain high-quality usable synthetic training data and identifying causal structures under the paradigm of supervised learning. Additionally, we advocate for joint causal inference as a fundamental framework for handling different intervention settings. Our numerical experiments demonstrate that our method outperforms existing intervention causal discovery methods on a wide range of datasets and metrics. However, despite the existence of supervisory identifiability, due to our adoption of traditional manual feature engineering, certain features may incur negative returns in specific situations for the final causal identifiability judgment. Therefore, future work will explore extending it to fully end-to-end neural network models, which may be more suitable by automatically capturing features instead of Go to Appendix Index45 manual operations. Finally, as a potential benefit of the design of joint causal reasoning, we can naturally identify unknown intervention targets while discovering causality, enabling scientists to apply these methods to automated experiment design and scientific discovery. Go to Appendix Index46