Exogenous Matching: Learning Good Proposals for Tractable Counterfactual Estimation

Authors: yikang  chen
dehui  du
lili  tian
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 10:02:58

Abstract: We propose an importance sampling method for tractable and efficient estimation of counterfactual expressions in general settings, named Exogenous Matching. By minimizing a common upper bound of counterfactual estimators, we transform the variance minimization problem into a conditional distribution learning problem, enabling its integration with existing conditional distribution modeling approaches. We validate the theoretical results through experiments under various types and settings of Structural Causal Models (SCMs) and demonstrate the outperformance on counterfactual estimation tasks compared to other existing importance sampling methods. We also explore the impact of injecting structural prior knowledge (counterfactual Markov boundaries) on the results. Finally, we apply this method to identifiable proxy SCMs and demonstrate the unbiasedness of the estimates, empirically illustrating the applicability of the method to practical scenarios.

Main paper: 
 exogenous matching: 
 In this section, we present only the final theoretical conclusions. Detailed derivations of all theories and equations can be found in App. A.Assumptions In addition to assuming that the SCM M = ⟨U, V, F, P U ⟩ is recursive, the proposed method also assumes: i) for any X ⊆ V, x ∈ Ω X , u ∈ Ω U , F x (u) is computable; ii) P U is sampleable and computable. These assumptions imply the feasibility of counterfactual generation, meaning we can draw arbitrary counterfactual samples, with F x and P U both treated as black boxes. For example, in the context of neural proxy SCMs, P U , as the latent distribution of the generative model, is typically modeled as a simple distribution, while F x is a neural network.Importance sampling for Monte Carlo integration A widely used method for estimating the Lebesgue integral is Monte Carlo integration. However, the indicator function 1 Ω U (Y * ) (u) in Eq. 2 limits the efficiency. We consider 1 Ω U (Y * ) (u) = 1 as a rare event. Since importance sampling is commonly used to handle Monte Carlo estimations involving rare events, we adopt the method of importance sampling on a parameterized proposal distribution Q U for tractable Monte Carlo integration, and the estimator can be derived as:EQUATIONwhere for each sampleu (i) ∼ Q U , it is required that if p(u (i) ) > 0, then q(u (i) ) > 0.The density ratio p(u)/q(u) is also termed as the importance weight and is usually assumed to be bounded. In this section, we assume that our derivations are based on an exogenous distribution defined over a continuous space. Naturally, similar conclusions hold for discrete exogenous distributions.In importance sampling, a proposal distribution with low variance regarding σ Y * is crucial, as it implies that fewer samples are required to achieve the same effect, leading to higher sampling efficiency. Previous works have attempted to constrain variance, such as providing upper bounds [18, 49] , introducing multiple proposals [16, 27] , and learning from sampled data [72, 83, 11, 65] . Given that the samples are i.i.d., with the same sample size n, a lower variance of σ Y * implies a lower estimator variance:EQUATIONIt is evident that this formulation can be directly optimized, and there are several methods employing similar approaches (e.g., transforming it into cross-entropy or χ 2 divergence), as shown in [32, 68] . This process is called proposal learning, which ultimately yields a proposal distribution Q u that is suitable for the estimator.The learned Q U is typically "one-off", rendering it unsuitable for other estimators. This necessitates multiple rounds of proposal learning to estimate an L 3 expression that involves various counterfactual events, which remains inefficient. To address this issue, we propose replacing the proposal distribution Q U with a conditional proposal distribution Q U|Y * . Intuitively, the proposal distribution Q U|y * corresponding to y * ∈ Ω Y * should be concentrated on the support Ω U (δ(y * )), thus allowing for reuse across different estimators.Building on this idea, we extend the methods from [32, 68] , as illustrated in App. B.1. Additionally, we find that when certain conditions are met in Y * , for any y * ∈ Y * , the proposal distribution Q U|y * shares a common variance upper bound: Theorem 1 (Variance Upper Bound). Let σ Y * (u) = (p(u)/q(u | y * )) 1 Ω U (Y * ) (u),EQUATIONwhere the constant c is solely dependent on κ and P U .It is noted that the expected term in Eq. 6 bears resemblance to cross-entropy, albeit not strictly. Intuitively, the original indicator function is now implied as potential responses within the condition, and optimizing this expected term can encourage the density of each proposal distribution at specific locations to match that of the exogenous distribution as closely as possible. Therefore, this approach is named Exogenous Matching (EXOM).The boundedness condition on the importance weights in Thm. 1 is likely to be violated for Y * with an infinite support set. We further discuss this issue in App. A.2, where we: i) provide two relaxed versions of Thm. 1 based on concentration inequalities; ii) demonstrate how to construct a guard proposal distribution that necessarily satisfies the boundedness condition.Generalize Thm. 1 Since any L 3 expression Q is composed of terms P (Y * ), and assuming that it is expressible by a finite set of P (Y * ), we can consider all Y * involved in Q as outcomes of a stochastic process. To allow different forms of Y * , we extend Thm. 1 to general cases, making it potentially applicable to the estimation of any counterfactual probability term P (Y * ) in an L 3 expression Q, with only one proposal learning process is required. EQUATIONwhere the constant c is the same as in Thm. 1.Disregarding the constant terms of the expected variance upper bound (Eq. 7) in Cor. 1, we obtain the optimization objective as follows:EQUATIONBased on our assumptions, P U and Y (s) * (u) are known, while the prior distribution of states Q S needs to be specifically designed for the particular expression in L 3 . The problem is now reformulated as a conditional distribution learning problem, and we can model this conditional distribution using neural networks. The complete learning algorithm is presented in Fig. 5 in App. A.3.Inference Given the conditional proposal Q U|Y * , we can take the expectation of the importance sampling results over Y * to intuitively further enhance robustness. This leads to another unbiased estimator in the form of multiple importance sampling [16, 27] :EQUATIONwhereEQUATIONy (i) * ∼ Q Y * , u (i) ∼ Q U|y (i) *. The distribution Q Y * could be any distribution with support covering Y * .Conditioning We attempt to find a function that maps y * to the vectorized parameters θ y * required for the proposal distribution Q U|y * . The information provided by y * (i.e., the counterfactual event Y * ∈ {y * }) can be interpreted as a set of 4-tuples, where each 4-tuple ⟨Y, y, X, x⟩ expresses y * Eq. 10= == ⇒        c ⟨Y 1 ,y 1 ,X 1 ,x 1 ⟩ c ⟨Y 2 ,y 2 ,X 2 ,x 2 ⟩ . . . c ⟨Y k ,y k ,X k ,x k ⟩        m = m 1 m 2 • • • m k ⇒ Thm. 3 h g θ y * = ⇒ Q U|y * Forward Mask Compute EncodingFigure 2 : Overview of the conditioning and masking process. y * serves as the input to the entire process, m represents the inferred mask, and the vectorized parameters θ y * of the proposal distribution Q U|y * are the output. Different colors represent information from different submodels. Both h and g represent neural networks.counterfactual information within the same submodel M x , including the observed variables Y and their values y, and the intervened variables X and their values x. We represent it as a vector:EQUATIONwhere ⊕ denotes vector concatenation. The functions π and ω map inputs to vectors, with each component of the vector corresponding to an exogenous variable X i ∈ V, defined as follows:EQUATIONTo encode the entire set{⟨Y 1 , y 1 , X 1 , x 1 ⟩, ⟨Y 2 , y 2 , X 2 , x 2 ⟩, .. . }, we introduce two functions, h and g. Here, h serves as an encoder, mapping the vector of a single 4-tuple c ⟨Y,y,X,x⟩ to a latent encoding, while g acts as an aggregator, capturing the permutation invariance of the set. The vectorized parameters θ y * are inferred through the interplay of these two functions:EQUATIONWe choose to model h using a multilayer perceptron (MLP), and model g as a function that ensures permutation invariance for encoding aggregation, such as summation or weighted summation, where we opt for attention [96] to enhance expressiveness.This section will explore how to inject available structural prior knowledge (specifically, counterfactual Markov boundaries) into neural networks used for conditional encoding, intuitively improving the quality and efficiency of learning.Counterfactual Markov boundary Markov boundaries [75] have been used in feature selection [111, 64, 104] and causal discovery [3, 87, 103] , since they reveal the local causal structure of variables, where all elements are strongly correlated with the variable [2] . Just as in feature selection, intuitively we can enhance the performance of distribution learning by masking redundant information.For modeling distributions defined over exogenous variables and conditioned on counterfactual variables, the structural prior information used here is named counterfactual Markov boundaries:Definition 2 (Counterfactual Markov Boundary). For an exogenous variable U j ∈ U and a set of counterfactual variables Y, along with their joint distributionP . If U j is independent of Y \ X given X under P , i.e. U j ⊥ ⊥ P (Y \ X) | X, then X is termed a Markov blanket of U j on Y. The collection of all Markov blankets of U j on Y is denoted as B j (Y). If X ∈ B j (Y) is a Markov blanket of U j on Y, and there exists no X ′ ⊊ X such that X ′ ∈ B j (Y), then X is termed a (counterfactual) Markov boundary of U j on Y, denoted as B j (Y).Methods for learning Markov boundaries based on independence in data distribution (e.g. [3, 87, 103] ) cannot be directly applied to learn counterfactual Markov boundaries unless exogenous variables are observable. Moreover, these methods may have low efficiency due to their reliance on independence tests. To obtain the counterfactual Markov boundary, in this work, we assume that another representation of SCM graphs called augmented graphs is known. Compared to causal graphs, augmented graphs incorporate exogenous variables U into the nodes, then remove bidirectional edges and add unidirectional edges(U j → V i ) for all V i ∈ U and U j ∈ U Vi .When the SCM is recursive, the augmented graph of any submodel is necessarily acyclic due to the closure of recursiveness under (perfect) intervention [8] , which implies that d-separation [75] always holds. We additionally assume faithfulness [2] , which means that only the variables that are d-separated in the graph are independent in the distribution, allowing us to directly compute counterfactual Markov boundary for any exogenous variable by graph algorithm.Under the above assumptions, the following theorem indicates that the counterfactual Markov boundary of U j is independent across different submodels. This not only reduces the search space of the Markov boundary, but also aligns with Eq. 12, where we encode counterfactual information according to each submodel.Theorem 2 (Counterfactual Markov Boundary Independence). If Y * = k i=1 Y i[xi]and each Y i[xi] corresponds to a different submodel M xi , then for each U j ∈ U, there exists a Markov boundaryB j (Y * ) = k i=1 B j (Y i[xi] ) on Y * , where B j (Y i[xi]) is a Markov boundary on M xi . Another theorem demonstrates how to obtain counterfactual Markov boundaries via d-separation, and indirectly proves their uniqueness: Theorem 3 (Counterfactual Markov Boundary on Graph). For an exogenous variable U j ∈ U and a counterfactual variable set Y x from the submodel M x , the counterfactual variableY x ∈ B j (Y x ) if and only if Y x ̸⊥ ⊥ G a x U j | Y x \ {Y x }, i.e., when given Y x \ {Y x }, Y x and U j are not d-separated on G ax , where G a x is the augmented graph induced from the submodel M x . The graph algorithm to determine counterfactual Markov boundaries B j (Y * ) for U j ∈ U is derived simply by combining Thms. 2 and 3 with the d-separation criterion as a subroutine, which can be found in Fig. 7 in App. A.4.We then inject the counterfactual Markov boundary into the neural network used for conditioning by vectorizing it as a weight mask, as briefly illustrated in Fig. 2 depicting the entire conditioning and masking process. Specifically, let θ y * ,j be any parameter corresponding to exogenous variable U j ∈ U in Q U|y * . According to the description in Thm. . This corresponds exactly to encoding the information for counterfactual events on each submodel M xi as depicted in Eq. 10. Therefore, we can modify Eq. 12 as follows:EQUATIONIn particular, m ij = ω(B i (Y j )), and B i (Y j ) are obtained through algorithm derived by Thm. 3. The correspondence between the parameters θ y * ,j and U j exists in the specific design of the model, such as the mean and covariance of each component in GMMs or the element-wise transformations in normalizing flows. As used in [31, 81, 13] , h is an MLP that allows weight masking, such that for all j, the i-th component of ∇ x h j (x, m ij ) ̸ = 0 if and only if m ij = 1. This ensures that θ y * ,j is only related to the Markov boundary of U j on y * . 
