Conditional Generative Models are Sufficient to Sample from Any Causal  Effect Estimand

Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 02:43:04

Abstract: Causal inference from observational data plays critical role in many applications in trustworthy machine learning.
While sound and complete algorithms exist to compute causal effects, many of them assume access to conditional likelihoods,
 which is difficult to estimate for high-dimensional (particularly image) data. Researchers have alleviated this issue by simulating causal relations with neural models. However, when we have high-dimensional variables in the causal graph along with some unobserved confounders, no existing work can effectively sample from the un/conditional interventional distributions. In this work, we show how to sample from any identifiable interventional distribution given an arbitrary causal graph through a sequence of push-forward computations of conditional generative models, such as diffusion models. Our proposed algorithm follows the recursive steps of the existing likelihood-based identification algorithms to train a set of feed-forward models, and connect them in a specific way to sample from the desired distribution. We conduct experiments on a Colored MNIST dataset having both the treatment ($X$) and the target variables ($Y$) as images and sample from $P(y|do(x))$. Our algorithm also enables us to conduct a causal analysis to evaluate spurious correlations among input features of generative models pre-trained on the CelebA dataset. Finally, we generate high-dimensional interventional samples from the MIMIC-CXR dataset involving text and image variables.

Main paper: 
 introduction: 
 Causal inference has been a central problem in many sciences and is also recently understood to be critical for developing more robust and reliable machine learning solutions (Xin et al., 2022; Zhang et al., 2020; Subbaswamy et al., 2021) . Although randomized controlled trials are known as the gold standard for estimating causal effects, modern formalisms such as the structural causal model (SCM) framework (Pearl, 2009) enable a data-driven approach to this problem: Given the qualitative causal relations, summarized by a causal graph, certain causal queries can be uniquely identified from the observational distribution. Today, we have a complete understanding of which causal queries can be uniquely identified, and which require further assumptions or experimental data, given this structure (Tian, 2002; Shpitser & Pearl, 2008; Bareinboim & Pearl, 2012) .For certain causal structures, sample-efficient ways to estimate (average) causal effect exist through, for example using propensity scores, or backdoor adjustments (Rosenbaum & Rubin, 1983; Pearl, 1993; Maathuis & Colombo, 2015) . However, the most general solutions assume that we have access to the joint observational probability distribution of data. For example, the sound and complete causal effect identification algorithm of Shpitser & Pearl (2008) writes interventional distributions as functions of the observational distribution. These functions can seemingly get arbitrarily complicated and in general there is no easy way to directly estimate them from data.This creates an important gap between causal inference and modern ML datasets: We typically observe high-dimensional variables, such as X-ray images of a patient, that need to be involved in causal effect computations. However, explicit likelihood-based models are impractical for such high-dimensional data. Instead, deep generative models have shown tremendous practical success in correctly sampling from such high-dimensional variables (Brock et al., 2018; Karras et al., 2019 ; Croitoru et al., 2023) . We are interested in the question: Can we leverage the representative power of deep conditional generative models for causal inference?As an example consider the causal graph in Figure 1 . We are interested in the distribution of variable Y after a hypothetical intervention on X. We use the do-operator of Pearl to represent the intervention, i.e., we are interested in p(Y |do(X)). The ID algorithm (Shpitser & Pearl, 2008) identifies this distribution in terms of the observational distribution as given on the right.This complicated expression can be evaluated if one has explicit access to the joint distribution P . However, it is not clear if one can sample from this probability distribution without an explicit likelihood-based model.In this paper, we show that this is possible. Specifically, we show that any identifiable interventional distribution can be sampled from by only training conditional generative models on the observational data. This enables us to leverage the state-of-the-art deep generative models, such as diffusion models, and even sample from p(Y |do(X)) when both X and Y are high dimensional. To the best of our knowledge, this is the first result that shows that conditional generative models are sufficient to sample from any identifiable interventional distribution. Our contributions are as follows:• We propose a training algorithm that learns a collection of conditional generative models using observational data and can sample from any identifiable interventional distribution. • We show that our algorithm is sound and complete, establishing that conditional generative models are sufficient to sample from any identifiable interventional query. The guarantees extend to sampling from any identifiable conditional interventional distribution. • We use diffusion models to demonstrate the performance of our algorithm on a synthetic image dataset, as well as a real-world COVIDx chest X-ray dataset. 
related work: 
 Over time, extensive literature has developed on the causal effect estimation problem. Given access to the causal graph and the probability distribution, Pearl proposed the do-calculus (Pearl, 1995) rules as a general solution for any identifiable causal effect. Shpitser & Pearl (2008) propose the popular sound and complete identification algorithm to express the causal effect of an arbitrary set of variables on other variables in terms of observational distributions or suggest a graphical criterion for non-identifiability. Bareinboim & Pearl (2012) ; Jaber et al. (2018) ; Lee et al. (2020) have proposed variants of the causal estimation problem based on different restrictions on the causal graph and the access to probability distributions. However, most of these studies have relied on explicit access to the observational probability distributions which limits their applicability to be employed in causal inference problems with high dimensional data.The use of deep neural networks for performing causal inference has been recently suggested by many researchers. Shalit et al. (2017) ; Louizos et al. (2017) proposed neural network-based approaches to estimate the causal effect. Sanchez & Tsaftaris (2022) employs energy-based generative models such as DDPMs (Ho et al., 2020) to generative high dimensional interventional samples. However, solutions of these methods do not generalize to arbitrary structures.The sampling-based approaches using deep generative models are limited in the literature. Kocaoglu et al. (2018) trains a collection of conditional generative models based on the causal graph and uses adversarial training. Pawlowski et al. (2020) employed a conditional normalizing flow-based approach to offer high dimensional interventional sampling as part of their solution. Chao et al. (2023) performs interventional sampling for arbitrary causal graphs employing diffusion-based causal mod-els with classifier-free guidance (Ho & Salimans, 2022) . However, all of these methods depend on the strong assumption that the system has no unobserved confounders. Xia et al. (2021; 2023) extends this idea to not only sample from interventional distribution but also test the identifiability of the interventional distribution. However, their approach does not handle high-dimensional variables. Perhaps more importantly, all these works train a forward model based on the causal graph structure. Conditional sampling then becomes tricky since it is not clear how to update the posterior of upstream variables using feedforward operations. For example, Kocaoglu et al. (2018) resorts to rejection sampling, which is slow and impractical for high-dimensional data.Perhaps the most conceptually related work is Jung et al. (2020) , where the authors identified an algorithm that can convert the expression returned by the ID algorithm into a form where it can be computed through a re-weighting function, similar to propensity score-based methods, to allow sample-efficient estimation. However, computing these reweighting functions from data is still highly nontrivial with high-dimensional variables in the system. 
background: 
 We first introduce the structural causal models (SCMs) and how they can model interventions. Definition 3.1 (Structural causal model, (SCM) (Pearl, 1980) ). An SCM is a tuple M = (G = (V, E), N , U, F, P (.)). V = {V 1 , V 2 , ..., V n } is a set of observed variables in the system. N is a set of independent exogenous random variables where N i ∈ N affects V i and U is a set of unobserved confounders each affecting any two observed variables. This refers to the semi-Markovian causal model. A set of deterministic functions F={f V1 , f V2 , .., f Vn } determines the value of each variable V i from other observed and unobserved variables asV i = f i (P a i , N i , U Si ), where P a i ⊂ V (parents), N i ∈ N (randomness) and U Si ⊂ U (common confounders) for some S i . P(.) is a product probability distribution over N and U and projects a joint distribution P V over the set of actions V representing their likelihood.An SCM M, induces an acyclic directed mixed graph (ADMG) G = (V, E) containing nodes for each variable V i ∈ V. For each V i = f i (P a i , N i , U Si ), P a i ⊂ V, we add an edge V j → V i ∈ E, ∀V j ∈ P a i . Thus, P a i (V i ) becomes the parent nodes in G. G has a bi-directed edge, V i ↔ V j ∈ E between V i and V j if and only if they share a latent confounder. If a path V i → . . . → V j exists, then V i is an ancestor of V j , i.e., V i = An G (V j ). An intervention do(v i ) replaces the structural function f i with V i = v i and in other structural functions where V i occurs. The distribution induced on the observed variables after such an intervention is represented as P vi (V). Graphically, it is represented by G Vi where incoming edges to V i are removed. Definition 3.2 (c-component). Given an ADMG a maximal subset of nodes where any two nodes are connected by bidirected paths is called a c-component. C(G) is the set of c-components of G.Pearl (1995) identified do-calculus rules, which relate different interventional and observational distributions to one another. Namely, rules identify conditions on the causal graph that allow i) removing/adding conditioning variables in the probability distribution (rule-1), ii) replacing dooperator with conditioning (rule-2) and iii) removing/adding do-operations (rule-3). These form the basis for any identification algorithm, and were shown to be complete to identify any identifiable interventional distribution (Shpitser & Pearl, 2008; Huang & Valtorta, 2012) . In other words, finitely many applications of these three rules is sufficient to convert an interventional distribution into a function of the observational distribution, if possible. Shpitser & Pearl (2008) ; Tian & Pearl (2002) ; Tian (2002) ; Huang & Valtorta (2012) provide systematic ways to apply them. Lemma 3.3 (c-component factorization (Tian & Pearl, 2002) ). Let M be an SCM that entails the causal graph G and P x (y) be the interventional distribution for arbitrary variables X and Y . Let C(G \ X) = {S 1 , . . . , S n }. Then we have P x (y) = v\(y∪x) P v\s1 (s 1 )P v\s2 (s 2 ) . . . P v\sn (s n ).Intuitively, P x (y) is factorized into c-factors P v\si (s i ) for each c-component S i in G \ X.Diffusion models. For the purposes of our algorithm, we aim to black-box the technical details of the generative models we learn. This is intentional: the field of conditional generative modeling is advancing rapidly and it is not at all clear that the current state-of-the-art frameworks will be used in the near future. For the purposes of our framework, we have the following simple requirement: given samples from a joint distribution, we need to be able to learn a model that provides approximate samples from any conditional distribution. In practice, we will use classifier-free diffusion guidance (Ho & Salimans, 2021) , but our framework is agnostic to the choice of generative model.Training Phase D[X,Z,Y] ~ P(V) Sample Phase MZ Z MY D[Z] D[X] MZ MY X=x Y Backdoor 2 X Z Y MZ MY Training Phase D[X,Z,Y] ~ P(X,Z,Y) Sample Phase D[X] Backdoor 1 MZ MY X=x X Z Y W Y Z Z Y Y Z Sampling graph Sampling graph (a) (b) D[Z] D[W] D[W]In this paper, we assume i) the causal model is semi-Markovian and ii) we have access to the ADMG, and iii) we assume that we can learn to sample high dimensional conditional distributions by training a classifier-free guidance diffusion model on samples from a joint distribution. 
