Consistency of Neural Causal Partial Identification

Authors: jiyuan  tan
jose  blanchet
vasilis  syrgkanis
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 20:02:09

Abstract: Recent progress in Neural Causal Models (NCMs) showcased how identification and partial identification of causal effects can be automatically carried out via training of neural generative models that respect the constraints encoded in a given causal graph [Xia et al. 2022, Balazadeh et al. 2022]. However, formal consistency of these methods has only been proven for the case of discrete variables or only for linear causal models. In this work, we prove the consistency of partial identification via NCMs in a general setting with both continuous and categorical variables. Further, our results highlight the impact of the design of the underlying neural network architecture in terms of depth and connectivity as well as the importance of applying Lipschitz regularization in the training phase. In particular, we provide a counterexample showing that without Lipschitz regularization this method may not be asymptotically consistent. Our results are enabled by new results on the approximability of Structural Causal Models (SCMs) via neural generative models, together with an analysis of the sample complexity of the resulting architectures and how that translates into an error in the constrained optimization problem that defines the partial identification bounds.

Main paper: 
 introduction: 
 Identifying causal quantities from observational data is an important problem in causal inference which has wide applications in economics [1] , social science [20] , health care [34, 18] , and recommendation systems [9] . One common approach is to transform causal quantities into statistical quantities using the ID algorithm [49] and deploy general purpose methods to estimate the statistical quantity. However, in the presence of unobserved confounding, typically the causal quantity of interest will not be point-identified by observational data, unless special mechanisms are present in the data generating process (e.g. instruments, unconfounded mediators, proxy controls). In the presence of unobserved confounding, the ID algorithm will typically fail to return a statistical quantity and declare the causal quantity as non-identifiable.One remedy to this problem, which we focus on in this paper, is partial identification, which aims to give informative bounds for causal quantities based on the available data. At a high level, partial identification bounds can be defined as follows: find the maximum and the minimum value that a target causal quantity can take, among all Structural Causal Models (SCMs) that give rise to the same observed data distribution and respect the given causal graph (as well as any other structural constraints that one is willing to impose). Note that in the presence of unobserved confounding, there will typically exist many structural mechanisms that could give rise to the same observational distribution but have vastly different counterfactual distributions. Hence, partial identification can be formulated as solving a max and a min optimization problem [3] maxEQUATIONsubject to P M (V ) = P M * (V ), andG M = G M * ,where θ(M) is the causal quantity of interest, M * is the true model, V is the set of observed nodes, P M (V ) is the distribution of V in SCM M, C is a collection of causal models and G M is the causal graph of M (see Section 2 for formal definitions). Two recent lines of work explore the optimization approach to partial identification. The first line deals with discrete Structure Causal Models (SCMs), where all observed variables are finitely supported. In this case, (1) becomes a Linear Programming (LP) or polynomial programming problem and tight bounds can be obtained [37, 4, 5, 27, 44, 45, 40, 8, 54, 56, 14] . The second line of work focuses on continuous models and explores ways of solving (1) in continuous settings using various techniques [24, 32, 29, 38, 3, 41] .Recently, Xia et al. [51] formalized the connection between SCMs and generative models (see also [33] for an earlier version of a special case of this connection). This work showcased that SCMs can be interpreted as neural generative models, namely Neural Causal Models (NCMs), that follow a particular architecture that respects the constraints encoded in a given causal graph. Hence, counterfactual quantities of SCMs can be learned by optimizing over the parameters of the underlying generative models. However, there could be multiple models that lead to the same observed data distribution, albeit have different counterfactual distributions. Xia et al. [51] first analyze the approximation power of NCMs for discrete SCMs and employ the max/min approach to verify the identifiability of causal quantities, without the need to employ the ID algorithm. Balazadeh et al. [3] and Hui et al. [30] , extend the method in [51] and re-purpose it to perform partial identification by solving the min and max problem in the partial identification formulation over neural causal models.However, for SCMs with general random variables and functional relationships, the approximation error and consistency of this optimization-based approach to partial identification via NCMs has not been established. In particular, two problems remain open. First, given an arbitrary SCM, it is not yet known if we can find an NCM that produces approximately the same intervention distribution as the original one. Although Xia et al. [51] show it is possible to represent any discrete SCM by an NCM, their construction highly relies on the discrete assumption and cannot be directly generalized to the general case. Moreover, Xia et al. [51] use step functions as the activation function in their construction, which may create difficulties in the training process since step functions are discontinuous with a zero gradient almost everywhere.Second, since we only have access to n samples from the true distribution, P M (V ), we need to replace the constraints in (1) with their empirical version that uses the empirical distribution of samples P M * n (V ) in place of the population distribution and looks for NCMs, whose implied distribution lies within a small distance from the empirical distribution. Moreover, even the NCM distribution is typically only accessible through sampling, hence we will need to generate m n samples from the NCM and use the empirical distribution of the m n samples from the NCM in place of the true distribution implied by the NCM. Thus, in practice, we will use a constraint of the form d(P M * n (V ), P M mn (V )) ⩽ α n , where d is some notion of distribution distance and α n accounts for the sampling error. It is not clear that this approach is consistent, converges to the correct partial identification bounds, when the sample size n grows to infinity. Balazadeh et al. [3] only show the consistency of this approach for linear SCMs. Consistency results concerning more general SCMs is still lacking in the neural causal literature.In this paper, we establish representation and consistency results for general SCMs. Our contributions are summarized as follows.• We show that under suitable regularity assumptions, given any Lipschitz SCM, we can approximate it using an NCM such that the Wasserstein distance between any interventional distribution of the NCM and the original SCM is small. Each random variable of the SCM is allowed to be continuous or categorical. We specify two architectures of the Neural Networks (NNs) that can be trained using common gradient-based optimization algorithms (Theorem 2, Theorem 3 and Corollary 1).• To construct the NCM approximation, we develop a novel representation theorem of probability measures (Proposition 1) that may be of independent interest. Proposition 1 implies that under certain assumptions, probability distributions supported on the unit cube can be simulated by pushing forward a multivariate uniform distribution. • We discover the importance of Lipschitz regularization by constructing a counterexample where the neural causal approach is not consistent without regularization (Proposition 2). • Using Lipschitz regularization, we prove the consistency of the neural causal approach (Theorem 4).Related Work There exists a rich literature on partial identification of average treatment effects (ATE) [37, 4, 5, 27, 44, 45, 40, 8, 54, 56, 14, 7, 39, 26] . Balke and Pearl [4, 5] first give an algorithm to calculate bounds on the ATE in the Instrumental Variable (IV) setting. They show that regardless of the exact distribution of the latent variables, it is sufficient to consider discrete latent variables as long as all endogenous variables are finitely supported. Moreover, they discover that ( 1) is an LP problem with a closed-form solution. This LP-based technique was generalized to several special classes of SCMs [44, 8, 27, 46] . For general discrete SCMs, [14, 56, 55] consider transforming the problem (1) into a polynomial programming problem. Xia et al. [51] discover the connection between generative models and SCMs. They show that NCMs are expressive enough to approximate discrete SCMs. By setting C in problem (1) to be the collection of NCMs, they apply NCMs for identification and estimation.For causal models with continuous random variables, the constraints in (1) become integral equations, which makes the problem more difficult. One approach is to discretize the constraints. [24] uses stochastic process representation of the causal model in the continuous IV setting and transforms the problem into a semi-infinite program. [32] relaxes the constraints to finite moment equations and solves the problem by the Augmented Lagrangian Method (ALM). The other approach is to use generative models to approximate the observational distribution and use some metric to measure the distance between distributions. [33] first propose to use GAN to generate images. Later, [29] uses Wasserstein distance in the constraint and transforms the optimization problem into a min and max problem. Similarly, [3] solves the optimization problem using Sinkhorn distance to avoid instability during training. They propose to estimate the Average Treatment Derivative (ATD) and use ATD to obtain a bound on the ATE. They also prove the consistency of the proposed estimator for linear SCMs. [41] uses a linear combination of basis functions to approximate response functions. [26] uses sieve to solve the resulting optimization problem in the IV setting. [19] proposes a neural network framework for sensitivity analysis under unobserved confounding.In Section 2, we introduce the notations and some basic concepts used throughout the paper. Next, in Section 3, we demonstrate how to construct an NCM so that they can approximate a given SCM arbitrarily well. Two kinds of architecture are given along with an approximation error analysis. In Section 4, we highlight the importance of Lipschitz regularization by giving a counterexample that is not consistent. Then, leveraging the previous approximation results, we are able to prove the consistency of this approach under regularization. Finally, we compare our method with the traditional polynomial programming method empirically in Section 4.1.-Ascoli theorem, F K L [-K, K] d V i,in , [-K, K] d V i,outare precompact set with respect to the infinity norm in space of continuous functions. And thus F Vi are compact sets.[-K, K] d U j are tight , P [-K, K] d U jare compact with respect to weak topology by Prokhorov's theorem. And the Wasserstein distance metricizes the weak topology. Thus, P [-K, K] d U j are compact and the space (M, d) is compact space. Closed set Θ ∞ ⊂ M is also compact. 