related works: 
 Importance sampling with normalizing flows In recent years, the combination of normalizing flows and importance sampling has been widely employed in tasks related to Monte Carlo integration. For instance, a plethora of literature [70, 99, 101, 53, 60, 10, 51, 67] utilizes this combination to efficiently compute the partition function of energy. Some studies [68, 30, 36] also employ this method to solve integration problems of complex functions in high-dimensional spaces. When the integrand involves rejection sampling [6, 89] , particle events [29, 9, 88] , or failure events [33, 20] , this method can improve the efficiency of sampling rare events. This combination is also applied for posterior inference [19, 1, 79, 21] as an alternative to integrating over the intractable distribution.Identifiable neural proxy SCMs Our approach is applicable to pre-trained neural proxy SCMs, where causal mechanisms are modeled as neural networks. Increasing attention has been paid to the identifiability of these models. BGM [69] combines bijections to demonstrate the identifiability of causal mechanisms in several special cases. CausalNF [42] builds upon the conclusions of [106] to prove the identifiability of its causal mechanisms up to invertible functions. NCM [108] extends the findings of their previous work [107] and proposes a sound and complete algorithm to identify counterfactual queries.Arbitrary conditioning In dealing with exponentially many conditional distributions, our work establishes a connection with another unsupervised learning task, known as arbitrary conditioning.The task of this paper can be viewed as extending arbitrary conditioning from a single observational distribution to multiple counterfactual distributions. Relevant works include VAE-AC [40] based on VAE, NC [7] based on GAN, ACE [90] based on energy models, and ACFlow [59] based on normalizing flows. The most similar to our work is [91] , which matches on the posterior distribution of VAE, akin to our matching on the exogenous distribution. 
experiments: 
 Experiment settings and metrics We first define the prior distribution Q S over the state space of a stochastic counterfactual process Y * and train according to Eq. 8. Subsequently, based on Q S , we construct a distribution P Y * concerning counterfactual events for the evaluation of the estimator (Eq. 4). Inspired by the metrics for rare event sampling effectiveness, we utilize Effective Sample Proportion (ESP) and Failure Rate (FR) with threshold m to measure the performance of the importance sampling, defined as:EQUATIONwhere Our experiments involve the following fully specified SCMs, which fall into three categories: i) Markovian diffeomorphic [42] , where there is no hidden confounder and causal mechanisms are diffeomorphic, including SIMPSON-NLIN, TRIANGLE-NLIN and LARGEBD-NLIN; ii) Semi-Markovian continuous [107] , where hidden confounders may exist and endogenous variables are continuous, including M and NAPKIN; iii) Regional canonical [108] , where hidden confounders may exist and endogenous variables are discrete and finite, including FAIRNESS and FAIRNESS-XW. The consistent performance observed in experiments across these different types of SCMs will serve as a reference for evaluating the robustness of the proposed method. In most cases, EXOM with Markov boundaries masked (orange bar) exhibits superior performance compared to when the Markov boundaries are not masked (blue bar).η(Y * ) = n i=1 1 Ω U (Y * ) (u (i) ) /nConvergence We integrated Exogenous Matching with various density estimation models (GMM, MAF [73] , NICE [24] , SOSPF [41] ) to conduct experiments on the stochastic counterfactual process Y B * (with |s| = 3) over different SCM settings. We present the results of ESP and FR during the training process. As shown in Fig. 3 (with all results detailed in App. C.5), both the ESP and FR metrics change as expected until convergence, along with the decrease in training loss Eq. 8 (equivalent to the increase in LL in Fig. 3 ). This validates the effectiveness of Thm. 1 and Cor. 1, and empirically suggests that the conditions we assumed are generally not violated under our settings.Comparison We compare our proposed method with three different but related sampling methods: i) Rejection Sampling (RS); ii) Cross-Entropy based Importance Sampling (CEIS, [32] ); iii) Neural Importance Sampling (NIS, [68] ). To apply the latter two methods under the same experimental setup, we extended them as detailed in App. B.1. As shown in Tab. 1, our method outperforms other approaches in both continuous and discrete cases. Specifically, since MAF exhibits stronger representational capabilities than GMM, it is expected that EXOM using MAF as the conditional distribution model generally outperforms EXOM using GMM. Here, "O" represents the original SCM, and "P" represents the proxy SCM. For SIMPSON-NLIN, the proxy SCM is CausalNF, and the metric used is FR ("-" indicates FR equals 1); whereas for FAIRNESS, the proxy SCM is NCM, and the metric used is the average bias w.r.t. the ground truth. The subscript denotes the 95% CI error bound over 5 trials. For more details, see App. C.9 4 , the use of Markov boundaries as masks significantly improves the performance of EXOM under various settings. We also observed that when masks are used, the width of the hidden layers of neural networks in conditioning affects practical performance, which will be discussed in App. C.8.We combine identifiable proxy SCMs with EXOM for counterfactual estimation. We employ CausalNF [42] for SIMPSON-NLIN (based on the stochastic process Y B * ) and NCM [108] for FAIRNESS (based on the stochastic process Y Q * ). For the combination of CausalNF and SIMPSON-NLIN, we estimate the counterfactual probability on cubes δ l (y * ) with side length l = 0.02 centered around 1024 randomly sampled y * from P Y * (this can be used further to estimate counterfactual density), measure the sampling FR, and employ a dimensionregularized 95% CI error bound. For the combination of NCM and FAIRNESS, we estimate 4 different counterfactual queries, measure the average bias of the query results relative to the ground truth, and use a 95% CI error bound to demonstrate the unbiasedness of the estimates.The results in Tab. 2 empirically demonstrate the effectiveness of EXOM in counterfactual density and effect estimation tasks compared to RS. Specifically, when estimating the cube δ l (y * ) (equivalent to counterfactual density), RS almost fails in high-dimensional settings, while EXOM exhibits good sampling efficiency and lower error compared to RS. When estimating counterfactual queries, the bias and error of EXOM are similar to those of RS (partly due to the lower sampling difficulty in discrete cases). Further experimental details and analysis of the conclusions can be found in App. C.9. 
conclusion: 
 We propose Exogenous Matching for tractable estimation of expressions in L 3 in general settings. Specifically, leveraging importance sampling, we find a variance upper bound for estimators potentially applicable to any relevant counterfactual probability (Cor. 1), and then introduce an optimization objective (Eq. 8) amenable to sampling, transforming the problem into a conditional distribution learning problem. Theoretical and empirical results demonstrate that this approach can be combined with identifiable proxy SCMs to address practical problems. We also explore injecting Markov boundaries as prior knowledge and empirically validate effectiveness in several scenarios.Limitations i) The proposed method does not directly estimate counterfactuals through the available distribution but requires a partially specified SCM. This work explores its effectiveness in combination with identifiable proxy SCMs, which are still an active area of research; ii) It is important to emphasize that the "general case" mentioned in this paper does not cover all scenarios. For instance, our estimand, an L 3 expression Q, requires the assumptions of finite submodels and approximability by a finite number of counterfactual probabilities; iii) Furthermore, the specific structures and tricks demonstrated in this paper (such as the injection of Markov boundaries) may not be applicable to SCMs with general parameter settings, which warrants further investigation. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments? Answer: [Yes] Justification: In the experimental figures and tables, we present error bars, with 95% CI and quartile outliers. Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments? Answer: [Yes] Justification: In App. C.4, we describe the computation time and computational resources. 