conditional generative models for causal sampling: 
 Given a causal graph G, dataset D ∼ P (V ), our objective is to generate samples from a causal query P (Y |do(X)). Before we describe our algorithm formally, we express the challenges and outline our main ideas which lead to the general algorithm over examples.As the first example, consider the backdoor graph in Figure 2 (a). Suppose we have a dataset D ∼ P (X, Y, Z). Consider the causal query P X (Y ). An application of ID algorithm shows that P X (Y ) = z P (Y |X, z)P (z). We now intuit how to sample from this distribution. Suppose, for example, that we could sample from P X (Y, Z); in other words suppose we had a mechanism that when provided a value for x would provide samples (y, z) ∼ P x (Y, Z). From this we can derive a way to sample from P X (Y ), by sampling (y, z) ∼ P x (Y, Z) and only keeping the y variable: Dropping Z from the joint sample is equivalent to sampling Y from the distribution where Z is marginalized out. Now focusing on P X (Y, Z), ID gives P X (Y, Z) = P (Y |X, Z)P (Z). Sampling from this distribution seems more attainable; indeed if we could sample from both P (Z) and P (Y |X, Z), we could use the following procedure: first sample z ∼ P (Z) and then use this z and the specified interventional x to sample from P (Y |x, z). Hence the only ingredients we need are sample access to P (Z) and P (Y |X, Z). From the joint dataset D, we can train conditional generative models M Z and M Y to approximate sampling from these distributions. With these models in hand, they can be wired together in a sequential structure, which we visualize as a DAG, where each variable corresponds to a node and its (conditional) generative model. Sampling from this DAG can be performed by sampling from each node in a topological order, and passing sampled values to descendant nodes.Next we consider the backdoor graph in 2(b), where we again want to sample from P X (Y ). Applying the ID algorithm yields P X (Y ) = z P (z)P (Y |X, W, z), which looks identical to the previous example, except for the added term W in the conditional distribution. This raises the important question of where this W comes from. Do-calculus ensures in this case that P X (Y ) = P X,W (Y ); in other words, the causal effect of W on Y is irrelevant assuming we also intervene on X. Hence we can pick any value of W to apply here. Then intuitively, sampling can proceed in a similar fashion to the previous example. First attain a sampling mechanism M Z which samples z ∼ P (z), and then attain a sampling mechanism M Y providing sample access P (Y |X, W, Z). Sampling can then proceed by using M Z to get z, picking any arbitrary value w, and sampling y from M Y . Again this can be arranged in a DAG structure, where sampling can be done according to the topological order. P X (Y ) = w1 P (X, Y | W 1 , W 2 )P (W 1 ) w1 P (X | W 1 , W 2 )P (W 1 ).While this expression seems more complicated, it is helpful to note that, with another call to the ID algorithm, we can show that P W2 (X, Y ) := w1 P (X, Y |W 2 , W 1 )P (W 1 ). Hence, if we define P ′ (V ) := P W2 (V ), sampling from P X (Y ) is equivalent to sampling from P ′ (X, Y )/P ′ (X) = P ′ (Y |X). In other words, if we could sample a dataset D ′ ∼ P ′ (V ), which implicitly depends on W 2 , then we could sample from P X (Y ) by learning a conditional model M Y using data from P ′ (V ) that emulates samples y ∼ P ′ (Y |X). We note that since P ′ depends on W 2 , M Y should be a function of W 2 , so we can pass it as an argument to M Y . This would be correct, only if we could generate a dataset D ′ ∼ P ′ (V ) = P W2 (V ). But sampling from an interventional distribution is precisely the problem we are seeking to solve: this suggests that a recursion to smaller subproblems must play a part in our algorithm, similar to the ID algorithm. On the other hand, we notice that P W2 (X, Y ) has a similar structure to the backdoor graph, which we have already established we can sample from. Thus a sampling procedure for P X (Y ) must involve several steps: generate a dataset D ′ ∼ P W2 (V ), use D ′ to train M Y , and then sample from M Y .We can summarize the insights from these examples as follows: for an identifiable query P X (Y ) and access only to joint samples, we were able to generate a directed acyclic graph where each node corresponds to a variable and all but the interventional variables are imbued with a conditional generative model that only depends on its parents in the DAG. Sampling from the interventional query can be performed by sampling each node in this graph according to its topological ordering.The only ingredient we have required throughout this procedure is the ability to learn models to approximate samples from a conditional distribution given only joint sample access. We have been intentionally informal during these examples, but we formalize this procedure in the next section.Now we are ready to present our algorithm in full generality. The key idea is that we will be given a causal query P X (Y ), a dataset D from the joint distribution P , and a causal graph G. Our algorithm will follow the same recursive structure of the ID algorithm and return a DAG where each node corresponds to a variable and is associated with a conditional generative model that depends only on its parent nodes in that DAG. Sampling can then be performed by providing values for the intervened variables and calling each generative model according to the topological ordering of the DAG.The first step of our algorithm is an identifiability check: we run ID on our provided causal query and if ID returns FAIL, then the query is not identifiable and we return FAIL. Otherwise, we proceed.Our key contribution is the ID-DAG algorithm in Algorithm 1. We will spend the next several paragraphs outlining its operation. We note that, identical to the ID algorithm, ID-DAG proceeds through a series of 7 steps, several of which are recursive. The non-recursive steps are called base cases. In the recursive steps, we point out several key differences from the ID algorithm. Where ID passes distributions to its recursive calls, ID-DAG passes datasets, assumed to be joint samples from the corresponding distribution in the ID algorithm. Further, where ID returns arithmetic expressions of likelihood functions, ID-DAG returns sampling networks, i.e., a collection of conditional generative models that form a DAG structure. These need to be incorporated into a larger global sampling for each Vi ∈ Y in topological order do 7:Let MV i be a model trained on {Vi,V (i-1) π } ∼ D such that MV i (V (i-1) π ) ∼ P (Vi|V (i-1) π ) 8: Add node (Vi, MV i ) to H 9: Add edge Vj → Vi to H for all Vj ∈ V (i-1) π 10: Return H 11: if V \ An(Y)G ̸ = ∅ then {Step 2} 12: Return ID-DAG(Y, X ∩ An(Y)G, D ′ = D[An(Y)G], G An(Y) ) 13: Let W = (V \ X) \ An(Y)G X {Step 3} 14: if W ̸ = ∅ then {} 15: Return ID-DAG(Y, X = X ∪ W, D, G) 16: if C(G \ X) = {S1, . . . , S k } then {Step 4} 17: for each Si ∈ C(G \ X) = {S1, . . . , S k } do 18: Hi=ID-DAG(Si, X = V \ Si, D, G) 19: Return ConstructDAG({Hi} ∀i ) 20: if C(G \ X) = {S} then 21: if C(G) = {G} then {Step 5} 22: throw FAIL 23: if S ∈ C(G) then {Step 6} 24: H = ∅ 25:for each Vi ∈ X do 26:Add node (Vi, ∅) to H 27:for each Vi ∈ S do 28:Let MV i be a model trained on {Vi,V (i-1) π } ∼ D such that MV i (V (i-1) π ) ∼ P (Vi|V (i-1) π ) 29: Add node (Vi, MV i ) to H 30: Add edge Vj → Vi to H for all Vj ∈ V (i-1) π 31: Return H. 32: if (∃S ′ ) such that S ⊂ S ′ ∈ C(G) then {Step 7} 33: XZ = X \ S ′ , H S ′ = ID-DAG(S ′ , XZ , D, G), D ′ ∼ H S ′ (XZ ) 34: Return ID-DAG(Y, X, D ′ , G S ′ ,X Z )network. Hence the operation of each recursive step can be summarized according to the following scheme: i) generate the appropriate dataset from joint samples to pass to the recursion, ii) make the recursive call and acquire the returned sampling networks, iii) update the global sampling network according to the output of the recursion. As we describe the algorithm, we will briefly motivate each case and describe how the operation fits into the above scheme.Base Cases: We will describe the base cases of ID-DAG first. These are step 1 and step 6.Step 1 occurs when the current intervention set is empty. Here, the desired distribution we wish to sample from is P (Y ), and it suffices to only return a generative model for each variable in Y from the dataset D. This can be done by learning a generative model for each Y i ∈ Y conditioned on its parents according to the topological order of G and wiring them to only depend on their parents.Step 6: When we enter step 6, we are seeking to sample from P X (Y ) where Y is entirely in a single c-component S, and X is disjoint from S. In this case, ID algorithm asserts that we can replace intervening on X by conditioning on X. Similar to step 1, then, we can leverage samples from the dataset D to train a model for each variable V i in the c-component, conditioned only on ancestor nodes according to the topological order. Edges are added according to these conditioning variables.Next we can consider the recursive cases, keeping in mind the schema that we need to i) alter the dataset, ii) acquire the sampling network from the recursion, iii) incorporate it into the global sampling network. We start with steps 2 and 3 as they only make a single recursive call to ID-DAG.Step 2: We enter step 2 when there are variables X i ∈ X that are not ancestor of any Y i ∈ Y .Intervening on such variables cannot affect Y and we can safely drop them from the intervention set. We restrict the dataset D by only considering the variables that are ancestors of Y according to the topological order and stop considering the irrelevant X ′ i s. While bottom-up, we simply return the sampling network returned by the recursive call.Step 3: In step 3, ID asserts that intervening on extra variables W will have no influence on Y assuming we already intervene on X. Here the dataset is the same, and we simply augment the interventional variables. We note that the choice of W can be arbitrary, as discussed in Section 4.1.Step 4: We enter this case when there are multiple c-components in the sub-graph G \ X and one or more c-components are affected by X. Multiple recursive calls are made here, one for each ccomponent, with no alterations made to the dataset D. Each call will return a sampling network for a c-component, but special care must be taken to merge these appropriately. The straightforward parallel to the ID algorithm would be to use Tian's factorization, where P X (V ) = i P V \Si (S i ). However each returned sampling network has one node per variable considered, depending only on its parents according to the global topological ordering. Here we make a call to the subroutine ConstructDAG, outlined in the appendix, to handle the correct edges to add to merge each of these graphs together. This step then returns a valid sampling network.Step 7: This case occurs when i) G \ X is a single c-component S, ii) all the variables in Y are contained within a single c-component S ′ ⊂ S, but iii) the variables in X can be partitioned into those that are contained within S ′ and those that are not in S ′ . Letting X Z be those that are not contained within S ′ , ID asserts that evaluating P X (Y ) is equivalent to evaluating P ′ X∩S ′ (Y ) where P ′ is defined as P X Z (V ). Hence we need to generate a dataset D ′ ∼ P X Z (V ). This is handled by a recursive call to ID-DAG. Then we make the same recursive call as ID with the modification that we need to keep the variables X Z in the graph as conditioning variables. Please see appendix for more details on why this modification is needed.Sampling from ID-DAG: Finally, after ID-DAG has returned a sampling network, we can sample from it as follows: specify values for any interventional variables in X, and choose arbitrary values for any variables W that may have been added to the interventional set during step 3. Then proceed through the topological ordering of the DAG and call each generative model in turn. Of this joint sample, only return the variables in Y . We can show that this procedure is sound and complete: Theorem 4.1. (informal) ID-DAG is sound and complete for sampling from any identifiable P X (Y ).We provide the full conditional interventional sampling algorithm in Algorithm 5 and its soundness and completeness proof in Appendix C. 
