{'fraction_non_alphanumeric': 0.03721406623420963, 'fraction_numerical': 0.00805735745988392, 'mean_word_length': 4.509969902182092, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Disentangling Observed Causal Effects from Latent Confounders using Method of Moments

Discovering the complete set of causal relations among a group of variables is a challenging unsupervised learning problem. Often, this challenge is compounded by the fact that there are latent or hidden confounders. When only observational data is available, the problem is ill-posed, i.e. the causal relationships are non-identifiable unless strong modeling assumptions are made. When interventions are available, we provide guarantees on identifiability and learnability under mild assumptions. We assume a linear structural equation model (SEM) with independent latent factors and directed acyclic graph (DAG) relationships among the observables. Since the latent variable inference is based on independent component analysis (ICA), we call this model SEM-ICA. We use the method of moments principle to establish model identifiability. We develop efficient algorithms based on coupled tensor decomposition with linear constraints to obtain scalable and guaranteed solutions. Thus, we provide a principled approach to tackling the joint problem of causal discovery and latent variable inference.


Introduction
Discovery of causal relationships among a group of variables lies at the heart of scientific discovery [1]. Inferring causes and effects, and separating them from spurious correlations, is central to scientific domains such as medical drug trials, high energy physics, and econometrics. From a machine learning perspective, discovery of causal effects can lead to more robust, reliable, and interpretable predictions.
Causal inference aims to learn strength of causal interactions on a known graph, while causal discovery also aims to learn the causal graph (and the parameters) from data. This problem is challenging especially in the presence of hidden or latent confounders, which make the problem ill posed. These are variables that are not directly observed but interact with the observables. These hidden variables cause spurious associations among the observed variables and make it challenging to differentiate causal relations from spurious correlations. For instance, in drug trials, there are many unobserved factors that affect drug efficacy in a patient (e.g. their lifestyle choices), and not all of that is recorded. As another example, to do system identification in physical systems, like control systems, not all the contributors to the dynamics are observed. This can be due to the difficulty in data collection or the limitations in data fieldity [2]. On the other hand, it is feasible to intervene or manipulate the physical system to record different kinds of data. Using such interventions to learn causal relationships is an interesting problem.
The problem we consider is two-fold. We aim to simultaneously discover the relationships among a group of observed variables, and in addition, discover latent factors and their relationships with the observables. We aim to solve both causal structure discovery and latent variable inference.
Causal Discovery with Interventions While observational data is easy to obtain, causal relationships cannot always be uniquely identified using only observational data. Interventions provide essential information for discovering the causal effects. For example, in structural equation models (SEMs) [3], the act of assigning a variable X to a specific value x allows us to dislodge X from its prior causes and evaluate the expected change of other variables. More generally, softly intervening on a subset of variables induces a new probability distribution over the variables [1]. This has been utilized previously for estimating causal effects when the confounders are also observed variables [1]. What we consider is the general interventional framework in this paper. Interventions are usually represented using Do operator. For example, in Figure 1, Do(X 1 ) represents intervention on variable X 1 . As a consequence, the incoming edge from H to X 1 will be removed after the intervention. Causal Discovery with hidden confounders Previous works that focus on learning causal models with hidden confounding variables generally consider only observational data [4][5][6]. However, these works do not consider information obtained through interventions and hence, require strong modeling assumptions. Hence, they may not be applicable in many settings. For instance, the linear, non-Gaussian, and acyclic model (LiNGAM) has been investigated for causal discovery using the independent component analysis (ICA) model under non-gaussian noise assumptions on observational data [7,8]. It is recently extended to the hidden confounder setting [6]. However, using observational data, it can only recover a class of causal models that are equivalent and compatible with observational data, but this does not guarantee unique identifiability of the causal relationships.
Latent variable models In terms of latent variable modeling without causal relationships, there is a rich body of work proposing different classes of models and algorithms. We consider the independent component analysis (ICA) which assumes the presence of independent latent factors that mix together to generate the observed data [9]. Learning the mixing process provides a natural disentangled representation of the observed domain [10]. The simplest setting is when the mixing process is linear. In this case, the mixing matrix can be recovered provably and efficiently under mild assumptions using tensor decomposition methods and higher order (3rd or 4th) moments [11].
In contrast to tensor methods, standard approaches for learning latent variable models such as expectation maximization (EM) are not guaranteed to converge to the globally-optimal solution and are prone to failure in high dimensions. We consider the method of moments approach that fits observed moments in order to estimate the model. For a wide class of latent variable models that includes ICA, [11] showed that low order moments (up to third or fourth) is sufficient to provably learn the model. Efficient tensor decomposition algorithms have been developed that guarantee learning with polynomial sample and computational complexity.
Our formulation In this paper, we propose a new model termed as SEM-ICA that combines the two frameworks, viz., linear SEMs and ICA, to incorporate both observables and latent confounders. We use a mixing matrix to represent the confounding effects between the latent variables and the observable variables. We use a DAG graph to represent the relationship among observable variables. Thus, we assume a latent generative process, given by the ICA model, and this acts as latent confounders to the causal discovery process among the observables. This is a reasonable assumption in many scenarios. For example, in power networks, our intervention on different distributors would not change the latent variables such as unrecorded weather or temperature.
We assume that (soft) interventions are available and this helps disentangle the causal and latent effects. A key observation we make is the fact that there is invariance of the latent effects across different interventions, since we assume that they are driving the generative process of observations. This insight allows us to relate the data obtained under different interventions. Our method generally allows for both hard interventions and soft interventions. But we focus our analysis and simulation on the hard interventions.
Algorithmic Contributions We establish a novel algorithm for learning linear SEM-ICA based on the method of moments. We consider moments computed both from observational and interventional data. In particular, each of these datasets can be expressed as another ICA model. Thus, we cannot hope to learn the correct underlying model by considering only observational data. In particular, we need to understand how the different datasets from different interventions are related to one another. We show that the mixing matrices of these models are related to the parameters of the underlying SEM-ICA model through a linear system of equations. We assume that the distribution of the generative latent factors is invariant under interventions. This holds in many cases when environmental conditions being invariant under controlled interventions. For example, in power networks, our intervention would not change the latent variables such as unrecorded weather or temperature. We then formulate a joint optimization method that incorporates both the tensor decomposition and the linear systems with invariance constraints to obtain recovery of the mixing and observed matrices.
Theoretical Guarantees We analyze a variant of our algorithm and prove identifiability for the mixing and observed causal matrices. We demonstrate that under mild non-degeneracy conditions our method can identify both the mixing and observed causal matrices up to permutation, sign and scale. This requires the mixing matrix of the ICA model to be full-rank and the adjacent graph among observables to be directed acyclic graph(DAG). Both are natural for identifiability, which implies that the latent factors are not (linearly) degenerate. We also establish the learning error guarantees for the whole algorithm and show that we achieve polynomial sample and computational complexity. This is the first work that provides guaranteed causal discovery with hidden confounders.

Problem Setup
Given a set of n observable variables X := {X 1 , . . . , X n } and a set of m confounding variables The edge set E consists of all directed edges between vertices, which represent the causal relations between two neighboring variables. For a vertex X ∈ X , the parent set of X is denoted by Pa(X) ⊆ X H such that X ∈ Pa(X) if and only if there is a direct edge from X ∈ X H to X. For modeling the causality among confounding variables and observable variables, we assume that the confounding variables are independent and there only exist paths from confounding variables to observable variables, but not the other way around. Moreover, we assume there is no self-loop and the partial graph with vertices in X is acyclic, whence can be regarded as a DAG. Therefore, the graph G is a quasi-bipartite graph, as shown in Figure 2 (a).
We fix the notation for this paper. We use boldface letters X and H for random vectors and denote the i-th entry of X as X i . Observed measurements are realizations of X, written as X [1] , . . . , X [N] and N is the number of samples. Matrices are denoted by uppercase boldface letters A, B, C and D. For a matrix A, the i-th row is written as a i and the i-th column is denoted by a (i) . The (i, j)-th entry of A is denoted by a i,j . The superscript T refers to the transpose of a matrix or vector.

Structural Equation Models with Latent Variables
The set of observable variables is represented by X as a vector in R n . The confounding (hidden) variables are denoted by a vector H in R m . We consider a linear model for the observable and confounding variables: wherein A ∈ R n×m denotes the mixing matrix. It represents the correlation between observable and confounding variables (and equivalently, the edges between confounding variables and observable variables in Figure 2). Similarly, the matrix B ∈ R n×n is a observed causal matrix, which encodes the causal relationship among observable variables (corresponding to directed edges in Figure 2). Additive noise is written as N. Rearranging the terms, we obtain the following alternative expressions: where C : Here, I is an n × n identity matrix. Below we present the assumption of our model. Assumption 1. We assume that A is a full-rank matrix. We also suppose that the faithfulness assumption holds so that a j,k = 0 if there is a causal path from H k to X j and b i,j = 0 if there is a causal path from X j to X i , for all i, j = 1, . . . , n and k = 1, . . . , m.
With the help of n interventions, we can discover the causal ordering of X 1 , X 2 , ..., X n (see (2.3)).
Hence after suitable permutation of the rows and columns of B, we can formulate it as strict lower triangular and the matrix C is full-rank. In this work, we consider a new approach for identifying causal effects, i.e., recovering the mixing matrix A and observed causal matrix B, by conducting interventions, defined in the next section.

Interventions
An intervention in causal discovery refers to the experiment that sets the values of the variables and isolates them from their causes. Definition 2.1. (Do(·) operator): We write Do(·) as the intervention operator. Graphically, each ideal intervention Do(X i = λ i ) corresponds to removing all the directed edges from Pa(X i ) to X i and set the coordinate X i in X to be λ i . We write X i ∈ R n as a vector with the i-th entry set to λ i . (1) still holds by modifying the mixing matrix and observed causal matrix as follow:

After an intervention
where A i and B i are formed by setting the i-th rows of A and B as the all-zero vector and the unit vector e i . Rearranging the terms, similarly we obtain n linear systems by intervening Do(X i = λ i ) for i = 1, . . . , n and obtain the following alternative expressions: where each matrix D i := (I−B i ) −1 A i ∈ R n×m is defined as the response matrix to the intervention The goal of this work is to use both the observational data and the experimental data obtained by interventions Do(X i = λ i ) to identify the mixing matrix A and the observed causal matrix B in (1), under the assumption presented in Assumption 1. In the following section, we discuss independent component analysis ( ICA) models and the associated tensor decomposition methods, that can be utilized for recovering the response matrices C, D 1 , . . . , D n , based on observed samples of X, Y 1 , ..., Y n .

Recovery of Causal Order
It is well known that if single-variable interventions are conducted on the variables, causal ordering can be recovered [12] [13]. Intuitively, by doing single-variable interventions, the causal ordering can be recovered by checking the variables that have changed after the interventions. Previous work provides method that can be used to recover the causal directions with n − 1 interventions [13]. In our framework, we can recover the causal ordering immediately after getting the experimental data by doing interventions Do(X i = λ i ) on the variables X 1 , X 2 , ..., X n . Furthermore, we can also determine matrices A and B using the same n single-variable interventions as shown below.
3 Tensor Decomposition for ICA

Independent Component Analysis
The classical independent component analysis (ICA) model is a well-known blind source separation method for separating a set of source signals from a set of mixed signals, with no or a little prior knowledge about the source signals or mixing process. ICA-based methods decompose a multivariate random signal into a set of independent non-Gaussian signals. Therefore, they can be applied to the linear structural equations in (2) and (4). Taking the equation in (2) as an example, our goal is to recover C whose columns are defined as independent components (ICs), provided with data samples of X. The corresponding graphical representation is depicted in Figure 2. The classical ICA problem is an ill-posed task since it has inevitable indeterminacy of permutation and dilation of the independent components. Therefore, without additional information, we can only recover columns of the matrix C in equation (2) up to permutation and scaling. Under Assumption 1, the matrix C in equation (2) is full-rank. Thus, it is identifiable if latent variables is non-Gaussian.

Tensor Decomposition
Various techniques have been proposed to recover the factor matrix C of an ICA model X = CH+N.
In this work, we use the method in [11] for guaranteed recovery of ICs. This procedure, which we refer to as the tensor method, uses higher-order moments of the random vector X to recover the columns of D. Here, we use the 4-th order cumulant of zero-mean vector X defined as: where T is the following 4-th order tensor, for 1 ≤ i 1 , i 2 , i 3 , i 4 ≤ m: Denote by κ i the excess kurtosis, a measure of non-Gaussianity defined as . The 4-th order moment M 4 can be estimated from observed measurements and tensor decomposition algorithms can be applied to recover the columns of C. In particular, the alternating power updates proposed in [11] can be applied to solve the tensor decomposition problem in (5).

Causal Effects Identification
Based on the causal ordering restored, the matrix B after suitable permutation of the rows and columns becomes a lower triangular matrix. Therefore, without loss of generality, we assume that B is strict lower triangular and thus matrix C is full-rank. The next step is to recover the mixing matrix A and observed causal matrix B, from the linear systems: (I − B i )D i = A i , for i = 0, . . . , n, which is derived from (4) and (2).

Recovery of Response Matrices
To obtain matrices D 0 , . . . , D n satisfying above, the labels of rows and columns of the matrices should be preserved, i.e., the rows and columns should correspond to the same confounding or observable variables. However, the tensor decomposition method reported in [11] leaves the permutation and scaling (PS) for columns of any response matrix unknown. We name this issue PS indeterminacy. Without considreation of time complexity, we can solve this issue in a brute force manner, by trying out different permutations of the hidden vectors since only the original mixing and observed causal matrices make the entire linear system solvable. Hence this issue will not affect the identification of our method. Also, heuristic ways exist for preserving the correspondence between columns of two response matrices D i and D j in practice. With prior knowledge, adding regularizers can also eliminate indeterminacy in an ICA model ( [14]).

Recovery of Mixing and Observed Causal Matrices
For robust identification of A and B, instead of solving the linear systems in section ( 4), we combine the optimization objective in the tensor methods solving ICA with the linear systems to formulate a joint optimization objective: (6c) Here, K 4 is the Kronecker product. In practice, we use Tensorly library [15] to jointly optimize D i together with A and B. Specifically, the library utilizes the alternative least square (ALS) method for solving the optimization problems. In practice, we can also add regularization to match columns of D i and C to eliminate the optimization difficulties brought by permutation. We outline the recovery process in Algorithm 1.
Algorithm 1: Algorithm for recovering mixing and observed causal matrices.
Data: Observed measurements (as copies of X) and model parameters Result: Estimated mixing and observed causal matrices A and B Step 1. Recover the matrix C (using tensor decomposition methods).
Step 2. Obtain response matrices via interventions: for i = 1 to n do Perform the interventions Do(X i = λ i ). Collect experimental measurements (as copies of Y i ). Recover the response matrix D i (using tensor decomposition methods). end Step 3. Recover causal ordering.

Recovery Guarantees
With the recovered response matrices, in this section, we analyze the correct recovery of the mixing matrix A and the observed causal matrix B, with interventions. All the proofs are in the Appendix.
Since the optimization in Algorithm 1 is non-convex, we analyze a variant of Algorithm 1. This variant involves recovering each of C and D i separately using tensor methods using Eq. 7, and finding a permutation on columns of C and D i such that is minimizes the linear system 8a -8b. min and using the recovered C and D i , we optimize the following min A,B,perm(C),perm(Di) where perm(·) represents permutation of the columns. Thus, we find an optimal permutation that minimizes (8b). The above optimization is computationally more expensive than Algorithm 1 since it requires search over all permutations. But it allows us to provide theoretical guarantees on recovery. We now provide theoretical guarantees for this variant.

Recovery with infinite samples
We first analyze the case when infinite number of samples are available to the algorithm. In this case, the following theorem holds. The mixing matrix A ∈ R n×m and the observed causal matrix B ∈ R n×n are uniquely recovered under the assumptions that A has full column rank and B is strict lower triangular.

Recovery with finite samples
We analyze the number of data samples, denoted by N , required in Algorithm 1. Suppose we have no prior information and the subroutine of tensor method is performed in an unsupervised manner 1 .
To be consistent with the settings of Algorithm 1 in [11], the 4-th order moments M 4 with respect to the response matrices C, D 1 , . . . , D n are estimated empirically from N observed measurements. We assume the following conditions for the theorem to hold: (1). Matrix A has full column rank and B is strict lower triangular. Consider the estimated matrices A and B. Let G := (I − B) −1 . The next theorem states that given the conditions in above, the ∞-norm error for the recovered rows of A and columns of G is bounded from above. Theorem 5.2 (Sample Complexity). Suppose the above conditions (1)-(3) hold. Then the rows of the estimated mixing matrix A and the columns of the estimated observed causal matrix G given by Algorithm 1 satisfy that for i = 1, . . . , n, with high probability, when sign ambiguity and column permutation are taken into account.

Necessary Conditions on Interventions
Consider interventions Do(X i = λ i ) on the n variables X 1 , . . . , X n denoted by {I i : i ∈ {0} ∪ {1, . . . , n}} where I 0 is the null intervention of passive observations. The following necessary conditions on the number of interventions hold. Theorem 5.3 (Necessary conditions; SEM). There exists a graph G = (X ∪ H, E) such that any collection interventions {I i : i ∈ S} with max{|S|, m} ≤ n − 1 cannot recover its edge set E. Remark 1. Our results indicate that even with an SEM model, if m ≤ n − 1, then the worst-case recovery of causal effects still requires n interventions. In contrast, for graphical causal models with perfect data and faithfulness assumption similar to our Assumption 1, the authors of [13] (see Proposition 2) have shown that n − 1 interventions, counting the null intervention of passive observations (with each randomizing a single variable) are necessary to determine the causal relations of n variables.

Simulation
We validate Algorithm 1 by showing the trade-offs of the number of samples used and the accuracy of recovery. We generate each H k i.i.d. according to a Laplacian distribution with mean 1 for k = 1, . . . , m and each X i i.i.d. according to a Gaussian distribution with mean 0 and variance 1 for i = 1, . . . , n. The noise vector N has i.i.d. Gaussian entries with mean 0 and variance 10 −3 . Figure 3 depicts the sample complexity results for recovering the observed causal matrix. In particular, we generate data from a linear system with different numbers of observables and hidden variables. We restart the optimization 20 times to find the best optimization error. We plot the best recovery error on B using Algorithm 1 on different number of variables when there are same number of hidden variables and the observables (3 variables: Figure 3 Left, 4 variables: Figure 3 Middle, 5 variables: Figure 3 Right, 10 variables: Figure 4 Left). We can see that the estimation error on B decrease fast with more samples. We also see that the recovery error increases with the number of the observables and hidden variables.

Ablation Study
We conduct an ablation study on misspecified number of hidden variables and less number of interventions to see how robust our method is. Figure   Middle: Performance with increasing number of interventions with a 4-variable model. We can observe that we can use less intervention to achieve similar empirical results with full interventions. Right: Performance with misspecified number of latent variables on a 3-variable model. We can see that the model is relatively robust to small perturbations on number of latent variables.

Related Work
Causal discovery without interventions: Traditionally, if the causal sufficiency assumption holds, i.e., no unobserved common causes exist for any pair of variables that are under consideration in the model, causal structure learning algorithms such as PC [17] and IC [1] can identify a class of causal models that are equivalent to each other in terms of the conditional independence. Given more background knowledge, the models that fit the observed data can be narrowed down. In more general structure learning, assumptions about structures such as trees and polytress are made to learn from observational data [18,19]. In this paper, we focus on modeling possible unobserved common causes, i.e., latent confounders, and establish a unique identification method using experimental data. Utilizing hidden variables not only leads to a more descriptive representation of the model but also relaxes the computational complexity of many algorithms through dimensionality reduction.
Causal discovery using interventions: Studying causal discovery using interventions has a long history. There is a large body of work on experimental design for casual discovery based on combinatorics and search [20][21][22]. Representing the confounding effects implicitly and allowing multi-variable interventions, a linear cyclic causal model can be identified with presence of latent confounders [12]. Recent work focuses on developing efficient randomized algorithm to learn a causal graph with confounding variables in a minimal or constrained budget [23][24][25][26][27][28][29]. Combined with active learning and optimization theory, interventions can be regarded as active data queries with different cost [30,31]. In our case, we assume the invariance of latent confounders to relate mixing matrices across interventions. We restrict to the case of single-variable intervention and learn the parameters of the mixing and causal matrices explicitly.
Causal inference and discovery with latent variable models: Previous research on causal learning and latent variable models focus on the recovering the causal order [32]. Lingam model has also been generalized to hidden variable cases [33] and latent Gaussian confounders [34]. Recent study of causal inference using latent variable models [35] assumes the conditional independence between observable variables given the latent confounders. It opens up a line of research about utilizing latent variables in regression problems for estimating the potential outcome. In contrast to that, our problem setting is a causal discovery one, where we recover both the causal relationships between observable variables together with the confounding effects.

Conclusion and Discussion
In this work, we propose a novel approach to discover causal relation from experimental data based on recovering mixing and causal matrices in the latent linear structural equation model. We provide theoretical analysis regarding the conditions of unique identification and sample complexity for the method. The effectiveness of our approach is validated by experiments. An interesting extension of the methods is to consider non-linear models including deep neural networks, where invariance of a mixing matrix A does not hold in general but need to be enforced by constraints.
A Proof of Theorem 5.1 Proof. Given (I − B)C = A and the linear equations obtained via n interventions Do(X i = λ i ) for i = 1, . . . , n: we show that the matrices A and B are unique up to permutation and scaling.
Subtracting both sides of (9) from (I − B)C = A, for i = 1, . . . , n, which leads to Rearranging the terms we obtain where we have used the triangle inequality to derive (14). Plugging in (13) and (11) into (14), it follows that for all i = 0, . . . , n with high probability, for all j = 1, . . . , m, Therefore, for each matrix F i , with high probability, we can bound the ∞ -norm by its Frobenius norm from above as Since F i is rank-1, there is only one positive singular value, denoted by σ i > 0 and F i can be written as F i = σ i u i v T i . Applying the perturbation bound for SVD in Theorem 1.1 of [37], for all i = 0, . . . , n with high probability, · max{n, m} σ i √ m ε (m, n, N ) , yielding the desired results with ε (m, n, N ) defined in (12).

C Proof of Theorem 5.3
Proof of Theorem 5.3. Suppose G = (X ∪ H, E) is an arbitrary graph with a mixing matrix A and a causal matrix B. For any i ∈ S with |S| ≤ n − 1, the intervened SEM model in (3) is We show that there exists another graph G = G with a mixing matrix A and a causal matrix B such that G is consistent with the SEM in (15). The two matrices B, B ∈ R n×n satisfy that supp(B) = supp (B ) and dim (kernel(B − B )) = n − 1. Moreover, X i = (I − B i ) −1 A i H implies dim (span({X i : i ∈ S})) ≤ max{|S|, m}. Therefore max{|S|, m} ≤ n − 1 ensures that X i ∈ kernel(B − B ) for all i ∈ S with |S| ≤ n − 1. Hence, yielding that at least n interventions (including the null intervention) are necessary.