: 
 Counterfactual reasoning, considered one of the advanced cognitive abilities of humans, aims to address questions about hypothetical worlds that have not been observed. Answering these questions, typically in the form of "what-if" or "why," is crucial for attribution and explanation, thus counterfactual reasoning is widely applied in generating explanations [48, 47, 94, 63] , making decisions [71, 37, 95, 58] , and evaluating fairness [56, 113, 105, 78, 117] . A class of generative models called Structural Causal Models (SCMs) provides a semantic characterization of counterfactuals, consisting of a set of endogenous mechanisms and an exogenous distribution, from which a collection of distributions described by specific formal languages is induced. These languages are organized into a progressively refined hierarchical structure, and each is associated with human activities: language L 1 is about seeing (observational), L 2 is about doing (interventional), and L 3 is about imagining (counterfactual). This hierarchy is referred to as the causal ladder [77] or the Pearl Causal Hierarchy (PCH) [5] . Fig. 1 provides a brief illustration of the relevant concepts. As the highest level of the PCH, counterfactuals entail the most subtle information in the hierarchy, making the evaluation of counterfactual languages an appealing objective that has received widespread attention in recent years.Counterfactual reasoning consists of two subtasks: counterfactual identification and counterfactual estimation. For an expression Q in L 3 , identification is to determine the uniqueness of the answer to Q [86, 17] or its bounds [114, 25, 66] under some assumptions, while estimation is to compute the specific value of Q. A series of literature [86, 17] transforms Q into a combination of expressions in L 1 and L 2 through identification. However, estimating the transformed expression may be intractable.Another series of literature focuses on identifying and estimating Q in specific settings, such as back door [82, 38, 115, 102] , front door [43, 109, 98] , instrumental variables [14, 100, 35] , representation decoupling [15, 97, 61] , categorical exogenous [22] , exponential family [85] , domain counterfactuals [116] . There are also some literature applicable to relaxed cases, such as [28, 44, 45, 46, 57] , but still less developed.In recent years, a class of methods has emerged that simulate the structure of the original SCM using neural network architectures and learn causal mechanisms based on observed or intervened distributions, which we refer to as neural proxy SCMs. These methods include the use of VAEs [110, 55, 54, 92, 23] , GANs [52, 108] , normalizing flows [74, 50, 4, 69, 42] , DDPMs [84, 12] . Here, we focus only on models that provide identifiability results. Using bijections, [69, 42] demonstrates the identifiability of learned causal mechanisms. However, constrained by bijections, they only provide tractable counterfactual estimation methods for fully observed endogenous variables. The closest method to performing tractable counterfactual estimation in general cases is NCM [107, 108] , which employs an optimization algorithm to concurrently achieve tractable counterfactual estimation and identification (or partial identifiability [93] ). However, the counterfactual estimation depends on rejection sampling, thereby lacking scalability.It is commonly recognized that performing tractable counterfactual estimation in general cases is challenging, even when the SCM is fully specified. From the perspective of definition [5] , one needs to integrate over the entire exogenous space, which involves determining whether the constraints of different hypothetical worlds are simultaneously satisfied. From the perspective of complexity, although [34] has demonstrated that counterfactual reasoning is tractable if observational and interventional reasoning are tractable, it relies on algorithms with exponential complexity. Furthermore, according to [112] , performing marginal inference on SCMs, including parameterized SCMs such as NCM, has been shown to be NP-hard.This paper presents a tractable and efficient counterfactual estimation method in general settings based on importance sampling. Here, by "general settings", we mean that the exogenous variables of the SCM can be discrete or continuous, the expressions to be estimated can involve an arbitrary finite number of hypothetical worlds, and the observations and interventions can also be arbitrary.Specifically, the contributions of this paper are as follows:1. We propose a tractable and efficient importance sampling method for counterfactual estimation in general settings, which is based on optimizing an upper bound on the variance of the estimators, and is formulated as a conditional distribution learning problem. 2. We explore the injection of structural prior knowledge (counterfactual Markov boundaries) into neural networks used for parameter optimization. 3 . The experiments conducted on different SCM settings demonstrate that the proposed method outperforms existing importance sampling methods. Ablation studies empirically highlight the effectiveness of injecting structural prior knowledge. The experiments on counterfactual estimation tasks conducted on two identifiable proxy SCMs illustrate the feasibility of applying the method in conjunction with proxy SCMs to real-world problems.probability when X take values X ⊆ Ω X . The lowercase p(x) represents the probability density function of P X if Ω X is continuous.Structural causal models An SCM is a 4-tuple M = ⟨U, V, F, P U ⟩, where U is a set of exogenous variables; V is a set of endogenous variables; F, which describes the causal mechanisms of M, is a collection of functions {f V1 , f V2 , . . . , f Vn } such that each endogenous variable V i is determined by a function f Vi . Each f Vi ∈ F is a mapping from (the domains of) U Vi ∪ Pa Vi to V i , where U Vi ⊆ U and Pa Vi ⊆ V \ V i . The entire set F forms a mapping from U to V, while the uncertainty comes from a probability distribution P U over exogenous variables U.An SCM yields a causal graph G, where each V i ∈ V is a vertex, and the edges come in two types: there is a directed edge (V i → V j ) between each pair V i and V j ∈ Pa Vi ; and there is a bidirected edge (V i ↔ V j ) between each pair V i and V j if U Vi ∪ U Vj ̸ = ∅. This resulting graph is also known as a directed mixed graph (DMG). Here, we assume that the SCM is recursive, which means that there exists an order in F such that for any pair f i , f j ∈ F, if f i < f j , then V j / ∈ Pa Vi . In particular, this implies that given values for exogenous variables (also known as a unit or an individual) u ∈ Ω U , all values of endogenous variables in V are also fixed (i.e., there exists a unique solution for endogenous variables w.r.t. u), and the resulting causal graph G is an acyclic directed mixed graph (ADMG). We denote Y(u) as the solution for Y ⊆ V given u.In the paradigm of SCM, an (perfect) intervention is described as replacing the mechanisms of some variables X with constants x. The model after this intervention is termed a submodel, defined asEQUATIONIn submodel M x , the solution of Y ⊆ U given u, denoted as Y Mx (u), is termed as the potential response, and is typically abbreviated as Y x (u).The random variable Y ∈ V corresponding to the submodel M x is denoted as Y x , termed a counterfactual variable, with its domain denoted as Ω Yx .A set of counterfactual variables from the same submodel M x is denoted as Y x , and the set of counterfactual variables from k different submodels M x1 , . . . , M x k is denoted asY * = k i=1 Y i[xi] ,with the corresponding domainΩ Y * = × k i=1 Ω Y i[x] . Let σ-algebra Σ Y * ⊆ 2 Ω Y *, where 2 Ω Y * is the power set of Ω Y * . The counterfactual probability space is then a triplet ⟨Ω Y * , Σ Y * , P Y * ⟩, where P Y * is a probability measure P Y * : Σ Y * → [0, 1], and for any measurable set Y * ∈ Σ Y * , its value equals the Lebesgue integral, Counterfactual probability is the cornerstone of the language L 3 in PCH. Specifically, all expressions in L 3 consist of inequalities between polynomials (and their Boolean combinations) over terms of the form P (Y * ). Here, we assume that the L 3 expressions to be estimated can be represented by a finite number of P (Y * ) terms (including those approximated by Monte Carlo methods). Therefore, the estimation of P (Y * ) will become a focal point of the subsequent discussion.EQUATIONThe normalizing flow T θ : R d → R d is a transformation parameterized by θ that maps observable samples x with d features to latent samples z distributed according to a simple distribution P Z with probability density function p(z). The mapping T θ is a diffeomorphism, which allows us to compute the density p(x) through the change-of-variables formula: log p(x) = log p(T θ (x)) + log |det(∇ x T θ (x))| .(EQUATIONWhen the Jacobian determinant of T θ (x) is tractable to compute, we can use samples from the observable distribution and maximum likelihood estimation to train the flow T θ . Additionally, we can easily generate i.i.d. samples conforming to the learned observable distribution by first sampling z (i) from the base distribution P Z , and then evaluating the inverse transformation T -1 θ (z (i) ). 
 
 Appendix: In this section, we demonstrate the validity of counterfactual probability-related concepts, providing intuitions based on derivations in discrete spaces, and proving lemmas to streamline subsequent derivations. Firstly, we establish the validity of counterfactual probabilities as defined in Eq. 2.Proposition 1. The probability measure defined by Eq. 2 (on recursive SCMs) adheres to the Kolmogorov axioms.Proof. (i) For any Y * ∈ Σ Y * , 0 ≤ P Y * (Y * ) ≤ 1. Since 0 ≤ 1 Ω U (Y * ) (u) ≤ 1, the proof followsfrom the Monotonicity of Lebesgue integration:EQUATIONwhereP Y * (Y * ) = Ω U 1 Ω U (Y * ) (u) dP U , thus 0 ≤ P Y * (Y * ) ≤ 1. (ii) P Y * (Ω Y * ) = 1.Since SCM is recursive, given any u ∈ Ω U , the values of exogenous variables are determined, and hence all values in Y * are determined, that is, the potential responses Y * (u).EQUATIONEQUATION=EQUATION=EQUATION= i≥1 P Y * (Y (i) * ) (23) If for all i, Y * (u) / ∈ Y (i) * , then i≥1 1 Ω U (Y (i) * ) (u) = 1 Ω U ( i≥1 Y (i) * ) (u) = 0.Therefore, both sides of the equation are zero, so it holds. Hence, we have provedP Y * ( i≥1 Y (i) * ) = i≥1 P Y * (Y (i) * ).Given that axioms (i), (ii) and (iii) hold, which are the Kolmogorov axioms, it follows that Eq. 2 satisfies the definition of a probability measure.From a measure-theoretic perspective, the indicator function in Eq. 2 can be transformed into the Dirac measure, and vice versa:EQUATIONProof. Given Y * ∈ Σ Y * and u ∈ Ω U , we can consider Y * as u being constant, thusEQUATIONEQUATION=EQUATIONwhere the integral Y * dδ Y * (u) remains a function of u.Then, substituting Lem. 1 into Eq. 2, we have:EQUATION=EQUATION=EQUATIONwhere P Y * ,U satisfying dP Y * ,U = dδ Y * (u) dP U is called the joint distribution of counterfactual variables and exogenous variables, which we abbreviate as the joint distribution.In words, the joint distribution is determined by the exogenous distribution and the Dirac distribution, while the counterfactual probability is the integral over the endogenous space after marginalizing out the exogenous variables from the joint distribution.Importance sampling estimator for counterfactual estimation Eq. 4 only illustrates the case where exogenous variables are continuous. Here, we will provide a complete derivation. For discrete cases, the Lebesgue integral can be expressed as a summation, thus,EQUATIONEQUATIONEQUATIONFor continuous cases, the Radon-Nikodym theorem is required. Suppose that the probability density functions of P U and Q U exist, denoted by p(u) and q(u), respectively. According to the definition of probability density functions, p(u) = dP U / du and q(u) = dQ U / du. ThenEQUATION=EQUATIONWe assume that the integrand is Lebesgue integrable. It can be observed that the difference between the discrete and continuous cases lies only in the importance weights. Therefore, in subsequent derivations, we will proceed in the continuous case, with the discrete case requiring only the substitution of the importance weights.The expectation term can be estimated using Monte Carlo methods, resulting in the basic form of the importance sampling estimator:EQUATIONwhereσ Y * (u) = (p(u)/q(u))1 Ω U (Y * ) (u), and each sample u (i) ∼ P U .Variance of the importance sampling estimator The variance of the estimator in Eq. 4 is:V 1 n n i=1 σ Y * (u (i) ) = 1 n 2 n i=1 V u (i) ∼Q U σ Y * (u (i) ) Linearity (40) = n n 2 V u (i) ∼Q U [σ Y * (u)] i.i.d (41) = 1 n V u∼Q U [σ Y * (u)] .(5, the first equation)V u∼Q U [σ Y * (u)] can be further expanded as:EQUATIONEQUATIONEQUATION=EQUATION=EQUATIONEQUATION= E u∼P U [σ Y * (u)] -P 2 (Y * ).(5, the second equation)The proposal distribution Q U only affects the termE u∼P U [σ Y * (u)]in the variance Eq. 5. The only difference between it and the estimator P(Y * ) = E u∼Q U [σ Y * (u)] is that E u∼P U [σ Y * (u)]is the expectation over the exogenous distribution P U .Optimization of Eq. 5 and the limitations Reducing variance implies improving estimation efficiency; in other words, it means achieving the same level of estimation accuracy with fewer samples. One direct method to optimize variance is by sampling from the exogenous distribution P U using Monte Carlo methods to estimate the term E u∼P U [σ Y * (u)] in the variance Eq. 5, and then minimizing the estimated value of this term using stochastic optimization techniques. Alternatively, as suggested by [32, 68] , this term can be transformed into an optimization problem involving the cross-entropy or χ 2 divergence between the proposal distribution and the integrand, and then optimizing the values of these mismatches estimated by sampling under the proposal distribution Q U can also indirectly optimize variance. This optimization process is referred to as proposal learning.However, a learned individual proposal distribution is typically only applicable to a single estimator, meaning it can only estimate the probability for a single counterfactual event Y * , as the indicatorfunction 1 Ω U (Y * ) (u) in σ Y * (u)operates solely for Y * . But L 3 expressions can typically contain multiple counterfactual probability terms, and this limitation implies that, in the worst case, we need to perform a separate optimization for each counterfactual probability, still affecting the efficiency of the counterfactual estimation task.To address this issue, our solution is to replace the sampling proposal distribution Q U with a conditional proposal distribution Q U|Y * , adding an extra degree of freedom conditioned on Y * . Our overall idea is to let these different proposal distributions participate in the same optimization (i.e., sharing the same optimization objective) so that a single optimization yields different proposal distributions applicable to different counterfactual events.Additionally, to limit the scope of responsibility of each proposal distribution Q U|y * -that is, to prevent too many proposal distributions from targeting the same counterfactual event Y * and consequently neglecting some counterfactual events -we further impose a one-to-one correspondence between conditions and responsibilities for inference. Intuitively, this would make the proposal distribution Q U|y * corresponding to y * ∈ Ω Y * focus on the support Ω U (δ(y * )). Given that for any counterfactual event there is an approximation Y * ≈ ∞ y * ∈Y * δ(y * ), we can use the proposal distributions Q U|y * corresponding to these y * ∈ Y * to serve the estimation of P (Y * ).Variance upper bound When certain conditions are met in Y * , for any y * ∈ Y * , the proposal distribution Q U|y * shares a common upper bound on variance, which is exactly the common optimization objective we seek:Theorem 1 (Variance Upper Bound). Let σ Y * (u) = (p(u)/q(u | y * )) 1 Ω U (Y * ) (u), where q(u | y * )denote the density of the proposal distribution Q U|y * , and let Y * (u) be the potential response w.r.t. u. If for any y * ∈ Y * , there exists κ ≥ 1 such that 1/κ ≤ p(u)/q(u | y * ) ≤ κ holds almost surely on the support Ω U (Y * ), then for any Q U|y * where y * ∈ Y * :EQUATIONwhere the constant c is solely dependent on κ and P U .Proof. For a Y * ∈ Σ Y * that satisfies bounded importance weights (i.e., there exists κ ≥ 1 such that 1/κ ≤ p(u)/q(u | y * ) ≤ κ holds almost surely on the support Ω U (Y * )), the following inequality holds for any y * , y ′ * ∈ Y * holds almost surely:EQUATION)In the following derivations, cases where Eq. 48 does not hold will be implicitly omitted. This is because, according to the assumption of bounded importance weights, the exogenous distribution's probability measure corresponding to the cases where Eq. 48 does not hold is zero (in other words, the probability that the inequality Eq. 48 holds is 1). Therefore, omitting the cases where it does not hold does not affect the expectation w.r.t. the exogenous distribution.For ease of implementation and computation, we focus on log-weights by constructing ξ Y * (u) to bridge the inequality conclusions over E u∼P U [ξ y * (u)] and the inequality relationship with the variance, whereEQUATIONNext, we look for a lower bound forE u∼P U [ξ Y * (u)].Firstly, according to the bounded density ratio condition, for any u ∈ Ω U (Y * ) with p(u) > 0 and y * ∈ Y * , κ -1 ≤ p(u)/q(u | y * ) ≤ κ holds. Therefore, based on the concavity of the logarithmic function, the following inequality holds (assuming κ > 1, as the conclusion is trivially true when κ = 1):EQUATION)Due to the boundedness of E u∼P U [ξ Y * (u)] in the logarithmic space (with ≤ log κ and ≥log κ), by employing the monotonicity and linearity of the Lebesgue integral, we obtain:EQUATIONEQUATIONEQUATIONEQUATION)Substitute the above inequalities into Eq. 5:EQUATIONEQUATIONEQUATIONEQUATION≤ - κ 2 -1 2κ log κ E u∼P U [log q(u | Y * (u))] + κ 2 -1 2κ log κ E u∼P U [log p(u)] + 3κ 2 -1 2κEq. 55 (65)EQUATIONThe constant c is solely dependent on κ and P U , whereas E u∼P U [log p(u)] equals the negative entropy of P U .The bounded importance weights assumption (i.e., there exists κ ≥ 1 such that 1/κ ≤ p(u)/q(u | y * ) ≤ κ holds almost surely on the support Ω U (Y * )) may be too strong for general cases, especially when the exogenous distribution has an infinite support set and Ω U (Y * ) is also improper. For instance, even for two Gaussian distributions, their density ratio cannot be guaranteed to be bounded ( [18] , E.g. 4.1).One remedy to mitigate this strong assumption is to weaken the statement of the theorem from "almost surely" to "with probability at least 1δ." In this way, if δ is small, the conclusion can still hold with high confidence. To reach such a conclusion, we typically need to use concentration inequalities and introduce other weaker assumptions. One of such weaker assumptions we introduce is that for any y * ∈ Ω Y * , the expected value of the square of the importance weights underP U is finite, that is, E u∼P U p 2 (u)/q 2 (u | y * ) < +∞.Another assumption is about the boundedness of the importance weights in the empirical distribution P U ; this is a discrete relaxation of the original boundedness assumption of importance weights in the exogenous distribution, which makes Thm. 1 trivially hold on the empirical distribution.Corollary 2 (Thm. 1 with Cantelli's Upper Bound). Assume that for any y * ∈ Y * , the expectationE u∼P U p 2 (u)/q 2 (u | y * ) < +∞. Let ν = sup y * ∈Y * E u∼P U p 2 (u)/q 2 (u | y * ). Denote by P U the empirical distribution formed by n i.i.d. samples u (i) ∼ P U , i = 1, . . . , n drawn from P U , and assume there exists a constant κ ≥ 1 such that 1/κ ≤ p(u (i) )/q(u (i) | y * ) ≤ κ for all i = 1, . . . , n.Then, with probability at least 1δ, the following inequality holds:EQUATIONwhere the constant c is the same as in Thm. 1.Proof. First, consider the empirical distribution P U formed by n i.i.d. samples u (i) ∼ P U . Given the existence of the bound κ for these samples, meaning that the importance weights exist on the support of P U , it is evident that Thm. 1 holds for P :EQUATIONwhere V denotes the empirical variance, which is calculated as:EQUATIONwhere we treat P (Y * ) as a constant.Next, we need to establish the inequality relationship between the empirical variance and the true variance. According to Cantelli's inequality, we have:EQUATION)Reversing the inequality direction, substituting the random variable X with E [X] = 1 n n i=1 X i (where each X i is i.i.d.), and letting λ =1 δ -1 V[X]n , we obtain:EQUATION)Let X i = σ Y * (u (i) ), and subtract P (Y * ) from both sides, then we have:EQUATION)whereEQUATIONEQUATION≤ supEQUATIONSubstituting Eq. 68 into this equation, we obtain the result:EQUATIONwhere the constant c is the same as in Thm. 1.The following obtains a tighter bound after introducing the constraint of sub-Gaussian conditions (specifically, assuming that the two-sided Bernstein's condition is satisfied):Corollary 3 (Thm. 1 with Bernstein's Upper Bound). Assume that for any y * ∈ Y * , the expectationE u∼P U p 2 (u)/q 2 (u | y * ) < +∞. Let ν = sup y * ∈Y * E u∼P U p 2 (u)/q 2 (u | y * ). Additionally, assume that the variable σ Y * (u) satisfies Bernstein's condition, meaning that there exists a parameter b > 0 such that for any λ ∈ (-1/b, 1/b), the following inequality holds:EQUATION)Denote by P U the empirical distribution formed by n i.i.d. samples u (i) ∼ P U , i = 1, . . . , n drawn from P U , and assume there exists a constant κ ≥ 1 such that 1/κ ≤ p(u (i) )/q(u (i) | y * ) ≤ κ for all i = 1, . . . , n. Then, with probability at least 1δ, the following inequality holds:EQUATION)where the constant c is the same as in Thm. 1.Proof. When Bernstein's condition is satisfied, the following inequality holds:EQUATION)The remaining derivation is similar to that in Cor. 2.It is trivially provable that 1 δ -1 ≫ 2 log 1 δ when δ approaches zero. This implies that, for the same δ close to zero, the remainder term in Cor. 3 is typically smaller than that in Cor. 2, meaning that the former provides a much tighter upper bound with a high probability.A guard proposal distribution to ensuring boundedness For a proposal distribution Q U|y * , define the bounded part asΩ [κ] U = {u | 1/κ ≤ p(u)/q(u | y * ) ≤ κ}.If the probabilities of the unbounded part under the exogenous distribution and the proposal distribution,P U Ω U \ Ω [κ] U = E u∼P U 1 Ω U \Ω [κ] U and Q U|y * Ω U \ Ω [κ] U = E u∼Q U|y * 1 Ω U \Ω [κ] U, are known, we can splice the exogenous distribution over the unbounded part to form a new proposal distribution that ensures the bounded importance weight condition: Proposition 2 (Importance Weight Bound Guard). For a proposal distribution Q U|y * , construct a new proposal distribution as follows:EQUATIONwhich we refer to as the guard of the proposal distribution Q U|y * . Then, for any u ∈ Ω U , there exists a constantEQUATIONsuch that 1/κ ′ ≤ q [κ] (u | y * )/p(u) ≤ κ ′ .Proof. First, we prove that the constructed proposal distribution is a valid distribution, mainly by showing that it satisfies the normalization condition:EQUATION=EQUATIONNext, we prove the boundedness: when u ∈ Ω[κ]U , according to the definition ofΩ [κ] U , we have 1/κ ≤ p(u)/q(u | y * ) ≤ κ. Since κ ′ ≥ κ, it follows that 1/κ ′ ≤ p(u)/q [κ] (u | y * ) ≤ κ ′ . When u / ∈ Ω [κ]U , it trivially holds that p(u)/q [κ] (u | y * ) = 1.In the subsequent experiments, we assume that the probability of the unbounded part could be non-zero but approximates zero. Formally,P U Ω [κ] U ≈ 1 and Q U|y * Ω [κ] U ≈ 1.This implies that the error terms in Cor. 2 and Cor. 3 would be small enough, or that when using the guard proposal distribution Prop. 2, we could obtain a relatively unbiased estimate without the need for re-weighting (where1 -P U Ω U \ Ω [κ] U / 1 -Q U|y * Ω U \ Ω [κ] U ≈ 1).The generalization of Thm. 1 Cor. 1 is generalized from Thm. 1 in a straightforward manner through monotonicity. In particular, the relaxation and guarantees discussed earlier also apply to Cor. 1 when the scope of assumptions extends from a single counterfactual variable Y * to multiple counterfactual variables of different forms, Y EQUATIONwhere the constant c is the same as in Thm. 1.Proof. For any state s ∈ S, and anyY (s) * ∈ Σ (s)Y * , according to Thm. 1, if Y (s) * satisfies the bounded density ratio condition, we have:EQUATIONAccording to the monotonicity of expectation, we take the expectation w.r.t. PY * on both sides, then:EQUATIONEQUATIONEQUATION= -EQUATIONwhere We then take the expectation w.r.t. P U on both sides of the inequality, producing:E u∼P U log q(u | Y (s) * (u)) is independent of Y (s) * ,EQUATIONEQUATION= -EQUATION= -EQUATIONwhere c is a constant independent of any state s ∈ S.A. for i ← 1 and i ≤ n do 3: s ∼ Q S , u ∼ P U 4: θ y * ← θ y * -η • ∇ θy * log q(u | Y (s) * (u)) 5:end for 6: end while Learning Cor. 1 provides an upper bound that directly gives an optimization objective:arg min q -E s∼QS E u∼PU log q(u | Y (s) * (u)) . (8)This is equivalent to learning the conditional distribution Q U|Y * using maximum likelihood.To obtain an estimate of this log-likelihood, we can directly use Monte Carlo methods, which simply involve two sampling processes and one computation: i) sampling state s from the prior state distribution Q S ; ii) sampling u from the exogenous distribution P U ; iii) computing the log density (or probability for discrete exogenous variables) log q(u | Y (s) * (u)). Assuming that we can compute the gradient of q(u | y * ) w.r.t. the parameters θ y * of Q U|y * , as modeled in the paper, then combining with the minibatch gradient descent algorithm, we can write the complete algorithm for the training process as in Fig. 5 . Inference If for any y * ∈ Y * , the support of the proposal distribution Q U|y * covers the support of P U , then it constitutes a valid proposal, and an unbiased estimate can be obtained using Eq. 4. More can be achieved through multiple importance sampling, as any estimate obtained from the proposal Q U|y * is unbiased for such y * ∈ Y * . Consequently, the expectation over these estimates remains unbiased, i.e.,Input: Prior distribution Q Y * , conditional pro- posal distribution Q U|Y * , sample size n. 1: S ← ∅ 2: for i ← 1 and i ≤ n do 3: y (i) * ∼ Q Y * , u (i) ∼ Q U|y (i) * 4: σ (i) ← p(u (i) )/q(u (i) | y * ) 1 ΩU(Y * ) (u (i) ) 5: S ← S ∪ {σ (i) } 6: end for 7: Output: P (Y * ) ← n i=1 σ (i) /nEQUATIONEQUATIONwhich corresponds to the estimator used in Eq. 9 and provides more robust results when an appropriate prior Q Y * is chosen. Our inference algorithm involves two sampling processes and one computation: i) sampling y * from the prior distribution Q Y * ; ii) sampling u from the proposal distribution Q U|y * ; iii) computing the importance weights and checking if Y * (u) ∈ Y * to obtain σ Y * [y * ] (u). This process is detailed in Fig. 6 . In our experiments, because the tasks are density estimation (for continuous) and point probability queries (for discrete), we actually still use a single proposal distribution at y * , which is a special case where Q Y * is a Dirac distribution.In order to express counterfactuals that span multiple submodels using a single model, the twin SCM concept is often used to articulate counterfactuals for a particular setting within a unified framework.Definition 3 (Twin SCM). For an SCM M = ⟨U, V, F, P U ⟩ and a set of its submodels{M 1[x1] , M 2[x2] , . . . , M k[x k ] }, we refer to M * = ⟨U, V * , F * , P U ⟩ as a Twin Structural Causal Model (Twin SCM), where V * = k i=1 V i[xi] and F * = k i=1 F i[xi] . V i[xi]represents the endogenous variables in the i-th submodel, and F i[xi] denotes the corresponding causal mechanisms.In general, the twin SCM amalgamates multiple distinct submodels into one SCM, thereby rendering conclusions drawn on a single SCM applicable to the twin SCM. It is trivial to demonstrate that if the original SCM is recursive, then the twin SCM remains recursive. Consequently, the causal graph derived from the twin SCM is also an ADMG, and the augmented graph is a DAG.Our proof is based on d-separation [75, 76] , which is a criterion for quickly determining conditional independence on a Directed Acyclic Graph (DAG). Suppose a distribution is denoted by P and its corresponding DAG is G. We say that a path p is blocked by a set of nodes Z if and only if: i) p contains a chain A → B → C or a fork A ← B → C and the intermediate node B is in Z; ii) p contains a collider A → B ← C,P . If U j is independent of Y \ X given X under P , i.e. U j ⊥ ⊥ P (Y \ X) | X, then X is termed a Markov blanket of U j on Y. The collection of all Markov blankets of U j on Y is denoted as B j (Y). If X ∈ B j (Y) is a Markov blanket of U jon Y, and there exists no X ′ ⊊ X such that X ′ ∈ B j (Y), then X is termed a (counterfactual) Markov boundary of U j on Y, denoted as B j (Y).This definition relies on conditional independence. With the aid of the augmented graph G a and the d-separation based on the twin SCM, it is possible to derive Thm. 2 and Thm. 3 from the properties on the graph.Theorem 2 (Counterfactual Markov Boundary Independence). If Y * = k i=1 Y i[xi] and each Y i[xi]corresponds to a different submodel M xi , then for each U j ∈ U, there exists a Markov boundaryB j (Y * ) = k i=1 B j (Y i[xi] ) on Y * , where B j (Y i[xi] ) is a Markov boundary on M xi .Proof. We aim to demonstrate thatk i=1 B j (Y i[xi] ) forms a Markov boundary for Y * , where B j (Y i[xi] ) denotes the Markov boundary over Y i[xi] . (i) k i=1 B j (Y i[xi]) is a Markov blanket of U j : For any submodel M j , according to the faithfulness assumption, sinceB j (Y i[xi] ) is a Markov blanket over Y i[xi], it blocks all paths between U j andY i[xi] \ B j (Y i[xi]) in the augmented graph of the twin SCM. We proceed to prove that, givenB j (Y i[xi] ), it does not alter the blocked paths of B j (Y i[xi]). For a blocked path p ∈ P, unblocking can occur only at the intermediate node and its descendants in the colliders. However, the nodes inB j (Y i[xi] ) are not in Y i[xi], thus demonstrating that its blocking status remains unchanged. Therefore,k i=1 B j (Y i[xi] ) is a Markov blanket of U j . (ii) Any proper subset of k i=1 B j (Y i[xi] ) is not a Markov blanket of U j : Suppose that a proper subset D ⊂ k i=1 B j (Y i[xi] ) is a Markov blanket. Let D j = D ∩ Y i[xi]be a subset in some submodel of this Markov blanket. Obviously, since any proper subset ofB j (Y i[xi] ) is not a Markov blanket in Y i[xi] , D j is not a Markov blanket on Y i[xi]. This implies that there exists a subset S j of Y i[xi] which is not independent of U j given D j . We now show that S j remains dependent on U j given D. By the faithfulness assumption, as S j is not independent of U j given D j , there still exists an unblocked path from S j to U j in the augmented graph of the twin SCM. Continuing with the conditionInput: An exogenous variable U j ∈ U, a set of counterfactual variables Y * , an augmented graph G a of SCM M = {V, U, F, P U }.1: B j (Y * ) ← ∅ 2: for Y i[xi] ⊆ Y * do 3: B j (Y i[xi] ) ← ∅ 4: for Y ∈ Y i[xi] do 5: if ¬ D-SEPERATION({Y xi }, {U j }, Y i[xi] \ {Y xi }, G a xi ) then 6: B j (Y i[xi] ) ← B j (Y i[xi] ) ∪ {Y } ▷ Thm. 3 7: end if 8:end for 9: on D \ D j , we discuss all the paths P from S j to U j . For an intermediate node of a collider on a path p ∈ P, since it is unblocked, its descendants must be conditioned, thus continuing to condition on D \ D j does not affect the colliders. For a non-collider on a path p ∈ P, for p to be blocked, there must exist an intermediate node of a non-collider in D \ D j , however, (D \ D j ) ∩ Y i[xi] = ∅, so this condition cannot be met. Hence, continuing to condition on D \ D j does not affect any unblocked paths. According to d-separation, given D, S j and U j remain dependent. Therefore, D is not a Markov blanket, which contradicts the assumption. This means that no proper subset ofB j (Y * ) ← B j (Y * ) ∪ B j (Y i[xi] ) ▷ Thm.k i=1 B j (Y i[xi] ) is a Markov blanket of U j .Theorem 3 (Counterfactual Markov Boundary on Graph). For an exogenous variable U j ∈ U and a counterfactual variable set Y x from the submodel M x , the counterfactual variableY x ∈ B j (Y x ) if and only if Y x ̸⊥ ⊥ G a x U j | Y x \ {Y x }, i.e., when given Y x \ {Y x }, Y x and U j are not d-separated on G ax , where G a x is the augmented graph induced from the submodel M x .Proof. (i) Y x ∈ B j (Y x ) implies Y x ̸⊥ ⊥ G a x U j | Y x \ {Y x }: Assuming Y x ⊥ ⊥ G a x U j | Y x \ {Yx }, according to the faithfulness assumption, all paths between U j and Y x are blocked given Y x \ {Y x }. If a path consists solely of colliders, then, since the intermediate node are given, the path remains unblocked, contradicting the assumption. This implies that on every path, there exists at least one non-collider intermediate node, denoted this set by R. Since each path has at least one non-collider intermediate node, given R, all paths between Y x and U j are blocked. LetR b = R ∩ B j (Y x ) and R c = R \ R b .For R c , the Markov boundary property holds, blocking all paths between U j and Y x given B j (Y x ). Thus, for every path between Y x and U j , it either passes through nodes in R b and is blocked given B j (Y x ) \ Y x , or it passes through nodes in R c . We need to prove that the paths that pass through R c are also blocked given B j (Y x ) \ Y x . Taking any R ∈ R c , for any path p between R and U j , if p does not pass through Y x , then p is blocked given B j (Y x ) \ Y x , which means that any path from Y * passing through R and then through p to reach U j is also blocked by B j (Y x ) \ Y x . If p passes through Y x , it must pass through another R ′ ∈ R and eventually reach U ′ through a path that does not include Y x , which is also blocked by B j (Y x ) \ Y x . Therefore, any path involving Y x is blocked by B j (Y x ) \ Y x , which means that any path from U j to Y x is also blocked. For colliders on ancestors of Y x , not conditioning on Y x adds blocking to the paths, not removing any existing blocks; for Y x being a non-collider, any path involving Y x is blocked by B j (Y x ) \ Y x , thus remaining blocked even if Y x is not conditioned on. Therefore, Y x does not influence the blocking status of any relevant paths. Consequently, given B j (Y x ) \ Y x , all the paths between any other nodes in Y x (including Y x ) and U j are blocked. This implies that B j (Y x ) \ Y x is also a Markov blanket. However, it is a proper subset of B j (Y x ), which contradicts B j (Y x ) being a Markov boundary. Therefore,Y x ∈ B j (Y x ) implies Y x ̸⊥ ⊥ G a x U j | Y x \ {Y x }. (ii) Y x ̸⊥ ⊥ G a x U j | Y x \ {Y x } implies Y x ∈ B j (Y x ):It can be proven that the contrapositive holds. If Y x is not in the Markov boundary, then it is independent of U j given the Markov boundary Y x ∈ B j (Y x ). Clearly, this conclusion follows directly from the definition of the Markov blanket. Therefore, Y x ̸⊥ ⊥ G ax U j | Y x \ {Y x } holds trivially. Our final algorithm for efficiently computing Markov boundaries is obtained directly from dseparation, Thms. 2 and 3 as shown in Fig. 7 . However, a limitation of these theorems is the reliance on the faithfulness assumption, which may not always hold, even in the case of a fully specified SCM, where numerical violations can potentially occur.To estimate Eq. 2, a straightforward approach is to start from the definition and employ rejection sampling, as implemented in [108] :EQUATION)This method operates effectively within a countable, finite exogenous space. However, as discussed in the main text, it exhibits low efficiency in general settings due to the infrequency with which the indicator function in the counterfactual probability expression is effective (i.e., non-zero). One solution to this problem is importance sampling, where selection of the proposal distribution is a critical issue. In this section, we discuss related works in detail.Cross-entropy based importance sampling It is straightforward to prove that when the density of the proposal distribution q(u) ∝ p(u)|f (u)| (where, in the context of counterfactual estimation, f (u)is 1 Ω U (Y * ) (u)), the corresponding importance sampling estimator possesses the minimum variance. This proposal distribution is known as the optimal proposal distribution, denoted as Q * U . The idea behind a class of methods is to minimize the cross-entropy (or equivalently, the KL divergence) between the conditional proposal distribution Q U and the optimal proposal distribution Q * U . The expression for this in the context of counterfactual estimation is:EQUATIONIn the task of sampling rare events, due to the difficulty of directly sampling from P u to obtain effective samples for which f (u) is valid, it is common to estimate the cross-entropy by sampling from a proposal distribution Q u . This method, known as importance sampling, is considered adaptive because it continuously adapts based on the samples already drawn from the proposal distribution. For instance, [32] modeled Q U as a Gaussian Mixture Model (GMM) and optimized an equivalent form of cross-entropy to improve sampling efficiency. When applied to counterfactual estimation, the corresponding optimization objective is:EQUATION)Neural importance sampling In addition to GMM, normalizing flows provide a more expressive capability as density estimators. Therefore, as enumerated in the main text, numerous recent methods have leveraged normalizing flows as proposal distributions for importance sampling. Many of these methods still rely on optimizing cross-entropy (and KL divergence), as introduced earlier. Another optimization method discussed in [68] is based on the χ 2 divergence between the optimal proposal distribution Q * U and the proposal distribution Q U :EQUATIONEQUATIONIt can be seen that the only difference between Eq. 104 and Eq. 5 lies in a constant term. Thus, directly optimizing this equation is equivalent to optimizing the variance itself. Naturally, considering the difficulty of sampling from P U as discussed in the paper, an adaptive approach is employed. This involves sampling from the proposal distribution Q U , with the corresponding optimization objective being:EQUATIONIn [68] , the proposal distribution Q U is modeled as NICE [24] .Extension and limitation The above method, without extensions, is only applicable to estimate a single counterfactual probability P (Y * ). We extend these methods similarly to our proposed approach by introducing a conditional proposal distribution Q U|Y * , and employing the same estimator as in Eq. 9 during inference as baselines. During the variance optimization phase, by introducing prior distributions Q S and P U and letting y * = Y (s) * (u ′ ) be sampled from stochastic counterfactual processes (u ′ ∼ P U and s ∼ Q S ), we reformulate the optimization expressions for cross-entropy based importance sampling (CEIS) as:EQUATIONand neural importance sampling (NIS) as:EQUATIONwhere δ(y * ) denotes a small region around y * . The above optimization objective tasks Q U|y * with sampling from Ω U (δ(y * )), where any Y * can be decomposed into a set of y An unavoidable drawback of the extended method is the presence of indicator functions in the optimization term. If the initial proposal distribution Q u is similar to P u in performance, the entire optimization process requires a prolonged warm-up period before effective samples start to emerge from the proposal distribution Q u . For example, if the entire space is the unit plane [0, 1] 2 and Ω U (δ(y * )) is a small region {(x, y) | x ≤ 0.001, y ≤ 0.001}, then if the chosen initial proposal distribution is uniform, it has only a probability 10 -6 of sampling an effective sample, resulting in a low efficiency for optimizing the loss functions.In contrast, the optimization term in Exogenous Matching according to Thm. 1 does not include the indicator function, ensuring that every sample is effective during the optimization process, thus significantly enhancing learning efficiency.In this section, we provide a detailed description of the construction of the proxy SCMs used in the experiments and their identifiability.The Neural Causal Model (NCM, [107, 108] ) is a recently proposed general approach to mimicking the original SCMs to construct proxy SCMs. For identifiability, its primary assumption is the causal graph of the original SCM. By replacing the causal mechanisms in the original SCM with neural networks that have a strong expressive power, it can effectively proxy the original SCM.Specifically, an NCM is a 4-tuple ⟨U, V, F, P U ⟩, which is almost identical to the SCM definition, except that each f Vi ∈ F is a neural network with sufficient expressive power, and P U is any well-defined distribution over U.Their method for inducing submodels is the same as for a regular SCM, i.e., replacing causal mechanisms under perfect interventions, which they term neural mutilation. For counterfactual estimation, they use the rejection sampling introduced earlier because the endogenous variables in their experimental scenarios are discrete and finite, making rejection sampling efficient enough to sample effective samples.Regarding identifiability, this work assumes causal graphs and has been shown to operate effectively when endogenous variables are discrete and finite. Specifically, they construct a class of models called G-constrained NCMs, which have the same causal graph G as the original SCM. Then, they prove that any NCM in the G-constrained NCM family has dual identifiability with the causal graph G, i.e., the NCMs in the G-constrained NCM family (given some observational or interventional distributions) output the same result when answering a Q ∈ L 3 , if and only if the graph algorithms for counterfactual identification output identifiable results given input Q and observational or interventional distributions.Furthermore, based on this duality, they develop a sound and complete identification algorithm, and since this identification algorithm is based on the learning process of NCMs, it only requires gradient optimization. Certainly, although this work guarantees theoretical results for identifiability in discrete and finite settings [108] as well as partial identifiability in Lipschitz continuous settings [93] , further improvements are necessary to apply its estimation methods in broader scenarios.Causal Normalizing Flow Causal Normalizing Flow (CausalNF, [42] ) is a recently summarized method to construct proxies of SCMs based on flow models. For identifiability, it assumes the causal graph or causal ordering of the true SCM, and the true SCM is Markovian and diffeomorphic. Specifically, a CausalNF is a 4-tuple ⟨U, V, T -1 , P U ⟩, where T = {f Vi (u i ; pa Vi ) | V i ∈ V} is a sequence-preserving Autoregressive Normalizing Flow (ANF), and P U is the base distribution of this flow model.Their approach to derive submodels differs from conventional SCMs because the modeled CausalNFs have multiple architectures, and substituting causal mechanisms is only applicable in the recursive case. Their solution is to alter exogenous distributions and leverage the Markov assumption to ensure that the density of the intervened exogenous distribution is:EQUATION)δ denotes the Dirac distribution, where a density only exists when the condition inside its parentheses is satisfied. In practice, the distribution of the intervention is determined by the following steps: i) Sample u fromP U ; ii) Compute v 1 = T -1 (u); iii) Replace v 1 , resulting in v 2 = {v i | v i ∈ v 1 , v i ̸ ∈ x} ∪ x; iv) Compute u 2 = T (v 2 ); v) Replace u 2 , yielding u 3 = {u i | u i ∈ u 2 , v i ̸ ∈ x} ∪ {u i | u i ∈ u 1 , v i ̸ ∈ x}; vi) Compute v x = T -1 (u 3 ).Here, step (v) corresponds to replacing relevant parts of the density of the intervened exogenous distribution with the Dirac distribution, as shown in Eq. 108. The last four steps among these six steps are employed for counterfactual reasoning, but this form of counterfactual reasoning is not for generalized cases. Specifically, the counterfactual query corresponding to these four steps isP (V x | V 1 ). Although the submodel cannot be directly obtained, the v x obtained from the above steps is indeed the potential response V x (u), so our method is still applicable to this model.As for identifiability, this work requires numerous assumptions, including Markovianity, diffeomorphism, causal order, and uses results from [106] to assist in proving the identifiability of the learned model up to invertible transformations.There is a challenge of modeling an exponentially large number of conditional distributions with a single model. Several existing works have focused on this problem, which is termed arbitrary conditioning. Arbitrary conditioning is an unsupervised learning task that, given a dataset over V, aims to answer any query of the form P (X | Y = y), where Y ⊆ V and y ∈ Ω Y are arbitrary.The initial step to address the problem of arbitrary conditioning is to encode the condition y into a fixed-length vector. Subsequently, this vector is learned through various modeling methods, and the training process incorporates various random masks to simulate different forms of Y. Specifically, related works introduce a binary mask b to indicate whether a variable has been observed. In this context, the conditional probability can be expressed asP (V 1-b | v b , b).Consequently, the learning objective can equivalently be viewed as maximizing the following log-likelihood term:EQUATIONwhere P b represents the prior distribution of the random mask, while P V denotes the observational distribution. The model q corresponds to the probabilistic model employed for this task. Specific methods employed in related works are enumerated in the main text, involving a variety of generative models that are suited to different task scenarios. For instance, some methods permit direct training with partially observed data, some can only output samples without density estimation, and others are incapable of performing efficient marginal probability computation.The similarity between this work and these methods lies in the encoding scheme and the training objectives. For example, the conditional encoding method in this paper also employs binary masks, introducing both the observation mask ω(X) and the intervention mask ω(Y). Furthermore, the training objective, as defined in Eq. 8, closely resembles Eq. 109.The distinction between our work and these methods lies in the available data and the target tasks. The problem of arbitrary conditioning is based on learning from the observed distribution and is typically applied to imputation tasks. In contrast, our conditional distribution modeling can be seen as an extension of arbitrary conditioning, generalizing the observed distribution to counterfactual distributions. This allows the conditions to originate from multiple hypothetical worlds and is applied in counterfactual reasoning.The work most similar to ours in terms of overall structure and optimization objectives is Posterior Matching [91] . It addresses the problem of arbitrary conditioning by matching the encoder of a given VAE. The optimization formulation is as follows:EQUATIONwhere P U|v represents the posterior distribution in a VAE, also referred to as the encoder. By comparing the optimization objective of our method, Eq. 8, with Eq. 110, we observe a notable similarity between the two. We can establish a one-to-one correspondence between the terms in the optimization formulations: the mask b can be analogized to the state s in the stochastic counterfactual process; the decoder of the VAE (mapping from latent variables to observable variables) corresponds to the SCM's potential responses (mapping from exogenous to endogenous variables); the encoder of the VAE (mapping from observable variables to latent variables) corresponds to the abduction process in the SCM (mapping from endogenous to exogenous variables). Consequently, the posterior distribution aligns with the abduction process in the SCM, represented by the product of the indicator function and the exogenous distribution.An intuitive explanation is that posterior matching attempts to align the conditional distribution with the posterior distribution, whereas our work can be interpreted as attempting to align the conditional distribution with the optimal proposal distribution. The interpretation of the optimal proposal distribution is detailed in App. B.1, where in the counterfactual estimation task, the optimal proposal distribution is proportional to the product of the indicator function and the exogenous distribution.This work involves injecting prior structure into the learning process. Specifically, such methods of injecting prior structure are achieved by directly influencing the internal connectivity of the model. [31] , Zuko [81] , and StrNN [13] . MADE interprets from the perspective of conditional probability, altering the connectivity between inputs and outputs in autoregressive models by masking, thereby enabling the model to output specific conditional distributions p(x i | y) only when there exists connectivity from inputs in y to output x i ; however, this method only respects autoregressive properties. The method in Zuko is primarily used for implementing masked autoregressive normalizing flows, but is not limited to autoregressive dependencies. Based on the interpretation of Jacobian matrices, it alters the connectivity between inputs and outputs in MLPs through masking, ensuring that the Jacobian matrix of MLP outputs is non-zero only when the mask is one, thus representing the dependency between input and output variables. However, the masking algorithm in Zuko only approximately preserves the maximum number of connections under dependency constraints in variable order. StrNN extensively investigates the issue of variable dependency and connectivity, formalizing it as an integer programming problem, and providing exact and greedy algorithms that outperform the previous two methods in terms of performance.The tasks typically addressed by the above methods include representing autoregression (autoregressive models), constraining the dependency of Bayesian networks (Bayesian network inference), and constraining the dependency of causal graphs (causal inference). Different from these tasks, this paper applies them to representing the dependency of Markov boundary (Markov boundary feature selection). The masking algorithm adopted in this paper is Zuko, as it can be directly combined with the normalizing flows it implements and used directly in subsequent experiments.We formally introduce all the metrics used in the subsequent experiments, which aim to reflect the quality of counterfactual estimation:Effective Sample Proportion (ESP) As mentioned earlier, effective samples here refer to samples that make the indicator function equal to one (i.e., the constraints of potential responses are satisfied).For a specific Y * , we defineEQUATIONwhich reflects the proportion of effective samples sampled in a counterfactual estimation for Y * . In rare event sampling, this metric essentially represents the success rate, indicating the proportion of samples in which the rare event occurs. ESP is about the expectation of this function over a prior distribution P Y * (related to the stochastic counterfactual process):EQUATIONSpecifically, ESP reflects the sampling efficiency when the importance weights are ignored. It focuses on the probability that the sampled u falls within the support Ω U (Y * ), for all Y * .Failure Rate (FR) This work claims that the proposed method holds the potential for counterfactual estimation in any form following a single training session. To verify the reliability of this claim, we introduce the FR metric, which reflects the efficacy of the proposal distribution post a single training across the entire state space of the relevant stochastic counterfactual processes, defined as follows:EQUATIONwhere 0 < m ≤ 1 is a constant representing the acceptance threshold. An η(Y * ) below this threshold is considered a failure. The closer the FR is to zero, the more it indicates that for the vast majority of Y * , there are sufficiently many samples of u falling within the support Ω U (Y * ), making the counterfactual estimation considered tractable.This paper introduces a stochastic counterfactual process (Def. 1), allowing the proposal distribution after a single training to potentially estimate counterfactual probabilities in multiple or even infinite forms in a tractable manner. This facilitates efficient estimation of a L 3 expression in real-world scenarios (assuming that the L 3 expression can be approximated by a finite set of counterfactual probabilities). Therefore, if the proposal distribution after a single training can tractably estimate these counterfactual probabilities, it implies that estimation of the L 3 expression only requires one training. Specifically, employing stochastic counterfactual processes, this paper conducts stochastic optimization on various forms of counterfactuals and demonstrates its feasibility (Cor. 1). In this section, we elaborate on its specific design in experiments. Q. Q ∈ {ATE, ETT, NDE, CtfDE}. Q Counterfactual Variables (Y * ) Expression ATE Y X=1 P (Y X=1 = 1) -P (Y X=0 = 1) Y X=0 ETT {Y X=1 , X} P (Y X=1 =1,X=1)-P (Y X=0 =1,X=1) P (X=1) {Y X=0 , X} X NDE {Y X=1,W =w , W X=0 }, w ∈ Ω W w P (Y X=1,W =w = 1, W X=0 = w) -P (Y X=0 = 1) Y X=0 CtfDE {Y X=1,W =w , W X=0 }, w ∈ Ω W w P (Y X=1,W =w = 1, W X=0 = w) -P (Y X=0 = 1) P (X = 1) Y X=0 X Construction of Y B * If Q ∈ L 3 is an arbitrary P (Y * ) such that Y * ⊆ Ω Y * isEQUATIONEQUATIONEQUATIONEQUATIONEQUATIONEQUATIONEQUATIONEQUATIONwhere Bernoulli(ρ 1 ) is a Bernoulli distribution with parameter ρ 1 , meaning it has a probability of ρ 1 to yield 1; P X is the endogenous distribution corresponding to X and for any x ∈ Ω X , p(x) > 0.This distribution, due to the constrained Bernoulli distribution parameters (0 < ρ 1 , ρ 2 < 1), implies that any X, Y ⊆ V may potentially be sampled; and for intervention values, assuming positivity (p(x (i) ) > 0), any possible intervention value can be sampled. Therefore, for any of these submodels, any counterfactual variable is covered, and hence any variable set is also covered.In subsequent experiments, we set ρ 1 = 0.2 and ρ 2 = 0.75, which are utilized during both the training and testing phases. The rationale behind employing Y B * for experiments lies in the large size of S, which provides more credible evidence to verify whether the conclusion of Cor. 1 satisfies the claim and whether the designs bridge the gap between theory and practice.To apply our method to real-world tasks, we construct stochastic counterfactual processes Y Q * for several common causal query tasks, where Q ∈ {ATE, ETT, NDE, CtfDE}. As shown in Tab. 3, the state spaces for these queries are finite or easily sampled, allowing the design of the prior distribution Q Q S to be a simple uniform distribution over these available states. We assume that the endogenous variable space is discrete and finite, which means that we only need to estimate the counterfactual probability for each of them. By computing according to their expressions, we can obtain the estimate for Q.In this section, we briefly introduce all fully specified SCMs used in our experiments.Markovian diffeomorphic SCMs These SCMs are sourced from [42] and are primarily used in the experiments for CausalNF, including: CHAIN-LIN-3, CHAIN-NLIN-3, CHAIN-LIN-4, CHAIN-LIN-5, COLLIDER-LIN, FORK-LIN, FORK-NLIN, LARGEBD-NLIN, SIMPSON-NLIN, SIMPSON-SYMPROD, TRIANGLE-LIN, and TRIANGLE-NLIN. Markovianity denotes the absence of hidden confounders in these SCMs, meaning that each exogenous variable is involved in the causal mechanisms of only one endogenous variable. The diffeomorphism is a parametric assumption, positing that the transformations between exogenous and endogenous variables satisfy the properties of a diffeomorphism. For detailed construction of these SCMs, see the Appendix of [42] .Semi-Markovian continuous SCMs These SCMs are described in [107] . We have reformulated their functional forms to adapt them to continuous scenarios. Specifically, this involves replacing their causal mechanisms with combinations of linear and nonlinear functions, all of which originate from the Markovian diffeomorphic SCMs already present above. These SCMs include: BACK-DOOR, FRONT-DOOR, M, and NAPKIN. Here, Semi-Markovianity indicates that there may be latent confounders within these SCMs, yet they still satisfy recursiveness. Continuity refers to the assumptions regarding their parameters.Regional canonical SCMs These SCMs are originated from [108] , where all endogenous variables are discrete and finite (binary), with complexity arises from the causal mechanisms being a type of stochastic mapping. These SCMs include: FAIRNESS, FAIRNESS-XW, FAIRNESS-XY, FAIRNESS-YW, where subscripts indicate the presence of hidden confounders between variables. Please refer to the Appendix of [108] for the specific constructions of these SCMs.In this section, we present details pertinent to reproducibility, encompassing the complete architecture of the models utilized, the computational resources employed, and the execution and random seeds of the reported experiments.Architecture details The entire model architecture involves three components: the sampler, the conditioning model, and the density estimation model.The sampler is used to sample s from the prior distribution Q S of the stochastic counterfactual process and u from the exogenous distribution P U , then compute the potential outcome Y (s) * (u).In our experiments, we selected a batch size of 256 and generated a training set of size 16,384 through the sampler. The validation set is also generated using this sampler. Specifically, we utilize the same sampler to construct P Y * , from which random y * is sampled. For discrete distributions, counterfactual events are assigned as Y * ∈ {y * }, and for continuous distributions, counterfactual events are assigned as Y * ∈ δ l (y * ), where δ l (y * ) is a cube centered at y * with side length l.A brief overview of the conditioning model is shown in Fig. 2 , where both g and h are neural networks. When masking with Markov boundaries, h is a neural network that allows dynamic influence on internal connectivity, as depicted in App. B.4. On the other hand, g is a permutation-invariant neural network used to aggregate hidden encodings from multiple different submodels. We conduct ablation studies related to the construction of these two neural networks in App. C.8.For the density estimation model, we chose normalizing flows and GMM. Specific models for normalizing flows include MAF [73] , NSF [26] , NCSF [80] , NICE [24] , NAF [39] , and SOSPF [41] . The implementations of these models are based on the Zuko library [81] , with some internal implementation details modified to fit our work. Regarding the hyperparameters of these models, the number of components in GMM is fixed at 10, the number of transformations for all flow models is set to 5, and all the neural networks involved contain 2 hidden layers with 64 neurons each (in some experiments, 256 neurons were used, which will be specifically indicated and discussed in subsequent sections).Hardware Each experiment was conducted on GPUs, and the primary computational overhead arising from the sampling process. All models, including the pre-trained SCM agents, were trained and tested on an NVIDIA RTX 4090 and Intel(R) Xeon(R) Gold 6430 platform. The training time for EXOM is primarily influenced by the number of variables in the SCM, the number of submodels, and the chosen density estimation model. For example, in the case of SIMPSON-NLIN with |s| = 5, one epoch (16, 384 samples, batch size of 256) including the time for the sampler to prepare data takes approximately 8 seconds. In all EXOM experimental setups, we fixed the training to 200 epochs. All experiments, including the unreported ones, consumed approximately 1000 GPU hours in total.Execution To reflect the error, each type of experiment was conducted five times with different random seeds, and the range of errors exhibited is reflected in the experimental results. During training, EXOM does not require sampling and is trained for a fixed 200 epochs. The sampling size for CEIS is set to 10 5 and trained for 200 epochs, whereas NIS has a sampling size of 10 4 and is trained for 50 epochs (we found that its training speed is excessively slow due to the time-consuming nature of sampling from the flow model). During validation and testing, the validation set queries 1024 counterfactual probabilities from the sampler, with sampling sizes for EXOM, CEIS, and NIS set to 10 3 , and RS set to 10 6 . For all EXOM, CEIS, and NIS training, we used AdamW [62] with an initial learning rate of 0.001 as the optimizer and ReduceOnPlateau with a patience of 5 and a factor of 0.5 as the learning rate scheduler. To improve training robustness, the conditional inputs y * and both the input and the output exogenous samples u were standardized. During inference, the probability densities and importance weights were processed in logarithmic form, with the terms where the indicator function is zero replaced byinf, and the log-sum-exp technique was used to estimate the probabilities.Here, we present additional results on convergence experiments. Our experiments span across combinations of 4 different types of SCMs (SIMPSON-NLIN, LARGEBD-NLIN, FAIRNESS-XW, NAPKIN) and 4 different density estimation models (GMM, MAF, NICE, SOSPF). All models were trained and validated based on the stochastic counterfactual process Y B * , which provides a sufficiently large state space to increase difficulty, where |s| = 3.As shown in Fig. 8 , for all combinations, as LL (log-likelihood, which is the negative form of the training objective Eq. 8, minimizing Eq. 8 is equivalent to maximizing LL) increases, ESP improves. In particular, there is a clear inflection point in the experiments with three continuous SCMs, after which the ESP rapidly improves until it stops increasing because of the convergence of LL. This series of experiments indicates that we can indeed optimize the upper bound (indirectly equivalent to LL) to optimize the variance (as discussed in App. C.1, ESP and variance are strongly related). Thus, empirically, these experiments support the correctness of Thm. 1, and also empirically indicate that our design potentially satisfies the bounded importance weight condition discussed in App. A.2. Cor. 1, and manifests the claim of simultaneously optimizing the variance of counterfactual estimators across the entire state space.We provide more comprehensive and detailed information regarding the comparison experiment in the main text. These comparisons are conducted on three different types of SCMs (SIMPSON-NLIN, NAPKIN, FAIRNESS-XW). EXOM employs two density estimation models (MAF, GMM), while NIS utilizes NICE, and CEIS uses GMM. All models are trained and validated based on a stochastic counterfactual process Y B * , where |s| = 1, 3, 5. As illustrated in Tab. 1, overall, our proposed method, EXOM, outperforms other importance sampling methods, while CEIS and NIS show comparable performance in the aforementioned experiments. The advantage of EXOM is particularly pronounced when dealing with multiple distinct submodels, significantly surpassing the other methods. All three importance sampling methods outperform RS, demonstrating that variance optimization indeed enhances estimation efficiency.In this section, we detail the performance of EXOM when using various density estimation models (GMM, MAF, NSF, NCSF, NICE, NAF and SOSPF) in different SCMs enumerated in App. C.3. These experiments are based on Y B * , with |s| = 3. In Figs. 10a and 11a , deeper colors indicate a higher number of effective samples, which implies greater sampling efficiency. For density estimation models, it is evident that MAF and NSF generally perform well, while other density estimation models are suitable for specific scenarios. Regarding SCM, we observe that FORK-LIN and TRIANGLE-NLIN perform worse than other SCMs, whereas some SCMs such as CHAIN-LIN-5, LARGEBD-NLIN, and NAPKIN show compatibility only with certain density estimation models. The performance improved significantly after increasing the hidden layer width from 64 to 256.In Figs. 10b and 11b , lighter colors indicate fewer failure cases, meaning that the method works effectively within the state space Y Q * . It can be noted that for any of the aforementioned SCMs, there is always a density estimation model that matches and keeps the FR at a low level. Furthermore, compared with Figs. 10a and 11a , there is a strong negative correlation between FR and ESP, consistent with the theoretical results. In regions where FR is particularly high, in addition to the poor performance of GMM due to its limited expressive capacity, there are some notable aspects. For example, FORK-NLIN and TRIANGLE-NLIN are the most challenging settings to learn. Furthermore, as the diameter of the SCM increases, the performance of EXOM gradually decreases, which is reflected in the performance under the CHAIN-LIN-N series settings. Finally, comparing Figs. 10b and 11b , some models with high parameter number requirements (such as NSF, NCSF, NAP, SOSPF) show significant improvements when the width of the hidden layers is increased. parameters, so setting the hidden layer width to 64 is sufficient. In contrast, models like NSF, NCSF, NAF, and SOSPF, which require a larger number of parameters, will need a wider hidden layer -256 is recommended -since, with a smaller Markov boundary, each parameter connection's hidden neurons are insufficient. Optimal width configuration should be fine-tuned according to the settings.Our ablation experiments were conducted on the architecture of the two neural networks g and h in the conditioning model.Let us first introduce how the mechanism of Markov boundaries operates on specific density estimation models. Firstly, we calculate the Markov boundary corresponding to each U i according to Theorem 3, and generate masks m i using the method described later. The subsequent processing is then model-specific.1. For GMM, h primarily outputs each component's mean, diagonal matrix, and anti-diagonal matrix.Here, we apply the mask m i only to the mean corresponding to the i-th exogenous variable in each component, ensuring that if the counterfactual variable is not on the Markov boundary, the inference of the mean is independent of it. 2. For Autoregressive models (including MAF, NSF, NCSF, UNF, and SOSPF), each layer transformation follows an autoregressive sequence, where all previous layers' hidden encodings may influence the subsequent layer. Consequently, each U i 's hidden encoding appears in each layer. Therefore, we only need to directly apply the mask to the inference of hidden encodings, ensuring that the inferred next layer hidden encodings are only relevant to the counterfactual variables when they are in the Markov boundary. 3. For Coupling models (including NICE), each layer transformation involves only a subset of hidden encodings from the previous layer; hence, each layer's output only corresponds to some U i . Thus, we create masks only for the U i that serve as the output for each layer transformation, ensuring that the inferred hidden encodings are only relevant to the counterfactual variables in the Markov boundary.Of course, the above initial design will raise some issues in practice:1. There are some inevitable discrepancies between theory and practice, which we have summarized in App. D.3. 2. Another issue worth discussing is how to handle the counterfactual Markov boundaries in submodels. In perfect interventions, all connections between intervened variables and their parent variables are cut (All Cut), which is also a way to represent causal graphs in submodels. However, theoretically, directly cutting connections would result in intervened variables never participating in the conditional reasoning of their parent exogenous variables. Hence, we have also devised two relaxed and compromised methods: one is to retain all connections (No Cut), pretending that the submodel has not been intervened upon, thus not reflecting the intervention through the mask of the Markov boundary but instead relying on the intervention information included in the conditional encoding; the second is to only retain connections with parent exogenous variables (Endo. Cut), while relationships between endogenous variables are still cut. 3. As discussed in App. C.7, the introduction of the Markov boundary theoretically increases the model's focus on specific variables. However, this comes at the cost of expressiveness, especially in models with a larger number of parameters. Therefore, in the ablation experiments, we use a hidden layer width of 256 for NICE and SOSPF, while for GMM and MAF, the hidden layer width is set to 64. This does not affect the fairness of the ablation, as the goal is to compare the performance difference with and without the Markov boundary. Under the same hidden layer width, the presence of a Markov boundary actually means fewer neurons are involved in reasoning. without a clear outperforming approach. However, we find that when the All Cut method performs weaker than that without Markov boundary masks, the latter two serve as compromises and indeed provide better results than All Cut. Given these observations, Endo. Cut was chosen in our previous experiments.The next problem to address is the choice of g. We have posited the use of permutation-invariant functions in this paper because of the permutation invariance of sets. Empirical evidence in our experiments supports this assertion.We compare 4 different choices of g: i) Concatenation: directly concatenate different encodings into a fixed-length vector, with zero-padding in empty positions. This is not a permutation-invariant function and is therefore included as a limiting case for comparison. ii) Summation: the simplest permutation-invariant method, which directly sums all encodings. iii) Weighted Summation: an extension of Summation that provides additional expressive power. Specifically, encodings are first mapped through a network w to obtain element-wise weights (via softmax), which are then multiplied by the encodings themselves to yield the result. iv) Attention: offering stronger expressive power, Attention differs from Weighted Summation in that it uses encodings before input h rather than after output h, and a new network a maps them to attention weights. Attention encoding is calculated through queries and keys softmax(h(x)a ⊤ (g(x))), where we set the values in the attention mechanism to be constantly 1.We no longer explore aggregation methods with stronger expressive power, as we found in our experiments that they are likely to converge prematurely to local optima. This causes the variance in the state space to increase as the optimization process progresses, which in turn prevents Eq. 8 from being uniformly applied across the entire state space.Fig. 14 and Fig. 15 respectively illustrate the impact of g selection on ESP and FR under different SCM and density estimation models. Overall, permutation-invariant methods (Summation, Weighted Summation, Attention) consistently outperform permutation-variant methods (Concatenation), aligning with theoretical expectations and intuition. Among all methods, Attention consistently outperforms others, thus serving as the choice for all preceding experiments. We combine identifiable proxy SCM with Exogenous Matching for counterfactual estimation. Specifically, we employ two distinct proxy SCM methods tailored for these scenarios to conduct counterfactual estimation on two types of SCMs.We conduct experiments on 3 different Markovian diffeomorphic SCMs. According to the previous arguments, CausalNF belongs to M IM(P V ,O) , indicating that they are L 3 -identifiable. In our experiments, we mainly measure the error in counterfactual probability queries P (δ l (y * )), comparing the SCM represented by CausalNF to the original SCM in 1024 different queries. These results are reflected in Tab. 4, with averages presented. For the original SCM, we train EXOM five times and calculate the error bound for the same query. For CausalNF, we first train 5 proxy SCMs using observational data, then train EXOM 5 times for each model, calculate the error bound for the same queries, and finally take the average error bound across proxy SCMs.Since in the continuous case it is challenging for RS to sample instances that satisfy the counterfactual event δ(y * ), obtaining ground truth in experiments becomes difficult. Therefore, we primarily use FR as the main evaluation metric, while the error bound is employed to reflect the bias. When measuring error bounds, it is imperative to regularize the dimensionality of the counterfactual variables. Failure to do so may result in significant discrepancies in the final outcomes (estimates in high-dimensional scenarios tend to be considerably lower than those in low-dimensional ones, leading to a severe underestimation of errors in high-dimensional cases). Therefore, we calculate the error bounds as follows: The results in Tab. 4 indicate that EXOM exhibits a relatively small error, comparable to that of RS. Furthermore, when EXOM is applied to the proxy SCM, it yields an even smaller error (in part because all failure cases are excluded from consideration).EQUATIONCounterfactual estimation with NCM We conducted experiments on 4 Regional canonical SCMs. According to their proof, NCM exhibits duality with respect to the identifiability of the counterfactual query Q and the graphical causal identifiability. In our experiments, we focus primarily on the unbiasedness of EXOM combined with NCM. The results are shown in Tabs 5 and 6.As RS is straightforward to sample in the discrete case, the ground truth can be computed. We used RS with 10 6 samples as ground truth (in the experiment, we compared against RS with 10 3 samples) and then computed the average bias. Additionally, since these queries only involve a single dimension (regarding solely to Y ), there is no need for dimensional regularization of the error bound.In estimating these causal queries, EXOM yielded relatively small errors and showed even smaller errors when applied to proxy SCMs. This indicates that integration of EXOM with proxy SCMs is practical to address real-world problems. In particular, the errors in ETT and CtfDE within the original SCM were substantially higher than those in RS, primarily because their expressions include the denominator P (X), and inaccuracies in P (X) markedly influence their results. This matter will be summarized in App. D.3.Our work aims to enhance the efficiency of counterfactual reasoning and make it applicable to theoretically broader contexts. On the positive side, current counterfactual reasoning has a positive impact on society in terms of counterfactual fairness and counterfactual decision making. However, on the flip side, if assumptions and premises are not applicable to real-world situations, counterfactual reasoning can still lead to potential biases, inequalities, and misuse. Therefore, careful consideration of the correctness of the assumptions and the appropriateness of usage scenarios is required before applying counterfactual reasoning.For continuous distributions, it is the density function rather than the probability measure that garners widespread attention and study within the relevant field. According to the definition of probability density, we can find its connection with the Lebesgue differentiation theorem, which directly provides a method for density estimation: This method is analogous to kernel density estimation (with the Parzen window as a cube). However, the precision of the estimates provided by this method is highly dependent on |δ(y * )|, further highlighting the difficulty of satisfying the indicator function in the counterfactual probability measure. We indirectly investigated the performance of EXOM on this issue in our related experiments on continuous counterfactual estimation, where we chose Y * ∈ δ l (y * ) as the counterfactual event for the estimation. The results of the related experiments confirmed the potential of EXOM in addressing this problem.Of course, kernel-based density estimation has many limitations. We hope that Thm. 4 will serve as a lemma for future work, leading to the derivation of more suitable density estimation methods for continuous counterfactual distributions.1. Dependency on counterfactual generation. In practical tasks, we typically only have access to observational or experimental data. In some settings, we also assume the acquisition of causal graphs or causal orderings through causal discovery and prior knowledge. However, the proposed method is not directly applicable to such real-world settings, as indicated by the assumption in Sec. 3, which requires the ability for counterfactual generation. In practical scenarios, one approach is to pre-train a neural proxy SCM to perform the counterfactual generation task. 2. Boundedness condition on importance weights. In the context of learning and optimization, we introduce additional designs: the bounded importance weight condition. We assumed conditional distribution model to satisfy the boundedness of importance weight as much as possible. However, in practical applications, this depends on the specific forms of the exogenous distribution and the conditional distribution model, necessitating case-by-case analysis. 3. Faithfulness. For injection of the Markov boundary, we also introduce the faithfulness assumption, which assists us in proving theorems and providing fast algorithms. However, this assumption might be violated, leading to more variables being included in the actual Markov blanket than our calculations predict. Additionally, modifications to autoregressive and coupling models do not fully align with the ideal effect. In fact, due to potential interactions within the hidden encodings of autoregressive and coupling models, even if each layer adheres to the Markov boundary dependencies, the aggregation of multiple layers may violate these dependencies. 4. Recursiveness. Another potential assumption in this work is the recursive SCM, which might be violated in more general settings, such as systems with cycles. All the theoretical foundations of this work are based on recursive SCMs, thus requiring these theories to be extended to more general cases. 5. Variance of L 3 expressions. In estimating counterfactual expressions Q ∈ L 3 , although the final result is unbiased, computing each term P (Y * ) individually introduces significant variance, as demonstrated in the experiments on counterfactual query estimation with real SCMs (see Tabs 5 and 6). Therefore, future work should focus on reducing the estimation variance of specific queries. 6. Potential sensitivity to the specific forms of SCM in practice. With respect to the parameter form of the SCM, although we impose no theoretical constraints on the form of SCM, experiments indicate that in some settings the form of SCM still significantly influences outcomes. For example, combinatorial experiments and ablation studies reveal unexpected scenarios, particularly with the combinations of long-chain SCMs (CHAIN-LIN-4, CHAIN-LIN-5, LARGEBD-NLIN). As the causal chain lengthens, the performance of these density estimation models gradually weakens, suggesting that long causal chains complicate counterfactual inference. This phenomenon indicates that the form of SCM in practice affects the specific effectiveness of our methods, highlighting the necessity of selecting appropriate density estimation models and making corresponding adjustments.including • The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 