preliminary: 
 First, we introduce the definition of an SCM. Throughout the paper, we use bold symbols to represent sets of random variables.Definition 1. (Structural Causal Model) A Structural Causal Model (SCM) is a tuple M = (V , U , F, P (U ), G 0 ), where V = {V i } n V i=1 is the set of observed variables; U = {U j } n U j=1is the set of latent variables; P (U ) is the distribution of latent variables; G 0 is an acyclic directed graph whose nodes are V . The values of each observed variable V i are generated byEQUATIONwhereF = (f 1 , • • • , f n V ), Pa(V ) is the set of parents of V in graph G 0 and U Vi is the set of latent variables that affect V i . V i takes either continuous values R di or categorical in [n i ].We extend graph G 0 by adding bi-directed arrows between any V i , V j ∈ G 0 if there exists a correlated latent variable pair(U k , U l ), U k ∈ U Vi , U l ∈ U Vj .We call the extended graph G M the causal graph of M.When we write Pa(V i ), we refer to the parents of V i in the G 0 . To connect with the literature, the causal graph we define is a kind of Acyclic Directed Mixed Graph (ADMG), which is often used to represent SCMs with unobserved variables [49] . Note that we allow one latent variable to enter several nodes, which differs from the common definition. We use n U and n V to denote the number of latent variables and observable variables. Let T ⊂ V be a set of treatment variables. The goal is to estimate causal quantities under a given intervention T = t. Formally, the structural equations of the intervened model areT i = t i , ∀T i ∈ T , V i (t) = f i (Pa(V i ), U Vi ) , ∀V i / ∈ T .We denote V i (t) to be the value of V i under the intervention T = t and P M (V (t)) to be the distribution ofV (t) in M. The notion of a C 2 component [51] is defined as follows. Definition 2 (C 2 -Component). For a causal graph G, a subset C ⊂ V is C 2 -component if each pair V i , V j ∈ C is connected by a bi-directed arrow in G and C is maximal.We provide a concrete example in Appendix A to explain all these notions. We will make the following standard assumption about the independence of latent variables. Note that since we allow latent variables to enter in multiple structural equations, this is more a notational convention and not an actual assumption. Also note that under this convention a bi-directed arrow essentially represents the existence of a common latent parental variable. Assumption 1. All the latent variables in U are independent.To deal with categorical variables, we assume that latent variables that influence categorical variables contain two parts: the shared confounding that influences the propensity functions and the independent noise that generates categorical distributions. Assumption 2. The set of latent variables consists of two partsU = {U 1 , • • • , U n U } ∪ {G Vi : V i is categorical}. Precisely, if V i ∈ V is a categorical variable, the data generation pro- cess of V i satisfies V i = arg max k∈[ni] g Vi k + log (f i (Pa(V i ), U Vi )) k ∼ Categorical(f i /∥f i ∥ 1 ), where G Vi = (g Vi 1 , • • • , g Vi ni ) are i.i.d. standard Gumbel variables, U Vi ⊂ {U 1 , • • • , U n U }.This convention is without loss of generality at this point, but will be useful when introducing Lispchitz restrictions on the structural equation functions. The definition of NCMs we use is slightly different from that in [52] because we need to deal with mixed variables in our models. In (1), we usually take C to be the set of all SCMs. However, it is difficult to search over all SCMs since (1) becomes an infinite-dimensional polynomial programming problem. As an alternative, we can search over all NCMs. One quantity of common interest in causal inference is the Average Treatment Effect (ATE). Definition 4. (Average Treatment Effect). For SCM M, the ATE at T = t with respect to T = t 0 is given by ATEM (t) = E u∼P (U ) [Y (t) -Y (t 0 )].Partial identification can be formulated as estimating the solution to the optimization problems (1) [3] . The max and min values F and F define the interval [F , F ] which is the smallest interval we can derive from the observed data without additional assumptions. In particular, if F = F , then the causal quantity is point-identified.Notations. We use ∥• ∥, ∥ • ∥ ∞ for the 1-norm and ∞-norm and [n] for {1, • • • , n}. Bold letters represent sets of random variables. We let F L (K 1 , K 2 ) be the class of Lipschitz L-continuous functions f :K 1 → K 2 .We may omit the domain and use F L when the domain is clear from context. Let H(K 1 , K 2 ) be the set of homeomorphisms from K 1 to K 2 , i.e., injective and continuous maps in both directions. We define ϵ(F 1 , F 2 ) = sup f2∈F2 inf f1∈F1 ∥f 2 -f 1 ∥ for function classes F 1 , F 2 .We use standard asymptotic notation O(•), Ω(•). Given a measure µ on R d1 and a measurable function f : R d1 → R d2 , the push-forward measure f # µ is defined as f # µ(B) = µ f -1 (B) for all measurable sets B. We use P (X) to represent the distribution of random variable X.Let ∆ n = {(p 1 , • • • , p n ) : n i=1 p i = 1, p i ⩾ 0}be the probability simplex. We use Categorical(p), p ∈ ∆ n to represent categorical distribution with event probability p. We let W (•, •) be the Wasserstein-1 distance and S λ (µ, ν) be the Sinkhorn distance [12] with regularization parameter λ > 0.(i) ∩ C On(i+1) ̸ = ∅. Furthermore, C On(i) ∩ C On(i+1) is a d -1 dimensional cube.Proof of Proposition 1. Without loss of generality, we can assume that P is supported on [0, 1] d , vanishes at the boundary and P(B) ⩾ C 1 λ(B) for some constant C 1 > 0 and all measurable sets B ∈ [0, 1] d . Suppose we have proven the conclusion for this case. By Assumption 5, there existsf ∈ H([0, 1] d , K)∩F L ([0, 1] d , K) such that f -1 # P satisfies these conditions. There exists continuous function γ such that f -1 # P = γ # λ, which implies P = (f • γ) # λ. In particular, if γ is 1/d-Holder continuous, f • γ is also 1/d-Holder continuous because f is Lipschitz continuous.When d = 0, the constant map satisfies the requirement. We focus on the case d ⩾ 1 in the following. The proof is to modify the construction of the Hilbert space-filling curve. By changing the velocity of the curve, the push-forward measure of the Hilbert space-filling curve can simulate a great number of distributions. We use λ to represent the Lebesgue measure [0, 1].Proof Sketch: The construction of f is inspired by the famous Hilbert space-filling curve [28] . To illustrate the idea, let us first assume that P is absolutely continuous with respect to the Lebesgue measure λ and consider d = 2. In the k-th step, we divide the unit cube into 2 2k evenly closed cubesC n = {C 1,n , • • • , C 2 2k ,n } such that C n is a refinement of C n-1 .The construction of a standard Hilbert curve is to find a sequence of curves γ n that go through all the cubes in C n . The curve γ n has one special property. If γ n-1 (t) ∈ C k,n-1 , then γ m (t) ∈ C k,n-1 , ∀m ⩾ n. For example, in Figure 9 , the points on the curve in the lower-left cubes will stay in the lower-left cubes. Note that(γ n ) # λ(C k,n ) = λ γ -1 n (C k,n ) is the time curve γ n stays in cubes C k,n .The idea is to change the speed of the curve so that P and (γ n ) # λ agree on cubes in C n , i.e.,P(C k,n ) = λ γ -1 n (C k,n) . Since we assume that P is absolutely continuous, P(∂C k,n ) = 0 and we don't need to worry about how to divide the mass on the boundary. For example, let the green cubes in the Figure 9 to be C 0 and suppose that γ 1 starts from C 0 . We will change the speed so that γ 1 spends P(C 0 ) time in this region. In the next step, we divide C 0 into four colored cubes C 1 , • • • , C 4 on the right. We change the speed again to let the time spent in each cube equal to P(C i ). Note that this construction preserves the aforementioned property, i.e., for t ∈ [0,P(C 0 )], γ n (t) ∈ C 0 , ∀n ⩾ 1.As n → ∞, it can be proven that γ n converges uniformly to a curve γ. We can also prove that P and γ # λ agree on ∪ ∞ n=1 C n . Given that ∪ ∞ n=1 C n generate the standard Borel algebra on [0, 1] 2 , we can conclude that P = γ # λ.In the general case, P may not be absolutely continuous. As a result, the boundary may not be P-measure zero. We will need to perturb the cutting planes to ensure their boundaries are measured zero. 2 . We also define the function R that measures the ratio between maximum edge and minimum edge. For a collection of cubesPreparation. For a n-dimensional cube C = n i=1 [a i , b i ], the diameter of C is diam(C) = n i=1 (b i -a i )C = {C 1 , • • • , C m }, C k = n i=1 a k i , b k i , R(C) is defined to be R(C) = max i,k |a k i -b k i | min i,k |a k i -b k i | .R(C) measures the shape of the cubes in a collection. We need to control R(C) in the construction to obtain Holder continuity in Step 2.Step 1: Define partition. First, we construct a sequence of partitions that has the property (*) recursively. Let C 0 = [0, 1] d . Suppose we have defined close cubes collectionC n = {C i1,i2••• ,in } ij =0,••• ,2 d -1 , , such that 1.C n is obtained by cutting the unit cube [0, 1] d by planes that are paralleled to the coordinate planes.2. C i1,i2••• ,in-1 = in C i1,i2••• ,in and diam(C i1,i2••• ,in-1 ) ⩽ 2/3 • diam(C i1,i2••• ,in ).n , |g n (θ) -g n (θ ′ )| = |S λn P θ * n (V ), P θ mn (V ) -S λn P θ * n (V ), P θ ′ mn (V ) | ⩽ |W P θ * n (V ), P θ mn (V ) -W P θ * n (V ), P θ ′ mn (V ) | + 4(log(mn) + 1)λ n ⩽ W P θ mn (V ), P θ ′ mn (V ) + 4(log(mn) + 1)λ n ⩽ O(d(θ, θ ′ )) + 4(1 + log(nm n ))λ n ,where we use Lemma 7 in the second inequality and triangle inequality and Theorem 1 in the third inequality. By the condition in Theorem 4, the second term (1 + log(nm n ))λ n → 0 as n → ∞. Therefore, ( 24) is verified with τ n = 4(1 + log(nm n ))λ n .We then verify the uniform convergence of g n . By triangle inequality,|g n (θ) -g(θ)| = |S λn P θ * n (V ), P θ mn (V ) -W P θ * (V ), P θ (V ) | ⩽ |S λn P θ * n (V ), P θ mn (V ) -W P θ * n (V ), P θ mn (V ) | + |W P θ * n (V ), P θ mn (V ) -W P θ * (V ), P θ (V ) |.By Lemma 7, |S λn P θ * n (V ), P θ mn (V ) -W P θ * n (V ), P θ mn (V ) | ⩽ 2(log(mn) + 1)λ n and we get|g n (θ) -g(θ)| ⩽ 2(log(mn) + 1)λ n + |W P θ * n (V ), P θ mn (V ) -W P θ * n (V ), P θ (V ) | + |W P θ * n (V ), P θ (V ) -W P θ * (V ), P θ (V ) | ⩽ 2(log(mn) + 1)λ n + W P θ mn (V ), P θ (V ) + W (P θ * n (V ), P θ * (V ))⩽ 2(log(mn) + 1)λ n + O W (P θ * n (U ), P θ * (U )) + W P θ (U ), P θ mn (U ) ,where we use Theorem 1 in the last inequality.We first bound sup θ W P θ (U ), P θ mn (U ) using standard VC dimension argument. Note that {U 1 , • • • , U n U } are independent. By Lemma 5,EQUATIONBy the construction in Theorem 3, P θ (U i ) ∼ j∈[N C,i ] p j P (f θi,j (Z i,j )), Z i,j ∼ Unif(0, 1), where f θi,j are neural networks with constant width and depth L 1,n + L 2,n , N C,i is the number of connected components of supp(P (U i )) and θ i,j are the parameters of the neural networks. By Lemma 4,sup θ W P θ (U i ), P θ mn (U i ) ⩽ O( max j∈[N C,j ]sup θi,j W (P (f θi,j (Z i,j )), P mn (f θi,j (Z i,j ))).By [6] , the pseudo-dimension of N N (Θ(1), L) is O L2 log L . By boundness of all the neural networks and Lemma 8 2 , with probability at least1-exp O (ϵ ′ n ) -d U max log(ϵ ′ n ) -1 + (L n,1 + L n,2 ) 2 log(L n,1 + L n,2 ) log ϵ -1 n -m n ϵ 2 n ,the following event happens.sup θi,j W (P (f θi,j (Z i,j )), P mn (f θi,j (Z i,j ))⩽ ϵ n + ϵ ′ n . Let ϵ ′ n = m -1/(d U max +2) n log m n , L n,i = m d U max /(2d U max +4) n log m n , ϵ n = Cm -1/(d U max +2) n log m n with constant C > 0 sufficiently large, we get P sup θi,j W (P (f θi,j (Z i,j )), P mn (f θi,j (Z i,j )))⩽ O m -1/(d U max +2) n log m n ⩾ 1 -O m -2 n .As long as m n = Ω(n), the Borel-Cantelli lemma implies that almost surely, there exists N > 0 such that when n > N , sup θi,j W (P (f θi,j (Z i,j )), P mn (f θi,j (Z i,j )))⩽ O m -1/(d U max +2) n log m n for all i, j. Therefore, almost surely, for sufficiently large n,EQUATION[17, Theorem 2] shows that for sufficiently large n, ifδ n = C 1 n -1/ max{d U max ,2} log 2 (n), where C 1 > 0 is a constant, P (W (P θ * (U i ), P θ * n (U i )) > δ n ) ⩽        exp(-CC 1 log 2 (n)) d U max = 1, exp -CC1 log 4 (n) log 2 (2+C -1 1 n 1/2 log -2 (n)) d U max = 2, exp(-CC 1 log 2/d U max (n)) d U max > 2,where C > 0 is a constant. Let C 1 sufficiently large such thatP (W (P θ * (U i ), P θ * n (U i )) > δ n ) ⩽ O n -2 .Therefore,∞ n=1 P (W (P θ * (U i ), P θ * n (U i )) > δ n ) < ∞.By Borel-Cantelli lemma, with probability 1, there exist N > 0 such that W (P θ * (U i ), P θ * n (U i )) ⩽ δ n , ∀n > N . By Equation (29) , with probability 1, there exist N > 0 such thatEQUATIONTherefore, by Equations ( 28), (30) and (31) , almost surely, the following inequality holds eventually,sup θ∈Θ∞ |g n (θ) -g(θ)| ⩽ O log(nm n )λ n + δ n + m -1/(d U max +2) n log m n = s n ,where s n is defined in Theorem 4. Proposition 1] shows that given a L-Lipschitz continuous function f : 3 For a multivariate function with output dimension m ′ , we can approximate each coordinate individually and get a neural network approximation f that is 
