Paper deep dive

Amortized Counterfactual Prediction with Graph Neural Networks

Hendrik Stumpf, Francesco Quinzan

Year: 2026Venue: arXiv preprintArea: cs.LGType: PreprintEmbeddings: 117

Intelligence

Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 92%

Last extracted: 10/13/2026, 2:58:23 AM

Summary

The paper introduces the Counterfactual Graph Neural Process (CFGNP), an end-to-end framework for counterfactual inference that operates directly on observational data without requiring prior knowledge of the underlying causal graph. By leveraging Graph Neural Networks (GNNs) and Neural Processes, CFGNP exploits the dual-graph structure to restrict message passing between factual and counterfactual variables. The model demonstrates superior predictive performance compared to existing methods on synthetic and real-world benchmarks, including the German Loan dataset, and exhibits strong generalization capabilities across varying graph sizes.

Entities (7)

CFGNP · model · 100%German Loan benchmark · dataset · 95%Graph Neural Networks · method · 95%Causal Graph · concept · 90%Dual Graph · concept · 90%Neural Processes · method · 90%Structural Causal Model · concept · 85%

Relation Signals (6)

CFGNP → eliminatesneedfor → Causal Graph

confidence 95% · answers counterfactual queries directly from observational data without first estimating a causal graph

CFGNP → evaluatedon → German Loan benchmark

confidence 95% · including on the German Loan benchmark, where it reduces MSE

CFGNP → uses → Graph Neural Networks

confidence 95% · CFGNP uses Graph Neural Networks (GNNs) to exploit the dual-graph structure

CFGNP → basedon → Neural Processes

confidence 90% · We approach this problem through meta-learning using Neural Processes (NPs)

CFGNP → outperforms → Existing approaches

confidence 90% · CFGNP achieves stronger predictive performance than existing approaches

Causal Graph → causeserrorpropagation → counterfactual predictions

confidence 85% · errors in the estimated causal graph can propagate to downstream counterfactual predictions

Cypher Suggestions (0)

No Cypher suggestions yet.

Abstract

Abstract:Counterfactual inference asks how an observed outcome might have differed under different circumstances, such as whether a patient would have recovered had they received a different treatment. Such questions arise whenever decisions depend on the outcomes of alternatives that were not observed. Answering counterfactual queries generally requires knowledge of the underlying causal structure, which is often unavailable in practice. A common approach is therefore to first estimate the causal graph from observational data and then perform counterfactual inference. However, errors in the estimated causal graph can propagate to downstream counterfactual predictions. To address this challenge, we introduce the Counterfactual Graph Neural Process (CFGNP), an end-to-end approach that answers counterfactual queries directly from observational data without first estimating a causal graph. CFGNP uses Graph Neural Networks (GNNs) to exploit the dual-graph structure by restricting cross-world message passing to corresponding factual--counterfactual variables. Across synthetic and image-based benchmarks, CFGNP achieves stronger predictive performance than existing approaches, including on the German Loan benchmark, where it reduces MSE from $4.51 \times 10^{-2}$ to $0.19 \times 10^{-2}$. Furthermore, on synthetic causal graphs, a model trained only on 10-node graphs maintains strong predictive performance on graphs with up to 30 nodes, without retraining.

PDF

Open source PDF →Open local PDF →

Trouble viewing inline? Open PDF directly →

Full Text

114,184 characters extracted from source content.

Expand or collapse full text