Decision: This paper aims to estimate counterfactual probabilities in Structural Causal Models, and the authors proposed EXOM, an importance sampling approach to solve this task. All reviewers are inclined towards a consensus of accepting the paper. During the rebuttal phase, the authors have added supporting arguments and results to justify the proposed method.

This is a particularly long paper for a conference: the paper ends on page 54. I encourage the authors to submit a cleaner, more concise version of the camera-ready version. In fact, the reviewers all raise concerns about the paper's clarity and comment that the paper is not easy to follow. The authors should also include the discussion about the assumptions of Theorem 1 in the camera-ready version.

Reviews: review: paper_summary: paper_summary: Based on the importance sampling methods, the authors propose an exogenous matching approach to estimate counterfactual probability in general settings. They derive the variance upper bound of counterfactual estimators and transform it into the conditional learning problem. They also employ the Markov boundaries information in the inference to improve the learning performances further. Extensive experiments validate the superiority and practicality of their method.
strength_weakness: strength_weakness: - This paper makes it clear in lines 140-144 about the assumptions needed for the proposed method. Regarding assumption ii), I think it is not mild, and I am wondering if the proposed method would be sensitive to the specified distribution $P_{\textbf{U}}$ of $\textbf{U}$. The authors might have performed such experiments, but it is not quite clear. 