experiments: 
 We apply our algorithms to two datasets. First, we consider a synthetic dataset using alterations to the thickness and color of MNIST images. Next we demonstrate that we can use Algorithm 1: ID-DAG to sample from an interventional distribution on a dataset involving real chest X-rays of COVID-19 patients. In both, we require the training of several neural networks, where full training details are deferred to the appendix. We remark that evaluation of correctness is challenging: The ground truth in these datasets is inaccessible and prior work was unable to compute interventional queries on high-dimensional data. Thus, our evaluations focus on the quality of neural network component and, in the case of MNIST, a surrogate ground-truth for a discrete version of the dataset.For conditional sampling with high-dimensional data, we train a diffusion model using classifier-free guidance (Ho & Salimans, 2022) . For conditional sampling of categorical data, we train a classifier using cross-entropy loss and apply temperature scaling as a means of calibration to ensure we are sampling in a calibrated fashion (Guo et al., 2017) . Training details and code are in the appendix.Data Generation: First we consider a synthetic dataset imbued over the napkin graph. Full data generation details are given in the appendix. We consider variables W 1 , W 2 , X, Y , where W 1 , X, Y are images derived from MNIST and W 2 is a discrete variable. We introduce latent confounders C, T , denoting color and thickness, where C can be any of {red, greed, blue, yellow, magenta, cyan}, and T can be any of {thin, regular, thick}. Data generation proceeds as follows: first we sample latent C, T from the uniform distribution. We color and reweight a random digit from MNIST to form W 1 . W 2 only keeps the digit value in {0 . . . 9} of W 1 and a restriction of its color: if the color of W 1 is red, green, or blue, W 2 's color is red, and it is green otherwise. X then picks a random MNIST image of the same digit as W 2 , is colored according to W 2 's color, and is reweighted according to the latent T . Then Y is the same original MNIST image as X, reweighted according to X's thickness and colored according to the latent C. Further, for every edge in the graph, we include a random noising process, where with p = 0.1, the information passed along is chosen uniformly randomly from the valid range. We are interested in sampling from the distribution P X (Y ).Component Models: The sampling net for the napkin graph is built by learning a model to sample from the distribution P ′ (Y |X, W 2 ), where P ′ (X, Y, W 2 ) := P W2 (X, Y ). Hence, given samples D ∼ P , we need to first generate a dataset D ′ for W 2 , X, Y from P W2 (X, Y ), or equivalently just the W 2 , X, Y samples from P (X, Y |W 1 , W 2 )P (W 1 ). Sampling from P (X, Y |W 1 , W 2 ) can be done by learning a conditional diffusion model trained on the observational distribution P . Samples from P ′ can be done by sampling W 1 ∼ P (W 1 ), choosing an arbitrary W 2 and sampling from the trained diffusion model. On the new data D ′ , we learn a diffusion model to sample from P ′ (Y |X, W 2 ).Evaluation: As the true ground truth interventional distribution is inaccessible, evaluation of our approach comes in two parts. First, we evaluate each of the trained neural nets for image fidelity and then we map our generated images to discrete variables where the ground truth is accessible.Examples from the trained diffusion models are given in Figure 4 (right). To evaluate the correctness of the sample, some preliminaries must be established. Note that each image in our data may be mapped to a discrete variable (Digit, Color, T hickness). Indeed, a discrete analogue of this dataset may be generated, for which exact likelihoods may be computed. Since we operate over images, we are only able to access these discrete properties through trained classifiers to identify digit, color, and thickness given an image. With classifiers in hand, we can estimate these properties of our sampled images from our learned P Y (X) and compare to the true (discrete) interventional and conditional distributions. We display these results for the color attribute in Table 1 , and see that our sampling much more closely emulates the interventional distribution than the true conditional.Data generation: Next we apply our algorithm to a real dataset using chest X-rays on COVID-19 patients. Specifically, we download a collection of chest X-Rays (X) where each image has binary labels for the presence/absence of COVID-19 (C), and pneumonia (N) (Wang et al., 2020) 1 . We imbue the causal structure of the backdoor graph, where C → X, X → N , and there is a latent confounder affecting both C and N but not X. This may capture patient location, which might Component Models: Applying ID-DAG to this graph requires access to two conditional distributions: P (X|C) and P (N |X, C). Since X is a high-dimensional image, we train a conditional diffusion model to approximate the former. Since N is a binary variable, we train a classifier that accepts X, C and returns a Bernoulli distribution for N . The generated sampling network operates by sampling an X given the interventional C, and then sampling an auxiliary C ′ ∼ P (C ′ ) and feeding X, C ′ to the classifier for P (N |X, C), finally sampling from this distribution.Evaluation: Again we do not have access to the ground truth. Instead, we focus on the evaluation of each component model, and we also perform an ablation on our diffusion model. We first evaluate the image quality of the diffusion model approximating P (X|C). We evaluate the FID of generated samples versus a held-out validation set of 10K X-ray images. When samples are generated with C taken from the training distribution, we attain an FID score of 16.17. We then evaluate the conditional generation by comparing class-separated FID evaluations and display these results in Table 2 (left). The classifier estimating P (N |X, C) has an accuracy of 91.9% over validation set.We note that we apply temperature scaling (Guo et al., 2017) to calibrate our classifier, where the temperature parameter is trained over a random half of the validation set. Temperature scaling does not change the accuracy, but it does vastly improve the reliability metrics; see Appendix.Finally we evaluate the query of interest P C (N ). Since we cannot evaluate the ground truth, we consider our evaluated P C (N ) versus an ablated version where we replace the diffusion sampling mechanism with P (X|C), where we randomly select an X-ray image from the held-out validation set. We also consider the query P C (N ) if there were no latent confounders in the graph, in which case, the interventional query P C (N ) is equal to P (N |C). We display the results in Table 2 (right). Return v\y P (v) 5: if V \ An(Y)G ̸ = ∅ then {Step:2} 6: Return ID(y, x ∩ An(Y)G, V\An(Y) G P, G An(Y) ) 7: Let W = (V \ X) \ An(Y)G X {Step:3} 8: if W ̸ = ∅ then 9:Return ID(y, x ∪ w, P, G) 3:10: if C(G \ X) = {S1, . . . , S k } then {Step:4} 11: Return v\(y∪x) i ID(si, v \ si, P, G 12: if C(G \ X) = {S} then 13: if C(G) = {G} then {Step:5} 14: Return FAIL(G, G ∩ S) 15: if S ∈ C(G) then {Step:6} 16: Return s\y {i|V i ∈S} P (vi|v i-1 π ) 17: if ∃S ′ s.t. S ⊂ S ′ ∈ C(G) then {Step:7} 18: Return ID(y, x ∩ S ′ , P = {i|V i ∈S ′ } P (Vi|V (i-1) π ∩ S ′ , v (i-1) π \ S ′ ), G S ′ ) Algorithm 3 ConstructDAG({Hi}if ∃α ∈ Z such that (Y ⊥ ⊥ α|X, Z \ {α})G X,α then 4: return IDC-DAG(Y, X ∪ {α}, Z \ {α}, D, G) 5: else 6: H1= ID-DAG(Y ∪ Z, X, D, G) 7: D ′ ∼ H1(X) 8: H2 = ∅ 9:Add node (X, ∅) and (Z, ∅) to H2 10:Let MY be a model trained on {Y, X, Z} ∼ D ′ such that MY (X, Z)∼ P ′ (Y |X, Z) i.e., MY (X, Z) ∼ PX (Y |Z) 11:Add node (Y, MY ) to H2 12:Add edge X → Y and Z → Y to H2 13: return H2B APPENDIX: THEORYHere we provide formal proofs for all theoretical claims made in the main paper, along with accompanying definitions and lemmas. Definition B.1. A conditional generative model for a random variable X ∈ V relative to distribution P (V ) is a function f : Pa X → X such that f (pa X ) ∼ P (X|pa X ), ∀pa ∈ Pa, where P a is a subset of observed variables in V . Definition B.2. A collection of conditional generative models for a set V of variables is said to form a sampling network if the directed graph obtained by connecting each X ∈ V to P a X via incoming edges is acyclic. This graph is called the sampling network. Assumption B.3. Each conditional generative model trained by ID-DAG correctly samples from the corresponding conditional distribution. Lemma B.4. Let H be a sampling network for random variables {V 1 , V 2 , . . . V n } formed by a collection of conditional generative models f i relative to P i for all V i . Then the tuple (V 1 , V 2 . . . V n ) obtained by sequentially evaluating each conditional generative model relative to the topological order of the sampling graph is a sample from the joint distributionΠ i P i (V i |P a i ).Proof. Without loss of generality, let (V 1 , V 2 , . . . , V N ) be a total order that is consistent with the topological ordering over the nodes in G. To attain a sample from the joint distribution, sample each V i in order. When sampling V j , we will have already sampled V i for all i < j, which is a superset of P a j by definition of topological orderings.Definition B.5. We say that a sampling network H is valid for an interventional distribution p X (Y ) if the following conditions hold:• Every variable y ∈ Y has a conditional generative model in H.• The only nodes in the sampling graph without conditional generative models are those in X.• If the values X = x are specified in H, then sampling Y from H is equivalent to sampling from P X (Y ). Lemma B.6. Termination: Let P X (Y ) be a query for causal graph G = (V, E) and D ∼ P (V ). Then the recursions induced by ID -DAG(Y, X, D, G) terminate in either step 1, 5, or 6.Proof. Since ID is sound and complete, it is guaranteed to terminate in its base cases of 1, 5, or 6. Since the steps of ID-DAG exactly mirrors the steps of ID, with respect to recursions and the arguments (Y, X, G) for every step except step 7, as long as ID-DAG does not enter step 7, it must terminate. In step 7, we make an extra call to ID-DAG of the form ID -DAG(S ′ , X Z , D, G). This follows exactly from Lemma 36 of Shpitser & Pearl (2008) : by definition, there are no bidirected arcs from X Z to S ′ , so we proceed immediately to step 6, which is a base case.Lemma B.7. If a non-identifiable query is passed to our algorithm, ID-DAG, it will return FAIL.Proof. We implicitly run an identifiability check before running the ID-DAG presented in the pseudocode. If the query is not identifiable, ID will return fail, and thus so will ID-DAG. Proof. Observe that P (Y ) may be factored asP (Y ) = Vi∈Y P (V i |Y ∩ V (i-1) π ).From Assumption 1, we may learn to sample from each conditional distribution in this product from joint samples. We can then add a node for each V i , with its associated conditional distribution as a sampling mechanism. This produces a sampling graph that is a DAG, where each variable V i ∈ Y has a sampling mechanism. Since every node has a sampling mechanism and X = ∅, the third case for correctness holds. By the factorization and Assumption 1, sampling from this graph produces samples from P (Y ), hence this sampling graph is correct for P (Y ).Lemma B.9. ID-DAG Base case (step 6): Let P X (Y ) be an identifiable interventional query over a causal graph G = (V, E), and that we have access to samples D ∼ P (V ) from the observational joint distribution over V . Suppose ID -DAG(Y, X, D, G) immediately enters step 6 and returns H. Then H is a valid sampling network for p X (Y ).Proof. In this case, both ID-DAG(Y, X, D, G) and ID(Y, X, P, G) enter step 6. By the condition of step 6, G\X has only one c-component {S}, where S ∈ C(G). Then the soundness of ID implies that P X (Y ) =S\Y {i|Vi∈S} P (V i |V (i-1) π )where π is a topological ordering for G. ID-DAG operates in this case by training, from joint samples D, a model to correctly sample eachP (V i | V (i-1) π) term, i.e., we learn a conditional generative model f i (V(i-1) π ) which produces samples from P (V i | V (i-1) π), which we can do according to Assumption 1. Then we construct a sampling network H by adding a node V i with sampling mechanism f i for each V i ∈ S. We add edges from V j → V i for each V j ∈ V (i-1) π . Since every vertex in G is either in S or in X, every edge either connects to a previously constructed node or a variable in X. When specified values for X and sampled according to topological order π, this sampling graph provides samples from the distribution {i|Vi∈S} P (V i |V (i-1) π ), i.e. P X (S). We assert the remaining conditions to show that this is correct for P X (Y ): certainly this graph is a DAG and every v ∈ Y has a conditional generative model in H. By the conditions to enter step 6, if C(G \ X) = {S}, then G = S ∪ X and S ∩ X = ∅. Then every node in H is either in S or is in X: hence the only nodes without sampling mechanisms are those in X as desired.Lemma B.10. Let H be a sampling network produced by ID-DAG from an identifiable query P X (Y ) over a graph G. If G has the topological ordering π, then every edge in the sampling graph of H adheres to the ordering π.Proof. We consider two factors: which edges are added, and with respect to which graphs. Since the only base cases ID-DAG enters are steps 1 and 6, the only edges added are consistent with the topological ordering π for the graph that was supplied as an argument to these base case calls. The only graph modifications occur in steps 2 and 7, and these yield subgraphs of G. Thus the original topological ordering π for graph G is a valid topological ordering for each restriction of G. Therefore any edge added to H is consistent with the global topological ordering Π.Theorem B.11. ID-DAG Soundness: Let P X (Y ) be an identifiable query given the causal graph G = (V, E) and that we have access to joint samples D ∼ P (V ). Then the sampling network returned by ID -DAG(Y, X, D, G) correctly samples from P X (Y ) under Assumption B.3.Proof. We proceed by structural induction. We start from the base cases, i.e., the steps that do not call ID-DAG again. ID-DAG only has three base cases: step 1 is the case when no variables are being intervened upon and is covered by Lemma B.8; step 6 is the other base case and is covered by Lemma B.9; step 5 is the non-identifiable case and since we assumed that P X (Y ) is identifiable, due to Lemma B.7, ID-DAG never enters this step. The structure of our proof is as follows. By the assumption that P X (Y ) is identifiable and due to Lemma B.6, its recursions must terminate in steps 1 or 6. Since we have already proven correctness for these cases, we use these as base cases for a structural induction. We prove that if ID-DAG enters any of step 2, 3, 4 or 7, under the inductive assumption that we have correct sampling graphs for the recursive calls, we can produce a correct overall sampling graph. The general flavor of these inductive steps adheres to the following recipe: i) determine the correct recursive call that ID algorithm makes; ii) argue that we can generate the correct dataset to be analogous to the distribution that ID uses in the recursion; iii) rely on the inductive assumption that the generated DAG from the recursion is correct.We consider each recursive case separately. We start with step 2. Suppose ID-DAG(Y, X, D, G) enters step 2, then by the same conditions ID(Y, X, P, G) enters step 2 as well. Hence the correct distribution to sample from is provided by ID step 2:P X (Y ) = ID(Y, X ∩ An(Y ) G , V \An(Y ) G P (V ), G An(Y ) ).Following our recipe, we need to generate the dataset sampled from V \An(Y ) G P , generated only with samples from D ∼ p(V ). This amounts to converting joint samples D to marginal samples D ′ . We do this by dropping all variables downstream of Y (in the graph G) from the dataset D, thereby attaining samples from the joint distribution V \An(Y ) G P (V ). Then we can generate the sampling network from ID-DAG(Y, X ∩ An(Y ) G , D ′ , G An(Y ) ) by the inductive assumption and simply return it.Next, we consider step 3. Suppose ID-DAG(Y, X, D, G) enters step 3. Then by the same conditions, ID(Y, X, P, G) enters step 3, and the correct distribution to sample from is provided from ID step 3 as p X (Y ) = ID(Y, X ∪ W, P, G)where W := (V \ X) \ An(Y ) G X . Since the distribution passed to the recursive call is P , we can simply return the sampling graph generated by ID-DAG(Y, X ∪W, D, G), which we know is correct for P X∪W (Y ) by the inductive assumption. While we do need to specify a sampling mechanism for W to satisfy our definition of a valid sampling network, this can be chosen arbitrarily, say W ∼ P (W ). Note that these freely-chosen distributions will not affect the final samples as they will cancel out due to conditioning on these variables.Next we consider step 4. Suppose ID-DAG(Y, X, D, G) enters step 4. Then by the same conditions, ID(Y, X, P, G) enters step 4 and the correct distribution to sample from is provided from ID step 4 as:V \(y∪x) i ID(s i , v \ s i , P, G)where S i are the c-components of G\X, i.e., elements of C(G\X). By the inductive assumption, we can sample from each term in the product by ID-DAG(s i , v \ s i , D, G). However, recall the output of ID-DAG: ID-DAG returns a 'headless' sampling network as follows: ID-DAG(Y, X, D, G) is a sampling network, i.e., a collection of conditional generative models where for each variable in G and every variable except those in X have a specified conditional generative model. To sample from this sampling network, values for X must first be specified. In the step 4 case, the values V \ s i need to be provided to sample values for s i , and similarly for i ̸ = j, values V \ s j are needed to sample values for s j . Since s i ⊆ (V \ s j ) and s j ⊆ (V \ s i ), this might lead to cycles in the corresponding directed graph of the sampling network obtained by combining each conditional generative model for each c-component. ) . The sum V \(y∪x) can be safely be ignored, because samples from the joint can be marginalized to attain samples from marginals. Hence H is correct for P X (Y ).Step 5 can never happen by the assumption that P X (Y ) is identifiable, and step 6 has already been covered as a base case. The only step remaining is step 7.Supposing ID-DAG(Y, X, D, G) enters step 7, then by the same conditions ID(Y, X, P, G) enters step 7. Then by assumption, C(G\X) = {S} and there exists a confounding component S ′ ∈ C(G) such that S ⊂ S ′ . The correct distribution to sample from is provided from ID step 7 asp X (Y ) = ID(Y, X ∪ S ′ , P ′ , G S ′ )whereP ′ := {i|Vi∈S ′ } p(V i |V (i-1) π ∪ S ′ , v (i-1) π \ S ′ ).Examining ID algorithm more closely, if we enter step 7 during ID, the interventional set X is partitioned into two components: X ∩S ′ and X Z := X \S ′ . From Lemmas 33 and 37 of Shpitser & Pearl (2008) , in the event we enter step 7, P X (Y ) is equivalent to P ′ X∩S ′ (Y ) where P ′ (V ) = P X Z (V ). Hence in order to sample correctly, we need to do two things: first we need to alter our samples D ∼ P to samples from D ′ ∼ P X Z (V ), and then we need to recurse on the query P ′ X∩S ′ (Y ) over the graph G S ′ . Generating samples D ′ ∼ P X Z (V ) can further be reduced to samples only from X ∪ Y , which we note is equivalent to S ′ ∪ X Z . Hence the dataset D ′ need only be supported over S ′ . Therefore we need to generate a dataset from D ′ via ID-DAG(S ′ , X Z , D, G). This is attainable via the inductive assumption and Lemma B.6. The only divergence from ID during generation of D ′ is that ID presumes pre-specified values for X, where we train a sampling mechanism that is agnostic a priori to the specific choice of X and hence X Z . To sidestep this issue, we generate a dataset with all possible values of X Z and be sure to record the values of X Z in the dataset D ′ . Next, we need to map the recursive call ID(Y, X ∩ S ′ , P ′ , G) to ID-DAG. Since instead of passing the distribution P ′ , we pass the dataset D ′ , which does not have specified values for X Z , we need to pass the interventional values for X Z into this recursive call. Hence, instead of intervening in the recursion on X ∩ S ′ , we intervene on X or equivalently (X ∩ S ′ ) ∪ X Z . Additionally, since we are intervening on X ∩ S ′ and X Z , we need to ensure X Z is included in the restricted graph, hence we replace G S ′ with G S ′ ,X Z . By the inductive assumption, we can generate a correct sampling graph from the call ID-DAG(Y, (X ∩ S ′ ) ∪ X Z , D ′ , G S ′ , XZ ), and hence the returned sampling graph is correct for P X (Y ).Since we have shown that every recursion of ID-DAG ultimately terminates in a base case, that all the base cases provide correct sampling graphs, and that correct sampling graphs can be constructed in each step assuming the recursive calls are correct, we conclude that ID-DAG returns the correct sampling graph for P X (Y ).Conditional sampling: Given a conditional causal query P X (Y |Z), we sample from this conditional interventional query by calling Algorithm 5: IDC-DAG. This function finds the maximal set α ⊂ Z such that we can apply rule-2 and move α from conditioning set Z and add it to intervention set X. Z\α) Px∪α(Z\α) . Next, Algorithm 1: ID -DAG(.) is called to obtain the sampling network that can sample from the interventional joint distribution P X∪α (Y, Z \ α). We use the sampling network to generate samples D ′ through feed-forward. A new conditional model M Y is trained on D ′ that takes Z \ α and X ∪ α as input and outputs Y . Finally, we generate new samples with M Y by feeding input values such that Y ∼ P X∪α (Y, Z \ α) i.e, Y ∼ P X (Y |Z). Theorem C.1 (Shpitser & Pearl (2008) ). For any G and any conditional effect P X (Y|W) there exists a unique maximal set Z = {Z ∈ W|P X (Y|W) = P X,Z (Y|W \ Z)} such that rule 2 applies to Z in G for P X (Y|W). In other words, P X (Y|W) = P X,Z (Y|W \ Z).Precisely, P X (Y |Z) = P x∪α (Y |Z \ α) = Px∪α(Y,Theorem C.2 (Shpitser & Pearl (2008) ). Let P X (Y|W) be such that every W ∈ W has a back- for Vj ∈ Hi do 5:door path to Y in G \ X given W \ {W }. Then P X (Y|W) is identifiable in G if and only if P X (Y, W) is identifiable in G.if Vj.f = ∅ and ∃V k ∈ Hr, ∀r such that Vj.name = V k .name and V k .f ̸ = ∅ then 6:Vj = V k 7: H = {Hi} ∀i 8: Return H D.1 SUB-PROCEDURE: ConstructDAG()Each sampling network H i is a set of conditional trained models connected with each other according to a directed acyclic graph. Any node V j in the sampling network, represents a variable name V j .name and its corresponding trained model V j .f . This sub-procedure takes multiple sampling networks {H i } ∀i as input and combines them to construct a larger consistent sampling network H. We iterate over the nodes V j of each sampling network H i (lines 3-4). If V j does not have any corresponding trained model associated with it, it must have its conditional model located in some other sampling network H r as node V k . When we find that node V k in another sampling network H r (line 5), we combine H i and H r by combining V j and V k since V j and V k are the same variable (line 6). In this manner, we combine all sampling network to construct one single network H.ALGORITHM ID algorithm Step 4 1: if C(G \ X) = {S 1 , . . . , S k } then {Step:4} 2: Return v\(y∪x) i ID(s i , v \ s i , P, G ID-DAG algorithm Step 4 1: if C(G \ X) = {S 1 , . . . , S k } then {Step 4} 2:for eachS i ∈ C(G \ X) = {S 1 , . . . , S k } do 3: H i =ID-DAG(S i , X = V \ S i , D, G) 4: Return ConstructDAG({H i } ∀i )Step 4 of the ID algorithm performs Tian's factorization, i.e, splits the variables in G\X into multiple c-components and estimates P v\si (s i ) recursively. Finally, multiplies these factors and marginalizes all variables except y and x. To sample from the corresponding causal query P x (Y), one might feel tempted to sample S i from each P v\si (s i ) in topological order and then combine them. However, this might lead to cyclic situation. For example, consider the causal graph in Figure 6 . Here, according to the ID algorithm step 4, P x (y) = w1,w2 P w1 (x, w 2 ) * P x,w2 (w 1 , y). To sample X, W 2 ∼ P w1 (x, w 2 ) we need W 1 as input which has to be sampled from P x,w2 (w 1 , y). But to sample W 1 ∼ P x,w2 (w 1 , y), we need X, W 2 as input which has to be sampled from P w1 (x, w 2 ). Therefore, no order helps to sample all X, W 1 , W 2 , Y consistently. ID-DAG solves this cyclic issue by avoiding direct sampling from the c-components. Rather ID-DAG builds a sampling network H i consists of trained models for each c-component S i . After that, sub-procedure ConstructDAG is called to merge the sampling network found from each c-component to build a one single sampling network. We use this single sampling network to sample [S 1 , . . . , S k ] ∼ P V \S1 (S 1 ) * . . . P V \S k (S k ). ID-DAG drops all variables from [S 1 , . . . , S k ] except Y and these will be the samples from P x (Y ). Return ID(y, x ∩ S ′ , P ={i|Vi∈S ′ } P (V i |V (i-1) π ∩ S ′ , v (i-1) π \ S ′ ), G S ′ ) ID-DAG algorithm Step 7 1: if ∃S ′ s.t. S ⊂ S ′ ∈ C(G) then {Step 7} 2: X Z = X \ S ′ , 3: H S ′ = ID-DAG(S ′ , X Z , D, G) 4: D ′ ∼ H S ′ (X Z ) 5: Return ID-DAG(Y, X, D ′ , G S ′ ,X Z )Here we show the differences between step 7 of the ID algorithm and the ID-DAG algorithm which represent the technical novelty of our method. This case occurs when i) G \ X is a single ccomponent S, ii) all the variables in Y are contained within a single c-component S ′ ⊂ S, but iii) the variables in X can be partitioned into those that are contained within S ′ and those that are not in S ′ . ID algorithm usesP = {i|Vi∈S ′ } P (V i |V (i-1) π ∩ S ′ , v (i-1) π \ S ′ ), x = x ∩ S ′ and G = G S ′as arguments values for the next recursive call.ID-DAG modification with respect to the dataset D: This new P is actually an interventional distribution, P X\S ′ (S ′ ) which is identifiable from obs data with the mentioned expression. It is not straightforward to see the role of the new distribution argument value used in the ID algorithm. Let X Z = X \ S ′ , i.e, the subset of intervention set X that is not contained within S ′ . ID-algorithm under the hood using the fact that evaluating P X (Y ) where P = P (V ) is equivalent to evaluating P ′ X∩S ′ (Y ) where P ′ = P X Z (V ). Therefore, to train the conditional models in the subsequent steps in the recursion, ID-DAG generates an interventional dataset from P X Z (V ). To do that ID-DAG first obtains the sampling network with a ID-DAG(S ′ , X Z , D, G) call. Then the algorithm uses this sampling network to obtain the interventional dataset D ′ ∼ P X Z (V ) which is used in the next recursive call.ID-DAG modification with respect to the graph G: The graph argument ID algorithm sends to the next recursive call is G S ′ , i.e, it removes the intervened variables X Z from the graph G. This works well for the ID algorithm because all variables values are assumed at the beginning and it estimates the corresponding probability. However, this does not work if our aim is sampling. Since what values of X Z we are using here depend on the values X Z has taken in another c-component, we consider all possible values of X Z while training all conditional models in any level of the recursion for a fixed c-component. Thus, we do not remove X Z from the causal graph but remove the incoming edges since we intervened on them. As a result, when ever we train any conditional model in a deeper level of the recursion, it will be used as an input condition. This is a very important step to sample from P X (Y ) consistently. As a warm-up, we outline the generation for the Napkin-MNIST dataset in a low-dimensional setting. When we consider in the next section the high-dimensional case, we simply replace some of these discrete variables with MNIST images which can be mapped back into this low-dimensional case.We start by enumerating the joint distribution and the support of each marginal variable. First lets define the sets • COLORS := {red, green, blue, yellow, magenta, cyan}.• RG COLORS := {red, green}.• THICKNESSES := {thin, regular, thick}.• DIGITS := {0, . . . , 9}.And then the definitions and support of each of the variables in our distribution:• (Latent) Color ∈ COLORS. • (Latent) Thickness ∈ THICKNESSES. • W 1 ∈ DIGITS × COLORS × THICKNESSES • W 2 ∈ DIGITS × RG COLORS. • X ∈ DIGITS × COLORS × THICKNESSES • Y ∈ DIGITS × COLORS × THICKNESSESNow we describe the full data generation procedure. A key hyperparameter is a noise-probability p. This defines the probability that any variable flips to a uniform probability. To ease notation, we define the function η p (v, S) defined as η p (v, S) := v with probability 1p U (S) otherwise and we define the mapping R : COLORS → RESTRICTED COLORS as R(c) := red if c ∈ {red, green, blue} green otherwise Where U (S) means a uniformly random choice of S. Then our data generation procedure follows the following steps:• Color := U (COLORS)• Thickness := U (THICKNESSES)• W 1 := U (DIGITS), η p (Color, COLORS), η p (Thickness, THICKNESSES)• W 2 := η p (W 1 . digit , DIGITS), η p (R(W 1 . color , RG COLORS) • X := η p (W 2 . digit , DIGITS), η p (W 2 . color , RG COLORS), η p (Thickness, THICKNESSES) • Y := η p (X. digit , DIGITS), η p (Color, COLORS), η p (X. thickness , THICKNESSES)It is easy to verify that this describes the Napkin graph, as each only Color, Thickness are latent and each variable only depends on its parents in the SCM.Secondly, observe that this structural causal model is separable with respect to digits, colors, and thicknesses. Since each digit only depends on parent digits, each color only depends on parent colors, and each thickness depends only on parent thicknesses, these can all be considered separately.Further, because this distribution is only supported over discrete variables, exact likelihoods can be computed for any conditional query. This is much more easily done programmatically, however, and we provide code in the attached codebase to do just that. We will claim without proof that in the case of thicknesses and digits, P Y (X) = P (Y |X). However in the case of colors, P Y (X) ̸ = P (Y |X). Hence we consider this case in the evaluations in the experiments section.The high-dimensional case follows the discrete case of the Napkin-MNIST dataset, with a few key changes. Namely, W 1 , X, and Y are MNIST images that have been colored and thickened. We explicitly outline these changes:• W 1 : A random MNIST image of the provided digit is used, then colored and thickened accordingly (noisy from latents).• W 2 : This is a discrete variable, only encoding the (noised) digit and (noised) restricted color of W 1 .• X: This is a random MNIST image of the (noised) digit obtained from W 2 , then colored with the (noised) restricted color from W 2 and thickened according to the (noised) latent thickness.• Y : This is the same base image of X, unless the noising procedure calls for a change in digit, then a random MNIST image of the specified image is used. The (noisy) color is obtained from the latent distribution, and the (noisy) thickness is obtained from X.To color the images, we convert each 1-channel MNIST image into a 3-channel MNIST image, and populate the necessary channels to generate these colors. Note that in RGB images: if only the RG channels are active, the image is yellow; if only the RB channels are active, the image is magenta; if only the BG channels are active, the image is cyan. To thicken the images, we use the MorphoMNIST package2 . Operationally, we generate a base dataset for our experiments of size equivalent to the original MNIST dataset. That is, the training set has a size of 60K, and the test set has a size of 10K. Because we have access to the latents during the data generation procedure, we are able to train classifiers for each variable to identify their digit, color and thickness. We use a simple convolutional network architecture for each of these cases and achieve accuracy upwards of 95% in each case.We train two diffusion models during our sampling procedure, and we discuss each of them in turn.To train a model to sample from P ( X, Y |W 1 , W 2 ), we train a single diffusion model over the joint (X, Y ) distribution, i.e., 6 channels. We train a standard UNet architecture where we follow the conditioning scheme of classifier-free guidance. That is, we insert at every layer an embedding of the W 1 (image) and W 2 (2-dimensional discrete variable). To embed the W 1 image, we use the base of a 2-layer convolutional neural network for MNIST images, and to embed the W 1 we use a standard one-hot embedding for each of the variables. All three embeddings are concatenated and mixed through a 2-layer fully connected network to reach a final embedding dimension of 64. Batch sizes of 256 are used everywhere. Training is performed for 1000 epochs, which takes roughly 9 hours on 2 A100 GPU's. Sampling is performed using DDIM over 100 timesteps, with a conditioning weight of w = 1 (true conditional sampling) and noise σ = 0.3.To train a model to sample Y from the generated dataset (W 2 , X, Y ), we follow an identical scheme.Perhaps the most correct thing is to train a single diffusion model for each choice of W 2 in our synthetic dataset, however we argue our schema still produces correct samples because: 1) W 2 can be arbitarily chosen, and thus should not affect Y , 2) we argue that the model fidelity benefits from weight sharing across multiple choices of W 2 , 3) the model is only ever called with a specified value of W 2 so we always condition on this W 2 .In addition to the evaluations presented in the main paper, we can further perform evaluations on the component models necessary to sample P X (Y ).P (X, Y |W 1 , W 2 ):We can evaluate the model approximating samples from P (X, Y |W 1 , W 2 ) on a deeper level than just visual inspection as provided in the main paper. In particular, assuming access to good classifiers that can predict the digit, color, and thickness of an MNIST image, we can compare properties of the generated images with respect to the ground truth in the discrete case.For example, assuming we have hyperparameter of random noise equal to p, we can compute the following quantities analytically on the discrete dataset as:• P [X d = W 2d ] = 1 -p + p 10 • P [X. c = W 2c ] = 1 -p + p 10 • P [X t = W 1t ] = (1 -p + p 3 ) 2 + p 3 2 • 2 • P [Y d = X d ] = 1 -p + p 10 • P [Y c = W 1c ] = (1 -p + p 6 ) 2 + big( p 6 2 • 5 • P [Y t = X t ] = 1 -p + p 3where V d , V c , V t refer to the digit, color, and thickness attributes respectively. These calculations follow from two formulas. In a discrete distribution with support S and |S| = K:• P [η p (z, S) = z] = 1 -p + p K • P [η p (z, S) = η p (z, S)] = (1 -p + p K ) 2 + p K 2 • (K -1)where in the second equation, it is assumed that η p (•, •) are two independent noising procedures.Then to evaluate, we can 1) consider a large corpus of joint data, 2) run each of W 1 , X, Y through a classifier for digit, color, and thickness, 3) evaluate the empirical estimate of each desired probability. We present these results for the synthetic dataset D synth sampled from the diffusion model approximating P (X, Y |W 1 , W 2 ), a dataset D orig generated according to the data generation procedure, and P true the true analytical probabilities. These results are displayed in the Table 3 .Here we compare our performance with two baselines and provide the results below: We provide a short description of the ground truth first. We have the napkin causal graph: < l a t e x i t s h a 1 _ b a s e 6 4 = " G W V t S K 6 e t T E 3 0 Algorithm execution: To easily understand the distinction among different methods, we chose two images i) of digit 3 and ii) of digit 5, both colored red as the intervention value for X. Then we used these images to sample from the corresponding distribution P (Y |X) and P (Y |do(X)) distributions using different methods. Although the implementation of these methods is different, we considered the performance of each method after running them for 300 epochs. The NCM algorithm took around 50 hours to complete while our algorithm took approximately 16 hours to complete.W 1 → W 2 → X → Y ; X ↔ W 1 ↔ Yv Z 2 V 1 b X 1 j s 7 B V 3 N 7 Z 3 d s v H R w 2 d Z w q h g 0 W i 1 i 1 A 6 p R c I k N w 4 3 A d q K Q R o H A V j C 6 m / q t J 1 S a x / L B j B P 0 I z q Q P O S M G i v VO i Z q g X v a n 4 n 9 d J T X j j Z 1 w m q U H J 5 o v C V B A T k + n X p M 8 V M i P G l l C m u L 2 V s C F V l B m b T d G G 4 C 2 + v E y a 5 x X v q n J Z v y h X b / M 4 C n A M J 3 A G H l x D F e 6 h BD i + G M d K A H h V 0 Y O 2 6 F Z h R C 3 q o s 3 g = " > A A A B 8 X i c b V D L T g J B E O z F F + I L 9 e h l I j G B C 9 k 1 v o 5 E L x 4 x k Y f C h s z O z s K E 2 Z n N z K w J Q f 7 C i w e N 8 e r f e P N v H G A P i l b S S a W q O 9 1 d Q c K Z N q 7 7 5 e S W l l d W 1 / L r h Y 3 N r e 2 d 4 u 5 e U 8 t U E d o g k k v V D r C m n A n a M M x w 2 k 4 U x X H A a S s Y X k 3 9 1 g N V m k l x a 0 Y J 9 W P c F y x i B B s r 3 d f L d 4 + h L L c r l V 6 x 5 F b d G d B f 4 m W k B B n q v e J n N 5 Q k j a k w h G O t O 5 6 b G H + M l W G E 0 0 m h m 2 q a Y D L E f d q x V O C Y a n 8 8 u 3 i C j q w S o k g q W 8 K g m f p z Y o x j r U d x Y D t j b A Z 6 0 Z u K / 3 m d 1 E Q X / p i J J D V U k P m i K O X I S D R 9 H 4 V M U W L 4 y B J M F L O 3 I j L A C h N j Q y r Y E L z F l / + S 5Z 2 5 A C n R Q A + 9 q q M = " > A A A B 9 H i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v g q s y I V Z d F N y 4 r 2 A e 0 Q 8 m k m T Y 0 y Y x J p l C G f o c b F 4 q 4 9 W P c + T d m 2 l l o 6 4 H A 4 Z x 7 u S c n i D n T x n W / n c L a + s b m V n G 7 t L O 7 t 3 9 Q P j x q 6 S h R h D Z J x C P V C b C m n E n a N M x w 2 o k V x S L g t B 2 M 7 z K / P a F K s 0 g + m m l M f Y G H k o W M Y G M l v y e w G S m R S i J m p X 6 5 4 l b d O d A q 8 X J S g R y N f v m r N 4 h I I q g 0 h G O t u 5 4 b G z / F y j D C 6 a z U S z S N M R n j I e 1 a K r G g 2 k / n o W f o z C o D F E b K P m n Q X P 2 9 k W K h 9 V Q E d j I L q Z e 9 T P z P 6 y Y m v P F T J u P E U E k W h 8 K E I x O h r A E 0 Y I o S w 6 e W Y K K Y z Y r I C C t M j O 0 p K 8 F b / v I q a V 1 U v a t q 7 e G y U r / N 6 y j C C Z z C O X h w D X W 4 h w Y 0 g c A T P M M r v D k T 5 8 V 5 d z 4 W o w U n 3 z m G P 3 A + f w D b 0 J I r < / l a t e x i t > ncm < l a t e x i t s h a 1 _ b a s e 6 4 = " Q I h K h L 5 / d f x 4 / A 0 A 4 k y Z U z m 6 r B w = " > A A A B 9 X i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v g q s y I V Z d F N y 4 r 2 A e 0 Y 8 m k m T Y 0 y Q x J R i l D / 8 O N C 0 X c + i / u / B s z 7 S y 0 9 U D g c M 6 9 y c k J Y s 6 0 c d 1 v p 7 C y u r a + U d w s b W 3 v 7 O 6 V 9 w 9 a O k o U o U 0 S 8 U h 1 A q w p Z 5 I 2 D T O c d m J F s Q g 4 b Q f j m 8 x v P 1 K l W S T v z S S m v s B D y U J G s L H S Q 0 9 g M 1 I i t d f p a a l f r r h V d w a 0 T L y c V C B H o 1 / + 6 g 0 i k g g q D e F Y 6 6 7 n x s Z P s T K M c D o t 9 R J N Y 0 z G e E i 7 l k o s q P b T W e o p O r H K A I W R s k c a N F N / b 6 R Y a D 0 R g Z 3 M U u p F L x P / 8 7 q J C a / 8 l M k 4 M V S S + U N h w p G J U F Y B G j B F i e E T S z B R z G Z F Z I Q V J s Y W l Z X g L X 5 5 m b T O q t 5 F t X Z 3 X q l f 5 3 U U 4 Q i O 4 R Q 8 u I Q 6 3 E I D m k B AQ 8 t E E V o n k k v V C r C m n A l a N 8 x w 2 o o V x V H A a T M Y X k / 8 5 i N V m k l x Z 0 Y x 9 S P c F y x k B B s r N W r l + 6 f W c b d Y c i v u F G i R e B k p Q Y Z a t / j V 6 U m S R F Q Y w r H W b c + N j Z 9 i Z R j h d F z o J J r G m A x x n 7 Y t F T i i 2 k + n 1 4 7 R k V V 6 K J T K l j B o q v 6 e SD i + G M d K A H h V 0 Y O 2 6 F Z h R C 3 q o s 3 g = " > A A A B 8 X i c b V D L T g J B E O z F F + I L 9 e h l I j G B C 9 k 1 v o 5 E L x 4 x k Y f C h s z O z s K E 2 Z n N z K w J Q f 7 C i w e N 8 e r f e P N v H G A P i l b S S a W q O 9 1 d Q c K Z N q 7 7 5 e S W l l d W 1 / L r h Y 3 N r e 2 d 4 u 5 e U 8 t U E d o g k k v V D r C m n A n a M M x w 2 k 4 U x X H A a S s Y X k 3 9 1 g N V m k l x a 0 Y J 9 W P c F y x i B B s r 3 d f L d 4 + h L L c r l V 6 x 5 F b d G d B f 4 m W k B B n q v e J n N 5 Q k j a k w h G O t O 5 6 b G H + M l W G E 0 0 m h m 2 q a Y D L E f d q x V O C Y a n 8 8 u 3 i C j q w S o k g q W 8 K g m f p z Y o x j r U d x Y D t j b A Z 6 0 Z u K / 3 m d 1 E Q X / p i J J D V U k P m i K O X I S D R 9 H 4 V M U W L 4 y B J M F L O 3 I j L A C h N j Q y r Y E L z F l / + S 5i) Image quality: First, we discuss the image quality of the generated samples from each approach that we provided in Figure 8 . The conditional model (row 1, row 4) and our algorithm (row 2, row 5) both generate high-quality images of digit 3 and digit 5 with a specific color from the six possible colors. Whereas, the colors in generated images from the NCM algorithm (row 3, row 6) can not remain at any specific value and get blended (such as green+yellow, etc).In Table 4 , we provide the Frechet Inception Distance (FID) of each method (the lower the better). We observe that our algorithm has the lowest FID score, i.e., our algorithm generates the most high-quality images from interventional distribution. ii) Correctness of the sampling distribution: We use a classifier to detect the color frequency from the generated images of the conditional model and our algorithm. We observe the probabilities in Table 5 .In the generated samples from the conditional model, the color of digit image X and digit image Y are correlated due to confounding through backdoor paths. For example, for a digit image X with color as red, Y takes a value from [Red, Green, Blue] with high probability. Thus, the conditional model does not have the ability to generate interventional samples. On the other hand, our algorithm generates samples from the interventional distribution P (Y |do(X)) and the generated samples chooses different colors:[Red, Green, Blue, Yellow, Magenta, Cyan] with almost the same probability.Therefore, we show our performance with a close baseline NCM and illustrate our capability of sampling from high-dimensional interventional distribution.Now we consider the Covid X-Ray Dataset. We first outline the steps we perform to preprocess our data and then discuss training details for each of the models we train.Note that a full pipeline of data preparation is contained in cxray/prep cxray dataset.sh in the provided codebase. We start by downloading the corpus approximately 30K Covid X-Ray images3 . Then we download Covid-19 labels and Pneumonia labels4 and attach labels to each image. Then we convert each image to black-white one-channel images and rescale each to a size of (128 × 128) pixels. Finally, a random split of the 30K images is performed: keeping 20K to be used during training, and 10K to be used as validation images. Note that the labels come with a set of 400 test images, but 400 images is too small to be an effective test set (say for FID computations). We will be explicit about where we use each data set.We train a diffusion model to approximate P (X|C). To do this, we use a standard UNet and classifier-free guidance scheme. We train for 100K training steps over the 20K-sized training set, using a batch size of 16. This takes roughly 10 hours on a single A100. The same classifier-free guidance parameters are used as in the NapkinMNIST diffusion training.To train a classifier to sample from P (N |C, X), we note that our inputs are a (1, 128, 128) image and a binary variable. Our architecture is as follows: we create an embedding of X by modifying the final linear layer of a ResNet18 to have output dimension 64 (vs 1000), and modify the input channels to be 1. We create an embedding for N by using a standard embedding layer for binary variables, with embedding dimension 64. These embeddings are then concatenated and pushed through 3 fullyLabels are from https://github.com/giocoal/CXR-ACGAN-chest-xray-generator-covid19-pneumonia/https://github.com/dccastro/Morpho-MNIST/https://www.kaggle.com/datasets/andyczhao/covidx-cxr2https://github.com/giocoal/CXR-ACGAN-chest-xray-generator-covid19-pneumonia 
 
 Appendix:  


Decision: The submission provides an algorithm for sampling from a causal interventional distribution using conditional generative models, building on Shpitser and Pearl's ID algorithm.

While most reviewers consider the contribution to be sound and valuable for the field, there has been concerns that the manuscript is hard to parse. Overall, the contribution seems technically sound, but highly inaccessible to a majority of the audience, in addition to the significance of the contribution being highly challenging to assess given the complexity of the assumptions and statement of the results. The AC recommends to make an effort towards the general readership and recommends acceptance.

Reviews: review: paper_summary: paper_summary: This paper provides an algorithm for sampling from a causal interventional distribution using conditional generative models, building on Shpitser and Pearl's ID algorithm. They discuss how their algorithm, ID-GEN, can sample from any identifiable distribution given a causal graph, and handles the presence of unobserved confounders (when identifiable). Empirically, they demonstrate their method can work for measurement, evaluation and interpretability purposes in the challenging setting where both the target and treatment are high dimensional e.g. images.
strength_weakness: strength_weakness: - I frankly got a bit lost in a few key parts of Section 3. I got the main ideas (I think) but missed a bunch of nuance. Some spots were: Example 4.1 (I don’t understand why ID fails in this case but ID-GEN succeeds - specifically the importance of merging is a bit lost on me), Step 4 of ID-GEN (again, merging), and Step 7 of ID-GEN (I think the logic around how training data is sampled, used and modified wrt the graph needs to be explained more clearly)
- Step 1 of ID-GEN confuses me - I don’t see why we can’t just learn a model of P(y) directly in this case? Also the 2nd equality in 203 doesn’t make sense to me - how is the sum over values of v equal to P(y)?
- in each experimental section I find I have at least a medium-sized point of confusion around the setup or evaluation - more care should be taken to explain empirical setup + results overall
- In 5.1, the authors state that U_color makes W1 and Y correlated by color - however, X contains color information and is a direct ancestor of Y, so this unobserved confounding seems trivial
- in 5.1, it seems like a better metric than a classifier (which may be unreliable and as you note isn’t useful for all possible images) would be something based specifically on the RGB values of the pixels themselves
- in 5.2, I don’t quite see why Young & Male have unobserved confounding - are they not fully determined by the shared + observed parent I_1?
- in 5.3, I don’t understand why the report is a causal ancestor in the graph - isn’t it generated upon viewing the X-ray?
- in 5.3, I think the setup with the labeller can be made clearer - how good is this labeller? How is it structured? Additionally, is the bottom row intended to be a success or failure examples? (label says it should be right lung base but all inferences name the left lung)


Smaller points:
- L115: are unobserved confounders only allowed to affect 2 variables in this framework? Is that more limiting than general SCMs?
questions: questions: - would be great to see clarification throughout Section 3, particularly in the highlighted areas and around merging and training data sampling
- experimental setups all need clarifications
- generally assuming I understood what's happening better I'd be happy to increase my score
limitations: limitations: n/a
paper_summary: paper_summary: This paper studies the problem of sampling from an interventional distribution of high-dimensional data, where the causal relationships are described by a acyclic directed mixed graph (ADMG). Motivated by the ID algorithm that provides a recursive way of identifying any causal effect from a conditional probability table, the authors propose ID-GEN that follows the recursive steps of ID but instead trains deep generative models to fit the conditional distributions. The final sampling model is then obtained by connecting all the trained networks together in some suitable way. The authors prove that ID-GEN is sound and complete, and run extensive experiments on both synthetic and real dataset to demonstrate the effectiveness of their approach.
strength_weakness: strength_weakness: I don't think this paper has obvious weaknesses. One thing that the authors may wish to improve is that the notations are a litle bit complex; and it would be better to more often remind the authors of their meanings.
questions: questions: 1. In some identification formula e.g. Eq.(1) in the paper, the probability on the denominator might be small. To what extent would this affect the stablity of the proposed algorithm.

2. Can your algorithm be straightforwardly adapted to the case of imperfect interventions i.e. some conditional probabilities in the structural causal model are modified, but no causal edges are removed?
limitations: limitations: N/A
paper_summary: paper_summary: This paper proposes an algorithm for sampling from intervention distributions under known causal models with hidden confounders, using conditional generative models.
This allows nonparametric distributional estimation of high-dimensional outcomes such as X-ray image generation, unlike existing methods.
The proposed method combines the existing ID algorithm with a generative model and inherits the ID algorithm's theoretical guarantees. That is, non-identifiable quantities are indicated to be non-identifiable, while all identifiable quantities can be estimated.
strength_weakness: strength_weakness: [W1] The motivation for high-dimensional distribution estimation is weak. For example, it does not seem very meaningful to me to generate synthetic X-ray images.

[W2] In particular, is it important in cases where there are bidirectional edges due to the presence of hidden confounding factors, but where the causal orientations are all identified among variables? A clear comparison with similar methods would be beneficial for readers, e.g., comparison in assumptions and targets (e.g., parametric/nonparametric, latent confounder, distributional estimation, etc.).

[W3] The base procedure seems to come from existing methods, such as the ID algorithm, and they just combine it with a generative model.
questions: questions: [Q1] Related to W2, is the proposed method a non-trivial novelty in situations where there are no unobserved confounding variables, i.e. no bidirectional edges?

[Q2] Related to W3, are there any non-trivial points in theoretical guarantees when the generative model is combined with the ID algorithm?
limitations: limitations: The method is limited to the cases where the causal direction is known.
paper_summary: paper_summary: The paper leverages state-of-the-art conditional generative models and algorithms from causal do calculus to perform "approximately correct" high-dimensional interventional sampling. Their contribution is ID-GEN, a recursive algorithm that uses diffusion models (among other generative models) to sample from any identifiable causal effect estimand in high-dimensional settings. The efficacy of the method is demonstrated on three diverse datasets: Colored MNIST, CelebA, and MIMIC-CXR.
strength_weakness: strength_weakness: - The paper is objectively hard to read. Many important graphical elements that should be paired with the text in the main paper are delayed to supplements. This is notably problematic for Example 4.1, where one would expect the full example to be self-contained in the main text. 
- The contributions are stated in the introduction, but it still seems hard to understand if the proposed method is "just" an implementation of the ID algorithm, replacing probability tables by samples from diffusion models. I appreciate that this is hard already, but it has much lower novelty then proposing a new recursive algorithm. This should be well explained in the manuscript. 
- The paper does not discuss the implications of the proposed algorithm. Is there any way to extend this to symbolic calculus, or to probabilistic programming? What are the obstacles for moving towards automatic causal inference with images (which would be a super exciting prospect).
questions: questions: - Can you clarify whether the proposed ID-GEN algorithm is primarily an implementation of the existing ID algorithm with diffusion models substituting for probability tables, or does it introduce fundamentally new recursive methodologies? How does this distinction affect the perceived novelty of your contribution?
- What are the potential extensions of your algorithm to areas like symbolic calculus or probabilistic programming?
- Are there any specific obstacles that need to be overcome to advance towards automatic causal inference using high-dimensional data such as images? How feasible is this prospect in the near future?
limitations: limitations: N/A

Comments: title: Thanks
comment: Dear Reviewer fUEc,

We cordially thank you for your response and highly appreciate that you raised your score. 

>For [W1], the lack of motivating examples has not been resolved.


We apologize that we misunderstood your concern about the motivational example. To address this concern, we will illustrate the scope of high-dimensional causal inference in another real-world scenario involving brain MRI [1] in our introduction section. Here the goal would be to generate samples from the interventional distribution involving MRI images with attributes such as age, sex, brain volume and ventricle volume [2]. 

Our algorithm, ID-GEN would play an important role in estimating causal effects when we consider the presence of confounders among variables, as existing works such as [2] does not consider confounders in the system. We will explain this motivational example in detail.


[1] Sudlow, Cathie, et al. "UK biobank: an open access resource for identifying the causes of a wide range of complex diseases of middle and old age." PLoS medicine 12.3 (2015): e1001779.\
[2] Ribeiro, Fabio De Sousa, et al. "High fidelity image counterfactuals with probabilistic causal models." arXiv preprint arXiv:2306.15764 (2023).
Thank you for your response.

For [W1], the lack of motivating examples has not been resolved.

For [W2], I understand the problem of unidentified causal direction would be another direction, and the proposed method can be extended.

For [W3], I understand that there would be some issues in a straightforward extension of the ID algorithm to generative models.

Based on the resolution of these concerns, the score is slightly increased.
We cordially thank you for reading our responses. We are also happy to know that  you would consider changing your score after the discussion with other reviewers.
I would like to thank the authors for their detailed response. I will consider changing my score while discussing with other reviewers.
Dear Reviewer 6o9W,

We thank you for your reply. We agree that some steps of our algorithm might seem a little hard to understand. We aim to add more details on those steps according to your feedback. If you prefer any further clarification, please let us know and we would be happy to elaborate.

Thank you.
Thanks for the rebuttal - I appreciate the corrections on some of the points. I think that given that my main issues were around general methodological confusion I probably will lean against increasing my score, but happy to consider this given discussion with other reviewers.
Dear Reviewer fUEc,

We thank you for your invested valuable time to provide constructive feedback for our paper.
We were wondering if we have addressed all your concerns including the theoretical and technical novelties of our algorithm.
If you have any additional questions in your mind, please let us know. We would be more than happy to answer them.

Also, we have used some latex notations in our responses. Please refresh the page if they do not render properly. 

Thank you.