Amortized Counterfactual Prediction with Graph Neural Networks Hendrik Stumpf† Affiliation: Department of Computer Science, University of Oxford Francesco Quinzan† Affiliation: † autohrs. Emails: hendrik.stumpf@cs.ox.ac.edu, francesco.quinzan@eng.ox.ac.uk Affiliation: Department of Engineering, University of Oxford Abstract Counterfactual inference asks how an observed outcome might have differed under different circumstances, such as whether a patient would have recovered had they received a different treatment. Such questions arise whenever decisions depend on the outcomes of alternatives that were not observed. Answering counterfactual queries generally requires knowledge of the underlying causal structure, which is often unavailable in practice. A common approach is therefore to first estimate the causal graph from observational data and then perform counterfactual inference. However, errors in the estimated causal graph can propagate to downstream counterfactual predictions. To address this challenge, we introduce the Counterfactual Graph Neural Process (CFGNP), an end-to-end approach that answers counterfactual queries directly from observational data without first estimating a causal graph. CFGNP uses Graph Neural Networks (GNNs) to exploit the dual-graph structure by restricting cross-world message passing to corresponding factual–counterfactual variables. Across synthetic and image-based benchmarks, CFGNP achieves stronger predictive performance than existing approaches, including on the German Loan benchmark, where it reduces MSE from 4.51×10−24.51× 10^-2 to 0.19×10−20.19× 10^-2. Furthermore, on synthetic causal graphs, a model trained only on 10-node graphs maintains strong predictive performance on graphs with up to 30 nodes, without retraining. 1 Introduction Understanding cause-and-effect relationships is a central goal in many scientific disciplines. In medicine, researchers may wish to determine whether a treatment improves patient outcomes [Stovitz and Shrier, 2019, Burotto et al., 2014]. In economics, policymakers often aim to understand how changes in policy affect employment or inflation [Heckman and Pinto, 2022]. Similar questions arise in biology, social sciences, and machine learning [Morgan and Winship, 2007, Sanderson et al., 2022, Schölkopf et al., 2021]. Especially with the increased deployment of machine learning systems in settings where interpretability is crucial, understanding cause-and-effect relationships becomes more important [Li and Han, 2023]. In all of these settings, the objective is not merely to observe statistical associations, but to understand how outcomes would change if different actions or conditions were imposed. Causal inference [Pearl, 2009] provides a framework for studying such questions. Unlike conventional statistical learning, causal inference seeks to uncover cause-effect relationships between variables, rather than statistical associations. The framework formalises this with the structural causal model (SCM), which specifies how variables are generated, and how they causally affect one another. Within the SCM framework, causal questions can be distinguished by the type of reasoning they require. Pearl [Pearl, 2009] differentiates between three levels: associational, interventional, and counterfactual queries. Associational queries concern relationships observed in data, while interventional queries ask how outcomes change when a variable is actively manipulated [Pearl, 1996]. Counterfactual queries go one step further by asking what would have happened in an alternative scenario, given what was observed in the factual world [Shpitser and Pearl, 2012]. Counterfactual queries are particularly important because they allow reasoning about individualized outcomes and explanations, rather than only population-level effects [Nguyen et al., 2020]. This makes them harder to solve than the other two query types as we additionally have to take the alternative hypothesis into account. Traditional approaches to counterfactual inference typically rely on substantial prior knowledge about the underlying cause-effect relationships between variables [Xia et al., 2022, Rasal et al., 2025, Ribeiro et al., 2025, Nasr-Esfahany et al., 2023]. In practice, this often means assuming access to the causal graph, which describes these relationships. Sometimes, less information, such as the topological node ordering of the causal graph, is required [Le et al., 2024]. However, in many real-world settings, such information is unavailable or highly uncertain. This is particularly the case in domains such as medicine, biology, and the social sciences, where causal relationships may be obscured by latent confounders, incomplete observations, and complex dependencies among variables [Weberpals et al., 2024, Ruiz-Perez et al., 2024, Loh, 2025]. Hence, if there is no prior knowledge about the causal graph, we must first use observational samples to estimate the causal graph [Spirtes et al., 2001] and then use the resulting causal graph with the existing methods for counterfactual inference. However, this has two drawbacks; first, the search space of causal graphs scales super-exponentially in the number of variables [Zheng et al., 2018] and second, choosing the wrong causal graph amplifies errors and may cause the counterfactual to be incorrect [Karimi et al., 2020b]. To avoid explicitly estimating the causal graph, we instead condition the counterfactual prediction directly on observational samples. This allows the model to extract information about the underlying causal structure that is relevant to the counterfactual query, without committing to a single estimated graph. We approach this problem through meta-learning using Neural Processes (NPs) [Garnelo et al., 2018], which learn to map a set of observations to a predictive distribution. Building on Graph Neural Processes [Carr and Wingate, 2019], we further incorporate the structural constraints imposed by the dual graph into the prediction. Based on this approach, we introduce the Counterfactual Graph Neural Process (CFGNP), an end-to-end framework for counterfactual inference under unknown causal structure. An overview of our setting is shown in Figure 1. Our contributions. • We introduce the Counterfactual Graph Neural Process (CFGNP), an end-to-end framework for learning counterfactual distributions without prior knowledge of the underlying causal graph (Section 4). • We develop a Graph Neural Process architecture that conditions counterfactual predictions on observational samples (Section 4.1). • We demonstrate that CFGNP outperforms existing counterfactual inference methods on both synthetic and real-world data (Sections 5.1 and 5.2). • We evaluate CFGNP on counterfactual queries involving features that are unseen during training and observe strong performance in this setting (Section 5.2). Figure 1: Overview of the task solved by CFGNP. Unknown variables are in light grey and known variables in dark grey. Blue edges represent missing knowledge about the causal graph. The counterfactual is predicted based on the counterfactual query and additional observational samples. 2 Related Work Causal Discovery. For causal inference tasks such as evaluating counterfactual queries, most existing methods rely on explicit knowledge of the causal graph. Estimating this graph from associational data is known as causal discovery [Verma and Pearl, 2022]. The primary goal of causal discovery is to uncover true cause-and-effect relationships by learning the underlying structural dependencies rather than relying on statistical correlations alone [Pearl, 2009, Peters et al., 2017]. Classical methodologies encompass constraint-based approaches, such as the PC algorithm [Spirtes et al., 2001], and score-based frameworks, such as Greedy Equivalence Search [Chickering, 2002]. More recently, modern techniques including LiNGAM [Shimizu et al., 2006], NOTEARS [Zheng et al., 2018], and DAG-GNN [Yu et al., 2019] have advanced the field by introducing functional form assumptions, continuous optimization, and deep learning architectures. Nevertheless, key challenges remain, such as unobserved confounding [Kaltenpoth and Vreeken, 2023, Diepen et al., 2023], complex nonlinearities [Li et al., 2024], and computational scalability [Amin and Wilson, 2024, Cheng et al., 2024]. Addressing these challenges continues to make causal inference an active area of research. Counterfactual Inference. Existing approaches attempt to learn counterfactual distributions under varying degrees of prior structural knowledge [Khemakhem et al., 2021, Karimi et al., 2020a, Caron et al., 2022]. For instance, Le et al. [2024] and Javaloy et al. [2023] use normalizing flows that require only a topological ordering of the graph nodes, whereas most alternatives assume the full causal graph. Within deep learning for counterfactual inference, Pawlowski et al. [2020] use normalizing flows to parameterize structural causal models, while recent flow-matching approaches similarly require full graph specification [Ribeiro et al., 2025]. Other methods explicitly exploit relational structure through Graph Neural Networks [Zečević et al., 2021, Sanchez-Martin et al., 2021]. This reliance on a known, estimated, or otherwise fixed causal structure can be restrictive [Wu et al., 2025, Caron et al., 2022]: the space of DAGs grows super-exponentially with the number of variables [Zheng et al., 2018], and errors in the assumed causal structure can propagate to downstream counterfactual predictions [Karimi et al., 2020b]. While end-to-end learning under uncertain causal structure has been explored for interventional inference [Geffner et al., 2022, Dhir et al., 2025], an analogous approach has not yet been established for counterfactual inference, where additional structural constraints from the dual graph can be exploited [von Kügelgen et al., 2022]. Image Counterfactuals. Several works focus specifically on image-based counterfactual generation and analysis [Pawlowski et al., 2020, Rasal et al., 2025, Luedtke and Fukumizu, 2025, Ma et al., 2026]. Among these, DeepBC [Kladny et al., 2024] performs counterfactual inference through latent exogenous variables under a known causal graph. Other recent approaches based on diffusion models [Rasal et al., 2025] and optimal transport [Ribeiro et al., 2025] similarly assume that the causal graph is known. In contrast, VCI [Wu et al., 2025] does not take an explicit causal graph as input, but instead assumes a prescribed treatment–covariate–outcome structure and uses a variational formulation to generate counterfactuals without requiring counterfactual training pairs. Neural Processes. To address these limitations, our approach builds on Neural Processes (NPs) [Garnelo et al., 2018], which have been successfully applied across a range of domains including Bayesian optimization and other scientific applications [Kim et al., 2019, Foong et al., 2020, Bruinsma et al., 2023]. NPs have recently been explored in causal settings, including treatment-effect estimation and Bayesian causal inference [Bynum et al., 2026, Dhir et al., 2025]. However, these works address related causal inference tasks rather than counterfactual inference. In contrast, our approach extends this framework to counterfactual inference, for which we additional leveraging GNNs to capture the structural relationships between factual and counterfactual variables [von Kügelgen et al., 2022]. Graph Neural Processes extend NPs to relational and graph-structured data [Carr and Wingate, 2019, Liang and Gao, 2021, Hu et al., 2023]. By combining these ideas, we directly map observational data to the target counterfactual distribution using GNNs over the dual-graph structure, without first estimating a causal graph. 3 Background We now introduce notions of causal inference, neural processes, and graph neural networks. Notation. Unstarred and starred variables denote the factual and counterfactual worlds, respectively. Let z∈ℝn×ddataz ^n× d_data be the factual observed unit, T∗=t∗T^*=t^* a counterfactual condition, and Y∗⊆VY^* V the target variables. The observational context is O=Oii=1NobsO=\O_i\_i=1^N_obs, consisting of i.i.d. samples from the same SCM ℳM, where n and NobsN_obs denote the number of endogenous variables and context samples, respectively. We write p⁡(Y∗∣O,z,T∗=t∗)p(Y^* O,z,T^*=t^*) when ℳM is informed by the observational context O, and p⁡(Y∗∣ℳ,z,T∗=t∗)p(Y^* ,z,T^*=t^*) when emphasizing the model. If ℳM is unspecified, we assume the true model. 3.1 Structural Causal Models and Causal Inference Structural Causal Model. A structural causal model (SCM) describes the causal relationships among a set of endogenous variables V=V1,…,VnV=\V_1,…,V_n\ through a collection of structural equations =1,…,nf=\f_1,…,f_n\. Each variable VkV_k is determined by its direct causes Pak⊆V∖VkPa_k V V_k and an exogenous noise variable UkU_k according to Vk=k(Pak,Uk).V_k=f_k(Pa_k,U_k). The exogenous variables U=U1,…,UnU=U_1,…,U_n capture factors outside the model and are assumed to be mutually independent. The causal relationships encoded by the structural equations can be represented by a directed and acyclic causal graph G, where an edge Vj→VkV_j→ V_k indicates that Vj∈PakV_j∈ Pa_k. In particular, with U coming from a probability distribution we can define associational queries (ℒ1L_1) to concern distributions p⁡(Y)p(Y) with Y⊂VY⊂ V. Observational samples OiO_i are sampled from p⁡(V)p(V). Interventional queries (ℒ2L_2) concern questions p⁡(Y∣do⁡(T=t))p(Y (T=t)) where we actively set nodes T⊂VT⊂ V to a specific value and consider the outcome Y⊂VY⊂ V. (Backtracking) counterfactual queries (ℒ3L_3) are probability distributions, with T∗,Y∗,Z⊂VT^*,Y^*,Z⊂ V and T∗∩Y∗=∅T^*∩ Y^*= , of the form P⁡(Y∗=y∗∣T∗=t∗,Z=z),P(Y^*=y^* T^*=t^*,Z=z), concerning questions such as: given that Z=zZ=z happened, what would have happened to variables Y∗Y^* if alternatively T∗=t∗T^*=t^*. In particular, we use notation X∗=x∗X^*=x^* instead of the dodo-operator to make clear that we are using backtracking counterfactuals [von Kügelgen et al., 2022] and not interventionist counterfactuals [Pearl, 2009]. While interventionist counterfactuals artificially modify the structural equations in the counterfactual world, they thereby locally violate the laws of nature. In contrast, backtracking counterfactuals keep the structural equations intact and instead consider alternative exogenous conditions that could have produced X∗=x∗X^*=x^*. Figure 2: Example of the dual graph for a counterfactual query with a causal graph of four nodes. Light grey nodes are unknown. Dual Graph. To answer backtracking counterfactuals, we need the dual graph. Figure 2 shows an example of the dual graph for backtracking counterfactuals. The dual graph consists of two copies of the same SCM. The left part, the observed factual unit, represents what actually happened, whereas the copy on the right, the counterfactual unit, represents the alternative circumstances. Both worlds are connected via the exogenous nodes, which share a probability distribution p⁡(U,U∗)p(U,U^*). We can solve the counterfactual query via the abduction–marginalisation–prediction framework [von Kügelgen et al., 2022], consisting of the following steps: 1.) Abduction: We infer the joint posterior over the exogenous nodes P(U,U∗∣Z=z,T∗=t∗)P(U,U^* Z=z,T^*=t^*). 2.) Marginalisation: We marginalise out U to obtain P⁡(U∗∣Z=z,T∗=t∗)P(U^* Z=z,T^*=t^*). 3.) Prediction: We use the SCM to predict P⁡(Y∗∣Z=z,T∗=t∗)P(Y^* Z=z,T^*=t^*). Identifiability. When the underlying causal model is only partially known, a causal query may not be uniquely determined by the available information. An ℒiL_i query is identifiable from ℒjL_j information if all models in the considered model class that agree on the available ℒjL_j information yield the same answer [Pearl, 2009]. This is relevant to our setting because counterfactual queries are not necessarily identifiable from observational data alone. 3.2 Neural Processes and Graph Neural Networks Neural Processes. Neural Processes (NPs) [Garnelo et al., 2018] are probabilistic models that learn distributions over functions from a collection of related tasks. Each task is represented by a context set C=(xi,yi)i=1NCC=\(x_i,y_i)\_i=1^N_C of observed input-output pairs and a target set G=(xj,yj)j=1NGG=\(x_j,y_j)\_j=1^N_G at which predictions are made. NPs learn the conditional predictive distribution p⁡(yG∣xG,C),p(y_G x_G,C), thereby amortizing inference across tasks by mapping context observations directly to predictive distributions. Extensions of NPs to graph-structured domains replace neural architectures with graph-based models to capture relational structure [Carr and Wingate, 2019, Ashman et al., 2024]. Graph Neural Networks. Graph Neural Networks (GNNs) learn representations of graph-structured data by iteratively propagating information between connected nodes. Let hi(k)h_i^(k) denote the representation of node viv_i at layer k, with hi(0)=xih_i^(0)=x_i. Following the message-passing framework [Gilmer et al., 2017], a GNN layer aggregates messages from the neighbourhood ⁡(i)N(i) of each node, mi(k)=AGG(k)⁡(ϕ(k)(hi(k),hj(k),eji)|vj∈⁡(i)),hi(k+1)=ψ(k)(hi(k),mi(k)),m_i^(k)=AGG^(k) ( \φ^(k)(h_i^(k),h_j^(k),e_ji)\, |\,v_j (i) \ ), h_i^(k+1)=ψ^(k)(h_i^(k),m_i^(k)), where ϕ(k)φ^(k) and ψ(k)ψ^(k) are learnable message and update functions, respectively, and AGG(k)AGG^(k) is permutation invariant. For directed graphs, the neighbourhood can be partitioned into incoming and outgoing neighbours, −(i)N^-(i) and +(i)N^+(i), allowing messages to be propagated separately along each direction. This distinction is particularly relevant for causal graphs, where incoming and outgoing edges represent causal relations in opposite directions. 4 A GNN Model for Counterfactual Inference without Causal Graph Knowledge This work focuses on predicting counterfactual distributions directly from observational data without prior knowledge about the causal graph. For each task, we observe a context set O=zii=1NobsO=\z_i\_i=1^N_obs of i.i.d. observational samples from an unknown SCM ℳ=(,)M=(G,f). Given a factual observation z and a counterfactual condition T∗=t∗T^*=t^*, our goal is to predict the counterfactual distribution p⁡(Y∗∣O,z,T∗=t∗)p(Y^* O,z,T^*=t^*) over a set of query variables Y∗Y^*. The observational context O provides evidence about ℳM. We aim to directly learn the model-averaged counterfactual predictive distribution p⁡(Y∗∣O,z,T∗=t∗)=ℳ∼p⁡(ℳ∣O,z)[p⁡(Y∗∣ℳ,z,T∗=t∗)],p(Y^* O,z,T^*=t^*)=E_M p(M O,z) [p(Y^* ,z,T^*=t^*) ], (1) thereby accounting for uncertainty over causal models without explicitly inferring or committing to a particular causal graph. Existing methods on the other hand try to estimate the causal graph first. Since we intend to approximate the counterfactual probability distribution, we use the Kullback-Leibler (KL) divergence [Kullback and Leibler, 1951] as the minimization objective. Suppose we have a list L of samples (τ,y∗)∼q(τ,y^*) q, where τ=((z,T∗,t∗,Y∗),O)τ=((z,T^*,t^*,Y^*),O) contains the factual representation z, counterfactual query, and observation pattern O, while y∗y^* contains the corresponding counterfactual target values. Let pθ(y∗∣τ)p_θ(y^* τ) denote the distribution predicted by the model with weights θ, where y∗y^* may contain multiple target variables. We then want to obtain the weights that minimize the KL divergence, hence θ∗ θ^* =argminθτ∼q[KL(p(Y∗=y∗|z,T∗=t∗)||pθ(Y∗=y∗|z,T∗=t∗,O))] = _θE_τ q [KL(p(Y^*=y^*|z,T^*=t^*)\ ||\ p_θ(Y^*=y^*|z,T^*=t^*,O)) ] =arg⁡maxθτ∼q[y∗|τ[log⁡pθ(Y∗=y∗|z,T∗=t∗,O)]]+C = _θE_τ q [E_y^*|τ[ p_θ(Y^*=y^*|z,T^*=t^*,O)] ]+C for some constant C independent of the model parameters. During training we assume access to counterfactual pairs, as we sample y∗|τ∼p⁡(Y∗=y∗|z,T∗=t∗,O)y^*|\ τ p(Y^*=y^*|z,T^*=t^*,O). Recovery of the Counterfactual Distribution. In the limit of infinite data and model capacity, minimizing the KL divergence recovers the conditional predictive distribution [Bruinsma, 2022] p⁡(Y∗∣O,z,T∗=t∗)p(Y^* O,z,T^*=t^*), with latent quantities, such as the exogenous variables, implicitly marginalized out [Gordon et al., 2019]. If the counterfactual distribution is identifiable from observational data, this predictive distribution coincides with the counterfactual distribution of the ground-truth causal model. Otherwise, the model instead recovers the model-averaged counterfactual distribution over SCMs compatible with the observations. 4.1 Approach Figure 3: Overview of CFGNP. Observational samples are first aggregated into node-wise context representations, which are combined with the observed factual unit. Message passing between factual and counterfactual unit then produces parameters of the counterfactual predictive distribution. Figure 3 illustrates the architecture of the Counterfactual Graph Neural Process (CFGNP). Given observational samples O, a factual observation z, and a counterfactual intervention x∗x^*, CFGNP predicts the corresponding counterfactual outcome distribution. The architecture consists of three stages: summarizing the observational context (middle), combining it with the factual observation (top right), and propagating the resulting representation to the counterfactual unit (bottom right). We describe each component below. Embedding Input. We project all nodes to the hidden dimension dhiddend_hidden using node-wise MLPs. For O, we apply the same MLP hOh_O across all nodes. For z and x∗x^*, we distinguish intervention nodes T∗T^* from the remaining nodes using separate MLPs: hz,inth_z,int and hz,commonh_z,common for z, and hx,inth_x,int and hx,commonh_x,common for x∗x^*. In x∗x^*, all entries outside T∗T^* are masked. Observational Context Encoding. We next incorporate the observational samples O by allowing nodes within each sample OjO_j to communicate and compressing their information into q. For each variable i, we introduce an aggregation node qiq_i, initialized to zero, and connect the corresponding node i from every observational sample to qiq_i. Since we have no prior knowledge of the causal graph, we additionally fully connect the nodes within each observational sample. The resulting graph consists of NobsN_obs disjoint fully connected n-node subgraphs and n aggregation nodes. A K-layer GAT [Brody et al., 2022] performs message passing, [O′,q′]=mθ(1)([O,q]),[O ,q ]=m_θ^(1)([O,q]), such that q′q captures the compressed observational context. We then integrate this context with the factual unit z. To keep the graph compact, we connect each qi′q_i only to its corresponding ziz_i and add a self-loop to each ziz_i. A 3-layer GAT propagates the contextual information to the factual unit, z′=mθ(2)(z,q′).z =m_θ^(2)(z,q ). Predicting Counterfactual Distribution. The final stage aims to approximate the abduction–marginalisation–prediction steps using the context-informed factual representation z′z and counterfactual condition x∗x^*. The dual graph enables these steps by coupling the factual and counterfactual worlds through connected exogenous nodes. Approximating these steps generally requires bidirectional message passing. Because we do not have prior knowledge of the causal graph, factual and counterfactual units are fully connected. Thus, each endogenous node has direct access to the information needed to reproduce its corresponding exogenous state, allowing us to remove explicit exogenous nodes and connect only corresponding factual and counterfactual variables across worlds. In this case, a unidirectional GNN can approximate the corresponding bidirectional GNN arbitrarily well (App. 1). As a design choice, we apply min⁡⌊n/2⌋,5 \ n/2 ,5\ GAT layers, scaling depth with graph size while capping it at five to limit over-smoothing [Li et al., 2018]. This yields [z′,h∗]=mθ(3)(z′,x∗)[z ,h^*]=m_θ^(3)(z ,x^*), where x∗x^* encodes t∗t^* at conditioned nodes and masks the rest. The resulting h∗h^* represents the counterfactual unit. To obtain the counterfactual predictive distribution, a node-wise MLP maps h∗h^* to its distributional parameters. We use a mixture of NcN_c Gaussians, allowing CFGNP to represent uncertainty arising from ambiguity over the underlying SCM. Training Objective. During training CFGNP minimizes the negative-log-likelihood objective, as this is equivalent to minimizing the KL divergence [Murphy, 2022]: ℒθ=−∑i=1|L|∑j=1dlog(∑k=1Ncwk(τi)j(yi,j∗|μk(τi)j,σk(τi)j)).L_θ=-Σ _i=1^|L|Σ _j=1^d ( _k=1^N_cw_k( _i)_jN (y^*_i,j| _k( _i)_j, _k( _i)_j ) ). Model Properties. The structure of CFGNP reflects two desirable properties. First, predictions should be independent of the ordering of variables. Permuting target nodes Y∗Y^* should induce the same permutation in the corresponding predictions, while permutations of the remaining nodes (Y∗)c(Y^*)^c should leave predictions for Y∗Y^* unchanged. This is provided by design of GNNs (App. .3). Second, a trained model should be applicable to graphs of different sizes and arbitrary numbers NobsN_obs of observational samples (App. .4). This enables CFGNP to generalize beyond the graph sizes and context sizes encountered during training. We note for completeness that, for n nodes, data/hidden dimensions ddata/dhiddend_data/d_hidden, NobsN_obs observations, K GNN iterations, and NCN_C mixture components, one forward pass has complexity (.2): Θ⁡(CLOSEOPENnddatadhiddenNobs+KNobsndhidden(dhidden+n)+ndhidden2+n2dhidden+ndhiddenNC). aligned \! (&n\,d_datad_hiddenN_obs+KN_obsnd_hidden(d_hidden+n)+nd_hidden^2+n^2d_hidden+nd_hiddenN_C ) aligned. 5 Results We evaluate CFGNP in terms of counterfactual prediction performance, generalization across causal graph sizes and unseen features, and computational efficiency, with additional experiments assessing its ability to recover analytically available counterfactual distributions with increasing context size (App. C.4). We also include an ablation study (App. C.5) where a normal transformer is run on the concatenated inputs. The results show that CFGNP performs much better than this. Baselines. We compare against the following baselines: ∙ VCI [Wu et al., 2025]: uses a conditional VAE with a counterfactual variational objective under a restricted causal structure. ∙ DeepBC [Kladny et al., 2024]: performs counterfactual inference through exogenous variables in a learned SCM, assuming a known causal graph (CG). ∙ VACA [Sanchez-Martin et al., 2021]: uses a variational graph autoencoder on a known CG. ∙ BGM [Nasr-Esfahany et al., 2023]: represents structural equations as invertible mappings between exogenous and observed variables given the causal structure. ∙ CFM [Le et al., 2024]: learns flow-based causal mechanisms assuming a given causal ordering. Unlike the baseline methods, CFGNP uses counterfactual pairs during training, but no causal graph. Metrics. We evaluate counterfactual predictions using the following metrics: ∙ MSE/MAE: Mean squared/absolute error (MSE/MAE) for numerical variables. ∙ Accuracy/ROC-AUC: Accuracy and ROC-AUC for categorical variables. For datasets containing both numerical and categorical variables, we report metrics separately for each variable. For synthetic experiments with multiple variables, we follow Sanchez-Martin et al. [2021] and report MSE aggregated across all variables and samples. Synthetic Datasets. For the synthetic experiments, we use the Triangle, M-shape, and German Loan benchmarks from Sanchez-Martin et al. [2021], together with randomly generated SCMs for the graph-size experiments. Triangle and M-shape have 3 and 5 nodes respectively and are evaluated with linear (LIN), nonlinear additive (NLIN), and nonlinear non-additive (NADD) mechanisms, while German Loan is a 7-node semi-synthetic seven-variable SCM inspired by the German Credit domain. Detailed definitions can be found in B.3. Further, we conduct experiments on causal graphs of increasing size where each sample comes from its own SCM. We follow Soleymani et al. [2020]. We generate the SCMs by randomly permuting the nodes into a topological order and adding each admissible edge Xi→XjX_i→ X_j with probability 12 12. Each variable X is generated from its direct causes using one of seven mechanisms f per sample, followed by posterior additive noise (App. B.4). We consider causal graphs of sizes =10,15,20,25,30S=10,15,20,25,30. Real-world Datasets. For the real-world experiments, we use the CheXpert dataset [Irvin et al., 2019], consisting of chest X-rays (App. B.2). We augment each image with four explicit variables: age, sex, disease, and AP/PA, where AP/PA denotes the radiographic viewing direction. Following Ribeiro et al. [2025], we supply the baselines with a causal structure in which age affects disease, while age, sex, disease, and AP/PA determine the image. Since ground-truth counterfactuals are unavailable [Ma et al., 2026], we evaluate whether generated images exhibit the intended counterfactual attributes using classifiers. Details on the classifiers and image-generating CFGNP variant are provided in App. B.2. We assess whether images resemble real-world images qualitatively. 5.1 Synthetic Experiments On synthetic experiments, we compare against VACA [Sanchez-Martin et al., 2021], BGM [Nasr-Esfahany et al., 2023], and CFM [Le et al., 2024], which have demonstrated strong performance in prior work. VACA and BGM require the full causal graph, whereas CFM only requires its topological ordering, allowing comparison across different levels of causal knowledge. For the graph-size experiments, we omit CFM as it was not computationally viable on larger causal graphs. Counterfactual Prediction. Table 1 shows the results for the synthetic datasets. CFGNP improves across almost all variants of the datasets and in many by a substantial amount. For instance, for German loan, CFGNP has an MSE of 0.19×10−20.19× 10^-2 compared to 4.51×10−24.51× 10^-2 of CFM. However, on the LIN variant of M-Shape, the MSE of 3.60×10−53.60× 10^-5 is higher than BGM’s with 1.46×10−51.46× 10^-5. Dataset Variant BGM VACA CFM CFGNP (Ours) MSE ↓ Triangle LIN (3.71±0.08)×10−3(3.71± 0.08)× 10^-3 (23.52±1.83)×10−2(23.52± 1.83)× 10^-2 (3.32±0.05)×10−3(3.32± 0.05)× 10^-3 (2.59±2.13)×−(2.59± 2.13)× 10^-7 NLIN (4.78±0.97)×10−6(4.78± 0.97)× 10^-6 (13.24±1.93)×10−2(13.24± 1.93)× 10^-2 (1.35±0.24)×10−5(1.35± 0.24)× 10^-5 (2.14±0.34)×−(2.14± 0.34)× 10^-7 NADD (1.80±1.80)×10−6(1.80± 1.80)× 10^-6 (22.50±2.61)×10−2(22.50± 2.61)× 10^-2 (1.54±0.40)×10−6(1.54± 0.40)× 10^-6 (2.81±1.11)×−(2.81± 1.11)× 10^-7 M-shape LIN (1.46±0.40)×−(1.46± 0.40)× 10^-5 (12.83±0.69)×10−2(12.83± 0.69)× 10^-2 (2.73±0.17)×10−5(2.73± 0.17)× 10^-5 (3.60±7.00)×10−5(3.60± 7.00)× 10^-5 NLIN (1.97±0.50)×10−5(1.97± 0.50)× 10^-5 (16.26±0.95)×10−2(16.26± 0.95)× 10^-2 (4.67±1.14)×10−5(4.67± 1.14)× 10^-5 (2.50±0.47)×−(2.50± 0.47)× 10^-6 NADD (1.97±0.07)×10−3(1.97± 0.07)× 10^-3 (41.06±0.31)×10−2(41.06± 0.31)× 10^-2 (1.89±0.03)×10−3(1.89± 0.03)× 10^-3 (1.08±0.35)×−(1.08± 0.35)× 10^-4 German Loan – (46.81±4.20)×10−2(46.81± 4.20)× 10^-2 (10.30±0.40)×10−2(10.30± 0.40)× 10^-2 (4.51±0.78)×10−2(4.51± 0.78)× 10^-2 (0.19±0.09)×−(0.19± 0.09)× 10^-2 Table 1: Mean MSE ± standard deviation achieved by BGM, VACA, CFM, and CFGNP (ours) across the Triangle, M-shape, and German Loan datasets. Lower values indicate better performance. Best scores are marked in bold. Graph-Size Experiments. For experiments with increasing graph size we use two variants of CFGNP. In one variant, CFGNP is retrained for each s∈s , while in the other, CFGNP is trained only on graphs with s=10s=10 nodes and applied to the other graph sizes without training. We refer to this model as 1010-CFGNP. Additionally, as the benchmarks require a fixed causal graph, we provide the fullest DAG, connecting each node to all topological successors. This highlights an advantage of CFGNP: when causal graphs vary across the dataset, baselines must still assume a fixed graph. Figure 4 shows MSE with increasing graph size (note the non-linear y-axis). 1010-CFGNP performs similarly to retrained CFGNP for Linear and Sum-sine, but worse for Log-sum-exp and Sqrt-sum. MSE increases for all methods, but both CFGNP variants outperform the baselines and degrade more slowly. Other mechanisms are reported in App. C.2. Figure 4: Performance (MSE) of models with regard to size of the causal graph. Models are compared on causal graphs with Linear (1.), Sum-sine (2.), Log-sum-exp (3.), Sqrt-sum (4.), and Sum-sqrt (5.) structural equations. MSE mean over distinct causal graphs. Vertical lines show 25th25^th and 75th75^th quantiles. 5.2 CheXpert Experiments We compare against VCI [Wu et al., 2025] and DeepBC [Kladny et al., 2024] which both have previously demonstrated strong performance. Again, VCI assumes treatment–covariate–outcome structure while DeepBC needs full graph knowledge. Counterfactual Prediction. Table 2 shows that CFGNP performs particularly well on age, reducing MAE from 10.8810.88 for DeepBC to 2.712.71, and consistently achieves high ROC-AUC on categorical variables. For disease, for example, accuracies are around 0.870.87, while CFGNP achieves a ROC-AUC of 0.770.77, compared with 0.750.75 for DeepBC and 0.450.45 for VCI. Generated counterfactual images are provided in Appendix C.2. Feature Metric CFGNP (Ours) DeepBC VCI Sex Accuracy ↑ 0.8452±0.06080.8452± 0.0608 0.8222±0.12200.8222± 0.1220 0.5173±0.02960.5173± 0.0296 ROC-AUC ↑ 0.9155±0.05400.9155± 0.0540 0.9204±0.07520.9204± 0.0752 0.5238±0.01220.5238± 0.0122 Age MAE ↓ 2.7136±0.73262.7136± 0.7326 10.8750±0.873210.8750± 0.8732 66.8172±62.176966.8172± 62.1769 AP/PA Accuracy ↑ 0.9664±0.01400.9664± 0.0140 0.9572±0.05520.9572± 0.0552 0.8597±0.00000.8597± 0.0000 ROC-AUC ↑ 0.9787±0.01360.9787± 0.0136 0.9658±0.05840.9658± 0.0584 0.5741±0.03420.5741± 0.0342 Disease Accuracy ↑ 0.8830±0.03680.8830± 0.0368 0.8720±0.02400.8720± 0.0240 0.8548±0.06670.8548± 0.0667 ROC-AUC ↑ 0.7710±0.05720.7710± 0.0572 0.7517±0.03270.7517± 0.0327 0.4516±0.02020.4516± 0.0202 Table 2: Performance of CFGNP, VCI, and DeepBC on the CheXpert dataset. MAE is reported for the continuous feature, while accuracy and ROC-AUC are reported for categorical features. Arrows indicate whether higher or lower values are preferred. Best values marked in bold. Out-of-Distribution Variable-Set Experiment. We evaluate CFGNP’s generalization to unseen features on CheXpert. For each OOD setting, we exclude one feature kOODk_OOD during training and evaluate on the full test set, leveraging CFGNP’s ability to handle graphs of different sizes. Table 3 shows strong generalization for sex and AP/PA, with ROC-AUC above 0.990.99, while age and disease are more challenging, with an MAE of up to 11.5311.53 and ROC-AUC between 0.610.61 and 0.640.64, respectively. Feature Metric OOD Sex OOD Age OOD AP/PA OOD Disease Sex Accuracy ↑ 0.9516±0.01080.9516± 0.0108 0.9634±0.01570.9634± 0.0157 0.9692±0.00840.9692± 0.0084 0.9627±0.00290.9627± 0.0029 ROC-AUC ↑ 0.9927±0.00370.9927± 0.0037 0.9941±0.00390.9941± 0.0039 0.9965±0.00170.9965± 0.0017 0.9964±0.00050.9964± 0.0005 Age MAE ↓ 10.7037±1.528310.7037± 1.5283 11.5291±0.095011.5291± 0.0950 9.2965±1.92059.2965± 1.9205 8.1875±0.42188.1875± 0.4218 AP/PA Accuracy ↑ 0.9917±0.00690.9917± 0.0069 0.9904±0.00850.9904± 0.0085 0.9995±0.00080.9995± 0.0008 0.9996±0.00050.9996± 0.0005 ROC-AUC ↑ 0.9998±0.00040.9998± 0.0004 0.9998±0.00030.9998± 0.0003 1.0000±0.00001.0000± 0.0000 1.0000±0.00001.0000± 0.0000 Disease Accuracy ↑ 0.6560±0.01760.6560± 0.0176 0.6678±0.01190.6678± 0.0119 0.6690±0.00730.6690± 0.0073 0.6579±0.02120.6579± 0.0212 ROC-AUC ↑ 0.6378±0.00680.6378± 0.0068 0.6296±0.01630.6296± 0.0163 0.6058±0.02640.6058± 0.0264 0.6155±0.01970.6155± 0.0197 Table 3: OOD generalization on CheXpert. Each column indicates the feature excluded during training; evaluation is performed on the full test set containing all features. Best results in bold. 6 Conclusion We propose CFGNP, an end-to-end architecture designed to learn counterfactual distributions directly from observational data without knowledge of the causal graph. CFGNP applies GNNs to enforce structural constraints on the dual graph implied by the notion of backtracking counterfactuals. By leveraging observational context samples alongside counterfactual queries, CFGNP directly approximates the Bayesian model-averaged counterfactual distribution. Evaluations on image and synthetic counterfactual queries have shown strong performance compared to baselines. Further experiments on synthetic datasets show that CFGNP retains strong performance on larger causal graphs compared with retraining from scratch, while outperforming baselines. Finally, we empirically validate that CFGNP performs adequately even on out-of-distribution features. In terms of limitations, the preceding capabilities come with a computational trade-off. The required observational context increases both training and inference costs. Moreover, if the counterfactual distribution is non-identifiable, CFGNP cannot recover the true distribution. However, this also applies to the baseline methods. In the future, one could incorporate priors about the causal graph. While the entire causal graph may not be known, an expert may have knowledge about the existence and absence of certain causal relationships. Taking this even further, one could provide simple distributions over edges, for instance obtained by amortized causal structure learning methods [Lorch et al., 2022]. AI use statement In this work, we used generative AI tools for refining writing, formulating mathematical claims, assisting in the writing of proofs, and feedback on research methodology or experiments. We have not used generative AI tools for generating synthetic data sets, helping develop theoretical models or conceptual frameworks, providing critical ingredients for proving mathematical claims, implementing methods, assisting with translation, cleaning and reformatting datasets, and interpreting results. Support for qualitative and thematic data analysis and proposing or refining hypotheses are not applicable to this work. We have reviewed all AI-assisted work. We checked that the generated formulation of the claim matches the intended claim. We take responsibility for the final content of this work, including text, claims or artifacts produced with the aid of generative AI. Ethics statement This work uses existing datasets and collects no new human-subject data. While counterfactual models may have applications in healthcare, their predictions can reflect biases and limitations in the data and should not be used for clinical decision-making without appropriate validation. Reproducibility statement To support reproducibility, we provide the source code and experimental setup. The repository can be found under https://github.com/Zelos1/CFGNP. Model architectures and hyperparameters are detailed in the appendix. References Amin and Wilson (2024) A. N. Amin and A. G. Wilson Scalable and flexible causal discovery with an efficient test for adjacency. In Proceedings of the 41st International Conference on Machine Learning, R. Salakhutdinov, Z. Kolter, K. Heller, A. Weller, N. Oliver, J. Scarlett, and F. Berkenkamp (Eds.), Proceedings of Machine Learning Research, Vol. 235, p. 1331–1358. Cited by: §2. Ashman et al. (2024) M. Ashman, C. Diaconu, J. Kim, L. Sivaraya, S. Markou, J. Requeima, W. P. Bruinsma, and R. E. Turner Translation equivariant transformer neural processes. arXiv preprint arXiv:2406.12409. Cited by: §3.2. Brody et al. (2021) S. Brody, U. Alon, and E. Yahav How attentive are graph attention networks?. arXiv preprint arXiv:2105.14491. Cited by: 2.§. Brody et al. (2022) S. Brody, U. Alon, and E. Yahav How attentive are graph attention networks?. External Links: 2105.14491, Link Cited by: §B.1, §4.1. Bruinsma et al. (2023) W. Bruinsma, S. Markou, J. Requeima, A. Y. K. Foong, T. Andersson, A. Vaughan, A. Buonomo, S. Hosking, and R. E. Turner Autoregressive Conditional Neural Processes. In International Conference on Learning Representations, External Links: Link Cited by: §2. Bruinsma (2022) W. Bruinsma Convolutional conditional neural processes. Ph.D. Thesis, Apollo - University of Cambridge Repository, (en). External Links: Document, Link Cited by: §4. Burotto et al. (2014) M. Burotto, V. L. Chiou, J. Lee, and E. C. Kohn The mapk pathway across different malignancies: a new perspective. Cancer 120 (22), p. 3446–3456. External Links: Document, Link, https://acsjournals.onlinelibrary.wiley.com/doi/pdf/10.1002/cncr.28864 Cited by: §1. Bynum et al. (2026) L. E. J. Bynum, R. Ranganath, and K. Cho Computational identifiability. External Links: 2606.19361, Link Cited by: §2. Caron et al. (2022) A. Caron, G. Baio, and I. Manolopoulou Counterfactual Learning with Multioutput Deep Kernels. arXiv. Note: arXiv:2211.11119 [cs] External Links: Link, Document Cited by: §2. Carr and Wingate (2019) A. Carr and D. Wingate Graph neural processes: towards bayesian graph neural networks. External Links: 1902.10042, Link Cited by: §1, §2, §3.2. Cheng et al. (2024) Y. Cheng, L. Li, T. Xiao, Z. Li, J. Suo, K. He, and Q. Dai Cuts+: high-dimensional causal discovery from irregular time-series. In Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 38, p. 11525–11533. Cited by: §2. Chickering (2002) D. M. Chickering Optimal structure identification with greedy search. Journal of machine learning research 3 (Nov), p. 507–554. Cited by: §2. Dhir et al. (2025) A. Dhir, C. Diaconu, V. M. Lungu, J. Requeima, R. E. Turner, and M. v. d. Wilk Estimating Interventional Distributions with Uncertain Causal Graphs through Meta-Learning. arXiv (en). Note: arXiv:2507.05526 [cs] External Links: Link, Document Cited by: §2, §2. Diepen et al. (2023) M. M. V. Diepen, I. G. Bucur, T. Heskes, and T. Claassen Beyond the markov equivalence class: extending causal discovery under latent confounding. In Proceedings of the Second Conference on Causal Learning and Reasoning, M. van der Schaar, C. Zhang, and D. Janzing (Eds.), Proceedings of Machine Learning Research, Vol. 213, p. 707–725. Cited by: §2. Foong et al. (2020) A. Foong, W. Bruinsma, J. Gordon, Y. Dubois, J. Requeima, and R. Turner Meta-learning stationary stochastic process prediction with convolutional neural processes. Advances in Neural Information Processing Systems 33, p. 8284–8295. Cited by: §2. Garnelo et al. (2018) M. Garnelo, J. Schwarz, D. Rosenbaum, F. Viola, D. J. Rezende, S. M. A. Eslami, and Y. W. Teh Neural Processes. arXiv (en). Note: arXiv:1807.01622 [cs] External Links: Link, Document Cited by: §1, §2, §3.2. Geffner et al. (2022) T. Geffner, J. Antoran, A. Foster, W. Gong, C. Ma, E. Kiciman, A. Sharma, A. Lamb, M. Kukla, N. Pawlowski, M. Allamanis, and C. Zhang Deep end-to-end causal inference. External Links: 2202.02195, Link Cited by: §2. Gichoya et al. (2022) J. W. Gichoya, I. Banerjee, A. R. Bhimireddy, J. L. Burns, L. A. Celi, L. Chen, R. Correa, N. Dullerud, M. Ghassemi, S. Huang, et al. AI recognition of patient race in medical imaging: a modelling study. The Lancet Digital Health 4 (6), p. e406–e414. Cited by: §B.2. Gilmer et al. (2017) J. Gilmer, S. S. Schoenholz, P. F. Riley, O. Vinyals, and G. E. Dahl Neural message passing for quantum chemistry. In International conference on machine learning, p. 1263–1272. Cited by: §3.2. Gordon et al. (2019) J. Gordon, J. Bronskill, M. Bauer, S. Nowozin, and R. E. Turner Meta-learning probabilistic inference for prediction. External Links: 1805.09921, Link Cited by: §4. He et al. (2016) K. He, X. Zhang, S. Ren, and J. Sun Deep residual learning for image recognition. In Proceedings of the IEEE conference on computer vision and pattern recognition, p. 770–778. Cited by: §B.2. Heckman and Pinto (2022) J. J. Heckman and R. Pinto Causality and econometrics. Working Paper Technical Report 29787, Working Paper Series, National Bureau of Economic Research. External Links: Document, Link Cited by: §1. Hornik et al. (1989) K. Hornik, M. Stinchcombe, and H. White Multilayer feedforward networks are universal approximators. Neural networks 2 (5), p. 359–366. Cited by: 1.§. Hu et al. (2023) J. Hu, Y. Liang, Z. Fan, H. Chen, Y. Zheng, and R. Zimmermann Graph neural processes for spatio-temporal extrapolation. In Proceedings of the 29th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, KDD ’23, p. 752–763. External Links: Link, Document Cited by: §2. Irvin et al. (2019) J. Irvin, P. Rajpurkar, M. Ko, Y. Yu, S. Ciurea-Ilcus, C. Chute, H. Marklund, B. Haghgoo, R. Ball, K. Shpanskaya, J. Seekins, D. A. Mong, S. S. Halabi, J. K. Sandberg, R. Jones, D. B. Larson, C. P. Langlotz, B. N. Patel, M. P. Lungren, and A. Y. Ng CheXpert: a large chest radiograph dataset with uncertainty labels and expert comparison. External Links: 1901.07031, Link Cited by: §5. Javaloy et al. (2023) A. Javaloy, P. Sánchez-Martín, and I. Valera Causal normalizing flows: from theory to practice. arXiv. Note: arXiv:2306.05415 [cs] External Links: Link, Document Cited by: §2. Javanmardi et al. (2018) M. Javanmardi, M. Sajjadi, T. Liu, and T. Tasdizen Unsupervised total variation loss for semi-supervised deep learning of semantic segmentation. External Links: 1605.01368, Link Cited by: §B.2. Kaltenpoth and Vreeken (2023) D. Kaltenpoth and J. Vreeken Nonlinear causal discovery with latent confounders. In Proceedings of the 40th International Conference on Machine Learning, A. Krause, E. Brunskill, K. Cho, B. Engelhardt, S. Sabato, and J. Scarlett (Eds.), Proceedings of Machine Learning Research, Vol. 202, p. 15639–15654. Cited by: §2. Karimi et al. (2020a) A. Karimi, J. v. Kügelgen, B. Schölkopf, and I. Valera Algorithmic recourse under imperfect causal knowledge: a probabilistic approach. arXiv (en). Note: arXiv:2006.06831 [cs] External Links: Link, Document Cited by: §2. Karimi et al. (2020b) A. Karimi, J. von Kügelgen, B. Schölkopf, and I. Valera Algorithmic recourse under imperfect causal knowledge: a probabilistic approach. External Links: 2006.06831, Link Cited by: §1, §2. Khemakhem et al. (2021) I. Khemakhem, R. P. Monti, R. Leech, and A. Hyvärinen Causal Autoregressive Flows. arXiv (en). Note: arXiv:2011.02268 [stat] External Links: Link, Document Cited by: §2. Kim et al. (2019) H. Kim, A. Mnih, J. Schwarz, M. Garnelo, A. Eslami, D. Rosenbaum, O. Vinyals, and Y. W. Teh Attentive neural processes. arXiv preprint arXiv:1901.05761. Cited by: §2. Kladny et al. (2024) K. Kladny, J. v. Kügelgen, B. Schölkopf, and M. Muehlebach Deep Backtracking Counterfactuals for Causally Compliant Explanations. arXiv (en). Note: arXiv:2310.07665 [cs] External Links: Link, Document Cited by: §2, item ∙ , §5.2. Kullback and Leibler (1951) S. Kullback and R. A. Leibler On information and sufficiency. The annals of mathematical statistics 22 (1), p. 79–86. Cited by: §4. Le et al. (2024) M. K. Le, K. Do, and T. Tran Learning Structural Causal Models from Ordering: Identifiable Flow Models. arXiv. Note: arXiv:2412.09843 [cs] External Links: Link, Document Cited by: §1, §2, item ∙ , §5.1. Li et al. (2024) C. Li, X. Shen, and W. Pan Nonlinear causal discovery with confounders. Journal of the American Statistical Association 119 (546), p. 1205–1214. Cited by: §2. Li et al. (2018) Q. Li, Z. Han, and X. Wu Deeper insights into graph convolutional networks for semi-supervised learning. In Proceedings of the AAAI conference on artificial intelligence, Vol. 32. Cited by: §4.1. Li and Han (2023) T. Li and L. Han Dealing with explainability requirements for machine learning systems. In 2023 IEEE 47th Annual Computers, Software, and Applications Conference (COMPSAC), Vol. , p. 1203–1208. External Links: Document Cited by: §1. Liang and Gao (2021) H. Liang and J. Gao How neural processes improve graph link prediction. External Links: 2109.14894, Link Cited by: §2. Loh (2025) W. W. Loh Causal inference with unobserved confounding: leveraging negative control outcomes using lavaan. Multivariate Behavioral Research 60 (5), p. 1029–1041. Cited by: §1. Lorch et al. (2022) L. Lorch, S. Sussex, J. Rothfuss, A. Krause, and B. Schölkopf Amortized inference for causal structure learning. External Links: 2205.12934, Link Cited by: §6. Loshchilov and Hutter (2019) I. Loshchilov and F. Hutter Decoupled weight decay regularization. External Links: 1711.05101, Link Cited by: §B.1. Luedtke and Fukumizu (2025) A. Luedtke and K. Fukumizu DoubleGen: Debiased Generative Modeling of Counterfactuals. arXiv (en). Note: arXiv:2509.16842 [stat] External Links: Link, Document Cited by: §2. Ma et al. (2026) H. Ma, F. Quinzan, T. Willem, and S. Bauer AI alignment in medical imaging: unveiling hidden biases through counterfactual analysis. External Links: 2504.19621, Link Cited by: §2, §5. Morgan and Winship (2007) S. L. Morgan and C. Winship Counterfactuals and causal inference. New York: Cambridge. Cited by: §1. Murphy (2022) K. P. Murphy Probabilistic machine learning: an introduction. MIT press. Cited by: §4.1. Nasr-Esfahany et al. (2023) A. Nasr-Esfahany, M. Alizadeh, and D. Shah Counterfactual Identifiability of Bijective Causal Models. arXiv (en). Note: arXiv:2302.02228 [stat] External Links: Link, Document Cited by: §1, item ∙ , §5.1. Nguyen et al. (2020) T. Nguyen, G. S. Collins, P. Landais, and Y. Le Manach Counterfactual clinical prediction models could help to infer individualized treatment effects in randomized controlled trials—an illustration with the international stroke trial. Journal of clinical epidemiology 125, p. 47–56. Cited by: §1. Pawlowski et al. (2020) N. Pawlowski, D. C. Castro, and B. Glocker Deep Structural Causal Models for Tractable Counterfactual Inference. arXiv (en). Note: arXiv:2006.06485 [stat] External Links: Link, Document Cited by: §2, §2. Pearl (1996) J. Pearl A causal calculus for statistical research. In Learning from Data: Artificial Intelligence and Statistics V, p. 23–33. Cited by: §1. Pearl (2009) J. Pearl Causality. 2 edition, Cambridge University Press. Cited by: §1, §2, §3.1, §3.1. Peters et al. (2017) J. Peters, D. Janzing, and B. Schölkopf Elements of causal inference: foundations and learning algorithms. The MIT press. Cited by: §2. Rasal et al. (2025) R. Rasal, A. Kori, F. D. S. Ribeiro, T. Xia, and B. Glocker Diffusion Counterfactual Generation with Semantic Abduction. arXiv (en). Note: arXiv:2506.07883 [cs] External Links: Link, Document Cited by: §1, §2. Ribeiro et al. (2025) F. D. S. Ribeiro, A. Santhirasekaram, and B. Glocker Counterfactual Identifiability via Dynamic Optimal Transport. arXiv (en). Note: arXiv:2510.08294 [cs] External Links: Link, Document Cited by: §1, §2, §2, §5. Ruiz-Perez et al. (2024) D. Ruiz-Perez, G. Isabella, S. Musfiqur, M. Kalai, and N. Giri Unfolding and de-confounding: biologically meaningful causal inference from longitudinal multi-omic networks using metalica. Msystems 9 (10). Cited by: §1. Sanchez-Martin et al. (2021) P. Sanchez-Martin, M. Rateike, and I. Valera VACA: Design of Variational Graph Autoencoders for Interventional and Counterfactual Queries. arXiv (en). Note: arXiv:2110.14690 [stat] External Links: Link, Document Cited by: §B.3, §2, item ∙ , §5.1, §5, §5. Sanderson et al. (2022) E. Sanderson, M. M. Glymour, M. V. Holmes, H. Kang, J. Morrison, M. R. Munafò, T. Palmer, C. M. Schooling, C. Wallace, Q. Zhao, et al. Mendelian randomization. Nature reviews Methods primers 2 (1), p. 6. Cited by: §1. Schölkopf et al. (2021) B. Schölkopf, F. Locatello, S. Bauer, N. R. Ke, N. Kalchbrenner, A. Goyal, and Y. Bengio Toward causal representation learning. Proceedings of the IEEE 109 (5), p. 612–634. Cited by: §1. Shimizu et al. (2006) S. Shimizu, P. O. Hoyer, A. Hyvarinen, and A. Kerminen A linear non-gaussian acyclic model for causal discovery. Journal of Machine Learning Research 7 (72), p. 2003–2030. External Links: Link Cited by: §2. Shpitser and Pearl (2012) I. Shpitser and J. Pearl What counterfactuals can be tested. arXiv preprint arXiv:1206.5294. Cited by: §1. Soleymani et al. (2020) A. Soleymani, A. Raj, S. Bauer, B. Schölkopf, and M. Besserve Causal feature selection via orthogonal search. arXiv preprint arXiv:2007.02938. Cited by: §5. Spirtes et al. (2001) P. Spirtes, C. Glymour, and R. Scheines Causation, prediction, and search. The MIT press. Cited by: §1, §2. Stovitz and Shrier (2019) S. D. Stovitz and I. Shrier Causal inference for clinicians. BMJ Evidence-Based Medicine 24 (3), p. 109–112. External Links: Document, ISSN 2515-446X, Link, https://ebm.bmj.com/content/24/3/109.full.pdf Cited by: §1. Verma and Pearl (2022) T. S. Verma and J. Pearl Equivalence and synthesis of causal models. In Probabilistic and causal inference: The works of Judea Pearl, p. 221–236. Cited by: §2. von Kügelgen et al. (2022) J. von Kügelgen, A. Mohamed, and S. Beckers Backtracking Counterfactuals. arXiv. Note: Version Number: 3 External Links: Link, Document Cited by: §2, §2, §3.1, §3.1. Weberpals et al. (2024) J. Weberpals, S. R. Raman, P. A. Shaw, H. Lee, M. Russo, B. G. Hammill, S. Toh, J. G. Connolly, K. J. Dandreo, F. Tian, et al. A principled approach to characterize and analyze partially observed confounder data from electronic health records. Clinical Epidemiology, p. 329–343. Cited by: §1. Wu et al. (2025) Y. Wu, L. McConnell, and C. Iriondo Counterfactual Generative Modeling with Variational Causal Inference. arXiv (en). Note: arXiv:2410.12730 [cs] External Links: Link, Document Cited by: §2, §2, item ∙ , §5.2. Xia et al. (2022) K. Xia, Y. Pan, and E. Bareinboim Neural Causal Models for Counterfactual Identification and Estimation. arXiv (en). Note: arXiv:2210.00035 [cs] External Links: Link, Document Cited by: §1. Yu et al. (2019) Y. Yu, J. Chen, T. Gao, and M. Yu DAG-GNN: DAG structure learning with graph neural networks. In Proceedings of the 36th International Conference on Machine Learning, K. Chaudhuri and R. Salakhutdinov (Eds.), Proceedings of Machine Learning Research, Vol. 97, p. 7154–7163. External Links: Link Cited by: §2. Zečević et al. (2021) M. Zečević, D. S. Dhami, P. Veličković, and K. Kersting Relating Graph Neural Networks to Structural Causal Models. arXiv (en). Note: arXiv:2109.04173 [cs] External Links: Link, Document Cited by: §2. Zhang et al. (2018) R. Zhang, P. Isola, A. A. Efros, E. Shechtman, and O. Wang The unreasonable effectiveness of deep features as a perceptual metric. External Links: 1801.03924, Link Cited by: §B.2. Zheng et al. (2018) X. Zheng, B. Aragam, P. Ravikumar, and E. P. Xing DAGs with no tears: continuous optimization for structure learning. External Links: 1803.01422, Link Cited by: §1, §2, §2. Appendix .1 Equivalence of a Direct MPNN to a Non-directional MPNN on the Dual Graph This section provides an expressivity result justifying the removal of explicit exogenous nodes from the dual-graph architecture. Proposition 1 (Expressive sufficiency after eliminating exogenous nodes). Let GexoG_exo denote the dual graph containing explicit exogenous nodes, and let GdirG_dir denote the reduced graph obtained by removing the exogenous nodes and directly connecting each factual–counterfactual pair (xi,xi∗)(x_i,x_i ). Assume that the endogenous subgraph of GdirG_dir is fully connected. Then, for every function of the endogenous nodes that can be approximated arbitrarily well by a bidirectional MPNN on GexoG_exo, there exists a direct MPNN on GdirG_dir, with sufficiently large hidden dimension, that approximates the same function arbitrarily well. Proof. A bidirectional MPNN layer on GexoG_exo has the form hv(ℓ)=Φℓ(hv(ℓ−1),∑(w,v)∈EMℓ→(hv(ℓ−1),hw(ℓ−1)),∑(v,w)∈EMℓ←(hv(ℓ−1),hw(ℓ−1)))h_v^( )= _ \! (h_v^( -1), _(w,v)∈ EM_ ^→(h_v^( -1),h_w^( -1)), _(v,w)∈ EM_ ^←(h_v^( -1),h_w^( -1)) ) where Mℓ→M_ ^→ and Mℓ←M_ ^← are neural networks handling incoming and outcoming edges respectively. Since the values of the exogenous variables are not observed, the initial exogenous embeddings are fixed constants, hui(0)=cU,hui∗(0)=cU∗.h_u_i^(0)=c_U, h_u_i^*^(0)=c_U^*. Consequently, the hidden state of every exogenous node is completely determined by the hidden states of neighbouring endogenous and exogenous nodes. GdirG_dir does not have exogenous nodes, so we model them by doubling the hidden dimension of the endogenous nodes. Instead of explicitly representing uiu_i and ui∗u_i^*, the direct MPNN stores their hidden states inside the corresponding endogenous representations, h~xi(ℓ)=(hxi(ℓ),hui(ℓ)),h~xi∗(ℓ)=(hxi∗(ℓ),hui∗(ℓ)). h_x_i^( )= (h_x_i^( ),h_u_i^( ) ), h_x_i^*^( )= (h_x_i^*^( ),h_u_i^*^( ) ). Assume inductively that these augmented representations reproduce the hidden states of the bidirectional MPNN after layer ℓ−1 -1. The bidirectional update of an exogenous node is a continuous function hui(ℓ)=Fℓ,i(hui(ℓ−1),hui∗(ℓ−1),hxi(ℓ−1)),h_u_i^( )=F_ ,i (h_u_i^( -1),h_u_i^*^( -1),h_x_i^( -1) ), and analogously for ui∗u_i^*. Since the endogenous graph is fully connected, node xix_i receives the representation of xi∗x_i^* through an ordinary directed message. Hence, all arguments of Fℓ,iF_ ,i are available locally in the augmented representations. By the universal approximation theorem [Hornik et al., 1989], an MLP can approximate Fℓ,iF_ ,i arbitrarily well, allowing the direct MPNN to update the virtual exogenous states. Next, consider the update of an endogenous node, hxi(ℓ)=Gℓ,i(hxi(ℓ−1),hy(ℓ−1):y∈Nend(xi),hui(ℓ−1)).h_x_i^( )=G_ ,i (h_x_i^( -1),\h_y^( -1):y∈ N_end(x_i)\,h_u_i^( -1) ). The virtual state hui(ℓ−1)h_u_i^( -1) is stored locally in h~xi(ℓ−1) h_x_i^( -1). Because the endogenous graph is fully connected, node xix_i receives the hidden state of every endogenous node. Hence, all arguments of Gℓ,iG_ ,i are available to the direct MPNN. Again, the update function can be approximated arbitrarily well by an MLP. After each layer, h~xi(ℓ)=(hxi(ℓ),hui(ℓ)), h_x_i^( )= (h_x_i^( ),h_u_i^( ) ), and similarly for xi∗x_i^*. By induction, the hidden representations of all endogenous nodes produced by the bidirectional MPNN can be reproduced arbitrarily closely by the direct MPNN. Since the readout is performed only on endogenous nodes, reproducing the exogenous hidden states explicitly is unnecessary. ∎ .2 CFGNP Complexity Analysis We calculate the complexity in each part of the module. We first note that the GAT layer has computational complexity Θ⁡(||dd′+|ℰ|d′) (|V|d +|E|d ) for a graph with |||V| nodes and |ℰ||E| edges and the hidden dimensions go from d to d′d [Brody et al., 2021]. Let |A||A| be the number of edges in the graph with adjacency matrix A. Changing Data Dimension Applying the MLP to each entry in O∈ℝ(Nobs⋅n)×ddataO ^(N_obs· n)× d_data and z,x∗∈ℝ(Ncount⋅n)×ddataz,x^* ^(N_count· n)× d_data takes (Nobsn+2Ncountn)ddatadhidden∈Θ⁡(nddatadhidden(Ncount+Nobs))(N_obsn+2N_countn)\,d_datad_hidden∈ \! (n\,d_datad_hidden(N_count+N_obs) ) operations. Since we use Ncount=1N_count=1 throughout this work we get complexity Θ⁡(nddatadhiddenNobs) \! (n\,d_datad_hiddenN_obs ) Propagation of Observational Information Each layer of the GAT uses adjacency matrix Gc=(Mnn⋯nInnMn⋱n⋯nMnInn⋯n).G_c= pmatrixM_n&0_n&·s&0_n&I_n\\ 0_n&M_n& & & \\ & & &0_n& \\ 0_n&·s&0_n&M_n&I_n\\ 0_n&·s&·s&0_n&0_n pmatrix. With |Mn|=n⁡(n−1)|M_n|=n(n-1) and |In|=n|I_n|=n we get |Gc|=Nobsn(n−1)+Nobsn|G_c|=N_obsn(n-1)+N_obsn. Further, we have ||=(Nobs+1)n|V|=(N_obs+1)n. For one layer we then get computational complexity Θ⁡((Nobs+1)ndhidden2+(Nobsn(n−1)+Nobsn)dhidden). ((N_obs+1)nd_hidden^2+(N_obsn(n-1)+N_obsn)d_hidden). For all K layers in the block, we have computational complexity Θ⁡(K⁡(Nobs+1)ndhidden2+K⁡(Nobsn2)dhidden)=Θ⁡(KNobsndhidden(dhidden+n)). (K(N_obs+1)nd_hidden^2+K(N_obsn^2)d_hidden)= (KN_obsnd_hidden(d_hidden+n)). Mixing Factual Information Each layer in the module for mixing factual information uses adjacency matrix GC′=(nInnIn).G_C = pmatrix0_n&I_n\\ 0_n&I_n pmatrix. With |GC′|=2n|G_C |=2n and ||=2n|V|=2n we get computational complexity Θ⁡(ndhidden2+ndhidden) (nd_hidden^2+nd_hidden). All three layers for this module thus have complexity Θ⁡(ndhidden2+ndhidden) (nd_hidden^2+nd_hidden). Estimating Counterfactual Graph Each layer in this block uses adjacency matrix GC′=(MnInInMn).G_C = pmatrixM_n&I_n\\ I_n&M_n pmatrix. We then get |GC′|=2n(n−1)+2n=2n2|G_C |=2n(n-1)+2n=2n^2 and ||=2n|V|=2n. This gives us computational complexity of Θ⁡(ndhidden2+n2dhidden) (nd_hidden^2+n^2d_hidden) in each layer. All min⁡⌊n2⌋,5 \ n2 ,5\ layers together give us computational complexity Θ⁡(ndhidden2+n2dhidden) (nd_hidden^2+n^2d_hidden) for the entire block. The MLP layer needs Θ⁡(ndhiddenNC) (n\,d_hiddenN_C) operations. Overall, the computational complexity of CFGNP is Θ⁡(CLOSE \! ( nddatadhiddenNobs+KNobsndhidden(dhidden+n) n\,d_datad_hiddenN_obs+KN_obsnd_hidden(d_hidden+n) OPEN+ndhidden2+n2dhidden+ndhiddenNC). +nd_hidden^2+n^2d_hidden+nd_hiddenN_C ). .3 Permutation Properties of CFGNP Proposition 2. Let F⁡(z,T∗,t∗,Y∗,O)F(z,T^*,t^*,Y^*,O) denote the output of the CFGNP model restricted to the target nodes Y∗Y^*. Then the CFGNP is 1. permutation equivariant with respect to permutations of the nodes in Y∗Y^*; and 2. permutation invariant with respect to permutations of the nodes in (Y∗)c(Y^*)^c. Proof. In the following, let MnM_n be the n×n× n matrix that has 00 on the diagonal and everywhere else 11. Let InI_n be the n×n× n identity matrix. Finally, let n0_n be the n×n× n matrix with only 00. Let P∈0,1n×nP∈\0,1\^n× n be the permutation matrix associated with a permutation π of the n variables. We use the convention Pei=eπ⁡(i),Pe_i=e_π(i), where eie_i denotes the i-th standard basis vector. Consequently, for every node-indexed matrix H∈ℝn×dH ^n× d, (PH)j=Hπ−1(j).(PH)_j=H_π^-1(j). A permutation represents only a relabelling of the variables and must therefore be applied consistently to every node-indexed input. We define zP:=Pz,xP∗:=Px∗,Ok,P:=POkz_P:=Pz, x_P^*:=Px^*, O_k,P:=PO_k for every observational sample OkO_k. We further define Tπ∗:=π⁡(T∗),Yπ∗:=π⁡(Y∗).T^*_π:=π(T^*), Y_π^*:=π(Y^*). Let τT∗∈0,1n _T^*∈\0,1\^n be the indicator vector of the intervention set T∗T^*, i.e., (τT∗)i=1,i∈T∗,0,i∉T∗.( _T^*)_i= cases1,&i∈ T^*,\\ 0,&i∉ T^*. cases The indicator vector of the permuted intervention set Tπ∗T^*_π is τT∗,P:=PτT∗. _T^*,P:=P _T^*. Indeed, for every j∈1,…,nj∈\1,…,n\, (τT∗,P)j=(PτT∗)j=(τT∗)π−1(j),( _T^*,P)_j=(P _T^*)_j=( _T^*)_π^-1(j), which is equal to 11 precisely when π−1(j)∈T∗π^-1(j)∈ T^*, or equivalently when j∈π⁡(T∗)j∈π(T^*). Recall that the counterfactual input x∗x^* contains the intervention values at the indices in T∗T^* and is zero elsewhere: xi∗=ti∗,i∈T∗,0,i∉T∗.x_i^*= casest_i^*,&i∈ T^*,\\ 0,&i∉ T^*. cases For every j, (xP∗)j=(Px∗)j=xπ−1(j)∗.(x_P^*)_j=(Px^*)_j=x^*_π^-1(j). (xP∗)j=tπ−1(j)∗,π−1(j)∈T∗,0,π−1(j)∉T∗.⇔(xP∗)j=(Pt∗)j,j∈Tπ∗,0,j∉Tπ∗.(x_P^*)_j= casest_π^-1(j)^*,&π^-1(j)∈ T^*,\\ 0,&π^-1(j)∉ T^*. cases (x_P^*)_j= cases(Pt^*)_j,&j∈ T^*_π,\\ 0,&j∉ T^*_π. cases Thus, the intervention values and the intervention indices are permuted consistently. We now show the characteristics for each of the model components. Equivariance of nodewise MLPs. We first consider the nodewise MLPs applied at the beginning of the model. For every observational sample OkO_k, the same MLP hOh_O is applied independently to every node. Therefore, hO(Ok,P)=PhO(Ok).h_O(O_k,P)=Ph_O(O_k). For the factual observation z, define the nodewise transformation hz(z,τT∗)i=(τT∗)ihz,int(zi)+(1−(τT∗)i)hz,common(zi).h_z(z, _T^*)_i=( _T^*)_ih_z,int(z_i)+ (1-( _T^*)_i )h_z,common(z_i). We now verify permutation equivariance coordinatewise. For an arbitrary j∈1,…,nj∈\1,…,n\: hz(zP,τT∗,P)j h_z(z_P, _T^*,P)_j =(τT∗,P)jhz,int((zP)j)+(1−(τT∗,P)j)hz,common((zP)j) =( _T^*,P)_jh_z,int ((z_P)_j )+ (1-( _T^*,P)_j )h_z,common ((z_P)_j ) =(τT∗)π−1(j)hz,int(zπ−1(j)) =( _T^*)_π^-1(j)h_z,int (z_π^-1(j) ) +(1−(τT∗)π−1(j))hz,common(zπ−1(j)) + (1-( _T^*)_π^-1(j) )h_z,common (z_π^-1(j) ) =hz(z,τT∗)π−1(j) =h_z(z, _T^*)_π^-1(j) =(Phz(z,τT∗))j. = (Ph_z(z, _T^*) )_j. Since this equality holds for every j, hz(Pz,PτT∗)=Phz(z,τT∗).h_z(Pz,P _T^*)=Ph_z(z, _T^*). The same argument applies to the counterfactual input. hx(Px∗,PτT∗)=Phx(x∗,τT∗).h_x(Px^*,P _T^*)=Ph_x(x^*, _T^*). Equivariance of the GATs. We now consider the GATs. A GAT mθm_θ is permutation equivariant in the sense that, for a node-feature matrix H, adjacency matrix G, and permutation matrix Q, mθ(QH,QGQ⊤)=Qmθ(H,G).m_θ(QH,QGQ )=Qm_θ(H,G). For the observational graph, the same permutation is applied to every observational sample and to the aggregation nodes q. The corresponding block-diagonal permutation matrix is QO=diag⁡(P,…,P⏟Nobs+1 times)=INobs+1⊗P.Q_O=diag ( P,…,P_N_obs+1 times )=I_N_obs+1 P. The adjacency matrix GcG_c consists only of blocks MnM_n, InI_n, and n0_n. Since PMnP⊤=Mn,PInP⊤=In,PnP⊤=n,PM_nP =M_n, PI_nP =I_n, P0_nP =0_n, we have QOGcQO⊤=Gc.Q_OG_cQ_O =G_c. This suppresses the unchanged adjacency matrix from the notation, mθ(1)(QO(O′q)) m_θ^(1) (Q_O pmatrixO \\ q pmatrix ) =QOmθ(1)((O′q))=QO(O~q′). =Q_Om_θ^(1) ( pmatrixO \\ q pmatrix )=Q_O pmatrix O\\ q pmatrix. It follows in particular that qP′=Pq′.q _P=Pq . For the graph containing the factual representation z and the aggregated observational representation q′q , define the permutation Qzq=(PnnP).Q_zq= pmatrixP&0_n\\ 0_n&P pmatrix. Because QzqGC′Qzq⊤=GC′,Q_zqG_C Q_zq =G_C , permutation equivariance gives mθ(2)(Qzq(zq′)) m_θ^(2) (Q_zq pmatrixz\\ q pmatrix ) =Qzqmθ(2)((zq′)) =Q_zqm_θ^(2) ( pmatrixz\\ q pmatrix ) =Qzq(z′q′) =Q_zq pmatrixz \\ q pmatrix =(Pz′Pq′). = pmatrixPz \\ Pq pmatrix. For the factual–counterfactual graph, define Qzx∗=(PnnP).Q_zx^*= pmatrixP&0_n\\ 0_n&P pmatrix. With Qzx∗GC′Qzx∗⊤=GC′,Q_zx^*G_C Q_zx^* =G_C , we similarly obtain mθcf(Qzx∗(z′x∗)) m_θ^cf (Q_zx^* pmatrixz \\ x^* pmatrix ) =Qzx∗mθcf((z′x∗)) =Q_zx^*m_θ^cf ( pmatrixz \\ x^* pmatrix ) =Qzx∗(z′x∗′) =Q_zx^* pmatrixz \\ x^* pmatrix =(Pz′Px∗′). = pmatrixPz \\ Px^* pmatrix. In particular, the counterfactual representation satisfies x∗P′=Px∗′.x^* _P=Px^* . The final MLP. The final MLP is applied independently to every node of x∗′x^* , with the same parameters at each node. Therefore, mθ(3)(Px∗′)=Pmθ(3)(x∗′).m_θ^(3)(Px^* )=Pm_θ^(3)(x^* ). Let R:=mθ(3)(x∗′)R:=m_θ^(3)(x^* ) denote the nodewise output containing the mixture parameters. The preceding arguments imply that the nodewise output of the model transforms as RP=PR.R_P=PR. The final nodewise output R∈ℝn×3NcR ^n× 3N_c contains one row for each node, and the CFGNP returns only the rows indexed by the target set: F⁡(z,T∗,t∗,Y∗,O)=R|Y∗.F(z,T^*,t^*,Y^*,O)=R |_Y^*. Moreover, the target set is relabelled according to Yπ∗=π⁡(Y∗)Y_π^*=π(Y^*). Since (PR)π⁡(i)=Rifor every i∈Y∗,(PR)_π(i)=R_i every i∈ Y^*, restricting PRPR to the indices in Yπ∗Y_π^* yields precisely the same nodewise outputs as restricting R to Y∗Y^*, but indexed according to the permuted target nodes. Hence, F⁡(zP,Tπ∗,tP∗,Yπ∗,OP)=(PR)|Yπ∗=PY∗(R|Y∗),F(z_P,T^*_π,t_P^*,Y_π^*,O_P)=(PR) |_Y_π^*=P_Y^* (R |_Y^* ), where PY∗P_Y^* denotes the permutation induced by π on the target outputs. Therefore, F⁡(zP,Tπ∗,tP∗,Yπ∗,OP)=PY∗F(z,T∗,t∗,Y∗,O).F(z_P,T^*_π,t_P^*,Y_π^*,O_P)=P_Y^*F(z,T^*,t^*,Y^*,O). This establishes permutation equivariance with respect to the target nodes. If the permutation acts only on nodes in (Y∗)c(Y^*)^c, then every target node is fixed and the induced permutation PY∗P_Y^* is the identity. The same identity therefore reduces to F⁡(zP,Tπ∗,tP∗,Y∗,OP)=F⁡(z,T∗,t∗,Y∗,O),F(z_P,T^*_π,t_P^*,Y^*,O_P)=F(z,T^*,t^*,Y^*,O), which establishes permutation invariance with respect to permutations of the non-target nodes. ∎ .4 Graph-Size and Observation-Count Transferability of CFGNP The CFGNP architecture can be applied to causal graphs with different numbers of nodes and to observational datasets containing different numbers of samples. In particular, the learnable parameters of the model do not depend on either the graph size n or the number of observational samples NobsN_obs. Consequently, a model trained using graphs with n nodes and NobsN_obs observational samples can be applied, without modifying its learned parameters, to graphs with m nodes and Nobs′N_obs observational samples. We formalize this property in the following proposition. Proposition 3 (Graph-size and observation-count transferability). Let Fθ(n,Nobs)F_θ^(n,N_obs) denote the CFGNP applied to a causal query with n endogenous variables and NobsN_obs observational samples. Suppose that the number of message-passing layers is fixed depending on n and independently of NobsN_obs, or, alternatively, repeated message-passing iterations reuse the same parameters. Then the parameter space Θ of Fθ(n,Nobs)F_θ^(n,N_obs) is independent of both n and NobsN_obs. Consequently, for arbitrary n,m,Nobs,Nobs′∈ℕ,n,m,N_obs,N_obs , the same parameter vector θ∈Θθ∈ defines both Fθ(n,Nobs)andFθ(m,Nobs′).F_θ^(n,N_obs) F_θ^(m,N_obs ). Hence, a CFGNP trained on one graph size and one number of observational samples can be evaluated on different graph sizes and different numbers of observational samples without introducing, removing, or resizing any learned parameters. Proof. We show that every learnable component of the CFGNP is independent of both the number of endogenous variables n and the number of observational samples NobsN_obs. First, consider the node-wise input transformations. Let hO:ℝddata→ℝdhiddenh_O:R^d_data ^d_hidden denote the MLP applied to the observational data. For O∈ℝNobs×n×ddata,O ^N_obs× n× d_data, the same map hOh_O is applied independently to every node in every observational sample, Ok,i↦hO(Ok,i),k=1,…,Nobs,i=1,…,n.O_k,i h_O(O_k,i), k=1,…,N_obs, i=1,…,n. Thus, changing either n or NobsN_obs changes only the number of evaluations of hOh_O, but not its input dimension, output dimension, or learned parameters. The same argument applies to the node-wise maps hz,int,h_z,int, hz,common,h_z,common, hx,int,h_x,int, and hx,commonh_x,common. Next, consider a GAT layer. For a node i, its update is of the form hi′=σ⁡(∑j∈⁡(i)αijWhj),h_i =σ ( _j (i) _ijWh_j ), where αij=exp⁡(eij)∑k∈⁡(i)exp⁡(eik), _ij= (e_ij) _k (i) (e_ik), and eij=LeakyReLU(a⊤[Whi∥Whj]).e_ij=LeakyReLU (a [Wh_i\|Wh_j ] ). The learnable quantities W and a depend only on the feature dimensions. They do not depend on the total number of nodes in the graph, the number of edges, or the number of observational samples represented in the graph. Changing n or NobsN_obs therefore changes only the set of nodes and neighborhoods on which the same local update rule is evaluated. We now consider the observational aggregation graph GcG_c. For n endogenous variables and NobsN_obs observational samples, its adjacency matrix is Gc(n,Nobs)=(Mnn⋯nInnMn⋱In⋱n⋯nMnInn⋯n).G_c^(n,N_obs)= pmatrixM_n&0_n&·s&0_n&I_n\\ 0_n&M_n& & &I_n\\ & & &0_n& \\ 0_n&·s&0_n&M_n&I_n\\ 0_n&·s&·s&0_n&0_n pmatrix. Thus, Gc(n,Nobs)G_c^(n,N_obs) contains NobsN_obs diagonal blocks MnM_n, one for each observational sample, and NobsN_obs corresponding blocks InI_n connecting the observational nodes to the aggregation nodes q. Changing NobsN_obs only changes the number of MnM_n- and InI_n-blocks, while changing n only changes their dimensions. The construction rule itself remains unchanged. Since the GAT mθ(1)m_θ^(1) applies the same attention and message-passing functions to every node and edge, its learned parameters are independent of both n and NobsN_obs. In particular, if [O′,q′]=mθ(1)([O′,q],Gc),[O ,q ]=m_θ^(1)([O ,q];G_c), then changing the number of observational samples only changes the size of O′O and the number of incoming messages available to the aggregation nodes q. It does not change the dimensionality of the representation qi′∈ℝdhiddenq_i ^d_hidden associated with each endogenous variable. Hence, the remainder of the architecture receives the same type of representation regardless of NobsN_obs. Next, consider the graph used to combine the observational representation with the factual observation. For n variables, GC′(n)=(nInnIn).G_C (n)= pmatrix0_n&I_n\\ 0_n&I_n pmatrix. For m variables, the corresponding graph is GC′(m)=(mImmIm).G_C (m)= pmatrix0_m&I_m\\ 0_m&I_m pmatrix. These graphs have different dimensions but are generated according to the same construction rule. Moreover, the number NobsN_obs no longer appears explicitly at this stage, because the observational samples have already been aggregated into the n representations contained in q′q . Therefore, the parameters of the GAT used in this stage are independent of both n and NobsN_obs. The same argument applies to the graph used for counterfactual inference. For n variables, GC′(n)=(MnInInMn),G_C (n)= pmatrixM_n&I_n\\ I_n&M_n pmatrix, whereas for m variables, GC′(m)=(MmImImMm).G_C (m)= pmatrixM_m&I_m\\ I_m&M_m pmatrix. Changing the graph size changes only the number of nodes and edges on which the shared message-passing functions are evaluated. The learned parameters of the GAT are unchanged. Furthermore, NobsN_obs does not affect the dimensionality of the node representations entering this stage. Finally, the output network is applied independently and identically to each target node, mθ(3):ℝdhidden→ℝ3Nc.m_θ^(3):R^d_hidden ^3N_c. Changing the number of endogenous variables changes only the number of node-wise predictions, while changing NobsN_obs has no effect on the dimensions of the output map. Hence, neither quantity changes the learned parameters of mθ(3)m_θ^(3). Consequently, every learned component of the CFGNP is parameterized only by quantities such as ddata,d_data, dhidden,d_hidden, Nc,N_c, and not by n or NobsN_obs. The matrices Mn,M_n, In,I_n, and n0_n and the repetition of the observational blocks determine only the structure and size of the computational graph; they contain no learned parameters. Therefore, Θn,Nobs=Θm,Nobs′=Θ _n,N_obs= _m,N_obs = for arbitrary n,m,Nobs,Nobs′∈ℕ.n,m,N_obs,N_obs . Thus, the same parameter vector θ can be applied to causal graphs with different numbers of nodes and to observational datasets containing different numbers of samples. ∎ This establishes architectural transferability across both graph size and observation count. It does not, by itself, guarantee identical predictive performance when n or NobsN_obs differs substantially from the values encountered during training, since such performance additionally depends on the distribution of training and test instances. Appendix A Extended Model Architecture In the following, let MnM_n be the n×n× n matrix that has 00 on the diagonal and everywhere else 11. Let InI_n be the n×n× n identity matrix. Finally, let n0_n be the n×n× n matrix with only 00. A.1 Model Input Our input is a tuple ((z,T∗,t∗,Y∗),O)((z,T^*,t^*,Y^*),O) with O∈ℝNobs×n×dataO ^N_obs× n× d_data, z∈ℝ1×n×dataz ^1× n× d_data, and x∗∈ℝNcount×n×datax^* ^N_count× n× d_data. A.2 Changing data dimension We apply MLPs on each node. On O we apply the same MLP hOh_O across all nodes. However, for the other parts we want the model to learn to distinguish between nodes we intervene on and the rest. Therefore, we apply an MLP hz,inth_z,int on the indices T∗T^* in z and an MLP hz,commonh_z,common on the other nodes. Similarly, we apply the MLP hx,inth_x,int on the indices T∗T^* of x∗x^* and an MLP hx,commonh_x,common on the other indices. A.3 Propagation of Observational Information We reshape O into O′∈ℝ(Nobs⋅n)×dhiddenO ^(N_obs· n)× d_hidden. Furthermore, we introduce a set of nodes q∈ℝn×dhiddenq ^n× d_hidden, initialized as q=q=0. For the GAT we use the adjacency matrix Gc=(Mnn⋯nInnMn⋱n⋯nMnInn⋯n).G_c= pmatrixM_n&0_n&·s&0_n&I_n\\ 0_n&M_n& & & \\ & & &0_n& \\ 0_n&·s&0_n&M_n&I_n\\ 0_n&·s&·s&0_n&0_n pmatrix. We then propagate information on GcG_c using a Graph Attention Network (GAT) mθ(1)m^(1)_θ consisting of K layers. A single GAT layer updates each node representation by computing attention coefficients over its neighbors. Let hi(t)h_i^(t) denote the representation of node i at layer k. The attention coefficients are computed as αij(k)=exp⁡(eij(k))∑k∈⁡(i)exp⁡(eik(k)),eij(k)=LeakyReLU(a⊤[Whi(k)∥Whj(k)]), _ij^(k)= (e_ij^(k)) _k (i) (e_ik^(k)), e_ij^(k)=LeakyReLU\! (a [Wh_i^(k)\|Wh_j^(k)] ), where W and a are learnable parameters. The node representations are then updated according to hi(k+1)=σ⁡(∑j∈⁡(i)αij(k)Whj(k)).h_i^(k+1)=σ\! ( _j (i) _ij^(k)Wh_j^(k) ). Applying the entire GAT to the initial node representations yields [O′,q′]=mθ(1)([O′,q]).[O ,q ]=m^(1)_θ([O ,q]). A.4 Mixing Factual Observation We construct a new graph GC′G_C and apply a GAT mθ(2)m^(2)_θ consisting of three layers. The resulting adjacency matrix is GC′=(nInnIn).G_C = pmatrix0_n&I_n\\ 0_n&I_n pmatrix. Applying the GAT three times yields [z′,q′]=mθ(2)([z,q′]).[z ,q ]=m^(2)_θ([z,q ]). A.5 Estimating Counterfactual Distribution We build a new graph GC′G_C with adjacency matrix GC′=(MnInInMn).G_C = pmatrixM_n&I_n\\ I_n&M_n pmatrix. To propagate information through the entire graph, we use min⁡⌊n2⌋,5 \ n2 ,5\ GAT layers, collected in the MPNN mθ(2)m_θ^(2), which yields [z′,x∗′]=mθ(2)([z′,x∗]).[z ,x^* ]=m_θ^(2)([z ,x^*]). Finally, we apply an MLP mθ(3)m_θ^(3) with hidden dimension 3Nc3N_c. After reshaping, the network outputs the parameters μk,σk,wkk=1Nc\ _k, _k,w_k\_k=1^N_c for each node. Appendix B Experimental Setup B.1 Model Hyperparameters For the (semi-)synthetic datasets, we use a hidden dimension of 80, 5 attention heads in the graph attention network (GAT) layers [Brody et al., 2022], and a batch size of 128. For the predicted MoG we use NC=3N_C=3 Gaussian mixture components. For the image-based counterfactual datasets, which consist of chest X-rays, we use a hidden dimension of 256, 9 attention heads in the GAT layers, and a batch size of 32. We choose a hidden dimension of 256 to match the dimensionality of the encoded image representations, providing a natural choice for the width of the network. We train each model for 200 epochs although we include early stopping after 2525 epochs. Further, we use the AdamW [Loshchilov and Hutter, 2019] with initial learning rate 10−410^-4 and a weight decay of 10−510^-5. B.2 CheXpert Experiment Setup We use Nobs=100N_obs=100 for the CheXpert samples. For each image, the disease information was represented by the 13 CheXpert findings. We combine them into one No Finding attribute. The Disease attribute is just the complement. The resulting dataset has some data imbalances. The Sex feature has 41.24%41.24\% female values and the rest is male. Only 17.24%17.24\% of samples have PA as value and only 9.05%9.05\% of the samples exhibit a disease. Finally, the average age in the samples is 60.6660.66 years. Counterfactual Evaluation. To evaluate which attributes a generated counterfactual image exhibits, we use classifiers/regression models. For the classifiers and regression models of attributes Sex, AP/PA, Age, and Disease we stick closely to Gichoya et al. [2022]. We apply the pretrained ResNet [He et al., 2016] and subsequently apply multiple MLP layers. For regression attributes we have final dimension 1 and for the classification attributes we apply the softmax activation function in the end to normalize the outputs to probabilities. We add here the performance of the classification/regression models of the features. The Sex classifier has a test accuracy of 0.99010.9901 and a ROC-AUC of 0.99950.9995. The Age regression model has MAE 6.14736.1473. Further, the AP/PA classifier has both accuracy and ROC-AUC of 1.01.0. Finally, the No Finding classifier has accuracy 0.90100.9010 and ROC-AUC 0.93140.9314. Figure 5 shows the training curves for each of the classifiers. Figure 5: The training curves of the classifiers and regression models used for the classification of the generated CheXpert images. CheXpert training objective. To the normal classification/regression losses we also add a reconstruction loss to ensure that the images look like actual images. We use MSE, LPIPS[Zhang et al., 2018], and the total variation [Javanmardi et al., 2018] of x^cf x_cf, to get ℒrec=λMMSE⁡(xreal,x^cf)+λLLPIPS⁡(xreal,x^cf)+λTVTV⁡(x^cf),L_rec= _MMSE (x_real, x_cf )+ _LLPIPS (x_real, x_cf )+ _TVTV ( x_cf ), (2) with λM=15 _M=15, λL=30 _L=30, and λTV=30 _TV=30. The resulting total loss combines the classification, regression, and reconstruction objectives: ℒ=∑i∈CE⁡(yi,y^i)+∑r∈ℛ1trMSE⁡(yr,y^r)+ℒrec.L= _i CE (y_i, y_i )+ _r 1t_rMSE (y_r, y_r )+L_rec. (3) where C and ℛR denote the sets of classification and regression features, respectively. B.3 Synthetic Experiments Setup This appendix gives the complete definitions of the SCMs used in the (semi-) synthetic experiments. Their definitions follow the benchmark specification in Appendix E of Sanchez-Martin et al. [2021]. They are reproduced here so that the experimental data-generating processes can be understood without the original paper. We also keep to Sanchez-Martin et al. [2021] for the computation of the standard deviation of the test MSE. We train the model with 5 different seeds and then compute mean and standard deviation of the runs. We use Nobs=100N_obs=100 in all synthetic experiments. B.3.1 Triangle Figure 6: The causal graph of the triangle dataset. The causal graph has three nodes and three edges. The Triangle SCM has the underlying causal graph as depicted in Figure 6. Hence X1X_1 affects X3X_3 through the direct path X1→X3X_1→ X_3 and through the mediated path X1→X2→X3X_1→ X_2→ X_3. Exogenous distributions. For the LIN and NLIN variants, U1 U_1 ∼12(−2,1.5)+12(1.5,1). 12N(-2,1.5)+ 12N(1.5,1). (4) For NADD, U1 U_1 ∼12(−2.5,1)+12(2.5,1). 12N(-2.5,1)+ 12N(2.5,1). (5) The remaining noises are LIN: U2∼⁡(0,1), U_2 (0,1), U3∼⁡(0,1), U_3 (0,1), (6) NLIN: U2∼⁡(0,0.1), U_2 (0,0.1), U3∼⁡(0,1), U_3 (0,1), (7) NADD: U2∼⁡(0,0.25), U_2 (0,0.25), U3∼⁡(0,0.0625). U_3 (0,0.0625). (8) Here the second parameter of ⁡(μ,⋅)N(μ,·) follows the convention used in the benchmark specification. Linear (LIN). X1 X_1 =U1, =U_1, (9) X2 X_2 =−X1+U2, =-X_1+U_2, (10) X3 X_3 =X1+0.25X2+U3. =X_1+0.25X_2+U_3. (11) Nonlinear additive (NLIN). X1 X_1 =U1, =U_1, (12) X2 X_2 =−1+31+exp⁡(−2X1)+U2, =-1+ 31+ (-2X_1)+U_2, (13) X3 X_3 =X1+0.25X22+U3. =X_1+0.25X_2^2+U_3. (14) The causal mechanisms are nonlinear, but each exogenous variable enters its structural equation additively. Nonlinear non-additive (NADD). X1 X_1 =U1, =U_1, (15) X2 X_2 =0.25sgn⁡(U2)X121+U22, =0.25\,sgn(U_2) X_1^21+U_2^2, (16) X3 X_3 =−1+0.1sgn⁡(U3)(X12+X22)+U3. =-1+0.1\,sgn(U_3) (X_1^2+X_2^2 )+U_3. (17) The function sgn⁡(⋅)sgn(·) returns the sign of its argument. B.3.2 M-graph Figure 7: The causal graph of the M-shape data. It has five nodes and four edges. The M-graph contains five endogenous variables as shown in Figure 7. All exogenous variables are independent standard normal variables: Ui∼(0,1),i=1,…,5.U_i (0,1), i=1,…,5. (18) Linear (LIN). X1 X_1 =U1, =U_1, (19) X2 X_2 =U2, =U_2, (20) X3 X_3 =X1+U3, =X_1+U_3, (21) X4 X_4 =−X2+0.5X1+U4, =-X_2+0.5X_1+U_4, (22) X5 X_5 =−1.5X2+U5. =-1.5X_2+U_5. (23) Nonlinear additive (NLIN). X1 X_1 =U1, =U_1, (24) X2 X_2 =U2, =U_2, (25) X3 X_3 =X1+0.5X12+U3, =X_1+0.5X_1^2+U_3, (26) X4 X_4 =−X2+0.5X12+U4, =-X_2+0.5X_1^2+U_4, (27) X5 X_5 =−1.5X22+U5. =-1.5X_2^2+U_5. (28) Nonlinear non-additive (NADD). X1 X_1 =U1, =U_1, (29) X2 X_2 =U2, =U_2, (30) X3 X_3 =X1U3, =X_1U_3, (31) X4 X_4 =(−X2+0.5X12)U4, = (-X_2+0.5X_1^2 )U_4, (32) X5 X_5 =−1.5X22U5. =-1.5X_2^2U_5. (33) B.3.3 German Loan Figure 8: The causal graph of the German loan data. The nodes are real attributes used for German credit rating system. they are connected based on the real-world connections. The Loan benchmark is a semi-synthetic SCM from the causal-recourse literature and is inspired by the German Credit domain. It contains seven endogenous variables as shown in Figure 8. In the figure, G is gender, A age, E education, L loan amount, D loan duration, I income, and S savings. The structural assignments are G G =UG, =U_G, (34) A A =−35+UA, =-35+U_A, (35) E E =−0.5+(1+exp⁡[1−0.5G−(1+exp⁡(−0.1A))−1−UE])−1, =-0.5+ (1+ [1-0.5G- (1+ (-0.1A) )^-1-U_E ] )^-1, (36) L L =1+0.01(A−5)(5−A)+G+UL, =1+0.01(A-5)(5-A)+G+U_L, (37) D D =−1+0.1A+2G+L+UD, =-1+0.1A+2G+L+U_D, (38) I I =−4+0.1(A+35)+2G+GE+UI, =-4+0.1(A+35)+2G+GE+U_I, (39) S S =−4+1.5I>0I+US. =-4+1.5\,I_\I>0\I+U_S. (40) Here I>0I_\I>0\ is one if I>0I>0 and zero otherwise. The exogenous distributions are UG U_G ∼Bernoulli⁡(0.5), (0.5), (41) UA U_A ∼Gamma⁡(10,3.5), (10,3.5), (42) UE U_E ∼⁡(0,0.25), (0,0.25), (43) UL U_L ∼⁡(0,4), (0,4), (44) UD U_D ∼⁡(0,9), (0,9), (45) UI U_I ∼⁡(0,4), (0,4), (46) US U_S ∼⁡(0,25). (0,25). (47) B.4 Structural Equations for Artificial Graphs Variables are generated as X=f⁡(Pa⁡(X))+ϵX=f(Pa(X))+ε, where f is sampled from the following set of transformations: Linear:f1(Pa(X)) : f_1(Pa(X)) =a∑X′∈Pa⁡(X)X′+b, =a _X (X)X +b, a=0.5,b=0, a=0.5,\;b=0, Sum-sqrt:f2(Pa(X)) -sqrt: f_2(Pa(X)) =a∑X′∈Pa⁡(X)|X′|+b, =a _X (X) |X |+b, a=0.5,b=0, a=0.5,\;b=0, Sum-sine:f3(Pa(X)) -sine: f_3(Pa(X)) =a∑X′∈Pa⁡(X)sin⁡(cX′)+b, =a _X (X) (cX )+b, a=1,b=0,c=0.5, a=1,\;b=0,\;c=0.5, Sum-tanh:f4(Pa(X)) -tanh: f_4(Pa(X)) =a∑X′∈Pa⁡(X)tanh⁡(cX′)+b, =a _X (X) (cX )+b, a=1,b=0,c=2, a=1,\;b=0,\;c=2, Geometric mean:f5(Pa(X)) mean: f_5(Pa(X)) =a(∏X′∈Pa⁡(X)|X′|)1|Pa⁡(X)|+b, =a ( _X (X)|X | ) 1|Pa(X)|+b, a=3,b=0.1, a=3,\;b=0.1, Log-sum-exp:f6(Pa(X)) -sum-exp: f_6(Pa(X)) =alog⁡(∑X′∈Pa⁡(X)eX′)+b, =a ( _X (X)e^X )+b, a=1,b=log⁡2, a=1,\;b= 2, Sqrt-sum:f7(Pa(X)) -sum: f_7(Pa(X)) =a|∑X′∈Pa⁡(X)X′|+b, =a | _X (X)X |+b, a=1,b=0. a=1,\;b=0. B.5 Identifiability in a Two-Node SCM Setup. Let =(X0,X1) X=(X_0,X_1) follow the SCM X0:=U0,X1:=aX0+⟨λ(X0),⟩,X_0:=U_0, X_1:=aX_0+ g_λ(X_0), W , (48) where U0∼⁡(0,1)U_0 (0,1), ∼⁡(,I2) W ( 0,I_2) are independent, and s⁡(x):=σu2+σv2x2,φ⁡(x):=arctan⁡(σvx/σu),s(x):= _u^2+ _v^2x^2, (x):= ( _vx/ _u), (49) λ(x):=s⁡(x)(cos⁡(λφ(x)),sin⁡(λφ(x))),λ∈[0,1]. g_λ(x):=s(x) ( (λ (x)), (λ (x)) ), λ∈[0,1]. (50) Here a controls the observational mechanism, whereas λ controls the cross-world coupling. Proposition 4 (Observational invariance). For every λ, X0∼(0,1),X1∣X0=x∼(ax,s(x)2).X_0 (0,1), X_1 X_0=x \! (ax,s(x)^2 ). (51) Hence the observational distribution is independent of λ. The same is true of the intervention Pdo⁡(X0=α)P_do(X_0=α). Proof. Since ‖λ(x)‖2=s(x)2(cos2⁡(λφ(x))+sin2⁡(λφ(x)))=s(x)2,\| g_λ(x)\|^2=s(x)^2 ( ^2(λ (x))+ ^2(λ (x)) )=s(x)^2, we have ⟨λ(x),⟩∼⁡(0,s(x)2) g_λ(x), W (0,s(x)^2), which gives (51). Under do⁡(X0=α)do(X_0=α), X1=aα+⟨λ(α),⟩∼⁡(aα,s(α)2),X_1=aα+ g_λ(α), W (aα,s(α)^2), again independently of λ. ∎ Backtracking semantics. Let the factual evidence be e=X0=x0,X1=x1,r:=x1−ax0.e=\X_0=x_0,\;X_1=x_1\, r:=x_1-ax_0. To construct a counterfactual world satisfying X0′=αX_0 =α, we use the minimum-change exogenous coupling (U0′,′)=(α,).(U_0 , W )=(α, W). (52) Equivalently, (52) is the solution of minu0′,′⁡c0|u0′−U0|2+cW‖′−‖2s.t. u0′=α, _u_0 , w c_0|u_0 -U_0|^2+c_W\| w - W\|^2 .t. u_0 =α, for any c0,cW>0c_0,c_W>0. Thus backtracking changes only the exogenous realization necessary to produce X0′=αX_0 =α while preserving the downstream noise W. Theorem 5 (Backtracking counterfactual). Under (48)–(52), X1∗(α)|e∼⁡(aα+kλr,τλ2),X_1^*(α) e (aα+k_λr, _λ^2 ), (53) where kλ:=s⁡(α)s⁡(x0)cos⁡(λΔφ),τλ2:=s(α)2sin2⁡(λΔφ),k_λ:= s(α)s(x_0) (λ ), _λ^2:=s(α)^2 ^2(λ ), (54) and Δφ:=φ⁡(α)−φ⁡(x0). := (α)- (x_0). Moreover, (53) is identical to the corresponding abduction–action–prediction counterfactual. Proof. The evidence fixes U0=x0U_0=x_0 and imposes ⟨λ(x0),⟩=r. g_λ(x_0), W =r. With ^0:=λ(x0)/s⁡(x0) g_0:= g_λ(x_0)/s(x_0), Gaussian conditioning gives |e∼⁡(rs⁡(x0)^0,I−^0^0⊤). W e ( rs(x_0) g_0,I- g_0 g_0 ). (55) Backtracking changes U0U_0 to α while preserving W, hence X1∗(α)=aα+⟨λ(α),⟩.X_1^*(α)=aα+ g_λ(α), W . Therefore ⁡[X1∗(α)∣e]=aα+rs⁡(x0)⟨λ(α),^0⟩,E[X_1^*(α) e]=aα+ rs(x_0) g_λ(α), g_0 , and Var⁡[X1∗(α)∣e]=s(α)2−⟨λ(α),λ(x0)⟩2s(x0)2.Var[X_1^*(α) e]=s(α)^2- g_λ(α), g_λ(x_0) ^2s(x_0)^2. Since ⟨λ(α),λ(x0)⟩=s⁡(α)s(x0)cos⁡(λΔφ), g_λ(α), g_λ(x_0) =s(α)s(x_0) (λ ), substitution yields (53)–(54). The interventional abduction–action–prediction construction evaluates the same quantity aα+⟨λ(α),⟩aα+ g_λ(α), W under the same conditional law (55), so the two counterfactual distributions coincide. ∎ Corollary 6 (Cross-world non-identifiability). Although λ is invisible to every observational and interventional distribution, it changes the backtracking counterfactual whenever λ↦(kλ,τλ)λ (k_λ, _λ) is non-constant. In particular, τλ2=0⇔λΔφ∈πℤ. _λ^2=0 λ ∈ . For λ=0λ=0, X1∗(α)|e=δaα+s⁡(α)s⁡(x0)r,X_1^*(α) e= _aα+ s(α)s(x_0)r, (56) whereas for λ>0λ>0 the counterfactual has positive conditional variance for almost every (x0,α)(x_0,α) with α≠x0α≠ x_0. For λ=1λ=1, 1(x)=(σu,σvx), g_1(x)=( _u, _vx), and k1=σu2+σv2x0αs(x0)2,τ12=σu2σv2(α−x0)2s(x0)2.k_1= _u^2+ _v^2x_0αs(x_0)^2, _1^2= _u^2 _v^2(α-x_0)^2s(x_0)^2. (57) Use the prior a∼⁡(0,σa2),λ∼πλ,a⟂λ,a (0, _a^2), λ _λ, a λ, and let Dobs=(x0(m),x1(m))m=1MD_obs=\(x_0^(m),x_1^(m))\_m=1^M be observational context data generated from (51). Writing ℰ:=Dobs∪(x0,x1),E:=D_obs∪\(x_0,x_1)\, we get for the optimal backtracking predictive pBCM(y∣x0,x1,X0′=α,Dobs)=∫pbt(y∣e,α,a,λ)p(a,λ∣ℰ)aλ.p_BCM(y x_0,x_1,X_0 =α,D_obs)= p_bt(y e,α,a,λ)p(a,λ )\,da\,dλ. (58) Proposition 7 (Exact posterior and model average). The likelihood is independent of λ, so p⁡(a,λ∣ℰ)=p⁡(a∣ℰ)πλ(λ).p(a,λ )=p(a ) _λ(λ). With wi:=s(x0(i))−2,w_i:=s(x_0^(i))^-2, we have a|ℰ∼⁡(μa,σa|ℰ2),a ( _a,σ^2_a ), (59) where σa|ℰ−2=σa−2+∑i∈ℰwi(x0(i))2,μa=σa|ℰ2∑i∈ℰwix0(i)x1(i).σ^-2_a = _a^-2+ _i w_i(x_0^(i))^2, _a=σ^2_a _i w_ix_0^(i)x_1^(i). (60) Consequently, pBCM(y∣⋅)=∫(y;μaβλ+kλx1,τλ2+σa|ℰ2βλ2)πλ(dλ),p_BCM(y ·)= (y; _a _λ+k_λx_1, _λ^2+σ^2_a _λ^2 ) _λ(dλ), (61) where βλ:=α−kλx0. _λ:=α-k_λx_0. Proof. Proposition 4 implies that the likelihood does not depend on λ, which yields the posterior factorization. (60) is the Gaussian posterior for weighted linear regression. Finally, substituting r=x1−ax0r=x_1-ax_0 into (53) gives aα+kλr=a⁡(α−kλx0)+kλx1=aβλ+kλx1.aα+k_λr=a(α-k_λx_0)+k_λx_1=a _λ+k_λx_1. Marginalizing the Gaussian posterior over a gives (61). ∎ Identifiable versus non-identifiable regimes. Let pBCM∗:=pbt(⋅∣e,α,a∗,λ∗)p_BCM^*:=p_bt(· e,α,a^*,λ^*) denote the backtracking counterfactual under the true mechanism. As M→∞M→∞, σa|ℰ2=Θ(M−1)→0,μa→a∗a.s.,σ^2_a = (M^-1)→ 0, _a→ a^* .s., and therefore pBCM⟹∫⁡(y,a∗α+kλr∗,τλ2)πλ(λ),r∗:=x1−a∗x0.p_BCM (y;a^*α+k_λr^*, _λ^2 ) _λ(dλ), r^*:=x_1-a^*x_0. (62) If πλ=δλ∗ _λ= _λ^*, the limit in (62) equals pBCM∗p_BCM^*. If instead πλ _λ is non-degenerate, then for almost every (x0,α)(x_0,α) with α≠x0α≠ x_0 the map λ↦(kλ,τλ)λ (k_λ, _λ) is non-constant, so the limiting Bayes predictor remains a genuine mixture over observationally indistinguishable cross-world mechanisms. Thus even infinite observational data identifies a but not λ: uncertainty about the observational mechanism vanishes, while uncertainty about the backtracking counterfactual need not. Appendix C Additional Results C.1 Additional CheXpert Results C.1.1 CheXpert Experiments We include some counterfactual images generated by the model. Figure 9: Example counterfactual of a CheXpert image. Observed factual unit (left) has age 8080 and the intervention value is 5858 in the counterfactual (middle). The difference between the pixels is on the right. Figure 10: Example counterfactual of a CheXpert image. Observed factual unit (left) has AP/PA as 11 and the intervention value is 00 in the counterfactual (middle). The difference between the pixels is on the right. Figure 11: An image where the observed factual unit of age is 7474 and the counterfactual intervention is 8585. C.1.2 CheXpert Out-of-Distribution Experiment This model was trained with sex chosen to be the out-of-distribution feature. However, during inference the sex feature is included again to test the performance. Out-of-Distribution Feature Sex. We here include the results when the model was trained with OOD feature sex. Figure 12: An image where the observed factual unit of age is 6565 and the counterfactual intervention is 7272. The model was trained with OOD index sex. Figure 13: An image where the observed factual unit of AP/PA is 0 and the counterfactual intervention is 1. The model was trained with OOD index sex. Out-of-Distribution Feature Age. We here include the results when the model was trained with OOD feature age. Figure 14: An image where the observed factual unit is age 3030 and the counterfactual intervention is age 5252. The model was trained with OOD index age. Figure 15: An image where the observed factual unit of sex is 1 and the counterfactual intervention is 0. The model was trained with OOD index age. Out-of-Distribution Feature AP/PA. We here include the results when the model was trained with OOD feature AP/PA. Figure 16: n image where the observed factual unit of sex is 0 and the counterfactual intervention is 1. The model was trained with OOD index AP/PA. Figure 17: An image where the observed factual unit is age 2626 and the counterfactual intervention is age 6464. The model was trained with OOD index AP/PA. Out-of-Distribution Feature Disease. We here include the results when the model was trained with OOD feature disease. Figure 18: An image where the observed factual unit of AP/PA is 0 and the counterfactual intervention is 1. The model was trained with OOD index disease. Figure 19: An image where the observed factual unit is age 6767 and the counterfactual intervention is age 3737. The model was trained with OOD index disease. C.2 Graph Size Results Figure 20: Performance (MSE) of models with regard to the size of the causal graph for sqrt-sum (left) and sum-sqrt (right) structural equations. The MSE is averaged over distinct causal graphs. Vertical lines show the 25th and 75th quantiles. Figure 21: Performance (MSE) of models with regard to the size of the causal graph for geometric-mean (left) and sum-tanh (right) structural equations. The MSE is averaged over distinct causal graphs. Vertical lines show the 25th and 75th quantiles. C.3 Runtime of Experiments As before, we use a batch size of 128128 for the (semi-)synthetic experiments and a batch size of 3232 for the CheXpert dataset. All experiments are run on an NVIDIA H100 GPU. Dataset CFGNP (Ours) BGM CFM VACA German Loan 31.5719s 5.2537s 115.1188s 5.3015s Triangle LIN 8.1853s 2.9878s 26.3716s 3.7258s Triangle NLIN 7.9175s 3.0323s 26.2195s 3.7365s Triangle NADD 7.9084s 3.0364s 26.3470s 3.7963s M-Shape LIN 12.9088s 3.3865s 41.5352s 3.8361s M-Shape NLIN 12.9274s 3.3941s 41.8969s 3.8385s M-Shape NADD 12.9784s 3.4396s 41.9043s 3.8209s Artificial Size 10 18.4545s 3.8820s – 2.4226s Artificial Size 15 27.6381s 5.3534s – 2.6127s Artificial Size 20 39.2469s 9.2543s – 2.6588s Artificial Size 25 54.7471s 8.5072s – 2.5591s Artificial Size 30 74.2878s 13.9569s – 2.7384s Table 4: Runtime comparison of CFGNP (Ours), BGM, CFM, and VACA across (semi-) synthetic datasets. The datasets were used throughout the results section. Dataset CFGNP (Ours) VCI DeepBC CheXpert 47.7369s 3.8506s 6032.4057s Table 5: Runtime comparison of CFGNP (Ours), VCI, and DeepBC on the CheXpert dataset. C.4 Convergence in (Non-)Identifiable Cases We empirically evaluate the convergence behaviour of CFGNP. To this end, we evaluate CFGNP on an SCM for which the ground-truth counterfactual distribution can be computed analytically. We consider a causal graph with two nodes and at most one edge. In Appendix B.5, we show that, under certain choices of parameters, such an SCM is identifiable. Furthermore, we analytically compute the ground-truth counterfactual distribution pT∗:=pT(Y∗=y∗∣z,T∗=t∗,,),p_T^*:=p_T(Y^*=y^* z,T^*=t^*,f,G), given the true causal graph and structural mechanisms. We compare the distribution predicted by CFGNP to both pT∗p_T^* and pT:=pT(Y∗=y∗∣O,z,T∗=t∗),p_T:=p_T(Y^*=y^* O,z,T^*=t^*), which conditions only on the observations rather than on knowledge of the SCM (,)(f,G). Subsequently, we repeat the experiment for a non-identifiable parameterization of the two-node SCM. In both cases, we use 4000040000 samples, split into 2100021000 train samples, 90009000 validation samples, and 1000010000 test samples. Figure 22: KL divergences for the two-node experiments across different numbers NobsN_obs of observations supplied to each sample. Vertical lines show the 10th and 90th quantiles. Results are shown for the identifiable (left) and non-identifiable (right) settings. C.5 Transformer Ablation To assess the importance of the CFGNP architecture, we consider a Transformer-based ablation that replaces the core CFGNP processing mechanism with a standard Transformer encoder while retaining the same input and output representations. The observation, factual unit, and counterfactual/intervention inputs are first embedded separately using the same embedding dimensionality as CFGNP. Again, eparate embeddings are retained for intervened and non-intervened variables. Then, we concatenate all these vectors and run a transformer architecture on them. The output of the transformer for x∗x^* is handled in the same way as for CFGNP Unless stated otherwise, all hyperparameters are kept identical to those used for CFGNP, including the base hidden dimension, number of attention heads, number of mixture components, and input/output dimensionality. As in CFGNP, the effective embedding dimension is given by the base hidden dimension multiplied by the number of attention heads. The Transformer-specific hyperparameters are five encoder layers, a feed-forward dimension of four times the embedding dimension, a dropout rate of 0.10.1, and GELU activations. Pre-layer normalization is used within each Transformer encoder layer, followed by a final layer normalization after the encoder stack. Table 6 compares the resulting Transformer ablation against the full CFGNP model. Dataset Variant CFGNP (Ours) Transformer MSE ↓ Triangle LIN (2.59±2.13)×−(2.59± 2.13)× 10^-7 (8.33±0.25)×10−4(8.33± 0.25)× 10^-4 NLIN (2.14±0.34)×−(2.14± 0.34)× 10^-7 (8.37±0.46)×10−4(8.37± 0.46)× 10^-4 NADD (2.81±1.11)×−(2.81± 1.11)× 10^-7 (5.02±0.03)×10−4(5.02± 0.03)× 10^-4 M-shape LIN (3.60±7.00)×−(3.60± 7.00)× 10^-5 (4.79±0.12)×10−3(4.79± 0.12)× 10^-3 NLIN (2.50±0.47)×−(2.50± 0.47)× 10^-6 (6.04±0.19)×10−3(6.04± 0.19)× 10^-3 NADD (1.08±0.35)×−(1.08± 0.35)× 10^-4 (3.04±0.17)×10−3(3.04± 0.17)× 10^-3 German Loan – (0.19±0.09)×−(0.19± 0.09)× 10^-2 (46.92±7.72)×10−2(46.92± 7.72)× 10^-2 Table 6: Mean MSE ± standard deviation of the full CFGNP model and the Transformer ablation across the Triangle, M-shape, and German Loan datasets. Lower values indicate better performance. Best scores are marked in bold.

Instant research discovery