- Are the Markov boundaries learned from the observational data via the d-separation, or they are given prior? The authors claimed that such Markov boundaries are structural prior knowledge in line 223, whereas they gave Theorem 3 to demonstrate how to obtain them. 

- It is suggested to offer the whole procedure or pseudo code of their proposed algorithm somewhere.

- In line 239, “augmentied” might be a typo error.
questions: questions: Please see my questions in Weaknesses.
limitations: limitations: Not applicable.
paper_summary: paper_summary: This paper presents Exogenous Matching (EXOM), a new importance sampling method for estimating counterfactual probabilities in Structural Causal Models (SCMs). EXOM transforms variance minimization into a conditional distribution learning problem, providing an upper bound on counterfactual estimator variance as per Theorem 1. It outperforms existing methods across various SCM settings and integrates well with identifiable neural proxy SCMs for practical applications. By incorporating prior knowledge through Markov boundaries, EXOM further enhances performance, demonstrating its potential as an efficient tool for counterfactual estimation in diverse scenarios.
strength_weakness: strength_weakness: 1. Theorem 1 relies on the assumption that the density ratio $q(\mathbf{u}|\mathbf{y}_ {\ast})/q(\mathbf{u}|\mathbf{y}_ {\ast}^\prime)\leq \kappa$ holds for all $\mathbf{u} \in \Omega_{\mathbf{u}}$ and $y_{\ast}$, $y_ {\ast} ^\prime \in \Omega_{\mathbf{Y}_{\ast}}$. This assumption may be overly stringent, as probability measures with infinite support sets might easily violate it. Could the authors elaborate on this assumption and provide examples of distributions that satisfy it?
  