approximation error of neural causal models: 
 In this section, we study the expressive power of NCMs, which serves as a key ingredient in proving the consistency result. In particular, given an SCM M * , we want to construct an NCM M such that the two causal models produce similar interventional results. Unlike in the discrete case [4, 51] , latent distributions can be extremely complicated in general cases. The main challenge is how to design the structure of NCMs to ensure strong approximation power.In the following, we first derive an upper bound on the Wasserstein distance between two causal models sharing the same causal graph. Using this result, we decompose the approximation error into two parts: the error caused by approximating structural functions via neural networks and the error of approximating the latent distributions. Then, we design different architectures for these two parts.Decomposing the Approximation Error First, we present a canonical representation of an SCM, which essentially states that we only need to consider the case where each latent variable U i corresponds to a C 2 component of G. Definition 5 (Canonical representation). A SCM M with causal graph G has canonical form if 1. The set of latent variables consists of two sets,U = U C : C is a C 2 -component of G ∪ G Vi = (g Vi 1 , • • • , g Vi ni ) : V i is categorical ,where U C and g Vi j are independent and g Vi j are standard Gumbel variables. 2. The structure equations have the formEQUATIONU Vi = {U C : V i ∈ C, C is a C 2 -component of G} and (x) k is the k-th coordinate of the vector x. We further assume that f i are normalized for categorical variables.Given a function class F, the SCM class M(G, F, U ) consists of all canonical SCM models with causal graph G such thatf i ∈ F, i ∈ [n V ].Proposition 4 in the appendix shows that any SCM satisfying Assumption 1,2 can be represented in this way and we provide an example in Appendix B.1 to illustrate how to obtain the canonical representation for a given SCM. Therefore, we restrict our attention to the class M(G, F, U ). For two SCM classes M(G, F, U ), M(G, F, Û ), we want to study how well we can represent the models in the first class by the second class. The Wasserstein distance between the intervention distributions is used to measure the quality of the approximation. To approximate the functions in the structural equations, we need to make the following regularity assumptions on the functions. Assumption 3.If V i is continuous, f i in (3) are L f -Lipschitz continuous. If V i is categorical, the propensity functions f i (Pa(V i ), U Vi ) ≜ P(V i = j|Pa(V i ), U Vi ), j ∈ [n i ] are L f -Lipschitz continuous. There exists a constant K > 0 such that max i∈[n V ],j∈[n U ] {∥V i ∥ ∞ , ∥U j ∥ ∞ } ⩽ K.The following theorem summarizes our approximation error decomposition. Theorem 1. Given any SCM model M ∈ M(G, F L , U ), let the treatment variable set be T = {T k } nt k=1 and suppose that Assumption 1, 2 and 3 hold for M with Lipschitz constant L, constant K. For any intervention T = t and M ∈ M(G, F, Û ) , we haveEQUATIONwhere C G (L, K) is a constant that only depends on L, K and the causal graph G and f i , fi are structural functions of M and M respectively.Theorem 1 separates the approximation error into two parts, which motivates us to construct the NCM in the following way. First, we approximate the functions f i in (3) by NNs fi . Then, we approximate the distribution of latent variables by pushing forward uniform and Gumbel variables using neural networks, i.e., ÛCj = ĝj (Z Cj ), where {C j } are C 2 components and Z Cj are multi-variate uniform and Gumbel random variables. The structural equations of the resulting approximated model M areEQUATIONwheres N C,j are constants to be specified later.For the first part, we need to approximate Lipschitz continuous functions. For simplicity, we assume that the domain of the functions are uniform cubes. Similar arguments hold for any bounded cubes. We denote N N k1,k2 (W, L) to be the set of ReLu NNs with input dimension k 1 , output dimension k 2 , width W and depth L. It has been shown that ϵ(N N k1,1 (2dEQUATION) [53] , where ϵ(•, •) denotes the approximation error defined in Section 2. For a vector valued function, we can use a wider NN to approximate each coordinate and get a similar rate.For the second part, we approximate each U i individually by pushing forward i.i.d. multivariate uniform and Gumbel variables ÛCi = ĝi (Z Ci ) since the latent variables are independent by Assumption 1. To do so, we examine under what assumptions on the measure P over R n we can find a NN ĝ such that W (ĝ # λ, P) is small, where λ(•) is some reference measure.In this subsection, we will extend the results in [43] to construct a wide neural network as the push-forward map. It turns out that to get a good approximation of the targeted distribution, the shape of the support is essential.Assumption 4 (Mixed Distribution). The support of measure P has finite connected componentsC 1 , C 2 • • • C N C , i.e., supp(P) = N C i=1 C i , and each component C i satisfies H([0, 1] d C i , C i ) ∩ F L ([0, 1] d C i , C i ) ̸ = ∅ for some d C i ⩾ 0. Recall that H(K 1 , K 2 )is the set of homeomorphisms from K 1 to K 2 defined at the end of Section 2.Assumption 4 encompasses almost all natural distributions. For example, distributions supported on [0, 1] d and closed balls, finitely supported distributions and mixtures of them all satisfy this assumption. Assumption 4 allows us to transform the support of the targeted distribution into unit cubes and the nice geometric properties of unit cubes facilitate our construction. Now, we briefly explain the construction of the push-forward maps. An example is provided in Appendix B.1 to illustrate the construction. The NN architecture consists of three parts (see Figure 1a ). The input dimension is the same as the number of connected components of the supportN C . For each component C i , let H i ∈ H([0, 1] d C i , C i ) ∩ F L ([0, 1] d C i , C i ), where d C i is the dimension of component C i in Assumption 4. By P = (H i ) # (H -1i ) # P on C i , we can approximate (H -1 i ) # P first, which is supported on a unit cube. [43] constructs a wide NN ĝ of width W and constant depth such that W (ĝ # λ,(H i ) -1 # P) ⩽ O(W -1/d C i )where λ is the uniform measure on [0, 1]. Then, we approximate the Lipschitz map H i to C i to pull the distribution back to C i . These are the first two parts (yellow and blue blocks in Figure 1a ) of the architecture. Suppose that the output of i-th coordinate in the first two parts is ⃗ v i , the Gumbel-softmax layer in the third part (green box) combines different components of the support. In particular, we want to output ⃗ v i with probabilityp i = P(C i ). Let V = [⃗ v 1 , • • • , ⃗ v i ],this can be achieved by outputting V X, whereX = (X 1 , • • • , X N C) T is a one-hot random vector with P(X i = 1) = p i . To use backpropagation in training, we use the Gumbel-Softmax trick [31] to (approximately) simulate such a random vector,Xτ i = exp((log pi+Gi)/τ ) N C k=1 exp((log p k +G k )/τ ), where τ > 0 is the temperature parameter (a hyperparameter) and G i ∼ Gumbel(0, 1) are i.i.d. standard Gumbel variables. As τ → 0, the distribution of Xτ converges almost surely to the categorical distribution [36, Proposition 1] . In particular, when τ = 0, we denote X0 i = I i=arg maxj {log pi+Gi} . The output of the last layer is V Xτ . Note that the Gumbel-softmax function is differentiable with respect to parameter {log(p i )} i=1,••• ,N C . Therefore, we can train the network with common gradient-based algorithms. Putting things together, we obtain the following theorem. Theorem 2. Given any probability measure P on R d that satisfies Assumption 4, let λ be the Lebesgue measure on [0, 1] N C , where N C is defined in Assumption 4. There exists a neural network ĝ = ĝτ 3 • ĝ2 • ĝ1 with the above architectsure (Figure 1a ) such thatW (ĝ # λ, P) ⩽ O W -1/ maxi{d C i } 1 + L -2/ maxi{d C i } 2 + (τ -τ log τ ) .Here, ĝi , i = 1, 2 has the separable form ĝ1i (x 1 ), • • • , ĝN C i (x N C ) and ĝj 1 ∈ N N 1,d C j (W 1 , Θ(d C j )), ĝj 2 ∈ N N d C j ,d (Θ(d • d C j ), L 2 ), j ∈ [N C ]. {d C j }are the dimension of cubes in Assumption 4. ĝτ 3 is the Gumbel-Softmax layer with temperature parameter τ > 0.3 (the Gumbel-softmax layer) is actually a random function since the coefficient vector Xτ is a random variable. In this sense, ĝ# λ can be viewed as pushing forward uniform and Gumbel variables using a neural net.In this subsection, we will show that under one additional assumption on the distribution, deep ReLu networks have a stronger approximation power in approximating distributions, which means we can use fewer computational units to achieve the same worst-case theoretical approximation error. Assumption 5 (Lower Bound). Suppose that P is supported on a compact set K ⊂ R D , there exists a constantC f > 0, f ∈ H([0, 1] d , K) ∩ F L ([0, 1] d , K), such that for any measurable set B ⊂ [0, 1] d , P (f (B)) ⩾ C f λ(B). Besides, if d > 0, f # P vanishes on the boundary ∂[0, 1] d .Assumption 5 implies that dλ/d(f -1 # P) exists and is lowered bounded by a constant C f . The next proposition extends the Skorohod representation theorem [48] . It shows that under Assumption 5, it is possible to simulate any distribution on the unit cubes with Hölder continuous curves and uniform distribution on [0, 1]. Proposition 1. Let λ be the Lebesgue measure on [0, 1]. Given any probability measure P that satisfies Assumption 5, there exists a continuous curve γ :[0, 1] → supp(P) such that γ # λ = P. Furthermore, if d ⩾ 1, γ is 1/d-Hölder continuous.Results from [53] show that we can approximate any Hölder continuous d-dimensional function with index α by a deep ReLu network with depth L and error O(L -2α/d ). Leveraging this result, we can replace the first part of the architecture in the previous subsection with deep ReLu networks (See Figure 8 in the appendix). The remaining two parts are the same as the construction in Figure 1a . The following theorem gives a sharper bound on the approximation error compared with Theorem 2. Theorem 3. Under the Assumption 4, if in addition, P constrained to each component satisfies Assumption 5, there exists a neural network ĝ = ĝτ 3 • ĝ2 • ĝ1 with the above architecture such thatW (ĝ # λ, P) ⩽ O L -2/ maxi{d C i } 1 + L -2/ maxi{d C i } 2 + (τ -τ log τ ) ,where d C i are the dimensions of the connected components in Assumption 4, ĝi , i = 1, 2 has the formĝ1 i (x 1 ), • • • , ĝN C i (x N C ) and ĝj 1 ∈ N N 1,d C j (Θ(d C j ), L 1 ). ĝj 2 , ĝτ 3 , τ are the same as in Theorem 2.Let N be the number of nonzero weights in a neural network, Theorem 3 shows that a deep NN can achieve O(N -2/d ) error while by Theorem 2 a wide network can only achieve O(N -1/d ) error. Now, we can put things together to construct NCM approximations. For simplicity, we omit the input and output dimensions of the neural network. As we mention in previous sections, our construction (5) consists of two parts, fi approximating the structural functions f i in (3), andĝj (Z j ) = ĝτ 3,j • ĝ2,j • ĝ1,j (Z j ) approximating the latent variables U j . Let NCM G (F 0 , F 1 , F 2 , τ) be a collection of NCMs with structural equation ( 5) andfi ∈ F 0 , ĝ1,j = (ĝ 1 1,j , • • • , ĝN C,j 1,j ), ĝi 1,j ∈ F 1 , ĝ2,j = (ĝ 1 2,j , • • • , ĝN C,j 2,j ), ĝi 2,j ∈ F 2 ,where N C,j is the number of connected components for U j and F i are function classes. Corollary 1. Given a causal model M * ∈ M(G, F L , U ), suppose that Assumptions 1-3 hold and the distributions of U C for all C 2 -component satisfy the assumptions of Theorem 3. Let d in max and d out max be the largest input and output dimension of f i in (3) and d U max be the largest dimension of all latent variables. There exists a neural causal modelM ∈ NCM G (N N (Θ(d in max d out max ), L 0 ), N N Θ d U max , L 1 , N N Θ (d U max ) 2 , L 2 , τ )with structure equations (5) . For any intervention T = t, M satisfiesW P M * (V (t)), P M(V (t)) ⩽ O(L -2/d in max 0 + L -2/d U max 1 + L -2/d U max 2 + (τ -τ log τ )).Similar approximation results also hold for wide NN approximation, as presented in Section 3.1. The proof can be easily adapted to the wide NNs architecture.i1,i2••• ,in ) = 0 and R(C n+1 ) ⩽ (1 + 2 -n ) R(C n ).Next, we construct C n+1 from C n that preserves these properties. We refine the divisionC n = {C i1,i2••• ,in } ij =0,••• ,2 d -1 by union of hyperplane P i,r = {y : y i = r} to get {C i1,i2••• ,in+1 } such that C i1,i2••• ,in = in+1 C i1,i2••• ,in+1 and diam(C i1,i2••• ,in-1 ) ⩽ 2/3 • diam(C i1,i2••• ,in ). By property 1, C n is obtained by dividing [0, 1] d using P 1,r 1 0 , • • • , P 1,r 1 2 n -1 , • • • , P i,r i j , • • • , P d,r d 2 n , r i 0 = 0, r i 2 n = 1,we can further cut the unit cube using hyperplane {P i,(r ij +r i j+1 )/2) } i=1,•••d,j=1,•••2 n -1 .However, this construction may not satisfy property 3, P(∂C i1,i2••• ,in+1 ) = 0. We need to perturb each hyperplane to ensure a measure-zero boundary. SinceZ i = r ∈ R : P x ∈ [0, 1] d , [x] i = r ̸ = 0is at most countable, we can perturb the hyperplanes a little bit, i.e., {P i,(r ij +r i j+1 )/2+ϵi,j) } i=1,•••d,j=1,•••2 n -1 to ensure property 3.In this way, we can choose |ϵ i,j | to be sufficiently small such thatdiam(C i1,i2••• ,in+1 ) ⩽ 2/3 • diam(C i1,i2••• ,in ),and max i,j r i j+1 -r i j /2 + |ϵ i,j | min i,j r i j+1 -r i j /2 -|ϵ i,j | ⩽ max i,j r i j+1 -r i j min i,j r i j+1 -r i j • 1 + 2 -n . Therefore, R(C n+1 ) ⩽ (1 + 2 -n ) R(C n )and it is easy to see that properties 1,2,3 are satisfied.Step 2: Construct Hilbert Curve. By Lemma 1, there exist an sequence of ordering{O n } of {C n } that satisfies • The first n digits of O n+1 (i) with O n (j), (j -1)• 2 d < i ⩽ j• 2 d . • C On+1(i) ∩ C On+1(i+1) are d -1 dimensional cubes.Let I 0 = {[0, 1]}, we define I n recursively. Suppose we have define interval collectionI n = {{I i1,i2••• ,in } ij =0,••• ,2 d -1 } such that I i1,i2••• ,in-1 = in I i1,i2••• ,in and P(C i1,i2••• ,in ) = |I i1,i2••• ,in |. Since P(∂C i1,i2••• ,in+1 ) = 0, we have P(C i1,i2••• ,in ) = in+1 P(C i1,i2••• ,in+1). With the ordering, we divide eachI On(j) , 0 ⩽ j ⩽ 2 nd -1 into 2 d closed sub-interval I On+1(i) , j2 d ⩽ i < (j + 1)2 d such that I On+1(i) = [ i-1 k=1 p k , i k=1 p k ] ,where p i = P(C On+1(i) ). This construction satisfies |I On+1(i) | = P(C On+1(i) ). Note that by condition 3 in step 1, P vanishes at the boundary of the cubes. Therefore,EQUATION)Now, we can construct a piecewise linear function γ n+1 (t) such thatEQUATION)The idea is to construct a piecewise linear curve going through all cubes in C n+1 exactly once and modify its speed. TakeEQUATION)The first two curves γ 1 , γ 2 are shown in Figure 9 . It is straightforward to verify that γ n+1 (i k=1 p k ) = v i and γ n+1 (I On+1(i) ) = γ n+1 ([ i-1 k=1 p k , i k=1 p k ]) ⊂ C On+1(i) .In fact, since v i-1 , v i are on different surfaces of the cube C On+1(i) , the line segment between v i-1 , v i lies inside the interior of C On+1(i) and γ n+1 goes through all the cubes C On+1(i) exactly once. We haveγ -1 n+1 (Int(C On+1(i) )) ⊂ [ i-1 k=1 p k , i k=1 p k ] and (γ n+1 ) # λ(C On+1(i) ) = P(C On+1(i) ), where Int(•) is the interior of a set. Besides, we also have(γ n+1 ) # λ(∂C On+1(i) ) ⩽ λ(γ -1 n+1 ({v j } j=1•••2 d(n+1) )) = 0. Finally, for any k 1 ⩾ k 2 ⩾ n, t ∈ [0, 1], by property 2 in step 1, ( 20) and ( 22), γ k1 (t),γ k2 (t) are in one cube C i1,i2••• ,in . and diam(C i1,i2••• ,in ) ⩽ 2 3 n √ d. Thus, |γ k1 (t) -γ k2 (t)| ⩽ 2 3n √ d which implies {γ n (t)} converges uniformly to one continuous function γ(t).P = γ # λ. By W ((γ k ) # λ, γ # λ) ⩽ ∥γ k -γ∥ ∞ → 0, (γ k ) # λ converges weakly to γ # λ. By construction, P(∂C i1,i2••• ,in ) = 0 and (γ n ) # λ(C i1,i2••• ,in ) = P(C i1,i2••• ,in ).Therefore, by condition 2 in step 1, for any k ⩾ n,EQUATIONwhere we useP(∂C i1,i2••• ,in ) = 0 and (γ k ) # λ(∂C i1,i2••• ,in ) = 0 in the second and the last equation. We claim that γ # λ(∂C i1,i2••• ,in ) = 0. For k > n, we define B k = ∂C i1,i2••• ,i k ∩ ∂C i1,i2••• ,in ̸ = ∅ C i1,i2••• ,i k . Let k > k 1 ⩾ n and B ϵ n = {x : d(x, ∂C i1,i2••• ,in ) < ϵ}, we have (γ k ) # λ (Int(B k1 )) ⩽ (γ k ) # λ(B k1 ) = P(B k1 ) ⩽ P B (2/3) k 1 √ d n , ∀k > k 1 ,where we use (23) and the fact that(γ k ) # λ vanishes on ∂C i1,i2••• ,i k in the first equality and diam(C i1,i2••• ,i k 1 ) ⩽ 2 3 k1 √ d in the last inequality. Let k → ∞, by Portmanteau Theorem, we get γ # λ(∂C i1,i2••• ,in ) ⩽ γ # λ (Int(B k1 )) ⩽ lim inf k→∞ (γ k ) # λ(B k1 ) ⩽ lim inf k→∞ P B (2/3) k 1 √ d n = P B (2/3) k 1 √ d n . Since k 1 is arbitrary, let k 1 → ∞, γ # λ(∂C i1,i2••• ,in ) ⩽ lim k1→∞ P B (2/3) k 1 √ d n = P(∂C i1,i2••• ,in ) = 0 and we conclude that γ # λ(∂C i1,i2••• ,in ) = 0. By Portmanteau Theorem, let k → ∞ in (23), we obtain γ # λ(C i1,i2••• ,in ) = P(C i1,i2••• ,in ). Let C * = ∞ i=1 C i . Notice that for any open set U ∈ [0, 1] d , U = C⊂U,C∈C * C,which means the σ-algebra generated by C * contains all Borel sets. Besides, C * is a π-system. Hence, λ γ -1 (U ) = P(U ) for all Lebesgue measurable set U .Step 4: Holder Continuity of γ. We verify the condition of [47, Theorem 1] . I n and C n in the above construction correspond to developmentsα k , β k in [47, Theorem 1]. Define f k (I O k (i) ) = C O k (i) , i = 1, • • • , 2 kd as a map from I n to C n . 1. By construction, if I 1 ∈ I k , I 2 ∈ I k-1 , I 1 ⊂ I 2 , f k (I 1 ) ⊂ f k-1 (I 2 ). 2. If I 1 , I 2 ∈ I k and I 1 ∩ I 2 ̸ = ∅, I 1 , I 2 are adjacent. By Lemma 1, f k (I 1 ), f k (I 2 ) share a (d -1)-dimensional boundary. Therefore, f k (I 1 ) ∩ f k (I 2 ) ̸ = ∅. 3. We have |C k-1 | ⩽ 2 d |C k |.4. By property 3,R(C n ) ⩽ n-1 k=1 1 + 2 -k < ∞ k=1 1 + 2 -k < ∞. Let A = ∞ k=1 1 + 2 -k .We have for any C ∈ C n , by Assumption 5diam d (C) ⩽ R(C n )λ(C) ⩽ AC -1 1 P(C) ⩽ AC -1 1 max I∈I (diam(I)) ,where C 1 is the constant in Assumption 5.[0, 1] m → R, there exist f ∈ N N m,1 (W, Θ(log(m)) such that ∥f -f ∥ ∞ ⩽ O W -1/d and f is √ mL- Lipschitz continuous.√ mm ′ L-Lipschitz continuous and ∥f -f ∥ ∞ ⩽ O W -1/m .Next, we define the truncate operator T K f = min{K, max{-K, f }}, which is a contraction mapping, and prove that applying the truncate operator will not increase approximation error. If |f | ⩽ K, we have∥f -T K f ∥ ∞ = ∥T K f -T K f ∥ ∞ ⩽ ∥f -f ∥ ∞ ⩽ O W -1/m .3 The original proof does not specify the depth of the neural network. The authors show that the network architectures can be chosen as consisting of O(W ) parallel blocks each having the same architecture. Each block is used to realize the function ϕ(x) = (min(mink̸ =s (x k -xs + 1), min k (1 + x k ), min k (1 -x k )))+, x ∈ R m .This function can be realized by feed-forward network with width O m 2 and depth O(log m).Therefore, we get thatEQUATIONTheorem 3 implies that for each latent variable U i , there exist a neural network g i with architecture in Corollary 1 such thatW (P (U i ), P (g i (Z i ))) ⩽ O( n i=1 L -2/dmax i,n + τ n (1 -log τ n )).By Assumption 3, we know ∥U i ∥ ∞ ⩽ K, which impliesEQUATIONCombining Equations ( 32) and (33) with the same proof as Corollary 1, it can be proven in the same way as Corollary 1 that there exists a θ n ∈ Θ n satisfyingd(θ * , θ n ) ⩽ ϵ n = O W -1/d in max 0,n + L -2/d U max 1,n + L -2/d U max 2,n + τ n (1 -log τ n ) .Let Θ∞ = {θ * }, we have verified the third assumption in Proposition 6.Take the Wasserstein radius to be α n = O(s n + ϵ n ), Proposition 6 implies the conclusion.C.3 Proof of Proposition 3 Lemma 2. Let ( T , Ŷ ) ∼ µ, (T, Y ) ∼ ν and suppose that f (t) = E ν [Y |T ], f (t) = E µ [ Ŷ | T ] are L-Lipschitz continuous and |f (t)| ⩽ K, | f (t)| ⩽ K, then we have (f (t) -f (t)) 2 dν(dt) ⩽ C W W (µ, ν),where C W = 4LK + 2K max{L, 1}.Proof of Lemma 2. By the duality formulation of Wasserstein-1 distance, we haveW (µ, ν) = sup g∈Lip(1) E µ [g( T , Ŷ )] -E ν [g(T, Y )].Let g 0 (t, y) = ( f (t)-f (t))y, we verify g 0 is a Lipschitz continuous function in {(t, y) : ∥(t, y)∥ ∞ ⩽ K}.|g 0 (t 1 , y 1 ) -g 0 (t 2 , y 2 )| ⩽ |g 0 (t 1 , y 1 ) -g 0 (t 1 , y 2 )| + |g 0 (t 1 , y 2 ) -g 0 (t 2 , y 2 )| = |( f (t 1 ) -f (t 1 ))(y 1 -y 2 )| + |y 2 ( f (t 1 ) -f (t 1 ) -( f (t 2 ) -f (t 2 )))| ⩽ 2K|y 1 -y 2 | + K(| f (t 1 ) -f (t 2 )| + |f (t 1 ) -f (t 2 )|) ⩽ 2K|y 1 -y 2 | + 2KL|t 1 -t 2 |.Let L g = 2K max{L, 1}, we have proven that g 0 is L g -Lipschitz continuous in {(t, y) : ∥(t, y)∥ ∞ ⩽ K}. Thus,EQUATIONNow, let h(t) = ( f (t) -f (t)) f (t). Following the same argument, it can be proven that h(t) is Lipschitz continuous with Lipschitz constant being L h = 4LK. Hence,E µ [h( T )] ⩾ E ν [h(T )] -L h W (µ, ν).Plug into (34) , and we get(L g + L h )W (µ, ν) ⩾ E ν [( f (T ) -f (T )) f (T ) -( f (T ) -f (T ))f (T )] = E ν ( f (T ) -f (T )) 2 . Proof of Proposition 3. If Y is not descendant of T , E[Y (t)] = E[Y ] and we have t2 t1 (E M [Y (t)] -E M[ Ŷ (t)]) 2 dt = t2 t1 (E M [Y ] -E M[ Ŷ ]) 2 dt = (t 2 -t 1 )(E M [Y ] -E M[ Ŷ ]) 2 ⩽ (t 2 -t 1 )W (P M (V ), P M( V )) 2 . Now, suppose that Y is a descendant of T . Let X = Pa(T ). Note that f y (x, t) = E M [Y |X = x, T = t] is L y -Lipschitz continuous with Lipschitz constant L y ⩽ (L + 1) n V. This is because Y can be written as Y = F 0 (X, T, U y ) where U y are latent variables independent of T, X and F 0 is composition of f i in structure equations. The composition of Lipschitz functions is still Lipschitz and F 0 is a composition of at most n V L-Lipschitz functions. F 0 is (L + 1) n V -Lipschitz and so isf y (x, t) = E U y [F 0 (X, T, U y )].Recall that ν, µ are the observation distributions of M and M respectively in Proposition 3. Similarly, fy (x, t)= E M[ Ŷ |X = x, T = t] is L y -Lipschitz continuous.By Lemma 2 and the overlap assumption, we have Therefore,C W W (µ, ν) ⩾ E ν (f y (X, T ) -fy (X, T )) 2 = (f y (x,f y (x, t)ν(dx) - fy (x, t)ν(dx) 2 ⩾ 1 2 f y (x, t)ν(dx) - fy (x, t)µ(dx) 2 - fy (x, t)µ(dx) - fy (x, t)ν(dx) 2 ⩾ 1 2 f y (x, t)ν(dx) - fy (x, t)µ(dx) 2 -L 2 y W 2 (µ, ν) = 1 2 (E M [Y (t)] -E M[ Ŷ (t)]) 2 -L 2 y W 2 (µ, ν).Combine all the equations, we get t2 t1(E M [Y (t)] -E M[ Ŷ (t)]) 2 P (dt) ⩽ 2C W δ W (µ, ν) + 2L 2 y W 2 (µ, ν)(t 2 -t 1 ).Proof of Corollary 2. In the proof of Theorem 4, we know with probability at least 1 -O(n -2 ), (6) has feasible solutions andW (P M * n (V ), P M * (V )) ⩽ O(α n ), W (P θn mn (V ), P θn (V )) ⩽ O(α n ),where notations θ and P θ are defined in the proof of Theorem 4 and θ n is one of the minimizers of (6) . By Lemma 7, we know thatW (P M * n (V ), P θn mn (V )) ⩽ S λn (P M * n (V ), P θn mn (V )) + 2(log(m n n) + 1)λ n ⩽ O(α n ).We get thatW (P M * (V ), P θn (V )) ⩽ W (P M * n (V ), P M * (V ))+W (P M * (V ), P θn (V ))+W (P θn mn (V ), P θn (V )) ⩽ O(α n ).By Proposition 3, we get|F n -F * | ⩽ O( √ α n ).The structure equations of the generative models are (5). We use three-layer feed-forward neural networks with width 128 for each fi and six-layer neural networks with width 128 for each ĝj . We use the Augmented Lagrangian Multiplier (ALM) method to solve the optimization problems as in [3] . We run 600 epochs and use a batch size of 2048 in each epoch. m n is set to be m n = n.The "geomloss" package [16] is used to calculate the Sinkhorn distance. To impose Lipschitz regularization, we use the technique from [23] to do layer-wise normalization to the weight matrices with respect to infinity norm. The upper bound of the Lipschitz constant in each layer is set to be 8.The τ in the Gumbel-softmax layer is set to be 0.For the choice of Wasserstein ball radius α n , we use the subsampling technique from [11, Section 3.4] . We can take the criterion function in [11] to be Q(θ) = W (P M * (V ), P M θ (V )) and its empirical estimation to be Q n (θ) = W (P M * n (V ), P M θ n (V )). In the first step, we minimize the Wasserstein distance θ * n = arg min θ W (P M * n (V ), P M θ n (V )). Then, [11] propose to refine the radius Q n ( θ * n ) by taking the quantile of subsample {sup θ:Qn( θ)⩽βQn( θ * n ) Q j,b (θ)}, where Q j,b denotes the criterion function estimated at j-th subsample of size b. However, it is time-consuming to solve this optimization problem many times. In practice, we set the radius α n to be the 95% quantile of {W ( P M * j }(V ), P M We consider the noncompliance binary IV example in [14, Section D.1] . The causal graph is shown in Figure 10 . We could see from 2 : The sample size is taken to be 5000. We average over 10 runs with different random seeds. SD is short for the standard derivation.Now, we turn to the continuous IV setting, where T is still binary, but Y and Z can be continuous.Let E λ i ∼ λZ 1 + (1 -λ)Unif(-1, 1), where Z 1 is the Gaussian variable conditioning on [-1, 1] and the structure equations of M λ to beE Y = E λ 1 , U = E λ 2 , Z = E λ 3 , P (T = 1|Z) = (1.5 + Z + 0.5U )/3, Y = 3T -1.5T U + U + E Y ,where E λ i are independent. It is easy to see the ATE of this model is 3 regardless of λ. In the experiment, we randomly choose ten λ ∼ Unif(0, 1). For each λ, we run the algorithms 5 times to get the bound. We choose different latent distributions (indexed by λ) in the experiments to create more difficulty.Since Autobounds can only deal with discrete variables, we discretize Z and Y . Suppose that Z ∈ [l, u], we map all points in intervals [l + i(u -l)/k, l + (i + 1)(u -l)/k], I = 0, 1 • • • , k -1 to the middle point l + (i + 1/2)(u -l)/k. We choose k = 8 in the following experiments. The choice will give rise to polynomial programming problems with 2 14 decision variables, which is quite large. 3 : We take a sample size of 5000. We randomly choose 10 λ ∼ Unif(0, 1) and get 10 models M λi . For each model M λi , we run the two algorithms for 5 times. The success rate is the number of times when the obtained bounds cover the true ATE divided by the total number of experiments. SD is short for the standard derivation.We test our neural causal method on the counterexample in Appendix C.1. We choose the noise U y to be the normal variable. The structure equations areP (U 1 = 1) = 3/5, P (U 1 = 0) = 2/5, P (U 2 = 1) = 1/3, P (U 1 = 0) = 2/3, X = -1 , w.p. 1/2, 0 , w.p. 1/2, P (T = 1|X = x) = 1/2 , x = -1, 5/8 , x = 0, Y =    T + U y , X = -1, U 1 + U y , X = 0, T = 1, U 2 + U y , X = 0, T = 0,We compare the confidence intervals of the unregularized neural casual method and the Lipschitz regularized one under different architectures. We choose the layerwise Lispchitz constants upper bound to be 5 and 1.2. As a benchmark, we also include the bound obtained by the Double Machine Learning estimator using the EconML package. The result is shown in Appendix D.3. For this identifiable case example, the double ML estimator produces better bounds for the ATE. The intervals given by regularized and unregularized NCM are similar to the regularized one slightly better in the left figure, where we use medium-sized NNs. However, in the right figure, the obtained intervals after regularization are tighter, although slightly biased, compared with the regularized setting. From these two experiments, we conclude that the architecture of NNs will also influence the results and adding regularization during the training process can prevent extreme confidence intervals or inconsistency. Lemma 3. Let µ, μ be two measures on compact set K ⊂ R d1 , for any measurable functions F, F :K → R d2 , if F is L-Lipschitz continuous, then W (F # µ, F# μ) ⩽ LW (µ, μ) + ∥F 1 -F 2 ∥ ∞ .Proof. For any 1-Lipschitz function g : R d2 → R,E X∼F # µ [g(X)] -E X∼ F# μ[g(X )] = E X∼µ [g • F (X)] -E X∼μ [g • F (X)] = (E X∼µ [g • F (X)] -E X∼μ [g • F (X)]) + (E X∼μ [g • F (X)] -E X∼μ [g • F (X)]).Note that for any x, y ∈ R d1 ,EQUATIONFor the second term,EQUATIONCombine ( 35) and ( 36), we have the conclusion.Lemma 4. Let µ 1 , • • • , µ n , μ1 , • • • , μn be measure on R d , for any p 1 , • • • , p n ∈ [0, 1] such that n k=1 p k = 1, we have W n k=1 p k µ k , n k=1 p k μk ⩽ n k=1 p k W (µ k , μk ).Proof. For any 1-Lipschitz function f : R d → R, we haveE X∼ n k=1 p k µ k [f (X)] -E X∼ n k=1 p k μk [f (X)] = n k=1 p k (E X k ∼µ k [f (X k )] -E X k ∼μ k [f (X k )]) ⩽ n k=1 p k W (µ k , μk ).By dual formulation of Wasserstein-1 distance, Proof. For any ϵ > 0, let π k be a coupling of µ k , μk such thatW n k=1 p k µ k , n k=1 p k µ k = sup f 1-Lipschitz E X∼ n k=1 p k µ k [f (X)] -E X∼ n k=1 p k μk [f (X)]E X k ,Y k ∼π k [∥X k -Y k ∥ 1 ] ⩽ W (µ k , μk ) + ϵ.Then,W (µ, μ) ⩽ E π1⊗•••⊗πn [∥X -Y ∥] = n k=1 E X k ,Y k ∼π k [∥X k -Y k ∥ 1 ]⩽ W (µ k , μk ) + ϵ.Let ϵ → 0, we get the result.Lemma 6 (Approximation error of Gumbel-Softmax layer). Given 1 > τ > 0 and p 1 ,EQUATION)where G k ∼ exp(-x -exp(-x)) are i.i.d. standard Gumbel variables. Let µ 0 be the distribution of X 0 = lim τ →0 X τ = X 0 1 , • • • , X 0 n . Then, W µ τ , µ 0 ⩽ 2(n -1)τ -2n(n -1)e -1 τ log τ.Proof. Without loss of generality, we may assume that n k=1 p k = 1. Otherwise, we can divide the denominator and numerator in (37) by n k=1 p k . Let ∆ n = {(x 1 , • • • , x n ) : n k=1 x k = 1}. We construct a transportation map T : ∆ n → ∆ n as follows. T (x 1 , • • • , x n ) is a n-dimensional vector with all zeros except that the i = arg max j x j -th entry being one. If there are multiple coordinates that is maximum, we choose the first coordinate and let the remaining coordinate be zero. It is easy to verify that T # µ τ and µ 0 have the same distribution [36] . For any δ > 0, we have E X τ ∼µ τ [∥T (X τ ) -X τ ∥ 1 ] ⩽ P (E δ )E X τ ∼µ τ [∥T (X τ ) -X τ ∥ 1 |E δ ] + P (E c δ ) E X τ ∼µ τ [∥T (X τ ) -X τ ∥ 1 |E c δ ] . Note that ∥T (X τ ) -X τ ∥ 1 ⩽ 2.For the first term, we haveP (E δ )E X τ ∼µ τ [∥T (X τ ) -X τ ∥ 1 |E δ ] ⩽ 2δn(n -1)e -1 .For the second term, under event E c δ , let k max = arg max i X i , if j = k max = arg max i (log p k +G k ) exp((log p j + G j )/τ ) Therefore, E X τ ∼µ τ [∥T (X τ ) -X τ ∥ 1 |E c δ ] ⩽ 2(n -1) exp(-δ/τ ). We get an upper bound of the transportation cost, W µ τ , µ 0 ⩽ E X τ ∼µ τ [∥T (X τ ) -X τ ∥ 1 ] ⩽ 2δn(n -1)e -1 + 2(n -1) exp(-δ/τ ).Let δ = -τ log τ , we get W µ τ , µ 0 ⩽ 2(n -1)τ -2n(n -1)e -1 τ log τ. + d F log ϵ -1 -nϵ 2 , sup f ∈F W (f # µ, f # µ n ) ⩽ δ + ϵfor all δ, ϵ > 0, where µ n is the empirical distribution of µ.Proof. By the dual formulation of the Wasserstein distance, sup where N (ϵ, S, ∥• ∥ 1 ) is the covering number of set S in ℓ 1 norm. Obviously, if h is 1-Lipschitz, N 1 (ϵ, h • F, n) ⩽ N 1 (ϵ, F, n). By Theorem 18.4 in [2] , for any fixed h,f ∈F W (f # µ, f # µ n ) = sup h∈F1(R d ,R d ),h(0)=0 sup f ∈F E X∼µn [h • f (X)] -E X∼µ [h • f (X)].N 1 (ϵ, h • F, n) ⩽ N 1 (ϵ, F, n) ⩽ e(d F + 1) 2e ϵ d F.By Theorem 17.1 in [2] ,P sup f ∈F E X∼µn [h • f (X)] -E X∼µ [h • f (X)] > ϵ ⩽ exp O d F log ϵ -1 -nϵ 2 .Let H δ be the δ-net of the setH = h : h ∈ F 1 [-K, K] d , [-K, K] d , h(0) = 0 in ℓ ∞ norm. By Lemma 6 in [21], |H δ | ⩽ exp O δ -d log δ -1 .Therefore, with probability no more than exp O δ -d log δ-1 + d F log ϵ -1 -nϵ 2 sup h∈H δ ,f ∈F E X∼µn [h • f (X)] -E X∼µ [h • f (X)] > ϵ.Notice that suph∈H,f ∈F E X∼µn [h • f (X)] -E X∼µ [h • f (X)] ⩽ 2δ + sup h∈H δ ,f ∈F E X∼µn [h • f (X)] -E X∼µ [h • f (X)],which implies that with probability no more than exp O δ -d log δ-1 + d F log ϵ -1 -nϵ 2 sup h∈H,f ∈F E X∼µn [h • f (X)] -E X∼µ [h • f (X)] > 2δ + ϵ.generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
consistency of neural causal partial identification: 
 In this section, we prove the consistency of the max/min optimization approach to partial identification.In the finite sample setting, we consider the following estimator F n of the optimal values of (1).EQUATIONs.t. S λn (P M mn (V ), P M * n (V )) ⩽ α n ,where P M * n , P M mn are the empirical distribution of P M * , P M with sample size n, m n , µ T is some given measure and F i,n will be specified later. For example, the counterfactual outcome E[Y (1)] is a special case of the objective. Our results can be easily generalized to any linear combination of objective functions of this form. We use the Sinkhorn distance because it can be computed efficiently in practice [15] . We want to study if F n gives a useful lower bound as n → ∞.To match the observational distribution, we need to increase the width or depth of the NNs we use. As the sample size increases, the number of parameters also increases to infinity, which creates difficulty in the analysis. To obtain consistency, we need to regularize the functions while preserving their approximation power. Surprisingly, if we do not use any regularization, the following proposition implies that consistency may not hold even if the SCM is identifiable. Proposition 2 (Informal, see Proposition 5 for a formal version). There exists a constant c > 0 and an identifiable SCM M * satisfying Assumptions 1-5 such that for any ϵ > 0, there exists an SCMM ϵ satisfying W (P M * (V ), P Mϵ (V )) ≤ ϵ and |ATE M * -ATE Mϵ | > c.Here, M * is the ground-truth model and M ϵ are the models we use to approximate M * . Proposition 2 implies that we need some regularization on M ϵ . Otherwise, even if the observation distributions are close, their ATEs can be far away. In particular, we may want to regularize the Lipschitz constant of the NN. Much work has been done to impose Lipschitz regularization during the training process [13, 50, 22, 42, 10] . We denote Lip(f ) to be the Lipschitz constant of a function f and define the truncated Lipschitz NN class,N N L f ,K d1,d2 (W, L) = {max{-K, min{f, K}} : f ∈ N N d1,d2 (W, L) , Lip(f ) ⩽ L f } .For simplicity, we omit the dimensions and use shorthand N N L f ,K (W, L). The next theorem gives the consistency result of the min estimator. To state the theorem, we defineF * = E t∼µ T E M * [F (V 1 (t), • • • , V n V (t))] to be the true value and F L to be the optimal value of the following optimization problem over SCMs with Lipschitz constant L.EQUATIONs.t.W (P M(V ), P M * (V )) = 0,where the minimum is taken overM G, F K L , U with F K L = {f : ∥f ∥ ⩽ K, f ∈ F L }and all latent distributions P (U ). Note that if L is the Lipschitz bound L f that we assume on our structural functions, then F L is the sharp lower bound. Theorem 4. Let M * be any SCM satisfying the assumptions of Corollary 1. Suppose that the Lipschitz constant of functions inM is L f , F : R n V → R in (6) is Lipschitz continuous and τ n > 0, let K > 0 be the constant in Assumption 3, Lf = d in max d out max L f , F 0,n =N N Lf ,K W 0,n , Θ log d in max , F 1,n = N N ∞,∞ (Θ(d U max ), L i,n ), F 2,n =N N ∞,K (Θ((d U max ) 2 ), L i,n), where d U max , d in max and d out max are defined in Corollary 1, take the radius to beα n = ϵ n + s n , ϵ n = O(W -1/d in max 0,n + 2 i=1 L -2/d U max i,n + τ n log τ n ), s n = O(m -1/(d U max +2) n log m n + δ n + log(nm n )λ n ), δ n = O(n -1/max{2,d U max } log 2 (n)). If m n = Ω(n), L i,n = Θ(m d U max /(2d U max +4) n ), i = 1, 2, lim n→∞ min{τ -1 n , W 0,n , (log(nm n )λ n ) -1 } = ∞ , then with probability 1, [lim inf n→∞ F n , lim sup n→∞ F n ] ⊂ [F Lf , F * ].Note that our theorem shows that the lower limit of the solution is large than the point F Lf , where Lf is slightly larger than the original constraint L f we impose on the structural functions. Hence, this point can potentially be slightly smaller than the sharp bound F L f . This worsening is due to the fact that we need to use NNs that satisfy a slightly worse Lipschitz property, to ensure that we have sufficient approximation power. Although Theorem 4 does not guarantee {F n } converges to a point, it states that {F n } may oscillate in the interval [F Lf , F * ], which will still give a useful lower bound to ground-truth value F * . In particular, if the graph is identifiable, we have F * = F Lf and {F n } converges to F * . Also, note that in Theorem 4, we use wide NNs rather than deep NN for F 0,n because results in [53] show that wide NNs can approximate Lipschitz functions while controlling their Lipschitz constants (a property that is not yet established for deep NNs). Similar results can be obtained if wide neural nets are used for all components, invoking Theorem 2.As a special case, we leverage Theorem 4 for a non-asymptotic rate for the ATE without confounding.Proposition 3 (Hölder continuity of ATE). Given two causal models M, M ∈ M(G, F L , U ) satisfying Assumption 1 and Assumption 3, let their observational distributions be ν, µ. Suppose the norms of all variables are bounded by K > 0. If (1) (Overlap) ν is absolutely continuous with respect to one probability measure P and the density p ν (t|Pa(T ) = x) ⩾ δ > 0 for x almost surely and t ∈ [t 1 , t 2 ] and (2) (No Confounding) there is no confounding in the causal graph G, we havet2 t1 (E M [Y (t)] -E M[ Ŷ (t)]) 2 P (dt) ⩽ 2C W δ W (µ, ν) + 2(L + 1) n V W 2 (µ, ν)(t 2 -t 1 ),whereC W = 4(L + 1) n V K + 2K max {(L + 1) n V , 1}. Corollary 2. Let F, µ T in Theorem 4 to be F (V ) = Y, µ T ∼ Unif([t 1 , t 2 ]), ϵ > 0. Suppose that the assumptions in Proposition 3 and Theorem 4 hold, with probability at least 1 -O(n -2 ), we have|F n -F * | ⩽ O( √ α n ), where F n , F * , α n are defined in Theorem 4.In this section, we examine the performance of our algorithm in two settings. We compare our algorithm with the Autobounds algorithm [14] in a binary IV example in [14] and in a continuous IV model 1 . Since Autobounds can only deal with discrete models, we discretize the continuous model for comparison purposes. The implementation details are provided in the Appendix D. The setting of the first experiment is taken from [14, Section D.1] . This is a binary IV problem and we can calculate the optimal bound using LP [4] . We find that our bound is close to the optimal bound. The second experiment is a general IV example where the treatment is binary but the rest of the variables are continuous. The program that Autobounds solves after discretization contains ≈ 2 14 variables. Even with such a fine discretization, the bound obtained by Autobounds is not tighter than our NCM approach. The details of the structural equations and analysis can be found in Appendix D. We also provide an extra experiment on the counterexample of Proposition 5 in the appendix to show the effect of Lipschitz regularization. 1 : Experiment results of 2 IV settings. The sample sizes are taken to be 5000 in each experiment. STD is the standard derivation. The experiments are repeated 10 times for binary IV and 50 times for continuous IV. In all experiments, the bounds given by both algorithms all cover the true values.In this paper, we provide theoretical justification for using NCMs for partial identification. We show that NCMs can be used to represent SCMs with complex unknown latent distributions under mild assumptions and prove the asymptotic consistency of the max/min estimator for partial identification of causal effects in general settings with both discrete and continuous variables. Our results also provide guidelines on the practical implementation of this method and on what hyperparameters are important, as well as recommendations on values that these hyperparameters should take for the consistency of the method. These practical guidelines were validated with a small set of targeted experiments, which also showcase superior performance of the neural-causal approach as compared to a prior main contender approach from econometrics and statistics, that involves discretization and polynomial programming.An obvious next step in the theoretical foundation of neural-causal models is providing finite sample guarantees for this method, which requires substantial further theoretical developments in the understanding of the geometry of the optimization program that defines the bounds on the causal effect of interest. We take a first step in that direction for the special case, when there are no unobserved confounders and view the general case as an exciting avenue for future work.A Illustration of Notions in Section 2V 2 V 1 V 3 V 4 V 5 U 1 U 2 U 3 U 4 Figure 2: An SCM example. V 2 V 1 V 3 V 4 V 5Figure 3 : The causal graph of this SCM.To further explain the notions in Section 2, we consider the following example. Let M be an SCM with the following structure equations.EQUATIONV 4 = f 4 (V 3 .U 3 , U 4 ), V 5 = f 5 (V 4 , U 3 ),where U 1 and U 2 are correlated and U 3 , U 4 and (U 1 , U 2 ) are independent. The causal model is shown in Figure 2 and its causal graph is shown in Figure 3 . In this example,U V1 = {U 1 }, U V2 = {U 1 , U 2 }, U V3 = {U 2 , U 4 }, U V4 = {U 3 , U 4 }, U V5 = {U 3 }. Since U 1 , U 2 are correlated, C 1 = {V 1 , V 2 , V 3 } is one C 2 component because all nodes in C 1 are connected by bi-directed edges. Note that {V 1 , V 2 , V 3 , V 4 } is not a C 2 component because V 4 and V 2 is not connected by any bi-directive edge. The rest of C 2 components are C 2 = {V 4 , V 5 }, C 3 = {V 3 , V 4 }.Now, we consider the intervention V 1 = t. Under this intervention, the structure equations can be obtained by setting V 1 = t while keeping all other equations unchanged, i.e.,V 1 (t) = t, V 2 (t) = f 2 (U 1 , U 2 ), V 3 = f 3 (V 1 , V 2 , U 2 , U 4 ), V 4 = f 4 (V 3 .U 3 , U 4 ), V 5 (t) = f 5 (V 4 , U 3 ).The figure of this model under intervention is shown in Figure 4 .V 2 (t) V 1 (t) V 3 (t) V 4 (t) V 5 (t) U 1 U 2 U 3 U 4 Figure 4: The SCM after intervention V 1 = t.In this section, we illustrate how to construct canonical representations and neural architectures given a causal graph via a simple example. We consider the example in the previous section. The causal graph is shown in Figure 3 . Following the construction in Proposition 4, we use one latent variable for each C 2 component. As we explain in Appendix A, this causal model has threeV 2 V 1 V 3 V 4 V 5 E 1 E 2 E 3V 2 V 1 V 3 V 4 V 5 Ê1 Ê2 Ê3 Z 2 Z 1 Z 3C 2 components, {V 1 , V 2 , V 3 }, {V 3 , V 4 }, {V 4 , V 5 }.In the canonical representation, exactly the latent variables enter their corresponding C 2 component. The structure equation of this SCM is as follows. 10) is equivalent to (8) . Therefore, we can see from this example that the canonical representation does not lose any information about the SCM. Now, we show how to construct the NCM architecture from a canonical representation. As we mentioned in Section 3, we approximate the latent distribution by pushing forward uniform and Gumbel variables. The structure equation of the NCM isEQUATIONV 4 = f 4 (V 3 .E 2 ), V 5 = f 5 (V 4 , E 3 ), If we set E 1 = (U 1 , U 2 ), E 2 = U 4 , E 3 = U 3 , (EQUATIONV 4 = f θ4 4 (V 3 .g θ1 2 (Z 2 )), V 5 = f θ5 5 (V 4 , g θ3 2 (Z 3 )),where f θi i , g θj j are neural networks, Z i are join distribution of independent uniform and Gumbel variables, and g θj j has the special architecture described in Section 3.1 and Section 3.2. Figure 6 shows the architecture of the NCM. Each in-edge represents an input of a neural net.The causal graph G does not specify how latent variables influence the observed variables. There could be many ways of recovering the latent variables from a graph G. Figure 7 shows an example where two different causal models have the same causal graph. The next proposition gives a canonical representation of a causal model. Proof of Proposition 4. Let G = {G Vi : V i is categorical} and the latent variables in the original modelV 1 V 2 V 3 U 1 V 1 V 2 V 3 U 2 U 1 U 3M be U = {U 1 , • • • , U n U } ∪ G.By Assumption 1 and Assumption 2, U are independent and the structure equations of M have the formV i = f i (Pa(V i ), U Vi ) , V i is continuous, arg max k∈[ni] g Vi k + log (f i (Pa(V i ), U Vi )) k , V i is categorical,where g Vi k are i.i.d. standard Gumbel variables andU Vi ⊂ {U 1 , • • • , U n U }contains the latent variables that affect V i . We regroup and merge the variables U i to make the model have a canonical form while not changing the functions in the structure equations. LetD 1 , • • • , D n C ⊂ V be the C 2 -component of M.For each U i , we define the vertices that are affected by U i asI(U i ) = {V j : U i ∈ U Vj }. We partition {U 1 , • • • , U n U } into n C sets Û 1 , • • • , Û n C in the following way. For each U k , U k is in the set Û i if I(U k ) ⊂ D i .If there are two components D i , D j satisfy the condition, we put U k into either of the sets Û i , Û j . Let Ûi have the same distribution as the joint distribution of the random variables in set Û i . Let M ∈ M(G, F, Û ) the SCM with structure equationsVi =    f i Pa( Vi ), Ûk1 , • • • , Ûkn i , V i is continuous, arg max k∈[ni] g Vi k + log f i Pa( Vi ), Ûk1 , • • • , Ûkn i k , V i is categorical, , Û ki ∩U Vi ̸ = ∅, i ∈ [n V ].Here, we slightly abuse the notationEQUATIONthe same distribution as U because we only merge some latent variables. In addition, the functions in the new model M has the same smoothness as the original model.We first verify that the causal graph of M is G. If there is a bi-directed arrow between nodes V i , V j in G, by the independence assumption, there must be one latent variableU k ∈ U Vi ∩ U Vj . There exist one D l such that I(U k ) ⊂ D l and Û l ∩ U Vi ̸ = ∅, Û l ∩ U Vj ̸ = ∅.Therefore, Ûl will affect both V i and V j and there is a bi-directed arrow between node V i , V j in G M. Suppose that there is a bi-directed arrow between node V i , V j in G M, it means there exist a Ûk such thatÛ k ∩ U Vi ̸ = ∅, Û k ∩ U Vj ̸ = ∅. Let U l1 ∈ Û k ∩ U Vi , U l2 ∈ Û k ∩ U Vj , then V i ∈ I(U l1 ) ⊂ D k , V j ∈ I(U l2 ) ⊂ D k . Since D k is a C 2 -component, there exist a bi-directed arrow between node V i , V j in G. Therefore, G = G M.Finally, we verify that M and M produce the same intervention results. The intervention distribution can be viewed as the push-forward of the latent distribution, i.e., V = F (U ), and the intervention operation only changes the function F . In our construction, we only merge the latent variables and leave the functions in structure equations being the same (except that they may ignore some coordinates in input). For any intervention T = t, suppose in M we have V (t) = F t (U ). Then, in M, we get V (t) = F t ( Û ). Since U and Û has the same distribution, the distribution of V (t), V (t) are the same.Proof of Theorem 1. The structure equations of M have the formV i = f i (Pa(V i ), U Vi ) , V i is continuous, arg max k∈[ni] g Vi k + log (f i (Pa(V i ), U Vi )) k , V i is categorical, ∥f i ∥ 1 = 1,where g Vi k are i.i.d. standard Gumbel variables and M ∈ M(G, F, Û ) has structure equationsVi =    fi Pa( Vi ), Û Vi , V i is continuous, arg max k∈[ni] ĝVi k + log fi Pa( Vi ), Û Vi k , V i is categorical, ∥ f1 ∥ 1 = 1,where ĝVi k are i.i.d. standard Gumbel variables. Let the treatment variables set beT = {T 1 , • • • , T n T }.We may give all vertices andU 0 = {U 1 , • • • , U n U }, an topology ordering (rearrange the subscript if necessary) U 1 , • • • , U n U , T 1 , • • • , T n T , V 1 (t), • • • , V n V -n T (t)such that for each directed edge (V i (t), V j (t)), vertex V i (t) lies before vertex V j (t) in the ordering, ensuring that all edges start from vertices that appear early in the order and end at vertices that appear later in the order. We put U 0 and T at the beginning because they are the root nodes of the intervened model. We denote µ T U 0,T ,1:k (resp. μT U 0,T ,1:k ) to be the distribution ofEQUATIONBy definition of Wasserstein-1 distance,EQUATIONEQUATION.For (1), if V k is a continuous variable, we haveEQUATIONwhere we have used the Lipschitz property of g in the first inequality.EQUATIONwhere we useEQUATION) Since g, f k are Lipschitz continuous functions, g (u, t, v 1 , • • • , v k-1 , f k (pa(v k ), u V k )) is (L + 1)- Lipschitz continuous with respect to (u, t, v 1 , • • • , v k-1). We haveEQUATION)EQUATIONCombine ( 12)-( 16), we prove Equation (11) . By induction, one can easily get W µ T U 0,T ,1:n V -n T , μT Û 0 ,T ,1:n V -n T ⩽ (L + 1) n V -n T W µ T U 0,T , μT Û 0 ,T + 2K 2 S/L -2K 2 S/L = (L + 1) n V -n T -1 L 2K 2 S + (L + 1) n V -n T W µ T U 0 ,T , μT Û 0,T Let C G (L, K) = 2K 2 •max{ (L+1) n V -n T -1 L , (L+1) n V -n T } and notice that W µ T U 0 ,T , μTEQUATIONB.3 Proof of results in Section 3.1Proof of Theorem 2. Recall that the neural network consists of three parts ĝ = ĝτ 3 • ĝ2 • ĝ1 and ĝ2 , ĝ1 have a separable form, dealing with each coordinate of the input individually. As mentioned before, each coordinate approximates the distribution of one connected component of the support. We first construct ĝ1 and ĝ2 , then use Gumbel-Softmax layer to combine each component together.To construct the first two parts, we only need to consider each coordinate individually. For the i-th component C i , let µ i = P Ci /P(C i ), where P Ci is measure P restricted to component C i . Under Assumption 4, there exists a Lipschitz mapg i 2 ∈ H([0, 1] d C i , C i )∩F L ([0, 1] d C i , C i ). By [43, Theorem VIII.1], there exists quantized ReLu network ĝi 1 of width W 1 and depth Θ(d C i ) (let s = 1 ⌈n(n-1)⌉ in [43, Theorem VIII.1]), such that W (g i 2 ) -1 # µ i , P (ĝ i 1 (U i )) ⩽ O W -1/d C i 1 ,where ∥ g i 2 j -ĝi,j 2 ∥ ∞ ⩽ O L -2/d C i 2 , Let ĝi 2 (x) = ĝi,1 2 (x), • • • , ĝi,d 2 (x) , we get ∥g i 2 -ĝi 2 ∥ ∞ ⩽ O L -2/d C i 2 .Thus, the width of ĝi2 is Θ d• d C i . Let ĝi (x) = ĝ1 i (x 1 ), • • • , ĝN C i (x N C ) , i = 1, 2, x k ∈ R di , d 1 = 1, d 2 = d, k = 1, • • • , N C ,where N C is the number of connected components. By Lemma 3, we getW µ i , P (ĝ i 2 • ĝi 1 (U i )) = W g i 2 # g i 2 -1 # µ i , ĝi 2 # P (ĝ i 1 (U i )) ⩽ LW g i 2 -1 # µ i , P (ĝ i 1 (U i )) + ∥g i 2 -ĝi 2 ∥ ∞ = O W -1/d C i 1 + L -2/d C i 2 .Next, let p i = P(C i ) and the distribution of ĝk 2 • ĝk 1 (U k ) be μk , then we have P = N C k=1 p k µ k . By Lemma 4, we haveW P, N C k=1 p k μk ⩽ N C k=1 p k W (µ i , μk ) = O W -1/ max{d C i } 1 + L -2/ max{d C i } 2 .Finally, we analyze the error caused by the Gumbel-Softmax layer. Note that as temperature parameter τ → 0,w τ = exp((log p i + G i )/τ ) N C k=1 exp((log p k + G k )/τ ) i=1,••• ,N C a.s.--→ One-hot(argmax i {G i + log p i })where One-hot(k) is N C -dimensional vector with the k-th coordinate equals to 1 and the remaining coordinates are 0 and G i ∼ Gumbel(0, 1) are i.i.d. Gumbel distribution. Let ν τ be the distribution of w τ . By Lemma 6, for 0 < τ < 1, we haveW ν τ , ν 0 ⩽ O(τ -τ log τ ).By Lemma 5,EQUATIONLetĝτ 3 (x 1 , • • • , x N C ) = N C k=1 w τ k • x k , x k ∈ R d , w τ = w τ 1 , • • • , w τ N C .We denote g 0 3 = lim τ →0 g τ 3 . By Assumption 4, the support is bounded. Thus, g i 2 and ĝi 2 are bounded, i.e., ∥ĝ i 2 (x)∥ ∞ ⩽ K. Since functions h(w, X) = w T X is Lipschitz continuous in the regionEQUATIONfrom (18) . Putting things together, we getW (P, P (ĝ τ 3 • ĝ2 • ĝ1 (U k ))) ⩽ W P, N C k=1 p k μk + W N C k=1 p k μk , P (ĝ τ 3 • ĝ2 • ĝ1 (U k )) = W P, N C k=1 p k μk + W P (ĝ 0 3 • ĝ2 • ĝ1 (U k )), P (ĝ τ 3 • ĝ2 • ĝ1 (U k )) ⩽ O τ -τ log τ + W -1/ max{d C i } 1 + L -2/ max{d C i } 2, where we use (19) . Before we prove Proposition 1, we need to introduce some important notions. Let {C n } be a sequence of partitions of the unit cube [0, 1] d . We say {C n } has property (*) if the following properties are satisfied.1. C 0 = {C -1 }, C -1 = [0, 1] d 2. C n+1 = {C i1,i2••• ,in+1 } ij =0,••• ,2 d -1 is a refinement of C n = {C i1,i2••• ,in } ij =0,••• ,2 d -1 , i.e., C i1,i2••• ,in = 2 d -1 in+1=0 C i1,i2••• ,in+1 .Besides, C n is obtained by cutting the unit cube [0, 1] d by planes that are paralleled to the coordinate planes.The following lemma is important in the construction of Hilbert curve. Lemma 1 (Algorithm 3 in [25] ). Given a sequence of partitions of the unit cube {C n } that has property (*), there exist a ordering for each partition C n , i.e., O n : 1,• • • , 2 nd → i 1 i 2 • • • i n : i j = 0, • • • , 2 d -1, such that the following two conclusions hold.1. If (j -1)• 2 d < i ⩽ j• 2 d , the first n digits of O n+1 (i) are the same as O n (j). 