2. In the Sampling and Optimization section, the distribution of the exogenous variable $\mathbf{U}$ is assumed to be known. However, in practical scenarios, $\mathbf{U}$ is often unknown, necessitating additional efforts to estimate $\mathbf{P_U}$ [A]​. Could the authors provide further clarification on this assumption and discuss potential methods for estimating $\mathbf{P_U}$​?
  
3. I understand the authors only consider models that provide identifiability results. However, it is encouraged to include neural proxy SCM methods based on VAE and DDPM as experimental baselines. While these may lack identifiability guarantees, comparing against them would further illustrate the superiority of the proposed method in relation to current state-of-the-art techniques.
  
4. While the method shows good performance on the tested SCMs, it's unclear how well it scales to larger, more complex causal models. The experiments are conducted on relatively small SCMs, and scalability to high-dimensional or densely connected causal graphs isn't thoroughly addressed.
  

[A] Ren, Shaogang, and Xiaoning Qian. "Causal Bayesian Optimization via Exogenous Distribution Learning." *arXiv preprint arXiv:2402.02277* (2024).
questions: questions: 1. In Table 1, the EXOM method with MAP shows significantly better performance on the SIMPSON-NLIN and NAPKIN datasets compared to EXOM with GMM, whereas the performance on the FAIRNESS-XW dataset is similar for both methods. Could the authors provide further explanation for this discrepancy? Why does the GMM approach underperform in these specific cases, and what factors contribute to the similar performance on FAIRNESS-XW?
  
2. In the ablation study investigating the impact of injecting Markov boundaries, could the authors please include the performance results of EXOM without Markov boundaries? This comparison would further illustrate the benefit of incorporating Markov boundaries.
limitations: limitations: The authors discuss the limitations of their work in Section 6 and Section D.3. Notably, this work does not have a negative societal impact.
paper_summary: paper_summary: This paper introduces an importance sampling method for efficient estimation of counterfactual expressions within general settings. It transforms the variance minimization problem into a conditional distribution learning issue, allowing integration with existing modeling approaches. The paper also explores the impact of incorporating structural prior knowledge, i.e. Markov boundaries, and applies the method to identifiable proxy SCMs, proving the unbiasedness of estimates and illustrating the method's practical applicability.
strength_weakness: strength_weakness: 1. Contributions are not disentangled well. All three points involve experimental or empirical findings.
2. Some results of the ablation study are abnormal. First, the results show that the approach proposed is not robust. Under setting SIMPSON-NLIN and M, the inclusion of Markov Boundary Mask significantly improves ESP. However, under setting LARGEBD-NLIN and NAPKIN, the Markov Boundary Mask harms the performance, especially with backbone SOSPF. Second, the ESPs under setting LARGEBD-NLIN with backbone SOSPF are almost 0, even when $|s|=1$, which is hard to explain if including Markov Boundary Mask is effective. Third, the variance under setting LARGEBD-NLIN with backbone NICE is extremely large. 
3. Insufficient explanation or legend for figures, making it difficult for readers to understand. For example, in Figure 1, $\mathcal{M}$ with a subscripted hammer is not explained. In Figure 3, the legend does not indicate what different colors mean.
4. Hard to follow. Some terminologies need explanation or reference. For example, in line 234, *faithfulness* is not defined.
questions: questions: I would like to know whether Theorem 2 and 3 are novel. Have the papers cited (like 81, 3, 94, 111)  and other papers proposed methods to obtain Markov boundaries? If so, what is the improvement of the method proposed in this paper?
limitations: limitations: The authors adequately discussed the limitations in the paper.
paper_summary: paper_summary: The manuscript introduces exogenous matching, an importance sampling method for efficient estimation of counterfactual expressions in various settings. This method transforms the variance minimization problem into a conditional distribution learning problem, allowing integration with existing modeling approaches. The authors validate their theoretical findings through experiments with different Structural Causal Models (SCMs), showing competetive performance in a range of counterfactual estimation tasks. They also examine the impact of structural prior knowledge and demonstrate the method's unbiased estimates and practical applicability in identifiable proxy SCMs.