since i k consists of one-dimensional intervals, we have diam(i: 
 k ) = gap(I k ).Therefore, γ is 1/d-Holder continuous.Proof of Theorem 3. The proof is similar to the proof of Theorem 2. The difference is that in the first part, we use a deep ReLU network to approximate the Holder continuous curve γ instead of wide ReLU network.Let µ i = P Ci /P(C i ). By Assumption 4 and Assumption 5, for each connected component C i of the support, there exists Lipschitz mapsg i 2 ∈ H([0, 1] d C i , C i ) ∩ F L ([0, 1] d C i , C i ) such that (g i 2 ) -1 # µ i (B) ⩾ C g λ(B)/P(C i ), for any measurable set B ⊂ [0, 1] d C i . By Proposition 1, there exist a 1/d C i -Holder continuous curve γ i : [0, 1] → [0, 1] d C i such that (γ i ) # λ = (g i 2 ) -1 # µ i .By [53, Theorem 2] (take ω(x) = C∥x∥ 1/d where C is the Holder continuity constant factor), there exists deep ReLu network ĝi,j 1 of width 12 and depth L 1 such that∥ĝ i,j 1 -(γ i ) j ∥ ∞ ⩽ O L -2/d C i 1 ,where(γ i ) j is the j-coordinate of γ i . Let ĝi 1 (x) = ĝi,1 1 (x), • • • , ĝi,d 1 (x) , we get ∥ĝ i 1 -γ i ∥ ∞ ⩽ O L -2/d C i 1 .Thus, the width of ĝi1 is Θ d C i . By Lemma 3 W g i 2 -1 # µ i , P (ĝ i 1 (U i )) = W (γ i ) # U i , P (ĝ i 1 (U i )) ⩽ ∥ĝ i 1 -γ i ∥ ∞ ⩽ O L -2/d C i 1 .The rest are the same as proof of Theorem 2.Proof of Corollary 1. Suppose that structure equations of M areV i = f i (Pa(V i ), U Vi ) , V i is continuous, arg max k∈[ni] g Vi k + log (f i (Pa(V i ), U Vi )) k , V i is categorical, ∥f i ∥ 1 = 1 , f i ∈ F L .Let d in i be the input dimension of f i , d out i be the output dimension,fi =    arg min fi∈N N d in i ,d out i (d out i (2d in i +10),L0) ∥ fi -f i ∥ ∞ , V i continuous, max{0, arg min fi∈N N d in i ,d out i (d out i (2d in i +10),L0) ∥ fi -f i ∥ ∞ }, V i categorical.We truncate the neural network at 0 for categorical variables because the propensity functions are required to be non-negative. This truncate operation will not influence the approximation error since f i are non-negative. According to Assumption 3,f i are Lipschitz continuous. By [53, Theorem 2], we have, if V i is continuous, ∥ fi -f i ∥ ∞ ⩽ O(L -1/d in i 0 ) ⩽ O(L -1/d in max 0).If V i is categorical, let S = ∥ f (pa(v k ), u V k ) -f (pa(v k ), u V k ) ∥ ∞ , p = f (pa(v k ), u V k ) , p = f (pa(v k ), u V k ) /∥ f ∥ 1 . If S > 1/(2K), ∥p -p∥ ∞ ⩽ 2KS. If S ⩽ 1/(2K), ∥p -p∥ ∞ ⩽ 1 ∥ f ∥ 1 ∥ f (pa(v k ), u V k ) -f (pa(v k ), u V k ) ∥ ∞ + 1 ∥ f ∥ 1 |∥ f ∥ 1 -1|∥f ∥ ∞ ⩽ (n i + 1) ∥ f ∥ 1 S ⩽ (n i + 1) 1 -n i S S ⩽ 2(K + 1)Swhere we use Assumption 3 that n i ⩽ K and S ⩽ 1/(2K). Thus,∥f i -fi /∥ fi ∥ 1 ∥ ∞ ⩽ O(∥f i -fi ∥ ∞ ) ⩽ O(L -1/d in max 0). By Theorem 3, there exists a neural network ĝj with architecture in Theorem 3 such thatW (P (U j ), P (ĝ j (Z j ))) ⩽ O L -2/d U j 1 + L -2/d U j 2 + (τ -τ log τ ) ,where Z j ∼ U ([0, 1] N C,j ) i.i.d., N C,j is the number of connected components of support of U j , d U j is the dimension of latent variables U j and L 1 , L 2 , τ are hyperparameters of the neural network defined in Theorem 3. Let Ûj = ĝj (Z j ), By Lemma 5,W (P (U 1 , • • • , U n U ), P ( Û1 , • • • , Ûn U )) ⩽ n U j=1 W (P (U j ), P (ĝ j (Z j ))) ⩽ O L -2/d U j 1 + L -2/d U j 2 + (τ -τ log τ ) .Let M be the casual model with the following structure equationsVi =    fi Pa( Vi ), (ĝ j (Z Cj )) U C j ∈U V i , V i is continuous, arg max k∈[ni] g k + log fi Pa( Vi ), (ĝ j (Z Cj )) U C j ∈U V i k , V i is categorical.By Theorem 1, we getW P M * (V (t)), P M(V (t)) ⩽ O( Vi continuous ∥f i -fi ∥ ∞ + Vi categorical ∥f i -fi /∥ fi ∥ 1 ∥ ∞ + W (P (U 1 , • • • , U n U ), P ( Û1 , • • • , Ûn U ))) ⩽ O(L -2/d in max 0 + L -2/d U max 1 + L -2/d U max 2 + (τ -τ log τ )), for any intervention T = t.Proposition 5. There exists a constant c > 0 and an SCM M * satisfying Assumptions 1-5 with a backdoor causal graph such that there is no unobserved confounding in M * . For any ϵ > 0, there exists an SCM M ϵ with the same causal graph such that W (P M * (V ), P Mϵ (V )) ≤ ϵ and|ATE M * -ATE Mϵ | > c.Proof. Let V = {X, T, Y } be the covariate, treatment and outcome respectively. T is a binary variable. Let structure equations of causal model M δ beX =    -1 , w.p. 1/2 0 , w.p. 1/4 δ , w.p. 1/4, P (T = 1|X = x) = 1/2 , x = -1 or 0, 3/4 , x = δ, Y = T + U y , X < 0 X/δ + U y , X ⩾ 0where all latent variables are independent and U y is mean-zero noise. The distribution of M δ is P X,T,Y (0, 1, y) = p Uy (y)/8, P X,T,Y (0, 0, y) = p Uy (y)/8 P X,T,Y (δ, 1, y) = 3p Uy (y -1)/16, P X,T,Y (δ, 0, y) = p Uy (y -1)/16.P X,T,Y (-1, 1, y) = p Uy (y -1)/4,P X,T,Y (-1, 0, y) = p Uy (y)/4. As δ → 0, distribution of M δ d -→ M * ,where structure equations of M * areP (U 1 = 1) = 3/5, P (U 1 = 0) = 2/5, P (U 2 = 1) = 1/3, P (U 1 = 0) = 2/3, X = -1 , w.p. 1/2, 0 , w.p. 1/2, P (T = 1|X = x) = 1/2 , x = -1, 5/8 , x = 0, Y =    T + U y , X = -1, U 1 + U y , X = 0, T = 1, U 2 + U y , X = 0, T = 0,where U 1 , U 2 , U y are independent. The distribution of M * is P X,T,Y (0, 1, y) = p Uy (y)/8 + 3p Uy (y -1)/16, P X,T,Y (0, 0, y) = p Uy (y)/8 + p Uy (y -1)/16,P X,T,Y (-1, 1, y) = p Uy (y -1)/4, P X,T,Y (-1, 0, y) = p Uy (y)/4.It is easy to see that M * satisfies Assumption 1-5. Some calculation gives W (P M * (V ), P M δ (V )) ⩽ δ/2. It is easy to calculate the ATE of the two models.ATE M * = 19/30, ATE M δ = 1/2.This example implies that even as the Wasserstein distance between P M * (V ), P M δ (V ) converges to zero, their ATEs do not change. As mentioned in the main body, this problem is caused by the violation of Lipschitz continuity assumption forM δ . Note that in M δ , E[Y |X, T ] = X/δ. As δ → 0, the Lipschitz constant explodes.This problem may also arise in (6) . If the distribution ball B n = { M : W (P M, P M * n ) ⩽ α n } includes a small ball around the true distribution S ϵ = { M : W (P M, P M * ) ⩽ ϵ} ⊂ B n and the NCM is expressive enough to approximate all the SCMs in the S ϵ , this example tells us the confidence interval may not shrink to one point as sample size increases to infinity even if the model is identifiable.To prove Theorem 4, we will need the following proposition. Proposition 6. Given a metric space (M, d) and setsEQUATIONand sup θ∈Θ∞ ∥g n (θ) -g(θ)∥ ⩽ δ n , g n ⩾ 0.3. There exists a compact subset Θ∞ ⊂ Θ ∞ such that for any θ 0 ∈ Θ∞ , there exists θ ∈ Θ n such that d(θ, θ 0 ) ⩽ ϵ n .Consider the following optimization problems:EQUATIONEQUATIONandEQUATIONWe assume that the feasible region of ( 26) and ( 27) are nonempty. Let f * n , f * , f * be the minimal of ( 25), ( 26) and ( 27) respectively. Then,[lim inf k→∞ f * n , lim sup k→∞ f * n ] ⊂ [f * , f * ].Proof of Proposition 6. By compactness of Θ ∞ and Θ∞ , the minimizers of ( 26) and ( 27) are achievable insider these two sets. Let θ * ∈ Θ ∞ , θ * ∈ Θ∞ be the minimizer of ( 26) and (27) . We first prove that lim sup n→∞ f * n ⩽ f * . Note that by (24) , the limiting function g is L g -Lipschitz. By condition 3, there exist θ n ∈ Θ n such that d θ n , θ * ⩽ ϵ n . Note thatg n (θ n ) ⩽ |g(θ n ) -g n (θ n )| + |g(θ n ) -g (θ * ) | + |g (θ * ) | ⩽ δ n + L g d (θ n , θ * ) ⩽ δ n + L g ϵ n .Therefore, θ n is a feasible point of (25) . We havelim sup n→∞ f * n ⩽ lim sup n→∞ f (θ n ) = f θ * = f * .Next, we argue that lim inf n→∞ f * n ⩾ f * . If this equation does not hold, there exists ϵ > 0 and subsequencef * n k k∈N such that f * n k < f * -ϵ, ∀k ∈ N. By compactness of Θ ∞ , for each k, there exists a subsequence of θ * n k ∈ Θ ∞ such that θ * n k satisfies constraint of (25) and f * n k = f θ * n k . By compactness, θ *n k k∈N has a converging subseqence. Without loss of generality, we may assume that θ * n k k∈N converges to θ * ∈ Θ ∞ . Since g k converge uniformly to g,lim k→∞ g k θ * n k ⩽ lim k→∞ |g k θ * n k -g θ * n k | + g θ * n k = g θ * ⩽ lim k→∞ α k = 0.θ * is a feasible point of Equation (26) . Since f is continuous on Θ,f θ * = lim k→∞ f θ * n k ⩽ f * -ϵ.which leads to contradiction.Proof of Theorem 4. We begin by defining a proper metric space. By Assumption 3, all random variables are bounded. Suppose that max i,j{∥V i ∥ ∞ , ∥U j ∥ ∞ } ⩽ K. A canonical causal model M ∈ M(G, F, U ) is decided by θ M = (f 1 , • • • , f n V , P (U 1 ), • • • , P (U n U )),where f i are functions in the structure equations ( 3) and U j are uniform parts of the latent variables. We denote M θ to be the SCM represented by θ, the underlying SCM be M θ * and P θ to be the distribution of M θ . We consider the spaceM = F V1 × • • • × F Vn V × P [-K, K] d U 1 • • • × P [-K, K] d U n U, whereF Vi =    F K L [-K, K] d V i,in , [-K, K] d V i,out , V i continuous, {f : ∥f ∥ 1 = 1, f ∈ F K L [-K, K] d V i,in , [-K, K] d V i,out } V i categorical, d V i,in , d V i,outare the input and output dimensions of f i ,F K L = {f : ∥f ∥ ∞ ⩽ K, f Lipschitz continuous} and P(K) is the probability space on K. For θ = (f 1 , • • • , f n V , P (U 1 ), • • • , P (U n U )) , θ ′ = f ′ 1 , • • • , f ′ n V , P (U ′ 1 ), • • • , P (U ′ n U ) , we define a metric on M d(θ, θ ′ ) = n V k=1 ∥f k -f ′ k ∥ ∞ + n U k=1 W (P (U k ), P (U ′ k )).Theorem 1 states that the Wasserstein distance between two causal models is Lipschitz with respect to metric d. Now, we define Θ n . LetP n = P (U ) : U is the latent distribution of M ∈ NCM G (F 0,n , F 1,n , F 2,n , τ n ) .In other words, P n contains all the push-forward measures of the uniform distribution by neural networks. We denoteΘ n = FV1,n × • • • × FVn V ,n × ×P n , where FVi,n = F 0,n , V i continuous, {f /∥f ∥ 1 : f ∈ F 0,n } V i categorical,and F 0,n is defined in Theorem 4. Note that by construction, latent variables are independent and P n can be decomposed into direct produceP n,1 × • • • × P n,n U . Let Θ ∞ = ∪ ∞ n=1 Θ n , g n (θ) = S λn P θ * n (V ), P θ mn (V ) , f (θ) = E t∼µ T E M θ [F (V 1 (t), • • • , V n V (t))]and Θ∞ = {θ * }. Note that f (θ) is a continuous function since by Theorem 1 and the fact that F is Lispchitz continuous,|f (θ) -f (θ ′ )| ⩽ E t∼µ T |W (P θ (V (t)), P θ ′ (V (t))| ⩽ E t∼µ T [O(d(θ, θ ′ ))] = O(d(θ, θ ′ )).Now, we verify the conditions in Proposition 6. 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [NA] Justification: This is a theoretical work. Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? Answer: [NA] Justification: This is a theoretical work. Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?38th Conference on Neural Information Processing Systems (NeurIPS 2024).The code can be found in https://github.com/Jiyuan-Tan/NeuralPartialIDNote that we truncate the neural network to ensure boundness in Theorem 4. This extra truncate operation f = max{-K, min{K, f }} can be viewed as two extra ReLu layers, so Lemma 8 is still applicable. 
 
 Appendix:  