Update: revising score upward following author rebuttal.
strength_weakness: strength_weakness: The manuscript is not always clear, probably because a great deal of material has been moved to the appendix to accommodate page count. The result is a somewhat disjointed text that would likely be better served by a journal publication than a conference paper. That said, I am generally supportive of this submission and would be willing to revise my score upward if my questions are adequately addressed (see below).
questions: questions: When we say the causal model is “not fully specified”, does that just refer to the structural equations or to the graphical structure as well? In general, I was not always certain just how much causal information is used as input to this method.

I’m a bit confused by Eq. 5. If this is meant to be a variance estimator, then presumably the RHS should be something like $\mathbb{E}[X^2] – (\mathbb{E}[X])^2$, where $X$ denotes the likelihood ratio $p(u) / q(u)$, correct? This is almost but not quite what we find here. Looking at the appendix, I don’t see why Eq. 48 follows from Eq. 47. Why do we get to drop the square from the first term?

What does the constant $c$ denote in Eq. 6? Is it just the entropy of $P$? A brief word on this would help with intuition.

Should there definitely be a negative both in front of the expectation *and* within it? I assume the first summand should just be $-\mathbb{E} [ \log Q(u \mid Y_* (u)) ] $? This would look more like the classic cross entropy formula, as indicated in the following line. Eq. 7 also suggests this.

On “augmented graphs” – does this “reverse projection” of the ADMG always work? Different DAGs can have the same ADMG, for instance if two latent variables have all the same endogenous children. Perhaps there’s an unstated minimality assumption at work here?

What is $m$ in Eq. 13?

It is not clear to me from Sect. 5 what the sample size and data dimensionality are for these tasks? In general, this section appears rushed. The performance metrics are also somewhat surprising. If data is simulated, then we presumably have ground truth with respect to counterfactual probabilities. If so, then why not just compute the mean square error of the proposed estimator, perhaps as a function of sample size?
limitations: limitations: Limitations are adequately addressed.

Comments: title: Re: Author rebuttal
comment: We are glad that the rebuttal addressed the reviewer's concerns, and we appreciate the reviewer for updating the rating accordingly.
Thanks for the rebuttal by the authors addresses my concerns. I will update my rating accordingly.
We thank the reviewer for reading our rebuttal and updating the score accordingly. If there are any further questions, we would be happy to provide additional information.
Thank you for the reviewer’s attention to this weakness. In fact, without concentration inequalities or approximations,  it is challenging to find concrete examples that fully satisfy this assumption on infinite support sets. Even the density ratios between Gaussian distributions often violate boundedness asymptotically, as illustrated in [1] (Eg. 4.1). The best example we could find is from [3] (Eq. 11), which provides upper and lower bounds between the binomial distribution and the Poisson distribution (with infinite support) under certain conditions.

Of course, once concentration inequalities or approximations are employed (as in relaxed formulations 2 and 3), the weakness regarding infinite support becomes less acute, provided we prove or assume that the probability measure of $P_\mathbf{U}$ on the unbounded portion is zero or very close to zero.

----
[3] [Dümbgen, L., Samworth, R., & Wellner, J. (2021). Bounding distributional errors via density ratios.](https://arxiv.org/abs/1905.03009)
Thank you to the authors for the detailed responses. Regarding Weakness 1, in addition to the relaxed assumptions, could the authors kindly provide some concrete examples of distributions with infinite support sets that satisfy these assumptions?
Many thanks to the authors for their detailed replies to my comments. I will be revising my score upward in light of these clarifications.