Decision: This work extends neural causal models (NCMs) to the continuous and mixed-type variables and develops consistency results for partial identification. 

Overall, all reviewers were positive about the contribution of the paper and agreed that the theoretical contribuitions of the work are novel and make NCM's more general. I agree with the assessemnt and recommend acceptance.

Reviews: review: paper_summary: paper_summary: This paper develop consistency results for partial identification via neural causal model with both continuous and categorical variables. Their results shed light on the impact of the neural network architecture and Lipschitz regularization during training. The resulting method can be trained via gradient-based optimization algorithms, and is validated on synthetic data.
strength_weakness: strength_weakness: - The empirical validation of the algorithm/result is not extensive, but I understand that the goal of the work is to establish identification results and a thorough empirical study might not be necessary.
- The paper would benefit from giving some proof sketches (or a brief overview of the proof strategy) in the main paper.
- Many of the assumptions and results are not well explained, making it not straightforward to understand the assumptions and implications of the results.
questions: questions: - A conclusion section does not seem to be provided.
- It would enrich the paper to briefly discuss the connection with optimization-based approaches in causal learning (many of which also used similar techniques including augmented Lagrangian and gumbel techniques), although the task is not completely the same (see e.g. https://arxiv.org/abs/2007.01754, https://arxiv.org/abs/1910.08527, https://openreview.net/forum?id=HsSLdHuAmnY).
- It would be helpful to elaborate more on the assumptions required and discuss some examples where the assumptions are satisfied/violated.
- Other related work that could be worth mentioning: https://arxiv.org/abs/2105.12891
limitations: limitations: The limitations have not been well discussed.
paper_summary: paper_summary: The paper extends neural causal models (NCMs) to the continuous and mixed-type variables. NCMs is a neural-networks based tool to perform an automated point/partial identification and estimation of causal queries. The authors provided new theoretic results (1) on how to construct a canonical representation of an SCM, (2) on how to approximate it with sufficiently expressive Lipshitz neural networks, and (3) on additional assumptions to perform the consistent identification/estimation of an average treatment effect (ATE) from an arbitrary causal digram and an observational distribution.  The paper provides experimental results for the partial identification of the ATE in the discrete and continuous outcomes settings.
strength_weakness: strength_weakness: I didn’t find any major or minor weaknesses in the paper. Yet, I have several suggestions on how to improve the paper’s understandability for the general public of the NeurIPS conference:
1. I encourage the authors to provide more examples of how to construct canonical representations for some simple SCMs and then the corresponding neural architectures of the NCMs.
2. I would provide a more precise explanation of Figures 1 and 5, e.g., the number of layers in the yellow and blue blocks. Also, multiple notation elements could potentially be added to the Figures, e.g., the number of connected components in the latent space or the number of confounded components.
3. I also encourage the authors to provide the code of the proposed method.
4. Seems like the conclusion is missing in the main part of the text. 
 
I am willing to increase my score if the authors implement the above-mentioned suggestions.
questions: questions: I appreciate the provided continuous IV experiment. I wonder, what is the performance of the method in other confounded ATE settings where the ground-truth is available, e.g., a no-assumptions bound for the ATE with hidden confounding?
limitations: limitations: The authors have fully stated the limitations of their theory.
paper_summary: paper_summary: This paper provides a novel perspective and solid contribution to the case where continuous and categorical variables both exist.
strength_weakness: strength_weakness: The assumption is somehow strong (e.g., assumption 2).
questions: questions: How do we argue the importance of solving the problem of ``the approximation error and consistency of the optimization-based approach to partial identification via NCMs''?
limitations: limitations: See above.
paper_summary: paper_summary: This paper studies theoretically the partial identification capabilities of neural causal models (NCMs), that is, the extent to which this type of models can approximate the interval on which a certain causal query falls on, $\theta(\mathcal{M}) \in [\underline{F}, \overline{F}]$. To this end, the authors develop a lot of results showing that under ideal scenarios with oracle knowledge on the exogenous distributions, NCMs can approximate any SCM if the proper architecture and regularization is used. Interestingly, the authors show that regularization is key in these settings, as otherwise the observational likelihood can be modelled arbitrarily well, while the interventional distribution do not look alike. It is worth-noting that the theory developed consider general scenarios with confounders and mixed-type data.
strength_weakness: strength_weakness: - **W1.** My personal main concern is that I do not immediately see any practical implication of the theory developed here. What I mean is that, while I find it interesting, if I have understood it correctly it works under the assumption that we know the exogenous distribution of the ground-truth model, which we never know, and which we could be arbitrarily different from whatever we model (we can always change the exogenous distribution and push the transformation to the generator $f$).
- **W2.** While trying to be general, I similarly feel that some of the assumptions can be quite restrictive in practice. E.g., global Lipschitz-continuity is quite a constraining assumption from my point of view.
- **W3.** I know it is a dense theory paper, but the presentation leaves a lot to be desired. There is a lot of similar notation, non-standard definitions (e.g. SCMs), and a lack of intuition/examples that makes the results hard to parse. The paper is simply not accessible to a large percentage of the community.
- **W4.** The experimental section is extremely short, and it would require of more settings, an ablation study of the different components of the theory, and a proper analysis of the results.
questions: questions: I have no specific questions, but I want to make clear that a dense paper like this would require many more hours for me to fully understand it and I may have missed some important points. I will reflect this in my confidence score.
limitations: limitations: Limitations are **not** discussed. In the checklist, it is said to be discussed on the Conclusion section, which does not exist.

Comments: title: Thanks
comment: Many thanks once again for the detailed response. I think I now get a bit better the intuition of the work.

Since the discussion period is over, I am going to increase my score to borderline accept, as I do not feel comfortable giving more score to a work for which I don't have a full grasp on the results. This is mostly on me, coming from a different area, and the "mathematical heaviness" of the work, which would require many more hours of me to fully grasp.

However, this should be not a problem, as my confidence score is already quite low.
Thank you for your question. 

We want to clarify that the approximation theorems (Theorem 1 and Corollary 1) are tools we use to prove the consistency theorem (Theorem 4), which is the central result of our paper. These approximation theorems state that given an SCM (with known latent distribution and structure functions) that satisfies some assumptions, we can approximate it using NCMs. However, in the consistency theorem, *we do not assume that we know the latent distribution or the structure functions (nor does the consistency theorem claim that we recovered the true latent distribution; just that the true latent distribution is part of the “search space”)*. The consistency theorem proved the consistency of the bounds for the target estimand. 

As we mentioned before, the goal of partial identification is to determine bounds for a causal quantity across all SCMs that have the same observation distribution as the ground-truth model. Since it is difficult to search over all SCMs,  we argue that it is sufficient—and more practical—to search over all NCMs. The approximation theorems are used to prove that the set of NCMs is ‘close’ to the set of SCMs so that the errors caused by this substitution in the partial identification problem can be controlled. 

Therefore, the short answer to your question is that we cannot rule out alternative causal models that result from transformations, as long as they satisfy all the assumptions (particularly Assumptions 3 and 4). Since we only observe a subset of the variables, multiple causal models could potentially lead to the same observational distribution. Without additional information about the true underlying model, it is often difficult to uniquely identify a single causal model from observations alone.

This is a typical case in causal inference. Even when the target estimand is point identified (e.g. ATE with no unobserved confounding), and hence our bounds converge to a singleton, that does not mean that the latent distribution of the structural causal model nor the structural functions are identifiable, due to the ambiguity you mention. Still, in such unique identification cases, all such transformations yield the same target estimand.

We hope this clarification helps you better understand our work, and we are happy to address any further questions you might have.
Thank you, now it is clear to me. I maintain my original positive score.
Thank you for your question. 

We would like to point out that the Manski bound [1] is derived specifically for binary treatment settings, whereas in our example, the treatment is continuous. Besides, the causal graph of [1] is different from the leaky mediation graph. For these reasons, it is not immediately clear whether the same bound can be applied directly in our context. A more comparable setting is the binary IV example in our paper, where both algorithms yield similar bounds. In the binary IV example, we use the algorithm from [2] to verify that the bounds we get are correct. 

Additionally, we impose Lipschitz constraints on the structure functions when solving the PI problem, which may result in tighter bounds. However, the analytic characterization of the optimal bounds under Lipschitz constraints has not been established and seems hard to establish in analytic form. 


We are happy to address any further questions you may have.


[1] Manski, Charles F. "Nonparametric bounds on treatment effects." The American Economic Review 80.2 (1990): 319-323.

[2] Balke, Alexander, and Judea Pearl. "Bounds on treatment effects from studies with imperfect compliance." Journal of the American statistical Association 92.439 (1997): 1171-1176.
Thank you for the clarifications.

I still have an important concern regarding the no-assumptions bound experiments (5). 

For the outcome $Y \in [a, b]$, the no-assumption (Manski) bounds on the ATE always have the width of $b - a$ [1]. In the example, provided by the authors, $Y \in [-11, 11]$ and, thus, the ground-truth width is $b - a = 22$. Therefore, it seems that the method proposed by the authors yields invalid bounds with a width of $8.32$ (whereas the Autobounds package provides a more valid width of $23.53$). Could the authors explain the source of the invalidity?   

**References**:
- [1] Manski, Charles F. "Nonparametric bounds on treatment effects." The American Economic Review 80.2 (1990): 319-323.
I thank the authors for their kind and detailed response, it should be clear by now that I am an outsider on this particular problem, and the initial barrier is a bit too tough to pass in the reviewing period. Clearly, other reviewers were much more positive than I initially was.

With that said, I understand the practical implications of this work, and I applaud it. However, I still have some questions for which I would appreciate an answer (I apologize I don't have currently the time to go through the proofs in the appendix, yet I'll try to go through the main paper one more time).

To my understanding, one key aspect of this work is splitting the estimation error in two terms, the error in the exogenous distribution and that of the causal mechanisms. From there, one can use the theory of NN approximation to estimate the error committed in each part, depending on depth and width. To achieve this result, it is necessary to assume a specific $\mathcal{M}^*$ in its canonical form, which would have a specific Lipschitz constant $L^*$ that the proofs should somehow use.

What I don't understand is that, without assuming that one knows the exogenous distribution of $\mathcal{M}^*$, it is possible to find alternative canonical representations by applying arbitrary transformations $\phi$ to the canonical form, i.e., having $(\phi \circ P_U, \phi^{-1} \circ f)$ instead of $(P_U, f)$, with arbitrary Lipschitz constants for the target SCM. As far as I can see, the definition of the canonical form does allow these cases. How does your work rule them out? I guess some of the assumptions do it?

Thanks in advance.
Thank you for the helpful replies. Most of my concerns have been addressed and I would like to maintain my rating.
Thank you so much for your feedback and for sharing the recent literature with us. We sincerely appreciate your contributions to improving our work.

Regarding the two papers you mentioned, they indeed focus on a specific type of causal graph, which differs from the general causal graph setting our algorithm addresses. However, we agree that these papers offer valuable insights for the partial identification literature. We will cite them in the camera-ready version of our paper.

As for your question about discretization, we appreciate the opportunity to clarify this point. As we mention in the paper, for a general discrete SCM, the Partial Identification (PI) problem becomes a Polynomial Programming (PP) problem [1], which is NP-hard in general. While discretizing continuous data and solving the PP problem to obtain bounds is possible, the problem's size grows *exponentially* with the cardinality of the support. This implies that the finer the discretization of the continuous data, the larger and more computationally expensive the resulting PP problem becomes.

In our experiments, we found that the PP problem quickly becomes intractable when we increase the support's cardinality of continuous variables. For instance, in a continuous IV setting, we compared our algorithm with Autobounds [1], the state-of-the-art algorithm for solving discrete PI problems, by discretizing all continuous variables such that their support's cardinality is 8 (details are in Section D.2). This produces a PP problem of size approximately $2^{14}$. However, even we solved such a large problem, the bound obtained by this approach was not as tight as the one obtained using the NCM approach. While finer discretization could potentially improve the bound, the resulting PP problem becomes so large that we were unable to obtain a solution after running the Autobounds algorithm for over a day. Therefore, we believe it is crucial to generalize the NCM approach to include continuous (and mixed) variables, as it proves to be more efficient and powerful compared to the discretization approach in these settings.

Thank you again for your detailed review and valuable questions. We are glad to take any further questions. 


[1] Guilherme Duarte, Noam Finkelstein, Dean Knox, Jonathan Mummolo, and Ilya Shpitser. An automated approach to causal inference in discrete settings, September 2021
It is a concise and high-quality response. Thank you! I raise my score to 6.

Besides, in fact, I'm quite familiar with Partial Identification, one of my main research topics. I also share two nice new papers and recommend you cite them as follows.

Of course, I acknowledge that this literature is somehow different from your settings, but they all make meaningful contributions to general PI. For instance, after my careful proofreading, Literature [1] 's main IFF theorem holds in discrete/continuous/mixed cases.

Finally, my last question is: Why is the extension from the discrete case to the continuous case highly non-trivial? For instance, when we get the continuous variable, why not conduct naive discretization, and the corresponding result may not change significantly? If this holds true, it will not be interesting. So, I wonder about some more insightful reasons/counterexamples.

Thank you! I will continuously consider raising my score.



[1] tight partial identification of causal effects with marginal distribution of unmeasured confounders. (Zhang, ICML2024 spotlight)

[2] Model-Agnostic Covariate-Assisted INference on Partially Identified Causal effects (Wen et al)