causal representation learning, causality, interventions

0.0.1 Overview↩︎

Causal representation learning (CRL) aims to form a causal understanding of the world by learning appropriate representations that support causal interventions, reasoning, and planning [1]. Specifically, CRL considers a data-generating process in which high-level latent causally-related variables are mapped to low-level, generally high-dimensional observed data through an unknown transformation. Formally, consider a causal Bayesian network [2] encoded by a directed acyclic graph (DAG) \(\mcG\) with \(n\) nodes and generating causal random variables \(\bZ \triangleq [Z_1,\dots,Z_n]^{\top}\). These random variables are transformed by an unknown function \(g: \R^n \to \R^d\) to generate the \(d\)-dimensional observed random variables \(\bX \triangleq [X_1,\dots,X_d]^{\top}\) according to: \[\bX = g(\bZ) \; . \label{eq:data-generation-process-intro}\tag{1}\] CRL is the process of using the observed data \(\bX\) and recovering (i) the causal structure \(\mcG\) and (ii) the latent causal variables \(\bZ\). Achieving these implicitly involves another objective of recovering the unknown transformation \(g\) as well. Addressing CRL consists of two central questions:

• Identifiability, which refers to determining the necessary and sufficient conditions under which \(\mcG\) and \(\bZ\) can be recovered. The scope of identifiability (e.g., perfect or partial) critically depends on the extent of information available about the data, the underlying causal structure, and the transformation. The nature of the identifiability results can be algorithm-agnostic and non-constructive without specifying how to recover \(\mcG\) and \(\bZ\). In particular, this is the case when considering CRL under a general transform \(g\) without parametric assumptions. Furthermore, the study of identifiability also investigates necessary conditions, e.g., which type of data is required for identifiability, regardless of the algorithmic approach.

• Achievability, which complements identifiability and pertains to designing algorithms that can recover \(\mcG\) and \(\bZ\) while maintaining identifiability guarantees. Achievability hinges on forming reliable estimates for the transformation \(g\).

0.0.1.0.0.1 CRL from interventions.

Identifiability is known to be impossible without additional supervision or sufficient statistical diversity among the samples of the observed data \(\bX\). As shown in [3], [4], this is the case even for the simpler settings in which the latent variables are statistically independent (i.e., graph \(\mcG\) has no edges). On the other hand, using data generated under interventions in addition to the observational data generated by the underlying latent model creates useful statistical diversity. Specifically, an intervention on a set of causal variables alters the causal mechanisms that generate those variables. Note that these causal mechanisms capture the effect of parents on the child variable. Such interventions, even when imposing sparse changes in the statistical models, create variations in the observed data sufficient for learning latent causal representations. This has led to CRL via intervention as an important class of CRL problems, which in its general form has remained an open problem [1]. This paper addresses this open problem by drawing a novel connection to score functions, i.e., the gradients of the logarithm of density functions.

We note that we use the interventions as a weak form of supervision via having access to only the pair of distributions before and after an intervention. Such supervision can be quite flexible and bodes well for practical applications in genomics [5] and robotics [6]. For instance, genomics experiments often involve interventions, which can be modeled as deterministic or stochastic interventions, depending on the experimental mechanism. In robotics, let us consider the causal variables to be the joint angles of a robotic arm. A stochastic intervention corresponds to setting the joint angle to take values within a suitable support. We note that the interventions here can be soft since the feasible support for the joint may still depend on the other joint angles after the intervention. In this setting, the observations are simply images of the entire arm captured by a camera positioned at a specific location.

0.0.1.0.0.2 Contributions.

This paper provides both identifiability and achievability results for CRL from interventions under general latent causal models and general transformations. We establish these results by uncovering hitherto unknown connections between score functions (i.e., the gradients of the logarithm of density functions) and CRL. We leverage these connections to design CRL algorithms that serve as constructive proofs for the identifiability and achievability results. We do not make any parametric assumption on the latent causal model, i.e., the relationships among elements of \(\bZ\) take any arbitrary form. For the transformation \(g\), we consider a diffeomorphism (i.e., bijective such that both \(g\) and \(g^{-1}\) are continuously differentiable) onto its image. Our results are categorized into two main groups based on the form of \(g\): (i) linear transformation, and (ii) general (nonparametric) transformation. We consider both stochastic hard and soft interventions, and our contributions are summarized below.

0.0.1.0.0.3 Linear transformations.

We first focus on linear transformations and investigate various extents of results (perfect and partial) given one intervention per node. In this setting, we consider both hard and soft interventions.

• On identifiability from hard interventions, we show that one hard intervention per node suffices to guarantee perfect identifiability (3).

• On identifiability from soft interventions, we show that one soft intervention per node suffices to guarantee identifiability up to ancestors – transitive closure of the latent DAG is recovered and latent variables are recovered up to a linear function of their ancestors (2).

• While establishing these results, we also establish partial identifiability of a latent variable given an intervention on single node (1).

• We further tighten the previous results and show that when the latent causal model is sufficiently nonlinear, perfect DAG recovery becomes possible using soft interventions. Furthermore, we recover a latent representation that is Markov with respect to the latent DAG, preserving the true conditional independence relationships (4).

• On achievability, we design an algorithm referred to as Linear Score-based Causal Latent Estimation via Interventions (LSCALE-I), which achieves the identifiability guarantees under both soft and hard interventions. LSCALE-I first forms an encoder estimate using correlation matrices of the observed score differences and then uses the learned encoder to estimate the transitive closure of the latent DAG. Then, in the case of hard interventions, it refines these outputs without requiring any conditional independence test.

0.0.1.0.0.4 General transformations.

In this setting, we do not have any restriction on the transformation from the latent space to the observed space. In this general setting, our contributions are as follows.

• On identifiability, we show that observational data and two distinct hard interventions per node suffice to guarantee perfect identifiability (8). This result generalizes the recent results in the literature in two ways. First, we do not require the commonly adopted faithfulness assumption on latent causal models. Secondly, we assume the learner does not know which pair of environments intervene on the same node. While proving these results, we also establish identifiability of a single variable given two hard interventions only for the said variable (5).

• More importantly, on achievability, we design the first provably correct algorithm that recovers \(\mcG\) and \(\bZ\) perfectly. This algorithm is referred to as Generalized Score-based Causal Latent Estimation via Interventions (GSCALE-I). We note that GSCALE-I requires only the score functions of observed variables as its inputs and computes those of the latent variables by leveraging the Jacobian of the decoders.

• We also establish new results that shed light on the role of observational data. Specifically, when two interventional environments per node are given in pairs, observational data is only needed for recovering the latent DAG (6). Furthermore, we show that observational data can be dispensed with under a weak faithfulness condition on the latent causal model (7). Remarkably, our results are tight in terms of the identifiability objectives, which we will discuss in 0.0.3.4.

• Finally, we show that extrapolation to unseen combinations of single-node interventions can be achieved without performing CRL, just by using score functions of the observed space (0.0.6.4).

Before describing our novel methodology, we recall the general approach to solving CRL. Recovering the latent causal variables hinges on finding the inverse of \(g\) based on the observed data \(\bX\), which facilitates recovering \(\bZ\) via \(\bZ = g^{-1}(\bX)\). In other words, referring to \((g^{-1},g)\) as the true encoder-decoder pair, we search for an encoder that takes observed variables back to the true latent space. A valid encoder should necessarily have an associated decoder to ensure a perfect reconstruction of observed variables from estimated latent variables. However, infinitely many valid encoder-decoder pairs satisfy the reconstruction property. The purpose of using interventions is to inject variations into the observed data, which can help us distinguish among these solutions.

0.0.1.0.0.5 Score-based methodology.

We start by showing that an intervention on a latent node induces changes only in the score function’s coordinates corresponding to the intervened node and its parents. This is because the intervention changes only the causal mechanism (i.e., conditional distribution) of the intervened node, which is a function of the latent node and its parents. Furthermore, when we consider two hard interventions on the same node, such changes will be limited only to the intervened node (parents intact). This implies that score changes of the latent variables \(\bZ\) are generally sparse across different environments. Furthermore, these score changes contain all the information about the latent causal structure. Motivated by these key properties, we formalize a score-based CRL framework, based on which we design provably correct distinct algorithms, LSCALE-I and GSCALE-I, for linear and general transformations, respectively, presented in 0.0.5 and 0.0.6. We briefly describe the key technique in GSCALE-I for general transformation via two interventions. LSCALE-I involves exploiting the linearity of the transformation to achieve identifiability using only one intervention per node.

0.0.1.0.0.6 Algorithm sketch for general transformations.

Consider two interventional environments in which the same node is intervened on. As described in the preceding paragraph, we show that the score functions of the latent variables under these two environments differ only at the coordinate of the intervened node. Subsequently, the key idea of the score-based framework is that tracing these sparse changes in the score functions of the latent variables can guide finding reliable estimates for the inverse of the transformation \(g\), which in turn facilitates estimating \(\bZ\). In particular, we will seek encoders that result in sparse score variations among the estimated latent variables, thereby matching the properties of the true latent variables. To this end, we consider a pair of interventional environments for each node and find an encoder that minimizes the variations in the score function across all pairs. We show that the encoder obtained via this procedure perfectly recovers \(\mcG\) and \(\bZ\). An important process in this methodology is projecting the score changes in observed data to the latent space so that we do not need to estimate latent scores for each possible encoder. We show that this can be done by multiplying the observed score difference by the Jacobian of the decoder associated with the encoder. Therefore, recovering \(Z_i\) is facilitated by solving the following problem: \[\label{eq:intro-OPT-X-verbal} For i-th int. env. pair : \;\; \min \Big\| \E\Big[\big | \text{Jac. of decoder} \times \big(\text{score difference of \bX} \big) \big\|\Big] - \be_i \Big \|^2 \;,\tag{2}\] in which \(\be_i\) denotes \(i\)-th standard basis vector, reflecting the sparsity property. Score differences of observed variables are computed across two environments with the same intervened node. When we have two interventions for each latent variable, we solve this problem for \(i\in[n]\) simultaneously to recover complete \(\bZ\), and subsequently, the graph \(\mcG\). Finally, the minimization is performed over the set of valid encoder-decoder pairs that ensure perfect reconstruction of \(\bX\).

0.0.1.0.0.7 Organization.

The rest of the paper is organized as follows. 0.0.2 provides an overview of the literature with the main focus on CRL via interventions. 0.0.3 provides the preliminaries for formulating the problem and specifies the notations and definitions used throughout the paper. 0.0.4 establishes the properties of score functions under interventions and presents the key lemmas and their implications, which will be used in our CRL framework. We present our CRL algorithms and results for linear transformations in 0.0.5, and general transformations in 0.0.6. Proofs of the building blocks of the results are presented in the main body of the paper, and the rest of the proofs are provided in the appendices. In 0.0.7, we empirically assess the performance of the proposed CRL algorithms for recovering the latent causal variables and the latent causal graph on both structured synthetic data, biological data, and image data. 0.0.8 concludes the paper with a discussion of the results and the future directions.

0.0.2 Related Work↩︎

In this paper, we study CRL in different interventional environments where the interventions act on the latent space. We first provide an overview of the literature that investigates CRL from interventional data, with the main results of the most closely related work summarized in 1 and 2. Then, we discuss the other relevant lines of work.

0.0.2.0.0.1 Interventional causal representation learning.

The majority of the rapidly growing literature on CRL from interventions focuses on parametric settings, i.e., a parametric form is assumed for the latent model, the transformation, or both. Among the related studies, that in [7] mainly considers polynomial transformations without restrictions on the latent causal model and shows identifiability under deterministic do interventions. It also establishes identifiability under general transformations, albeit requiring a combinatorial number of do interventions. Finally, it shows identifiability under soft interventions with independent support assumptions on latent variables. The study in [8] considers a linear latent model with a linear mapping to observations and proves identifiability under hard interventions. It also shows the impossibility of perfect identifiability under soft interventions and proves identifiability up to ancestors. The study in [9] focuses on linear Gaussian latent models and extends the results in [8] to prove identifiability for general transformations. [10] consider polynomial transformations under nonlinearity assumptions on latent models and prove identifiability up to ancestors under soft interventions. If the latent graph is restricted to polytrees, they further prove perfect identifiability. Identifying the non-intervened variables from the intervened variables by using single-node and multi-node soft interventions is studied in [11]. It also considers a new setting in which the latent DAG can change across data points and relies on the support invariance of non-intervened variables to identify them from the rest. The study in [12] considers a nonlinear latent model under linear transformation and uses multi-target do interventions to prove identifiability under certain sparsity assumptions. Linear transformation and linear non-Gaussian latent models are studied in [13], establishing identifiability up to surrounding parents using soft interventions. It also establishes sufficient conditions for multi-target soft interventions to ensure identifiability. The study in [14] takes a different approach and considers the task of intervention extrapolation. In this formulation, interventions are applied to exogenous action variables (e.g., instrumental variables) that linearly affect the latent variables.

On the fully nonparametric setting, [15] provide the most closely related identifiability results to ours. Specifically, they show that two coupled hard interventions per node suffice for identifiability under the faithfulness assumption on latent causal models. Our results have two significant differences: 1) We address achievability via a provably correct algorithm whereas [15] focus mainly on identifiability (e.g., no algorithm for recovery of the latent variables), 2) we dispense with the restrictive assumptions on identifiability results, namely, we do not require to know which two environments share the same intervention target (hence, uncoupled interventions), and do not require faithfulness on the latent models. While the importance of this relaxation is not fully apparent for single-node interventions, being able to work without knowing intervention targets is critical for investigating CRL under more realistic cases of multi-node interventions. Among the other studies on the nonparametric setting, [13] provide analogous results to [15] by considering two coupled soft interventions and identifying latent variables up to mixing with surrounding variables. [16] consider identifying the latent DAG without recovering latent variables, where it is shown that a restricted class of DAGs can be recovered. Finally, some studies employ stronger supervision signals, such as annotations of the ground-truth causal variables [17], or knowledge of the causal graph to recover the latent variables under hard interventions [18].

0.0.2.0.0.2 Multi-view causal representation learning.

A commonly used form of weak supervision in CRL involves multi-view data, non-i.i.d. samples (or views) are observed, typically generated by the same or closely related realizations of underlying latent variables. Multi-view scenarios can involve various data settings. Earlier studies have mostly considered paired data, in which pre- and post-intervention observations are generated from the same set of latent variables [19]–[21]. The more generalized and relaxed formulation typically involves partial observability, where multiple views are concurrently generated by an overlapping subset of latent variables, such as observing different camera angles of the same scene [22]. [23] generalizes the multi-view approach via a unified framework, which also allows partial observability with nonlinear transforms. We also note that [24] presents an even more general unifying view to interpret different CRL approaches, including interventional and multi-view, as special ways of aligning the representations to known data symmetries via leveraging the invariance principle.

0.0.2.0.0.3 Temporal causal representation learning.

Temporal CRL is particularly motivated by applications in domains where time-series data is readily available, such as robotics and control systems, which has led to a growing interest in this area. Some studies in this category assume that the latent causal variables are mutually independent, meaning there are no instantaneous causal effects between variables at the same time step [25]–[27]. However, more recent studies have relaxed this assumption, allowing instantaneous causal effects when sufficient diversity is present in the observed or interventional data [28]–[30].

0.0.2.0.0.4 Other approaches.

In some other approaches to CRL, [31] consider general mixing functions and general SCMs under sparsity constraints; [32] leverage nonlinear ICA principles to work with nonlinear mixing functions and polynomial latent SCMs; [33] consider disjoint groups of observational variables, with general mixing functions and pairwise latent causal relationships; [34] use score functions for linear mixing functions and nonlinear Gaussian SCMs to derive coarser identifiability results; and [35] study “domains” and interventions in a combined way for non-Markovian causal systems when the graph is known.

0.0.2.0.0.5 Identifiable representation learning.

As a special case of CRL, where the latent variables are independent, there is extensive literature on identifying latent representations. Some representative approaches include leveraging knowledge of the mechanisms that govern the system’s evolution [36] and using weak supervision with auxiliary information [37]. Nonlinear independent component analysis (ICA) also uses side information, in the form of structured time series, to exploit temporal information [38], [39] or knowledge of auxiliary variables that make latent variables conditionally independent [40]–[42]. [33] imposes additional constraints on observational mixing and the causal model to prove identifiability. [43] studies the identifiability of deep generative models without auxiliary information.

0.0.2.0.0.6 Score functions for causal discovery within observed variables.

Score matching has recently gained attention in the causal discovery of observed variables. [44] use score matching to recover nonlinear additive Gaussian noise models. [45] focus on the same setting, recover the full graph from Jacobian scores, and dispense with the computationally expensive pruning stage of the algorithm in [44]. [46] empirically demonstrate the robustness of score matching approaches against the assumption violations in causal discovery. [47] establish bounds on the error rate of score matching-based causal discovery methods. All of these studies are limited to observed causal variables, whereas in our case, we have a causal model in the latent space.

0.0.3 Preliminaries and Definitions↩︎

0.0.3.0.0.1 Notations.

For a vector \(\ba\in\R^{m}\), the \(i\)-the entry is denoted by \(a_i\). Random vectors are denoted by bold upper-case letters and their realizations are denoted by bold lower-case letters, e.g., \(\bX\) and \(\bx\). Matrices are denoted by bold upper-case letters, e.g., \(\bA\), where \(\bA_i\) denotes the \(i\)-th row of \(\bA\) and \(\bA_{i,j}\) denotes the entry at row \(i\) and column \(j\). For matrices \(\bA\) and \(\bB\) with the same shapes, \(\bA \preccurlyeq \bB\) denotes component-wise inequality. We denote the indicator function by \(\mathbb{1}\), and for a matrix \(\bA\in\R^{m\times n}\), we use the convention that \(\mathbb{1}\{\bA\}\in\{0,1\}^{m\times n}\), where the entries are specified by \([\mathbb{1}\{\bA\}]_{i,j}=\mathbb{1}\{\bA_{i,j}\neq 0\}\). For a positive integer \(n\in\N\), we define \([n] \triangleq \{1,\dots,n\}\). The permutation matrix associated with any permutation \(\pi\) of \([n]\) is denoted by \(\bP_{\pi}\), i.e., \([\pi_1 \;\; \pi_2 \; \dots \; \pi_n]^{\top}=\bP_{\pi}\cdot[1 \;\; 2 \; \dots \; n]^{\top}\). The \(n\)-dimensional identity matrix is denoted by \(\bI_{n \times n}\), and the Hadamard product is denoted by \(\odot\). We use \(\image(f)\) to denote the image of the function \(f\). We let \(\be_i\) denote the \(i\)-th standard basis vector of \(\R^n\). Given a function \(f: \R^s \to \R^r\) that has first-order partial derivatives on \(\R^s\), we denote the Jacobian of \(f\) at \(\bz\in\R^s\) by \(J_{f}(\bz) \in \R^{r \times s}\) with entries \([J_{f}(\bz)]_{i,j} = \partial f(\bz)_i / \partial \bz_j\).

0.0.3.1 Latent Causal Structure

Consider latent causal random variables \(\bZ \triangleq [Z_1,\dots,Z_n]^{\top}\). An unknown transformation \(g: \R^n \to \R^d\) generates the observable random variables \(\bX \triangleq [X_1,\dots,X_d]^{\top}\) from the latent variables \(\bZ\) according to: \[\bX = g(\bZ) \; . \label{eq:data-generation-process}\tag{3}\] We assume that \(d \geq n\), and the transformation \(g\) is continuously differentiable and a diffeomorphism onto its image (otherwise, identifiability is ill-posed). We denote the image of \(g\) by \(\mcX \triangleq \image(g)\subseteq \R^d\). The probability density functions (pdfs) of \(\bZ\) and \(\bX\) are denoted by \(p\) and \(p_{\bX}\), respectively. We assume that \(p\) is absolutely continuous with respect to the \(n\)-dimensional Lebesgue measure. Subsequently, \(p_{\bX}\), which is defined on the image manifold \(\image(g)\), is absolutely continuous with respect to the \(n\)-dimensional Hausdorff measure rather than \(d\)-dimensional Lebesgue measure⁴. The distribution of latent variables \(\bZ\) factorizes with respect to a DAG that consists of \(n\) nodes and is denoted by \(\mcG\). Node \(i\in[n]\) of \(\mcG\) represents \(Z_i\) and \(p\) factorizes according to: \[p(\bz) = \prod_{i=1}^{n} p_i(z_i \med \bz_{\Pa(i)}) \;, \label{eq:pz95factorized}\tag{4}\] where \(\Pa(i)\) denotes the set of parents of node \(i\) and \(p_i(z_i \med \bz_{\Pa(i)})\) is the conditional pdf of \(z_i\) given the variables of its parents. We use \(\Ch(i)\), \(\An(i)\), and \(\De(i)\) to denote the children, ancestors, and descendants of node \(i\), respectively. Accordingly, for each node \(i\in[n]\) we also define \[\begin{align} \label{eq:graph-notations} \overline{\Pa}(i) \triangleq \Pa(i) \cup \{i\} , \;\; \overline{\Ch}(i)\triangleq \Ch(i) \cup \{i\} , \;\; \overline{\An}(i)\triangleq \An(i) \cup \{i\} , \; and \; \overline{\De}(i)\triangleq \De(i) \cup \{i\}\;. \end{align}\tag{5}\] We denote the transitive closure and transitive reduction of \(\mcG\) by \(\mcG_{\rm tc}\) and \(\mcG_{\rm tr}\), respectively⁵. The parental relationships in these graphs are denoted by \(\Pa_{\rm tc}\) and \(\Pa_{\rm tr}\), and other graphical relationships are denoted similarly. Based on the modularity property, a change in the causal mechanism of node \(i\) does not affect those of the other nodes. We also assume that all conditional pdfs \(\{p_i(z_i \mid \bz_{\Pa(i)}) : i \in[n]\}\) are continuously differentiable with respect to all variables and \(p(\bz) \neq 0\) for all \(\bz \in \R^n\). We consider the general structural causal models (SCMs) based on which for each \(i\in[n]\), \[\begin{align} \label{eq:general-SCM} Z_i = f_{i}(\bZ_{\Pa(i)},N_{i}) \;, \end{align}\tag{6}\] where \(\{f_{i} : i \in [n]\}\) are general functions that capture the dependence of node \(i\) on its parents and \(\{N_{i} : i \in [n]\}\) account for the exogenous noise terms that we assume to have pdfs with full support. We specialize some of the results to additive noise SCMs, in which 6 becomes \[\begin{align} \label{eq:additive-SCM} Z_i = f_{i}(\bZ_{\Pa(i)}) + N_{i}\;. \end{align}\tag{7}\] Next, we provide several definitions that we will use frequently throughout the paper to formalize the framework and analyze it.

In this paper, without loss of generality, we assume that \((1,\dots, n)\) is a valid causal order. We also define a graphical notion that will be useful for presenting our results and analysis on CRL under a linear transformation.

The intuition behind surrounded nodes is that, the effect of \(i\) on its children \(\Ch(i)\) can be dominated by the effect of its surrounding node \(j\).⁷ Specifically, any effect of \(i\) on a node \(k \in \Ch(i)\) can also be interpreted as the effect of node \(j\) since \(k \in \Ch(j)\) as well. This effect causes ambiguities in the recovery of \(Z_i\) in the case of soft interventions.

0.0.3.2 Score Functions

The score function associated with a pdf is defined as the gradient of its logarithm. The score function associated with \(p\) is denoted by \[\bss(\bz) \triangleq \nabla_{\bz} \log p(\bz) \;. \label{eq:sz-def}\tag{9}\] Noting the connection \(\bX = g(\bZ)\), the density of \(\bX\) under \(\mcE^0\), denoted by \(p_{\bX}\), is supported on an \(n\)-dimensional manifold \(\mcX\) embedded in \(\R^{d}\). Hence, specifying the score function of \(\bX\) requires notions from differential geometry. For this purpose, we denote the tangent space of the manifold \(\mcX\) at point \(\bx \in \mcX\) by \(T_{\bx} \mcX\). Tangent vectors \(\bv \in T_x \mcX\) are equivalence classes of continuously differentiable curves \(\gamma \colon (-1,1) \to \mcX \subseteq \R^d\) with \(\gamma(0) = \bx\) and \(\gamma'(0) = \bv\). Furthermore, given a function \(f \colon \mcX \to \R\), denote its directional derivative at point \(\bx \in \mcX\) along a tangent vector \(\bv \in T_{\bx} \mcX\) by \(D_{\bv} f(\bx)\), which is defined as \[\require{physics} D_{\bv} f(\bx) \triangleq \left. \frac{\dd}{\dd t} (f \circ \gamma)(t) \right|_{t = 0} \;,\] for any curve \(\gamma\) in equivalence class \(v\). The differential of \(f\) at point \(\bx \in \mcX\), denoted by \(\require{physics} \dd f_{\bx}\), is the linear operator mapping tangent vector \(\bv \in T_{\bx} \mcX\) to \(D_{\bv} f_{\bx}\) [49], i.e., \[\require{physics} \dd f_{x} \colon T_{\bx} \mcX \ni v \mapsto D_{\bv} f(\bx) \in \R \;. \label{eq:differential-def}\tag{10}\] Let \(\bB \in \R^{d \times n}\) be a matrix for which the columns of \(\bB\) form an orthonormal basis for \(T_{\bx} \mcX\). Denote the directional derivative of \(f\) along the \(i\)-th column of \(\bB\) by \(D_i f\) for all \(i \in [n]\). Then, the differential operator can be expressed by the vector \[D f_{\bx} \triangleq \bB \cdot \big[ D_1 f_{\bx} \, \ldots \, D_n f_{\bx} \big]^{\top} \in \R^{d} \;, \label{eq:Df-defn}\tag{11}\] such that \[\require{physics} \dd f_{\bx}(\bv) = \bv^{\top} \cdot D f_{\bx} \;, \qquad \forall \bx \in \mcX \;, \quad \forall \bv \in T_{\bx} \mcX \;.\] Note that the differential operator \(\require{physics} \dd f\) is a generalization of the gradient. Hence, we can generalize the definition of the score function using the differential operator by setting \(f\) to the logarithm of pdf. Therefore, the score function of \(\bX\) under \(\mcE^0\) is specified as follows: \[\bss_{\bX}(\bx) \triangleq D \log p_{\bX}(\bx) \;, \qquad \forall \bx \in \mcX \;. \label{eq:sx-def}\tag{12}\]

0.0.3.3 Intervention Mechanisms

We consider two types of interventions. A soft intervention on node \(i\) (also referred to as imperfect intervention in literature), changes the conditional distribution \(p_i(z_i \med \bz_{\Pa(i)})\) to a distinct conditional distribution, which we denote by \(q_i(z_i \med \bz_{\Pa(i)})\). A soft intervention does not necessarily remove the functional dependence of an intervened node on its parents but rather alters it to a different mechanism. A stochastic hard intervention on node \(i\) (also referred to as perfect intervention) is stricter than a soft intervention and removes the edges incident on \(i\). A hard intervention on node changes \(p_i(z_i \med \bz_{\Pa(i)})\) to \(q_i(z_i)\) that emphasizes the lack of dependence of \(z_i\) on \(\bz_{\Pa(i)}\). Finally, we note that in some settings, we assume two hard interventions per node, in which case the two hard interventional mechanisms for node \(i\) are denoted by two distinct pdfs \(q_i(z_i)\) and \(\tilde{q}_i(z_i)\).

0.0.3.3.0.1 Interventional environments.

We consider atomic interventional environments in which each environment one node is intervened in, as is customary in the closely related CRL literature [7]–[9]. In some settings (linear transformation), we will have one interventional environment per node and denote the interventional environments by \(\mcE\triangleq \{\mcE^1,\dots,\mcE^n\}\), where we call \(\mcE\) the atomic environment set. We denote the node intervened in environment \(\mcE^m\) by \(I^m \in [n]\). For other settings (general transformation), we will have two interventional environments per node and denote the second atomic environment set by \(\tilde{\mcE}=\{\tilde{\mcE}^1,\dots,\tilde{\mcE}^n\}\). Similarly, we denote the intervened node in \(\tilde{\mcE}^m\) by \(\tilde{I}^m\) for each \(m\in[n]\). We assume that node-environment pairs are unspecified, i.e., the ordered intervention sets \(\mcI\triangleq (I^1,\dots,I^n)\) and \(\tilde{\mcI} \triangleq (\tilde{I}^1,\dots,\tilde{I}^n)\) are two unknown permutations of \([n]\). We also adopt the convention that \(\mcE^{0}\) is the observational environment and \(I^{0} \triangleq \emptyset\). Next, we define the notion of coupling between the environment sets \(\mcE\) and \(\tilde{\mcE}\).

Next, we define \(p^m\) as the pdf of \(\bZ\) in environment \(\mcE^m\). Hence, under soft and hard intervention for each \(m \in [n]\), \(p^m\) can be factorized as follows. \[\begin{align} soft int. in\mcE^m: & \quad p^m(\bz) = q_\ell(z_{\ell} \med\bz_{\Pa(\ell)}) \prod_{i \neq \ell} p_i(z_i \med \bz_{\Pa(i)}) \;, && \;\; where \quad \ell=I^m \;,\tag{13} \\ hard int. in\mcE^m: & \quad p^m(\bz) = q_{\ell}(z_{\ell}) \prod_{i \neq \ell} p_i(z_i \med \bz_{\Pa(i)})\;, && \;\; where \quad \ell=I^m \;. \tag{14} \end{align}\] Similarly, we define \(\tilde{p}^m\) as the pdf of \(\bZ\) in \(\tilde{\mcE}^m\), which can be factorized similarly to 14 with \(\tilde{q}_{\ell}\) replaced with \(q_{\ell}\). Hence, the score functions associated with \(p^m\) and \(\tilde{p}^m\) are specified as follows. \[\bss^{m}(\bz) \triangleq \nabla_{\bz} \log p^{m}(\bz) \;, \quad and \quad \tilde{\bss}^{m}(\bz) \triangleq \nabla_{\bz} \log \tilde{p}^{m}(\bz) \;. \label{eq:sz-defs}\tag{15}\] We denote the score functions of the observed variables \(\bX\) under \(\mcE^m\) and \(\tilde{\mcE}^m\) by \(\bss^{m}_{\bX}\) and \(\tilde{\bss}^{m}_{\bX}\), respectively. Note that the score functions change across different environments, which is induced by the changes in the distribution of \(\bZ\). Specifically, following 4 and 13 , the latent scores \(\bss(\bz)\) and \(\bss^{m}(\bz)\) are decomposed as \[\begin{align} \bss(\bz) &= \nabla_{\bz} \log p_{\ell}(z_{\ell} \med\bz_{\Pa(\ell)}) + \sum_{i \neq \ell} \nabla_{\bz} \log p_i(z_i \med \bz_{\Pa(i)}) \;, \tag{16} \\ and \quad \bss^{m}(\bz) &= \nabla_{\bz} \log q_{\ell}(z_{\ell} \med\bz_{\Pa(\ell)}) + \sum_{i \neq \ell} \nabla_{\bz} \log p_i(z_i \med \bz_{\Pa(i)}) \;. \tag{17} \end{align}\] where \(\ell = I^m\). Hence, \(\bss(\bz)\) and \(\bss^{m}(\bz)\) differ in only the causal mechanism of the intervened node in environment \(\mcE^m\). In 0.0.4, we investigate these discrepancies between \(\bss\) and \(\bss^{m}\) (or \(\tilde{\bss}^{m}\)) and characterize the relationship between the scores in the observational and interventional environments.

0.0.3.4 Identifiability and Achievability Objectives

The objective of CRL is to use observations \(\bX\) generated by the observational and interventional environments and estimate the true latent variables \(\bZ\) and causal relations among them captured by \(\mcG\). The first objective is identifiability, which pertains to determining algorithm-agnostic sufficient conditions under which \(\bZ\) and \(\mcG\) can be recovered uniquely up to a permutation and element-wise transform, which is the strongest form of recovery in CRL from interventions as shown in [15]. The second objective is achievability, which refers to designing algorithms that are amenable to practical implementation and generate provably correct estimates for \(\bZ\) and \(\mcG\), foreseen by the identifiability guarantees. In this subsection, we provide the definitions needed for formalizing these objectives.

We denote a generic estimator of \(\bZ\) given \(\bX\) by \(\hat{\bZ}(\bX):\R^d\to\R^n\). We also consider a generic estimate of \(\mcG\) denoted by \(\hat{\mcG}\). To assess the fidelity of the estimates \(\hat{\bZ}(\bX)\) and \(\hat{\mcG}\) with respect to the ground truth \(\bZ\) and \(\mcG\), we provide the following identifiability measures, which will be achieved when we have a complete set of atomic interventions.

We note that scaling consistency is a special case of component-wise latent recovery in which the diffeomorphism is restricted to scaling. Next, we provide node-level partial identifiability definitions which measure the recovery level of a single latent variable. These will be useful to assess the identifiability guarantees of the algorithms under an incomplete set of interventions.

0.0.3.4.0.1 Tightness of the identifiability definitions.

For general transformations, [15] show that identifiability up to component-wise diffeomorphisms is the best possible result under interventions, without additional assumptions. For linear transformations, [8] and [9] show that under hard interventions, identifiability up to scaling consistency and permutation ambiguity is the best achievable outcome without further information on \(\bZ\). Furthermore, for soft interventions, transitive closure recovery is the best latent graph identifiability result without additional assumptions on the causal model  [8]. Finally, [13] show that, under some non-degeneracy assumptions, consistency up to mixing with surrounding parents is the optimal result for soft interventions. Therefore, given a complete set of atomic interventions, the identifiability definitions presented in this section represent the ultimate objectives when using interventions.

0.0.3.5 Algorithm-related Definitions

For formalizing the achievability results and designing the associated algorithms, generating the estimates \(\hat{\bZ}(\bX)\) and \(\hat{\mcG}\) is facilitated by estimating the inverse of \(g\) based on the observed data \(\bX\). Specifically, an estimate of \(g^{-1}\), where \(g^{-1}\) denotes the inverse of \(g\), facilitates recovering \(\bZ\) via \(\bZ = g^{-1}(\bX)\) . Throughout the rest of this paper, we refer to \(g^{-1}\) as the true encoder. To formalize the procedures of estimating \(g^{-1}\), we define \(\mcH\) as the set of possible valid encoders, i.e., candidates for \(g^{-1}\). A function \(h\in\mcH\) can be such a candidate if it is invertible, that is, there exists an associated decoder \(h^{-1}\) such that \((h^{-1} \circ h)(\bX) = \bX\). Hence, the set of valid encoders is specified by \[\label{eq:mcH} \mcH \triangleq \{ h: \mcX \to \R^n \;\colon\; \exists h^{-1}: \R^n \to \R^d \;\; such that \;\; (h^{-1} \circ h)(\bx) = \bx \;, \;\; \forall \bx \in \mcX\} \;.\tag{21}\] Next, corresponding to any pair of observation \(\bX\) and valid encoder \(h\in\mcH\), we define \(\hat{\bZ}(\bX;h)\) as an auxiliary estimate of \(\bZ\) generated by applying the valid encoder \(h\) on \(\bX\), i.e., \[\label{eq:encoder} \hat{\bZ}(\bX;h) \triangleq h(\bX) = (h \circ g) (\bZ) \;, \qquad \forall \; h\in\mcH \;.\tag{22}\] The estimate \(\hat{\bZ}(\bX;h)\) inherits its randomness from \(\bX\), and its statistical model is governed by that of \(\bZ\) under the chosen \(h\). To emphasize the dependence on \(h\), we denote the score functions associated with the pdfs of \(\hat{\bZ}(\bX;h)\) under environments \(\mcE^{0}\), \(\mcE^m\), and \(\tilde{\mcE}^m\), respectively, by \[\bss_{\hat{\bZ}}(\cdot ;h)\;, \quad \bss_{\hat{\bZ}}^{m}(\cdot ;h), \quad and \quad \tilde{\bss}_{\hat{\bZ}}^{m}(\cdot;h)\;.\] We will be addressing both general and linear transformations \(g\). In the linear transformation setting, the true linear transformation \(g\) is denoted by the matrix \(\bG \in \R^{d \times n}\). Accordingly, we denote a valid linear encoder by \(\bH \in \R^{n \times d}\). For a given valid encoder \(\bH\), the associated valid decoder is given by its Moore-Penrose inverse, i.e., \(\bH^{\dagger} \triangleq \bH^{\top} \cdot (\bH \cdot \bH^{\top})^{-1}\).

0.0.4 Properties of Score Functions under Interventions↩︎

Score functions and their variations across different interventional environments play pivotal roles in our approach to identifying latent representations. In this section, we present the key properties of the score functions that will be leveraged in 0.0.5 and 0.0.6 to construct identifiability results, along with the corresponding algorithms.

We first investigate score variations across pairs of environments, such as the observational environment and an interventional one (under both soft and hard atomic interventions) or two interventional environments, either coupled or uncoupled. The key insight is that an intervention causes changes in only certain coordinates of the score function, as indicated by the sparse changes in the decompositions of the score functions in 16 and 17 . These sparse changes further reflect the graph structure in the score differences. For instance, for a single-node intervention \(I^m = \{i\}\) in \(m\)-th environment, \(j\)-th coordinate of the score difference \(\bss - \bss^{m}\) becomes \[\label{eq:def-score-diff-z-int-obs-single} \big[\bss(\bz) - \bss^{m}(\bz)\big]_j = \frac{\partial}{\partial z_j} \log p_{i}(z_i \mid \bz_{\Pa(i))}) - \frac{\partial}{\partial z_j} \log q_{i}(z_i \mid \bz_{\Pa(i)}) \;,\tag{23}\] which is zero if \(j \notin \Paplus(i)\). The following lemma formalizes this property for all relevant cases.

See 1.0.1.1.

1 provides the necessary and sufficient conditions for the invariance of the coordinates of the score functions of latent variables. The core idea of the score-based framework is that tracing these sparse changes in the score functions of the latent variables guides finding reliable estimates for the inverse of the transformation \(g\), which in turn facilitates estimating \(\bZ\). Intuitively, we will look for the encoders \(h \in \mcH\) such that the variations between the score estimates \(\bss_{\hat{\bZ}}(\hat{\bz};h)\), \(\bss_{\hat{\bZ}}^{m}(\hat{\bz};h)\), and \(\tilde{\bss}_{\hat{\bZ}}^{m}(\hat{\bz};h)\) will be similar to the true score variations given by 1. However, the scores of the latent variables are not directly accessible. To circumvent this, we need to understand the connection between score functions of \(\bX\) and \(\bZ\). In the following lemma, we leverage the change of variables formula for injective mappings and establish this relationship for any injective mapping \(f\) from latent to observed space.

See 1.0.1.2.
For CRL under linear transformations in 0.0.5, we use the following corollary of 2.

For CRL under general transformations in 0.0.6, we customize 2 as follows. Consider a candidate encoder \(h \in \mcH\). Recall that \(\hat{\bZ}(\bX;h) = h(\bX) = (h \circ g)(\bZ)\). Then, by letting \(f = (h \circ g)\), 2 gives \[\begin{align} between \mcE^{0}and\mcE^m&: \quad \bss_{\hat{\bZ}}(\hat{\bz};h) - \bss_{\hat{\bZ}}^{m}(\hat{\bz};h) = J_{f}^{-\top}(\bz) \cdot \big[\bss(\bz) - \bss^{m}(\bz)\big] \;, \tag{28} \\ between \mcE^{0}and\tilde{\mcE}^m&: \quad \bss_{\hat{\bZ}}(\hat{\bz};h) - \tilde{\bss}_{\hat{\bZ}}^{m}(\hat{\bz};h) = J_{f}^{-\top}(\bz) \cdot \big[\bss(\bz) - \tilde{\bss}^{m}(\bz)\big] \;, \tag{29} \\ between \mcE^mand\tilde{\mcE}^m&: \quad \bss_{\hat{\bZ}}^{m}(\hat{\bz};h) - \tilde{\bss}_{\hat{\bZ}}^{m}(\hat{\bz};h) = J_{f}^{-\top}(\bz) \cdot \big[\bss^{m}(\bz) - \tilde{\bss}^{m}(\bz)\big] \;. \tag{30} \end{align}\]

0.0.5 CRL under Linear Transformations↩︎

In this section, we consider CRL under linear transformation, in which the general transformation model in 3 becomes: \[\label{eq:data-generation-process-linear} \bX = \bG \cdot \bZ \;,\tag{31}\] where \(\bG \in \R^{d \times n}\) is an unknown full-rank matrix mapping the latent variables to the observed ones. We will present steps for leveraging the properties of the score functions presented in 0.0.4 to design an algorithm that identifies the true encoder and recovers the true causal representations. The algorithm is referred to as Linear Score-based Causal Latent Estimation via Interventions (LSCALE-I). The theoretical guarantees associated with the algorithm steps will also be presented, serving as constructive proof of identifiability.

LSCALE-I consists of three stages outlined in 3. In the first stage, we will learn an encoder \(\hat{\bH}\), row by row, where each row is an estimate of a row of the true encoder \(\bGinv\). This encoder estimate will recover each \(Z_i\) up to mixing with \(\bZ_{\Pa(i)}\). In the second stage, we will use the estimated latent variables \(\hat{\bH}\) to obtain the transitive closure of the latent graph \(\mcG\). Finally, as an optional third stage for hard interventions, we will leverage the independence statements implied by the hard interventions to refine our encoder estimate. We will demonstrate that this refinement yields scaling consistency and perfect DAG recovery. We will also discuss the relevance and distinctions of our results vis-à-vis the results in the existing literature.

0.0.5.0.0.1 Statistical diversity.

Similarly to the existing identifiability results, it is necessary to have some regularity conditions ⁸ on the probability distributions of observational and interventional environments. We adopt the following assumption for the case of linear transformations.

Essentially, \(\frac{\partial}{\partial z_k} \big(\log\frac{p_i(z_i \mid \bz_{\Pa(i)})}{q_i(z_i \mid \bz_{\Pa(i)})} \big)\) captures the effect of intervening on node \(i\) on the score associated with node \(k\). Therefore, 1 ensures that an intervention sufficiently differentiates the target variable and its parents in the score function. We note that this is a very mild assumption and holds for a wide range of commonly used models (including additive noise models under hard interventions) and is discussed in more detail in 0.0.5.5.

Before we explain the details of LSCALE-I and establish identifiability results in the following subsections, we first present the rationale and derive the partial identifiability result, which will serve as the building block. In the rest of the paper, for brevity, we denote the score difference functions for each \(m\in[n]\) by \[\label{eq:def-score-diff-z-int-obs} \bsd_{\bZ}^{m}(\bz) \triangleq \bss^{m}(\bz) - \bss(\bz) \quad and \quad \bsd_{\bX}^{m}(\bx) \triangleq \bss_{\bX}^{m}(\bx) - \bss_{\bX}(\bx) \;.\tag{32}\]

0.0.5.1 Rationale for LSCALE-I and Partial Identifiability

To construct an algorithm that achieves identifiability under linear transformations, we first analyze the score differences for the linear transformation. Note that for linear transformation \(\bX = \bG \cdot \bZ\), Jacobian \(J_{\bG}(\bz)\) is independent of \(\bz\) and equal to matrix \(\bG\). Also, we have \(\bz = \bGinv \cdot \bx\) for all \(\bx \in \colspace(\bG)\). Then, using 1 we have \[\begin{align} \label{eq:sx-from-sz-linear} \bsd_{\bX}^{m}(\bx) = [\bGinv]^{\top} \cdot \bsd_{\bZ}^m(\bz) \;. \end{align}\tag{33}\] Let us denote the correlation matrices of \(\bsd_{\bX}^{m}(\bx)\) and \(\bsd_{\bZ}^{m}(\bz)\), respectively, by \[\bR_{\bX}^{m} \triangleq \E \big[\bsd_{\bX}^m(\bx) \cdot [\bsd_{\bX}^m(\bx)]^{\top} \big] \; \quad and \quad \bR_{\bZ}^{m} \triangleq \E \big[\bsd_{\bZ}^m(\bz) \cdot [\bsd_{\bZ}^m(\bz)]^{\top} \big] \;. \label{eq:rx-rz-defn}\tag{34}\] Note that 33 implies that \(\bR_{\bX}^{m} = [\bGinv]^{\top} \cdot \bR_{\bZ}^{m} \cdot \bGinv\). The next lemma specifies that the structure of the column space of correlation matrix \(\bR_{\bX}^{m}\) is heavily constrained by the true encoder \(\bGinv\), the graph \(\mcG\), and the intervention target \(I^m\) in \(m\)-th environment.

3 crucially implies that we can achieve partial identifiability of a single latent variable given that we have an environment in which the said variable is intervened.

0.0.5.2 Identifiability via Soft Interventions

After obtaining the partial identifiability result from a single interventional environment, our goal is to construct an algorithm that uses the complete set of interventional environments \(\mcE = \{\mcE^1,\dots,\mcE^n\}\) with targets \(\cup_{m \in [n]} I^m = [n]\) to form estimates of the true encoder \(\bGinv\) and the latent graph \(\mcG\).

0.0.5.2.0.1 Encoder estimation.

In Stage L1 of 3, we form an encoder estimate \(\hat{\bH}\) via using the column spaces of \(\{\bR_{\bX}^{m} : m \in [n]\}\). Specifically, for each \(m \in [n]\) we select a vector randomly from \(\colspace(\bR_{\bX}^{m})\) and assign its transpose to \(m\)-th row of \(\hat{\bH}\).

0.0.5.2.0.2 Latent graph estimation.

Next, we use the encoder estimate \(\hat{\bH}\) to recover the transitive closure of the latent graph \(\mcG\). Motivated by 1(ii), we form \(\hat{\mcG}\) with parent sets \[\label{eq:linear-soft-parent-set} \hat{\Pa}(m) \triangleq \Big\{i \neq m: \E\Big[\big| \bss_{\hat{\bZ}}(\hat{\bZ}; \hat{\bH}) - \bss_{\hat{\bZ}}^{m}(\hat{\bZ}; \hat{\bH}) \big|_i \Big] \neq 0 \Big\} \;.\tag{37}\] Using the result that \(\hat{\bZ}(\bX; \hat{\bH})\) recovers \(\bZ\) up to mixing with parents, we show that the transitive closures of the estimate \(\hat{\mcG}\) and the true graph \(\mcG\) are the same under a graph isomorphism.

We note that the existing literature on CRL with linear transformations and one soft intervention requires the latent causal model to be either linear Gaussian [8], [9] or satisfy nonlinearity conditions [10]⁹. In contrast, 2 achieves identifiability guarantees for soft interventions without imposing any restrictions (distributional or structural) on the latent causal model. Furthermore, recovering \(\bZ\) up to mixing with parents improves upon the existing guarantees of recovering up to mixing with ancestors [8], [9]. Finally, we note that \(\mcG\) cannot be identified beyond its transitive closure when using soft interventions without making additional assumptions [8]. Therefore, the graph identifiability result in 2 is tight.

0.0.5.3 Identifiability via Hard Interventions

Next, we investigate hard interventions (i.e., perfect interventions) for general latent causal models. Hard interventions are special cases of soft interventions, in which the intervened node becomes functionally independent of its parents. Consequently, the identifiability guarantees for hard interventions are usually stronger, as we show in this section. Our analysis builds on leveraging the following property, which is exclusive to hard interventions.

The Markov property readily implies this property, that is, each variable in a DAG is independent of its non-descendants, given its parents. When node \(\ell\) is hard-intervened, it has no parents, and the statement follows directly. Motivated by this property, we want to ensure that the estimated latent variables conform to [fact:hard-indep]. To this end, we aim to update rows of the encoder estimate \(\hat{\bH}\) to remove the effects of \(\bZ_{\Pa(i)}\) from the estimate of \(Z_i\). We achieve this objective as follows.

Note that [fact:hard-indep] provides us with random variable pairs that are supposed to be independent, and the covariance of two independent random variables is necessarily zero. In Stage L3 of 3, for each row \(\hat{\bH}_m\), we consider \(\hat{\bZ}(\bX; \hat{\bH})\) in the corresponding environment in which \(\hat{Z}_m\) is intervened, and compute an unmixing vector \(\bu\) via a linear minimum mean square error (MMSE) estimator, \[\bu = \Cov(\hat{Z}_m, \hat{\bZ}_{\hat{\Pa}(m)}) \cdot \big[\Cov(\hat{\bZ}_{\hat{\Pa}(m)}) \big]^{-1} \;.\] Then, we update the corresponding row of the encoder as \(\hat{\bH}_m \gets \hat{\bH}_m - \bu \cdot \hat{\bH}_{\hat{\Pa}(m)}\). Finally, for the perfect DAG recovery, we compute the latent score differences under the final encoder estimate and reconstruct the graph \(\hat{\mcG}\) with parent sets \[\label{eq:linear-hard-parent-set} \hat{\Pa}(m) \triangleq \Big\{ i \neq m : \E\Big[\big| \bsd_{\hat{\bZ}}^{m}(\hat{\bZ}; \hat{\bH})\big|_i \Big] \neq 0 \Big\} \;.\tag{40}\]

See 1.0.2.2.

Similar to the restrictions in the existing results for soft interventions, the identifiability results for hard interventions in the existing literature restrict the latent causal model, e.g., to linear causal models [8], [9]. In contrast, 3 does not impose any restriction on the latent causal model and shows that one stochastic hard intervention per node is sufficient for the identifiability of general latent causal models.

Finally, we comment on the differences between the identifiability results for soft and hard interventions under missing intervention targets and discuss the identifiability of latent subgraphs under special cases.

0.0.5.4 Identifiability up to Surrounding Parents

For linear transformations, we finally investigate the conditions under which soft interventions are guaranteed to achieve identifiability results stronger than transitive closure and mixing up to parents. In particular, we specify one condition on the rank of the score function differences \(\{\bsd_{\bZ}^{m}: m \in [n]\}\), which is formalized next.

For insight into this assumption, it can be readily verified that for linear Gaussian latent models \(\rank(\bR_{\bZ}^{m}) \leq 2\) and on the other hand, for sufficiently nonlinear causal models, \(\rank(\bR_{\bZ}^{m})\) is \(\Paplus(I^m)\). This assumption is stronger than 1 since it implies that the effects of an intervention on all parents of the target variable are different. We will provide more discussions on this assumption in 0.0.5.5. Under this condition, 4 identifies the latent graph and the latent variables with the following key intuitions.

0.0.5.4.0.1 Latent graph estimation.

Recall that \(\bR_{\bX}^{m} = [\bGinv]^{\top} \cdot \bR_{\bZ}^{m} \cdot \bGinv\). Since \(\bGinv\) has full row-rank \(n\), 2 implies that \(\rank(\bR_{\bX}^{m}) = |\Paplus(I^m)|\) for all \(m\in[n]\), which further implies that 3 becomes an equality as well, i.e., \[\colspace(\bR_{\bX}^{m}) = \spann\big\{[\bGinv_i]^{\top} \,:\, i \in \Paplus(I^m)\big\} \;.\] Then, for any \(t, k \in [n]\), we have \[\begin{align} \colspace(\bR_{\bX}^{t}) \cap \colspace(\bR_{\bX}^{k}) &= \spann\big\{[\bGinv_i]^{\top} \,:\, i \in \Paplus(I^t) \cap \Paplus(I^k) \big\} \;, \tag{41} \\ \dim\big( \colspace(\bR_{\bX}^{t}) \cap \colspace(\bR_{\bX}^{k})\big) &= |\Paplus(I^t) \cap \Paplus(I^k)| \\ &= |\Pa(I^t) \cap \Pa(I^k)| + \mathbb{1}\big\{\{I^t \in \Pa(I^k)\} \lor \{I^k \in \Pa(I^t) \} \big\} \;. \tag{42} \end{align}\] In 4, we first obtain a topological order of the true DAG using the graph estimate from 3 and initialize the new estimate \(\hat{\mcG}\) with the empty graph. Consider a pair \((t,k)\) such that \(I^t \in \Pa(I^k)\). If \(\hat{\Pa}(k)\) does not contain node \(t\) in \(\hat{\mcG}\), then by 42 we have \(\dim\big( \colspace(\bR_{\bX}^{t}) \cap \colspace(\bR_{\bX}^{k})\big) > |\hat{\Pa}(t) \cap \hat{\Pa}(k)|\). We use this observation to sequentially identify the parents of each node, starting from the root node(s) and gradually advancing to the leaf node(s). In the end, resulting \(\hat{\mcG}\) equals \(\mcG\) under permutation \(\mcI\).

0.0.5.4.0.2 Encoder estimation.

In 4, we refine the encoder estimation step of 3 as follows. Extending 41 and using the graph recovery result, we have \[\begin{align} \textstyle\bigcap\limits_{k \, \in \, \hat{\Ch}(m) \, \cup \, m } \colspace(\bR_{\bX}^{k})&= \spann\big\{[\bGinv_i]^{\top} \,:\, i \in \textstyle\bigcap\limits_{k \, : \, I^k \in \Chplus(I^m)} \Paplus(I^k) \big\} \;. \label{eq:dx-col-intersection-full-rank} \end{align}\tag{43}\] Now, consider \(I^t \notin \sur(I^m)\). By the definition of surrounding parents, there exists \(I^j \in \Chplus(I^m)\) such that \(I^j \notin \Ch(I^t)\). Then, \(I^j\) is not in the intersection of the parent sets in the right-hand-side of 43 , which means \([\bGinv_j]^{\top}\) is not in the column space of the left-hand-side of 43 . Therefore, by choosing \(\hat{\bH}_m\) from this column space, we obtain for all \(m \in [n]\), \[\begin{align} [\hat{\bZ}(\bX; \hat{\bH})]_m = \hat{\bH}_{m} \cdot \bX &= c_{\ell} \cdot Z_{\ell} + \sum_{i \in \sur(\ell)} c_i \cdot Z_i \;, \qquad where \; \ell = I^m \;, \\ \hat{\bZ}(\bX; \hat{\bH}) &= \bP_{\mcI} \cdot \bC_{\rm pa} \cdot \bZ \;, \end{align}\] where \(\bP_{\mcI}\) is the permutation matrix of the intervention order \((I^1,\dots,I^n)\), \(\bC_{\rm sur}\) has nonzero diagonal entries and satisfies \(\big[\bC_{\rm sur}\big]_{i,j} = 0\) for all \(j \notin \overline{\sur}(i)\).

4 has two important implications. First, the latent DAG can be identified using only soft interventions under mild nonlinearity assumptions on the latent causal model. To our knowledge, this is the first result in the literature for fully recovering latent DAG with soft interventions without restricting the graphical structure, e.g., [10] require linear faithfulness assumption to achieve similar results, which is only shown to hold for nonlinear latent models with polytree structure. Secondly, the estimated latent variables reveal the true conditional independence relationships since they satisfy the Markov property with respect to the estimated latent graph, which is isomorphic to the latent graph. Recalling that the motivation of CRL is to learn useful representations that preserve causal relationships, our result shows that this can be achieved without perfect identifiability for a large class of models. Furthermore, [13] establish that under some nondegeneracy assumptions, mixing consistency up to surrounding parents is the best possible result when using single-node soft interventions. Hence, the results in 4 are tight for the considered setting.

0.0.5.5 Discussion on Assumptions 1 and 2

In this subsection, we elaborate on 1 and 2, which are relevant to 2 and 4, respectively. 1 essentially states that score changes in the coordinates of the intervened node and a parent of the intervened node are linearly independent. This property holds for (but is not limited to) the widely adopted additive noise models specified in 7 when we apply hard interventions.

When soft interventions are applied, the transitive closure of \(\mcG\) cannot be identified without making assumptions about the effect of the interventions. Specifically, for linear latent causal models, [9] prove impossibility results for pure shift interventions and [8] show that a genericity condition is necessary for identifying the transitive closure of the latent DAG. Therefore, 1 can be interpreted as the counterpart of the commonly adopted assumptions in the literature on soft interventions adapted to the setting of general latent causal models. Finally, the next example demonstrates the working of 1 on a linear Gaussian latent model.

Given the known result that perfect identifiability is impossible for linear Gaussian models given soft interventions [8], the purpose of 2 is to get more insight into the extent of identifiability guarantees under soft interventions. Intuitively, the mentioned impossibility results for linear Gaussian models are due to the rank deficiency of score differences for linear models – specifically, we know that \(\rank(\bR_{\bZ}^{m}) \leq 2\) for linear Gaussian models. In contrast, for sufficiently nonlinear causal models, \(\rank(\bR_{\bZ}^{m})\) can be as high as \(\Paplus(I^m)\). 2 ensures that this upper bound is satisfied with equality for all nodes. This condition holds for the class of sufficiently nonlinear models, such as quadratic causal models. In particular, we show that this condition holds for the two-layer neural networks (NNs) as a function class that can effectively approximate any continuous function. This result is formalized in the next lemma.

We discuss the nonlinearity and the proof of 7 in 1.0.4.2.

0.0.6 CRL under General Transformations↩︎

In this setting, we consider general transformations without any parametric assumption for the transformation \(g\). In the previous section, we exploited the transformation’s linearity and recovered the true encoder’s parameters sequentially using one intervention per node. For general transformations (parametric or nonparametric), however, we cannot use the same parametric approach and rely on the properties of linear transforms. To rectify these and design the general CRL algorithm, we use more information in the form of two interventions per node. Specifically, we will present the steps of leveraging the score functions under two interventions and design an algorithm that identifies the true encoder and recovers the true causal representations. The algorithm is referred to as General Score-based Causal Latent Estimation via Interventions (GSCALE-I), which will be summarized in 5. We also present additional results and discuss the role of various inputs and the distinction of our results compared to the existing literature.

0.0.6.0.0.1 Inputs.

The inputs of GSCALE-I are the observed data from the observational environment, the data from two interventional environments per node, whether environments are coupled/uncoupled, and a set of valid encoders \(\mcH\). Two sets of interventional environments are denoted by \(\mcE = \{\mcE^1,\dots,\mcE^n\}\) and \(\tilde{\mcE} = \{\tilde{\mcE}^1,\dots, \tilde{\mcE}^n\}\), as defined in 0.0.3.3. For these inputs, we compute the score functions \(\bss_{\bX}\), \(\{\bss_{\bX}^1,\dots,\bss_{\bX}^n\}\) and \(\{\tilde{\bss}_{\bX}^1, \dots, \tilde{\bss}_{\bX}^n\}\). Note that we will use 2 again to ensure access to the latent score differences by using these observed score functions.

0.0.6.0.0.2 Statistical diversity.

For being more informative than a single intervention mechanism, we assume that the two intervention mechanisms per node are sufficiently distinct. This is formalized by defining interventional discrepancy [18] among the causal mechanisms of a latent variable.

This condition ensures that the two distributions are sufficiently distinct, formally expressed as the partial derivative of their ratio with respect to the intervened variable is nonzero almost everywhere. For instance, two univariate, distinct Gaussians trivially satisfy this discrepancy condition. As shown by [18], even when the latent graph \(\mcG\) is known, for identifiability via one intervention per node, it is necessary to have interventional discrepancy between observational distribution \(p_i\) and interventional distribution \(q_i\), for all \(\bz_{\Pa(i)}\in\R^{|\Pa(i)|}\).

0.0.6.1 Rationale of GSCALE-I and Partial Identifiability

Similarly to LSCALE-I, analyzing GSCALE-I involves leveraging score differences. However, in contrast to linear transformations, estimated and true latent score differences are not always related by a constant matrix for general transformations. To circumvent this issue, we rely on 1(iii), i.e., the score difference of coupled interventional environments is one-sparse. As such, the key idea is to find an encoder \(h \in \mcH\) that adheres to this sparsity structure, and we will show that such an encoder ensures component-wise identifiability of the latent variables. To formalize these, corresponding to each valid encoder \(h \in \mcH\), we define the score change matrix \(\bD_{\rm t}(h)\) with entries: \[\begin{align} [\bD_{\rm t}(h)]_{i,m} &\triangleq \E\Big[\, \big|\bss_{\hat{\bZ}}^{m}(\hat{\bZ};h) - \tilde{\bss}_{\hat{\bZ}}^{m}(\hat{\bZ};h)\big|_i \Big] \;, \label{eq:Dh-entry} \end{align}\tag{45}\] where expectations are under the measures of latent score functions induced by the probability measure of observational data. We also denote the true score change matrix under true encoder \(g^{-1}\) by \(\bD_{\rm t} \triangleq \bD_{\rm t}(g^{-1})\) with entries \[\label{eq:true-score-change-matrix} [\bD_{\rm t}]_{i,m} \triangleq \E\Big[\big|\bss^{m}(\bZ)- \tilde{\bss}^{m}(\bZ)\big|_i\Big] \;, \quad \forall i,m \in [n] \;.\tag{46}\] For clarity in the exposition of the ideas, we consider coupled environments here, i.e., \(I^m = \tilde{I}^m\) for all \(m\in[n]\). Then, using 27 in 1 we have \[[\bD_{\rm t}]_{i,m} \neq 0 \;\; \iff \;\; i = I^m \;.\] This implies that \(\mathbb{1}\{\bD_{\rm t}\}\) is a permutation matrix, \(\mathbb{1}\{\bD_{\rm t}\} = \bP_{\mcI}^{\top}\) and \[\label{eq:column-Dt-coupled} \mathbb{1}\{[{\bD_{\rm t}}]_{:,m}\} = \be_{\ell} \;, \quad where \; \ell = I^m \;,\tag{47}\] and \(\be_{\ell}\) denotes the \(\ell\)-th standard basis vector. Subsequently, the key idea for identifying the encoder is that the changes between score estimates \(\bss_{\hat{\bZ}}^{m}(\hat{\bz};h)\) and \(\tilde{\bss}_{\hat{\bZ}}^{m}(\hat{\bz};h)\) should match the sparsity structure for those under the true encoder \(g^{-1}\). In the next result, we leverage the sparsity structure in 47 for two interventions on the same node to identify an intervened variable, formalized as follows.

We will show shortly that given a complete set of coupled interventions, 5 will readily imply the identifiability of all variables (in a similar spirit to going from partial to complete identifiability in 0.0.5). Furthermore, the results for uncoupled interventions will also leverage this key result. Hence, we provide the detailed proof of 5 as follows.

0.0.6.2 Complete Identifiability under Coupled Interventions

After achieving the node-level latent variable identifiability from one pair of coupled interventions, we extend the objective function in 5 for a complete set of coupled interventions.

0.0.6.2.0.1 Encoder estimation.

In Stage G1 of 5, we combine the objectives in ?? for all \(m\in[n]\) and form our encoder estimate by solving the following optimization problem \[\label{eq:general-OPT1-procedure} \min\limits_{h \in \mcH} \big\|\bD_{\rm t}(h) - \bI_{n \times n}\big\|_{\rm F}^{2} \;.\tag{51}\] Recall that the valid encoders set \(\mcH\) consists of invertible encoders over the input space \(\mcX\). Hence, 51 is readily equivalent to solving the following problem for any \(\lambda \in \R_+\): \[\label{eq:general-OPT1} \min\limits_{h} \big\|\bD_{\rm t}(h) - \bI_{n \times n}\big\|_{\rm F}^{2} \; + \; \lambda \; \E\Big[ \big\| h^{-1}(h(\bX)) - \bX \big\|_2^2\Big] \;.\tag{52}\]

0.0.6.2.0.2 Latent graph estimation.

Next, we use the encoder estimate \(\hat{h}\) for recovering the latent graph \(\mcG\). In Stage G2 of 5, we use the score differences between observational and interventional environments and form \(\hat{\mcG}\) with parent sets \[\label{eq:general-coupled-parent-set} \hat{\Pa}(m) \triangleq \Big\{i \neq m: \E\Big[\big| \bss_{\hat{\bZ}}(\hat{\bZ}; \hat{h}) - \bss_{\hat{\bZ}}^{m}(\hat{\bZ}; \hat{h}) \big|_i \Big] \neq 0 \Big\} \;, \quad \forall m \in [n] \;.\tag{53}\]

We emphasize that a key contribution of 6 beyond the identifiability result is that it provides a well-defined objective function, a minimizer of which provably recovers the latent variables and the graph up to identifiability guarantees. To our knowledge, this is the first such result for interventional CRL with general transformations. In comparison, [15] also use two hard interventions per node to prove identifiability, albeit without providing a provably correct algorithm. Furthermore, the result in [15] requires that the estimated latent distribution is faithful to the associated candidate graph for all \(h \in \mcH\). Although a faithfulness assumption does not compromise the identifiability result, it is a strong requirement to verify and poses challenges to devise recovery algorithms. In contrast, we only require observational data, which is generally accessible in practice.

Next, we shed light on the role of observational data. 8 requires only interventional data for identifying the encoder, whereas 9 uses observational data. We further tighten this result by showing that for DAG recovery, the observational data becomes unnecessary when we have additive noise models and a weak faithfulness condition holds.

Proof sketch: The recovery of latent variables follow from 8. For the recovery of latent DAG, we leverage 1(iv) under the additive noise assumption. Specifically, using the score difference between environments \(\mcE^i\) and \(\tilde{\mcE}^j\), we obtain \(\overline{\Pa}(i,j)\) for all \(i, j \in [n], i\neq j\). Finally, we perform at most \(n\) conditional independence tests to identify all parent sets \(\{\Pa(i) : i \in [n]\}\). See 1.0.3.1 for the complete proof.

0.0.6.3 Complete Identifiability under Uncoupled Interventions

In this setting, we consider the case where the interventional environments corresponding to the same nodes are not specified in pairs. That is, not only is it unknown which node is intervened in an environment, but the learner also does not know which two environments intervene on the same node. Hence, additionally, we need to determine the correct coupling between the interventional environment sets \(\mcE\) and \(\tilde{\mcE}\) as well. However, this is not a straightforward objective since we cannot readily determine the intervention target in an environment. To address this issue, we solve the problems of finding the correct coupling and optimizing the estimated score variations together. Let \(\sigma\) denote the permutation that takes \(\tilde{\mcI} = (\tilde{I}^1,\dots,\tilde{I}^n)\) to \(\mcI = (I^1,\dots,I^n)\), i.e., \(\sigma(\tilde{I}^m) = I^m\) for all \(m\in[n]\). We will modify the objective function in 52 to solve for permutation \(\pi\) as an estimate of \(\sigma\), in addition to solving for encoder \(h\). To this end, we first modify the score change matrix in 45 and define \(\bD_{\rm t}(h;\pi)\) as \[[\bD_{\rm t}(h;\pi)]_{i,m} \triangleq \E\Big[\, \big|\bss_{\hat{\bZ}}^{m}(\hat{\bZ};h) - \tilde{\bss}_{\hat{\bZ}}^{\pi_m}(\hat{\bZ};h)\big|_i \Big] \;, \label{eq:Dh-pi-entry}\tag{54}\] Unlike the case of coupled interventions, we will also leverage the observational data to identify the encoder. As such, we also define the score change matrices for \(\mcE\) and \(\tilde{\mcE}\) with respect to \(\mcE\) with entries \[\begin{align} [\bD(h)]_{i,m} &\triangleq \E\Big[\big|\bss_{\hat{\bZ}}(\hat{\bZ};h) - \bss_{\hat{\bZ}}^{m}(\hat{\bZ};h)\big|_i \Big] \;, \tag{55} \\ [\tilde{\bD}(h)]_{i,m} &\triangleq \E\Big[\big|\bss_{\hat{\bZ}}(\hat{\bZ};h) - \tilde{\bss}_{\hat{\bZ}}^{m}(\hat{\bZ};h)\big|_i \Big] \;. \tag{56} \end{align}\]

0.0.6.3.0.1 Encoder and coupling estimation.

We modify the objective in 52 and solve the following constrained optimization problem for any \(\lambda \in \R_+\): \[\label{eq:general-OPT2} \begin{align} \min\limits_{h, \pi} \big\|\bD_{\rm t}(h;\pi) &- \bI_{n \times n}\big\|_{\rm F}^{2} \; + \; \lambda \; \E\Big[ \big\| h^{-1}(h(\bX)) - \bX \big\|_2^2\Big] \\ such that \quad & \mathbb{1}\{\bD(h)\}= \mathbb{1}\{\tilde{\bD}(h)\} \cdot \bP_{\pi}^{\top} \\ & \mathbb{1}\{\bD(h)\} \odot \mathbb{1}\{[\bD(h)]^{\top}\} = \bI_{n \times n} \;. \end{align}\tag{57}\] Note that for the true encoder \(g^{-1}\), both \(\bD(g^{-1})\) and \(\tilde{\bD}(g^{-1})\) capture the graph structure via 1(i), albeit under different permutations. Hence, this equivalence of graph structures for the learned encoder is ensured by the first constraint of \(\mathbb{1}\{\bD(h)\}= \bP_{\pi} \cdot \mathbb{1}\{\tilde{\bD}(h)\}\). Also, the acyclicity of the graphs implied by the score differences is ensured by the second constraint, that is \(\bD(h)\) does not contain 2-cycles. To prove that solving this problem leads to the identifiability of the latent variables \(\bZ\), we first show that there exists a global minimizer of this problem.

Proof sketch: See 1.0.3.2 for the complete proof. The main intuition is that \(\pi^* = \sigma\) makes the problem more similar to the case of coupled interventions. The proof then follows by using the results of the coupled interventions to show that \(\bD_{\rm t}(g^{-1}; \sigma)\) is a scaled permutation matrix, using 48 to derive the relationship between the score differences under \(\hat{h}^*\) and the true encoder \(g^{-1}\), and using the sparsity structure of \(\bD(g^{-1})\) to show that the constraints are satisfied.

Next, we show that the objective in 57 can attain its minimum value zero only for the correct coupling.

Proof sketch: See 1.0.3.3 for the complete proof. The main intuition is that the constraints make it impossible to achieve \(\bD_{\rm t}(h;\pi)=\bI_{n \times n}\) for an incorrect coupling \(\pi \neq \sigma\). We prove it by contradiction. We assume that \(h^*\) is a solution, hence, \(\bD_{\rm t}(h^*) = \bI_{n \times n}\), and \(\mathbb{1}\{\bD(h)\}=\mathbb{1}\{\tilde{\bD}(h)\}\). Then, by scrutinizing the eldest node in \(\mcG\) with mismatched interventional environments, we show that \(\bD(h^*) \odot [\bD(h^*)]^{\top}\) cannot be a diagonal matrix, which contradicts the premise that \(h^*\) is a feasible solution of 52 .

Finally, we denote the minimizer encoder of 57 by \(\hat{h}\) and form the graph estimate similarly to the coupled intervention case via 53 . Combining these results, we present our strongest result for general transformations.

8 shows that using observational data enables us to resolve any mismatch between the uncoupled environment sets and shows identifiability in the setting of uncoupled environments. This generalizes the identifiability result of [15], which requires coupled environments.

0.0.6.3.0.2 Comparison to 6.

We note that 8 requires slightly stronger interventional discrepancy conditions than 6. In particular, when environments are coupled, we only need \((q_i,\tilde{q}_i)\) to satisfy the interventional discrepancy. On the other hand, to find the correct coupling while performing CRL, 8 requires each of the pairs \(\{(q_i,\tilde{q}_i),(p_i,q_i), (p_i,\tilde{q}_i)\}\) to satisfy interventional discrepancy.

0.0.6.4 Intervention Extrapolation without CRL

So far, we have studied learning the latent causal representations using the data from single-node interventional environments. A related research problem is extrapolating to unseen combinations of interventions. Namely, given a set of interventions \(\mcI = \{I^1,\dots,I^k\}\), it is desired to emulate the data from an unseen combination of these interventions, e.g., accessing multi-node interventional data using only the single-node interventional data. This is especially important in domains where interventions can be costly or not viable, e.g., not all combinations of different drugs can be clinically tested. Causal representation learning literature has also taken an interest in this problem. Specifically, [10] show for polynomial transformations that given interventions \(\mcI = \{I^1,\dots,I^k\}\), one can sample from any intervention \(I \subseteq \mcI\) after learning the latent representations. We argue that extrapolation to unseen combinations of interventions can be achieved on the observed space without performing CRL for the general transformations.

Consider two single-node interventional environments, \(\mcE^1\) and \(\mcE^2\). Without loss of generality, suppose that \(I^1 = \{1\}\) and \(I^2 = \{2\}\). Also consider the observational environment \(\mcE^0\) with \(I^0 = \emptyset\) and the unseen double-node interventional environment \(\mcE^{m}\) with \(I^{m} = \{1,2\}\). First, using the score function decompositions in 16 and 17 , we have \[\begin{align} \bss(\bz) &= \nabla_{\bz} \log p_1(z_1 \mid \bz_{\Pa(1)}) + \nabla_{\bz} \log p_2(z_2 \mid \bz_{\Pa(2)}) + \sum_{k=3}^{n} \nabla_{\bz}\log p_k(z_k \mid \bz_{\Pa(k)}) \;, \\ \bss^1(\bz) &= \nabla_{\bz} \log q_1(z_1 \mid \bz_{\Pa(1)}) + \nabla_{\bz} \log p_2(z_2 \mid \bz_{\Pa(2)}) + \sum_{k=3}^{n} \nabla_{\bz}\log p_k(z_k \mid \bz_{\Pa(k)}) \;, \\ \bss^2(\bz) &= \nabla_{\bz} \log p_1(z_1 \mid \bz_{\Pa(1)}) + \nabla_{\bz}\log q_2(z_2 \mid \bz_{\Pa(2)}) + \sum_{k=3}^{n} \nabla_{\bz}\log p_k(z_k \mid \bz_{\Pa(k)}) \;, \\ \bss^m(\bz) &= \nabla_{\bz} \log q_1(z_1 \mid \bz_{\Pa(1)}) + \nabla_{\bz} \log q_2(z_2 \mid \bz_{\Pa(2)}) + \sum_{k=3}^{n} \nabla_{\bz} \log p_k(z_k \mid \bz_{\Pa(k)}) \;. \end{align}\] Then, we have \[\begin{align} \label{eq:dsz-extrapolation} \big(\bss^1(\bz) - \bss(\bz)\big) + \big(\bss^2(\bz) - \bss(\bz)\big) = \bss^m(\bz) - \bss(\bz) \;. \end{align}\tag{58}\] Next, recall that 2 gives us the relationship for going from the latent score differences to observed score differences. Applying it to observed \(\bX\) in environment pairs \((\mcE^1,\mcE^0)\), \((\mcE^2,\mcE^0)\), and \((\mcE^m,\mcE^0)\), and using 58 we obtain \[\begin{align} \big(\bss^{1}_{\bX}(\bx) - \bss_{\bX}(\bx)\big) + \big(\bss^{2}_{\bX}(\bx) - \bss_{\bX}(\bx)\big) &= \Big[[J_g(\bz)]^\dagger\Big]^{\top} \cdot \Big[ \big(\bss^1(\bz) - \bss(\bz)\big) + \big(\bss^2(\bz) - \bss(\bz)\big) \Big] \\ &= \Big[[J_g(\bz)]^\dagger\Big]^{\top} \big( \bss^m(\bz) - \bss(\bz) \big) \\ \overset{\eqref{eq:score-latent-to-observed}}&{=} \bss^m_{\bX}(\bx) - \bss_{\bX}(\bx) \;. \label{eq:dsx-extrapolation} \end{align}\tag{59}\] This means that, given score functions of observed data from environments \(\mcE^0\), \(\mcE^1\), and \(\mcE^2\), we can obtain the score function of the unseen interventional environment with \(I^m = \{I^1, I^2\}\), and subsequently generate data from this synthetic environment (using techniques like Langevin sampling [50] that uses score function as the drift vector field). Note that this process does not require learning the latent causal representations and can be applied to obtain the observed score function \(\bss_{\bX}^m\) for any intervention \(I^m \subset \mcI\).

Finally, we note that concurrent work by [51] makes a similar observation for the extrapolation of two single-node interventions without performing CRL. In particular, [51] focuses on detecting pairwise interactions between two biological perturbations (i.e., interventions), investigates the conditions under which two perturbations are separable, e.g., distinct single-node interventions, and designs statistical tests for the separability of two perturbations. If the two perturbations are confirmed to be separable, our results indicate that we can use their score functions to generate synthetic environments that emulate combined perturbation effects.

0.0.7 Empirical Studies↩︎

In this section, we provide empirical assessments of our theoretical guarantees. Specifically, we empirically evaluate the performance of the LSCALE-I (0.0.7.1) and GSCALE-I (0.0.7.2) algorithms for recovering the latent causal variables \(\bZ\) and the latent DAG \(\mcG\) on synthetic data. In 0.0.7.3, we compare the performance of LSCALE-I to that of the existing algorithms in the closely related literature on both synthetic and biological data. Next, we also apply GSCALE-I on image data to demonstrate the potential of our approach in realistic high-dimensional datasets (0.0.7.5). Finally, we note that any desired score estimator can be modularly incorporated into our algorithms. In 0.0.7.6, we assess our performance’s sensitivity to the estimators’ quality. ¹⁰ Additional results and further implementation details are deferred to 1.0.5.

0.0.7.0.0.1 Evaluation metrics.

The objectives are recovering the graph \(\mcG\) and the latent variables \(\bZ\). We use the following metrics to evaluate the accuracy of LSCALE-I and GSCALE-I in recovering these (depending on the specifics of the transformations and interventions, we will also have more specific metrics). For each metric, we will report the mean and standard error over multiple runs.

• Structural Hamming distance: For assessing the recovery of the latent DAG, we report structural Hamming distance (SHD) between the estimate \(\hat{\mcG}\) and true DAG \(\mcG\). This captures the number of edge operations (add, delete, flip) needed to transform \(\hat{\mcG}\) to \(\mcG\).

• Mean correlation coefficient: For the recovery of the latent variables, we use the mean correlation coefficient (MCC), which was introduced in [41] and commonly used as a standard metric in CRL. Specifically, MCC measures the linear correlation between the estimated and ground-truth latent variables. Since the recovery of latent variables is up to permutations, it is reported for the best matching permutation between the components of \(\hat{\bZ}\) and \(\bZ\), i.e., \[\label{eq:MCC} {\rm MCC}(\bZ,\hat{\bZ}) \triangleq \max_{\pi} \frac{1}{n}\sum_{i \in [n]} {\rm corr}(Z_i, \hat{Z}_{\pi(i)}) \;.\tag{60}\]

0.0.7.0.0.2 Score functions.

LSCALE-I and GSCALE-I algorithms, in their first steps, compute estimates of the score differences in the observational environment. The designs of our algorithms are agnostic to how this is performed, i.e., any reliable method for estimating the score differences can be adopted and incorporated into our algorithms in a modular way. In our experiments, we adopt two score estimators necessary for describing the different aspects of our results.

• Perfect score oracle for identifiability: Identifiability, by definition, refers to the possibility of recovering the causal graph and latent variables under all idealized assumptions for the data. Assessing the identifiability guarantees formalized in Theorems 2–8 requires using perfect estimates for the score differences. Hence, we adopt a perfect score oracle for evaluating identifiability. Specifically, we use a perfect score oracle that computes the score difference inputs in LSCALE-I and GSCALE-I by leveraging 2 and using the ground truth score functions \(\bss,\bss^m\) and \(\tilde{\bss}^m\) (see 1.0.5.1 for details).

• Data-driven score estimates: For evaluating the accuracy of our algorithms in practice, we need real score estimates, which are inevitably noisy. For this purpose, when the pdf \(p_\bX\) has a parametric form and the score function \(\bss_{\bX}\) has a closed-form expression, then we can estimate the parameters to form an estimated score function. For instance, when \(\bX\) follows a linear Gaussian distribution, then we have \(\bss_{\bX}(\bx)= - \Theta \cdot \bx\) in which \(\Theta\) is the precision matrix of \(\bX\) and can be estimated from samples of \(\bX\). In other cases in which \(\bss_{\bX}\) does not have a known closed-form, we adopt nonparametric score estimators. In particular, we use sliced score matching with variance reduction (SSM-VR) for score estimation due to its efficiency and accuracy for downstream tasks [52]. We also introduce a classification-based score difference estimation method, inspired by [53]. The key observation is that given two distributions \(p\) and \(q\), the optimal minimum cross-entropy classifier for distinguishing the samples from two distributions is the log density ratio function \(\log p/q\). The difference between the score functions of these distributions, \(\nabla \log p - \nabla \log q = \nabla \log p/q\) is exactly the gradient of the learned function, which enables us to estimate score differences using a classifier directly.

0.0.7.1 LSCALE-I Algorithm for Linear Transformations

0.0.7.1.0.1 Data generation.

To generate \(\mcG\), we use Erdős-Rényi model with density \(0.5\) and \(n \in \{5,8\}\) nodes, which is generally the size of the latent graphs considered in CRL literature. We consider the observed dimension \(d=100\) and generate 100 latent graphs. For the causal mechanisms, we adopt both linear and nonlinear models:

1. Linear causal model: We adopt the linear Gaussian model with \[Z_i = \bA_i \cdot \bZ + N_i \;, \qquad \forall i \in [n] \;, \label{eq:linear-SEM--experiments}\tag{61}\] where \(\bA_i \in \R^{1 \times n}\) are the rows of the weight matrix \(\bA\) in which \(\bA_{i,j}\neq 0\) if and only \(j \in \Pa(i)\). The nonzero edge weights are sampled from \({\rm Unif}(\pm[0.5, 1.5])\), and the noise terms are zero-mean Gaussian variables with variances \(\sigma_{i}^2\) sampled from \({\rm Unif}([0.5, 1.5])\). For node \(i\), a hard intervention is given by \(Z_i = \bar {N}_i\) where \(\bar N_i \sim \mcN(0,\frac{\sigma_i^2}{4})\), and a soft intervention is given by \(Z_i =\frac{\bar \bA_{i}}{2} \cdot \bZ + N_{i}\).

2. Quadratic causal model: We adopt an additive noise model with \[\begin{align} Z_i =\sqrt{\bZ_{\Pa(i)}^{\top} \cdot \bQ_{i} \cdot \bZ_{\Pa(i)}} + N_{i}\;, \label{eq:quadratic-model--general-experiments} \end{align}\tag{62}\] where \(\{\bQ_{i}:i\in[n]\}\) are positive-definite matrices, and the noise terms are zero-mean Gaussian variables with variances \(\sigma_{i}^2\) sampled randomly from \({\rm Unif}([0.5,1.5])\). For node \(i\), a hard intervention is given by \(Z_i = \bar {N}_i\) where \(\bar N_i \sim \mcN(0,5 \cdot \sigma_i^2)\) to ensure a non-negligible change. A soft intervention is given by \[Z_i =\frac{1}{2}\sqrt{\bZ_{\Pa(i)}^{\top} \cdot \bQ_{i} \cdot \bZ_{\Pa(i)}} + \bar N_{i}\;.\]

3. Multi-layer perceptron (MLP) causal model: We adopt an additive noise model with \(Z_i = f_i(\bZ_{\Pa(i)}) + N_i\) where \(f_i\) is parameterized as a randomly initialized two-layer MLP, with hidden dimension 32 and ReLU activation function. For node \(i\), a hard intervention is given by \(Z_i = \bar N_i\) where \(\mcN(0,5 \cdot \sigma_i^2)\), and a soft intervention is given by \(Z_i =\frac{1}{2} f_i(\bZ_{\Pa(i)}) + \bar N_{i}\).

For each graph, we sample \(n_{\rm s}\) independent and identically distributed (i.i.d.) samples of \(\bZ\) from each environment. We consider \(n_{\rm s} \in \{5000, 10000, 50000 \}\) to investigate the effect of the number of samples on the performance of LSCALE-I. The observed variables \(\bX\) are generated according to \(\bX = \bG \cdot \bZ\), in which \(\bG \in \R^{d \times n}\) is a randomly sampled full-rank matrix.

0.0.7.1.1 Hard Interventions

3 ensures scaling consistency and perfect DAG recovery under hard interventions for linear transformations. As such, we assess the recovery of the latent DAG by the SHD between the estimate \(\hat{\mcG}\) and true graph \(\mcG\). For latent variable recovery, we report MCC between \(\hat{\bZ}\) and \(\bZ\). Note that \(\hat{\bZ} = (\bH \cdot \bG) \cdot \bZ\) implies that the effective transform recovery can also be measured by the closeness of \(\bH \cdot \bG\) to the identity matrix. Therefore, in addition to MCC, we also report the normalized effective transform error, defined as \[\label{eq:norm95scale95error} \ell_{\rm scale} \triangleq \| \bH \cdot \bG - \bI_{n \times n} \|_2 \;.\tag{63}\] 3 shows the performance of the LSCALE-I algorithm using perfect scores and noisy scores under hard interventions on linear causal models. The first observation is that we achieve excellent performance in latent variable recovery, as demonstrated by a perfect MCC and nearly zero effective transform error \(\ell_{\rm scale}\), even when using noisy score estimates. For graph recovery, SHD between the estimated and true latent graphs reduces to less than 0.1 even when using noisy scores given enough samples, e.g., \(n_s=50000\) in 3. Another key observation is that the results remain consistent while the dimension of observed variables increases from \(d=25\) to \(d=200\), as shown in 4. This confirms our analysis that the performance of LSCALE-I is agnostic to the dimension of the observations. Hence, we suffice by using \(d=100\) for the rest of the experiments.

0.0.7.1.1.1 Causal models with more complex score functions.

Next, we test the performance of LSCALE-I on quadratic and MLP causal models specified earlier. A difference in this setting compared to linear causal models is that we estimate the score functions of the observed variables using SSM-VR [52] as \(p_\bX\) is not amenable to parameter estimation. For quadratic causal models, 5 shows that using perfect scores, LSCALE-I performs nearly perfectly. Under noisy scores, an MCC of 0.93 indicates a strong performance for the latent variable recovery. The graph recovery performance suffers more from noisy score estimates, yet it remains reasonable with an approximate SHD of \(2.62\) (where the expected number of true edges is \(5\)). 6 shows a similar trend, where perfect scores ensure perfect performance, whereas noisy score estimates lead to a degradation in graph recovery performance.

0.0.7.1.1.2 Effect of noisy score estimates on graph recovery under noisy scores.

5 and 6 show that graph recovery performance suffers more than the latent variable recovery under noisy score estimates. This discrepancy is due to the difficulty of applying 40 under noisy scores. Specifically, we observe that the nonzero entries in true score differences (entries of \(\bsd^m\)) can be small, especially for MLP causal models. As such, under noisy score estimates, it becomes more difficult to successfully threshold the empirical score difference quantities to identify the edges of the latent graph. In 0.0.7.6, we provide further empirical evaluations on the effect of score estimation quality on the performance of LSCALE-I.

0.0.7.1.2 Soft Interventions

To evaluate soft interventions, we first consider linear causal models. In this setting, 2 ensures the recovery of the latent variables and the latent DAG up to ancestors when using soft interventions. So, we report the SHD between the transitive closure of the estimate \(\hat{\mcG}_{\rm tc}\) and that of the true graph \(\mcG_{\rm tc}\). Recall that recovery of latent variables for generic soft interventions is guaranteed up to mixing with parents. To measure the accuracy of \(\hat{\bZ} = \bH \cdot \bG \cdot \bZ\) with respect to this guarantee, we define \(\bL_{\rm pa}\) with entries \([\bL_{\rm pa}]_{i,j} = \mathbb{1}\{j \in \Paplus(i)\}\), and define \(\ell_{\rm pa}\) by \[\label{eq:norm95parents95error} \ell_{\rm pa} \triangleq \| \bH \cdot \bG \odot (\boldsymbol{1}_{n \times n} - \bL_{\rm pa}) \|_2 \;,\tag{64}\] which effectively measures the effect of incorrect mixing in estimated latent variables. 7 shows that by using perfect scores from as few as \(n_{\rm s}=5000\) samples, we meet the theoretical identifiability guarantees, implied by zero values of SHD and incorrect mixing norm \(\ell_{\rm pa}\). Similar to the case of hard interventions, we observe that increasing the number of samples improves the performance under noisy scores. For instance, given \(n_{\rm s}=50000\) samples, the average SHD between the transitive closures \(\hat{\mcG}_{\rm tc}\) and \(\mcG_{\rm tc}\) is approximately \(0.5\) for \(n=8\), where graph density being \(0.5\) implies that the expected number of edges is \(14\). Also, the latent variable recovery becomes near perfect, indicated by MCC of \(0.98\) and near-zero incorrect mixing norm \(\ell_{\rm pa}\).

Next, we consider quadratic causal models. In this case, 2 is satisfied, and 4 ensures the perfect recovery of the latent DAG and the recovery of the latent variables up to surrounding variables. Hence, we report SHD between true \(\mcG\) and estimate \(\hat{G}\), and the incorrect mixing norm for this setting given by \[\label{eq:norm95sur95error} \ell_{\rm sur} \triangleq \| \bH \cdot \bG \odot (\boldsymbol{1}_{n \times n} - \bL_{\rm sur}) \|_2 \;,\tag{65}\] where \(\bL_{\rm sur}\) is defined with entries \([\bL_{\rm pa}]_{i,j} = \mathbb{1}\{j \in \overline{\sur}(i)\}\). 8 shows that when using perfect scores, LSCALE-I performs nearly perfectly, verifying the results of 4. Similarly to the case of hard interventions, the performance of LSCALE-I suffers under noisy score estimations, while it remains reasonable, e.g., an SHD of \(2.79\) with respect to the true graph when using only soft interventions.

0.0.7.2 GSCALE-I Algorithm for General Transformations

Next, we focus on nonlinear transformations to showcase the settings for which the existing interventional CRL literature lacks provably correct algorithms and provides only identifiability results under nonlinear transformations.

0.0.7.2.0.1 Choice of nonlinearity.

In this section, we consider M-layer perceptrons with \(\tanh\) activation as our class of nonlinear transformations. Specifically, we consider functions of the form \[\label{eq:transf95exp-2} \bX = g(\bZ) = (\tanh \circ \bA^{M} \circ \cdots \tanh \circ \bA^{1}) (\bZ) \;,\tag{66}\] in which \(\tanh\) is applied element-wise, and parameters \(\{\bA^{1}, \dots, \bA^{M}\}\) are compatible, randomly sampled full column rank matrices. We consider two settings of increasing complexity: 1-layer MLP, and 3-layer MLP (for nonparametric transformations, see experiments on images in 0.0.7.5).

0.0.7.2.1 Single-layer MLP

First, consider the single-layer perceptron as the nonlinear transformation \(g\), given by \[\label{eq:transf95exp-1} \bX = g(\bZ) = \tanh (\bG \cdot \bZ) \;,\tag{67}\] where \(\bG \in \R^{d \times n}\). Because it is a relatively simple nonlinear transformation, this setting enables us to evaluate the performance of GSCALE-I more effectively. Specifically, leveraging 67 , we parameterize valid encoders \(h\) with parameter \(\bH \in \R^{n \times d}\), which gives \[\begin{align} \label{eq:tanh-GLM-parameterization} \hat{\bZ}(\bX;h) &= h(\bX) = \bH \cdot \arctanh (\bX) \;, \\ \hat{\bX} &= h^{-1}(\hat{\bZ}(\bX;h)) = \tanh (\bH^{\dagger} \cdot \hat{\bZ}(\bX;h)) \;, \\ \hat{\bZ} &= \bH \cdot \bG \cdot \bZ \;. \end{align}\tag{68}\] This is equivalent to reducing the original nonparametric functional estimation problem to a finite-dimensional parameter estimation problem. This enables stable training and large-scale testing. Furthermore, this parametrization enables us to compute the ground truth score functions of the observed variables, allowing us to directly assess the effect of score estimation errors on performance.

0.0.7.2.1.1 Data generation.

In this setting, we focus on quadratic causal models. To generate \(\mcG\) we use the Erdős-Rényi model with density \(0.5\) and \(n \in \{5,8\}\) nodes. For observational causal mechanisms, we use Equation 62 and the corresponding parameterization. For the two hard interventions on node \(i\), \(Z_i\) is set to \(N_{q,i} \sim \mcN(0,\sigma_{q,i}^2)\) and \(N_{\tilde{q},i} \sim \mcN(0,\sigma_{\tilde{q},i}^2)\), where we set \(\sigma_{q,i}^2 = \sigma_{i}^2 / 4\) and \(\sigma_{\tilde{q},i}^2 = 4 \cdot \sigma_{i}^2\). Similarly to LSCALE-I experiments, we consider a target dimension of \(d=100\) and generate 20 latent graphs, with \(n_{\rm s} = 30000\) samples per environment for each graph.

0.0.7.2.1.2 Results.

Recall that for general transformations, we can guarantee \(\hat{Z}_i = \phi_i(Z_i)\) for a diffeomorphism \(\phi_i\). However, MCC\((Z_i, \phi(Z_i))\) can largely deviate from 1 for a nonlinear \(\phi_i\). For instance, for \(Z_i \sim \mcN(0,1)\) and \(\phi(Z_i) = Z_i^3 + 0.1 Z_i\), we have MCC\((Z_i, \phi(Z_i)) \approx 0.786\). On the other hand, when we use the parameterization in 68 , the only element-wise diffeomorphism between \(\bZ\) and \(\hat{\bZ}\) is an element-wise scaling. Therefore, MCC remains a perfectly informative metric. 9 shows that we can almost perfectly recover the latent variables (indicated by MCC of \(1\)) latent DAG for \(n=5\) nodes by using perfect scores. Furthermore, increasing the latent dimension to \(n=8\) has no significant impact on performance.

0.0.7.2.1.3 The effect of the quality of score estimation.

When using noisy scores, the performance of GSCALE-I degrades significantly. For instance, MCC goes down to approximately \(0.85\) for \(n=5\) and \(0.75\) for \(n=8\). This degradation is similar to what happens when using LSCALE-I on quadratic causal models in 8. It is noteworthy that the transition from perfect to noisy scores is remarkably smoother when using LSCALE-I on linear causal models (3 and 7). This discrepancy is attributed to the distinct score estimation procedures adopted in the two experimental settings. Specifically, linear Gaussian latent models allow us to directly estimate the parameters of the closed-form score function \(\bss_{\bX}\). However, when using a quadratic latent model, we rely on a nonparametric score estimation via SSM-VR [52]. This comparison between two experiment settings and results underscores that the performance gap between the theoretical guarantees and practical results can be significantly mitigated through advances in general score estimation techniques. In 0.0.7.6, we provide a further empirical evaluation of how the quality of the score estimation affects the final performance.

0.0.7.2.2 Multi-layer nonlinearity

Next, we consider general MLPs by setting the number of layers to 3, which is in line with the MLP depths considered in CRL literature [18]. In this setting, we do not use any knowledge of the transformation parametrization and aim to learn a generic NN-based encoder-decoder pair \(h, h^{-1}\) parameterized by \(\theta_{\rm enc}, \theta_{\rm dec}\). In essence, this enables us to test our algorithm’s performance on a standard procedure of function approximation via neural networks.

0.0.7.2.2.1 Data generation.

In this setting, we focus on linear causal models. To generate \(\mcG\), we use the Erdős-Rényi model with density \(0.5\) and \(n=5\) nodes. For observational causal mechanisms, we use Equation 61 and its associated parameterization. For the two hard interventions on node \(i\), \(Z_i\) is set to \(N_{q,i} \sim \mcN(0,\sigma_{q,i}^2)\) and \(N_{\tilde{q},i} \sim \mcN(0,\sigma_{\tilde{q},i}^2)\), where we set \(\sigma_{q,i}^2 = \sigma_{i}^2/4\) and \(\sigma_{\tilde{q},i}^2 = \sigma_{i}^2 \cdot 4\). Due to the increased computational complexity of training, we consider target dimension value \(d=5\), and we generate 5 latent graphs and \(n_{\rm s} = 10000\) samples per graph.

0.0.7.2.2.2 Results.

GSCALE-I algorithm ensures latent variables recovery up to element-wise transformations. As a proxy, we report its linear analogue, i.e., MCC. We summarize the results of these experiments in 10.

In this setting, we obtain a smaller average MCC rate (\(0.54\)) compared to the 1-layer nonlinearity experiments in 9. Performance degradation can be attributed to two factors: First, making the transformation more complex and making it more difficult to construct any kind of autoencoder. In contrast to the parametric single-layer MLP setting in the previous section, fitting a generic NN has considerably higher sample complexity. Second, our methodology requires the estimation of score functions with relatively good accuracy. However, complex distributions can push forward distributions into very complex and, similarly, difficult-to-estimate landscapes, which can limit the accuracy of our algorithms.

0.0.7.3 Comparison with the Existing CRL Studies

In this section, we compare the performance of LSCALE-I with those of the approaches designed for comparable settings in  [8], [10].

0.0.7.3.1 Comparison with linear causal models on synthetic data

A closely related study in the linear transformation setting is [8]. As discussed in 0.0.5, [8] presents identifiability guarantees and provides algorithms for linear transformations with linear causal models. For a fair comparison, we perform comparisons for linear Gaussian SEMs under the parameterization detailed in Section 7.1. Experiments with synthetic linear causal models in [8] mostly focus on recovering the parameters of the true encoder \(\bG^\dagger\). Hence, in 11, we report performance comparisons in terms of the latent variable recovery metric \(\ell_{\rm scale}\). 11 shows that LSCALE-I achieves close to perfect accuracy with as few as \(n_s=5000\) samples, whereas the algorithm of [8] requires significantly more samples (more than \(n_s=50,000\)) to achieve a reasonable performance.

0.0.7.3.2 Comparison with DiscrepancyVAE on synthetic data

In this section, we compare LSCALE-I with the DiscrepancyVAE algorithm of [10] on the setting described in 0.0.7.1. We consider 10 runs of the quadratic causal model setting with soft interventions, latent dimension \(n = 5\), observed dimension \(d = 100\), and \(50000\) samples per environment. We observe that the average SHD of DiscrepancyVAE on this dataset is \(3.1 \pm 0.7\), which is comparable to LSCALE-I performance of \(2.79 \pm 0.18\) reported in 8.

To provide a visual comparison, we consider a single run and plot the ground truth latent graph (6 (a)), the estimate generated by LSCALE-I (6 (b)), and the estimate generated by DiscrepancyVAE (6 (c)). The differences between the latter and the ground truth graph are shown in these figures by red (added) and dashed (deleted) edges. We refer further comparison between these methods to 1.0.5.3.

0.0.7.3.3 Comparison with DiscrepancyVAE on biological data

In this section, we apply our algorithms to the Perturb-seq dataset of [54], which is also used in [10]. Following the same pre-processing steps, there are 8,907 unperturbed and 99,590 perturbed cells. Each cell (sample) is denoted by a sparse vector of dimension 5000, which represents the expression levels of select genes. The perturbed cells are generated by CRISPR activation [55] on one or two target genes out of 105 genes. Hence, this is modeled by a 5000-dimensional observable variable \(\bX\) and a 105-dimensional latent variable \(\bZ\) representing the perturbed genes. In this dataset, the ground truth latent graph and variables are unknown, with the exception of limited gene interactions that have been experimentally verified [54]. For graph estimation, instead of working with a latent dimension of 105, [10] first learns groups of perturbation targets subject to a regularity constraint and then learns a latent graph among these groups. Since the methodology of LSCALE-I does not naturally extend to learning such groups of latent nodes, we use the perturbation groupings reported by [10] as our super-nodes and run LSCALE-I to learn a latent graph over them. We list representative perturbation targets from each super-node in  7 (a).

The most important observation is that LSCALE-I and DiscrepancyVAE both estimate a causal edge from DUSP9 to the MAPK1/ETS2 group, which is a gene activation relation that is demonstrated experimentally [54]. Due to the lack of ground truth for other graph edges, comparing the algorithms for recovering other edges is not informative. Nevertheless, we provide the estimated latent graph from LSCALE-I ([fig:perturb-seq-graph-est:lscalei]) and DiscrepancyVAE ([fig:perturb-seq-graph-est:discrvae]) for completeness. This verifies that our score-based methodology can recover causal genomics relations.

0.0.7.4 Intervention Extrapolation on Biological Data

In 0.0.6.4, we have demonstrated that it is possible to sample from unseen combinations of given interventions via only learning the score functions of observed variables. We verify this on the Perturb-seq data discussed in 0.0.7.3 through the following steps:

1. Collect observed variables’ data from observational distribution \(p_\bX\), two single-node interventions \(p_\bX^1\) and \(p_\bX^2\), and double-node intervention \(p_\bx^{\{1,2\}}\).

2. Compute score functions \(\bss_{\bX}\), \(\bss_{\bX}^1\), and \(\bss_{\bX}^2\).

3. Perform score function extrapolation for double-node intervention via 59 : \(\tilde{\bss}_{\bX}^{\{1,2\}} = \bss_{\bX}^1 + \bss_{\bX}^2 - \bss_{\bX}\).

4. Single-node setting: Generate samples from \(\bss_{\bX}^1\) and \(\bss_{\bX}^2\) using Langevin dynamics [50], compare them to the original samples from \(p_\bX^1\) and \(p_\bX^2\).

5. Double-node intervention extrapolation: Generate samples from \(\tilde{\bss}_{\bX}^{\{1,2\}}\) using Langevin dynamics, compare them to the original samples from \(p_\bX^{\{1,2\}}\).

We note that [10] also report experiments on the same dataset for the single-node data generation and double-node intervention extrapolation settings. However, as discussed in 0.0.6.4, they perform these tasks using the VAE model for CRL, whereas our approach shows that learning the latent causal representations is unnecessary. Similarly to [10], we use the maximum mean discrepancy (MMD) [56] to measure performance quantitatively. In 12, we show that our approach performs at least as well as DiscrepancyVAE on the Perturb-seq dataset. Finally, in 1.0.5.4, we illustrate the simulated interventional distributions via UMAP clustering, which visually supports the performance implied by 12. These results validate the theoretical analysis that we do not need to learn latent variables to perform intervention extrapolation.

0.0.7.5 Experiments on Image Data

In this section, we perform experiments on CRL where the general transformation is image rendering, a highly nonlinear transformation, by applying the GSCALE-I algorithm to synthetic image data.

0.0.7.5.0.1 Image data generation.

For image-based experiments, we follow the setup of the closely related studies in [7], [9]. Specifically, we consider images of the form in 9, which are generated as follows. The pairs of latent variables \((Z_{2i-1},Z_{2i})\) describe the coordinates of the \(i\)-th ball’s center in a \(64 \times 64 \times 3\) RGB image. We use two and three balls in our experiments, which corresponds to \(n \in \{4, 6\}\) latent variables. We sample the latent graph \(\mcG\) from Erdős–Rényi model for \(n\) nodes and set the expected number of edges to \(2n\). Given this graph \(\mcG\), we adopt a truncated linear Gaussian latent causal model. The details of the linear model are the same as the 0.0.7.1, with the addition of a second set of interventions represented by the causal mechanisms \(\tilde{q}_{i}\), by setting \(Z_{i} = \tilde{N}_{i}\), where \(\tilde{N}_{i} \sim \mcN(0, \frac{\sigma_{i}^{2}}{4})\). Latent variables \(\bZ\) are sampled from the distribution defined by this linear model, but samples with any coordinate absolute value greater than \(1\) are discarded to truncate the distribution in the box \([-1, 1]^{n}\). For each latent sample \(\bZ\), we render an image of the balls centered at the specified coordinates (shifted and scaled to fit entirely within rendered image boundaries) with a radius of 8 pixels and use these as our observed data \(\bX\). Balls are color-coded, filled with three different colors, and the background is white. For both two and three balls, we generate 5 graphs and \(10000\) samples from each graph under each environment.

0.0.7.5.0.2 Candidate encoder and training.

For learning the latent variables, we construct a two-step autoencoder, depicted in 8. Specifically, in step one, we train an autoencoder on the observational image dataset with a bottleneck dimension of 64 using only a reconstruction objective. In step two, we train another autoencoder on the 64-dimensional output of the first encoder with the bottleneck dimension of \(n\) (the latent dimension) using both reconstruction loss and the score-based loss described earlier in 0.0.7.2. Architecture details and training schedule are given in 16. Finally, we repeat the training of the second autoencoder 5 times for each dataset.

0.0.7.5.0.3 Score difference estimation.

The main loss function for the second autoencoder of our setup is \(\norm{\bD_{\rm t}(h)-\bI_{n \times n}}_{1,1}\), which is equivalent to the score-based loss in 51 . For computing \(\bD_{\rm t}(h)\) for an encoder \(h\), we need to compute the score differences of image datasets. To this end, we adopt a binary classifier-based log density ratio (LDR) estimator and use the gradients of these learned LDRs as our score difference functions, as described earlier. Specifically, we parameterize the LDR at any image through a CNN-based model, details of which are given in 17. The output of this model is used to compute class probabilities, which are used to minimize cross-entropy to train the model. We train one LDR model for each pair of hard interventions on the same node and one for the observational-interventional environment pairs for one set of hard interventions.

0.0.7.5.0.4 Results.

GSCALE-I ensures perfect graph recovery and latent variables recovery up to permutation and element-wise scaling. Similarly to related work [7], [9], we report MCC between \(\hat{\bZ}\) and \(\bZ\) as a metric of latent variables recovery, shown in 13.

During training, we observed that MCC highly correlates with the reconstruction performance, i.e., when reconstruction errors in both autoencoders converge to small values, latent variable recovery is also successful. We note that the reconstruction does not have to be perfect. For instance, the samples in 9 are from a run with 3 balls with an MCC of 0.98. Conversely, for the runs with non-convergent reconstruction loss in either of the autoencoder training steps, latent variable recovery also struggled. Hence, as long as we can ensure reasonable reconstruction performance when we are minimizing the score-based loss, our algorithm manages to disentangle the true causal variables embedded in the latent domain. We report the results from related work in 13 and discuss them as follows. Overall, we achieve a better performance than that of [7], and slightly weaker than that of [9] for linear SCMs. This can perhaps be explained by the versatility of our approach, which works for general SCMs, whereas the algorithm of [9] is designed to work exclusively for linear SCMs, both theoretically and due to its loss function design. We demonstrate in 13 that our algorithm performs even better when the latent SCM is nonlinear (specifically, quadratic) than when it is linear, thereby verifying its versatility. Finally, we also test our algorithm when using one intervention per node, using the score differences between interventional and observational environments instead of two interventional environments. While two interventions per node are required for theoretical guarantees, our experiments show that using one intervention per node generally leads to competitive performance, as illustrated in 13.

0.0.7.5.0.5 Gap between the theory and practice in CRL.

We highlight that the nature of the results presented in this paper is theoretical. We have established unknown identifiability guarantees in several settings. Furthermore, we have designed algorithms by defining a differentiable loss function whose global optima achieve identifiability for general CRL and established that the algorithms generate provably correct representations. The experiments presented on various synthesized and real datasets demonstrate the potential of the score-based framework for complex nonlinear transforms, such as image rendering. We achieve this goal by showing reasonable performance in latent variable recovery. That being said, the considered image dataset and other synthetic datasets are still relatively simple compared to real-world problems, and a gap remains between theoretical guarantees and practical applications. In other words, the theory developed is only a necessary step for guaranteed practical performance. From our perspective, closing this gap and obtaining sufficient conditions for practical performance requires improvements in two aspects. First, using better score difference estimators (which we empirically investigate in 0.0.7.6) would enable a more effective application of our score-based CRL framework. Second, designing more suitable architectures for leveraging the invariance properties of causal models can help scale up CRL applications.

0.0.7.6 Sensitivity to Score Estimation Noise

As discussed earlier, score estimators can be modularly incorporated into our algorithm. For cases where \(p_{\bX}\) is not amenable to parameter estimation, we can use the noisy score estimates generated by SSM-VR as achievable baselines, as presented in this section. In this subsection, we evaluate the potential performance improvement that can be achieved as the score estimates become more accurate. To this end, we test the LSCALE-I algorithm under varying levels of score estimation noise. Specifically, we run the LSCALE-I algorithm using the scores generated by the following model: \[\hat{\bss}_{\bX}(\bx; \sigma^2) = \bss_{\bX}(\bx) \cdot \big( 1 + \Xi \big) \;, \quad where \quad \Xi \sim \mcN(0, \sigma^2 \cdot \bI_{d \times d}) \;.\] We consider hard interventions on graph size of \(n = 5\), set \(d = 25\), and vary the value of \(\sigma^2\) within \([10^{-4}, 10^{-2}]\). We repeat the experiments 100 times with \(n_{\rm s}=10000\) samples from each environment. We plot mean normalized \(\bZ\) error, defined as \[\label{eq:normalized95loss} \ell_{\rm norm}(\bZ,\hat{\bZ}) \triangleq \frac{\|\bZ-\hat{\bZ}\|_2}{\norm{\bZ}_2}\;.\tag{69}\] and mean SHD in 10 with respect to the signal-to-noise ratio (SNR). It is clear that when the score estimation error is small, indicated by a high SNR value, the LSCALE-I algorithm demonstrates a strong performance in recovering both latent causal variables and the latent graph. We also note that the baseline results with noisy scores computed via SSM-VR in 5 yield similar success at graph recovery at the SNR of approximately \(25\) dB and latent recovery at the SRM of approximately \(8\) dB. The curves also confirm our observations in 0.0.7.1 and 0.0.7.2 that graph recovery is more sensitive to score estimation errors. The trend of the curves in 10 indicates that our algorithm would greatly benefit from a better score estimator.

0.0.8 Discussion and Concluding Remarks↩︎

In this paper, we have proposed a score function-based CRL framework that uses stochastic interventions to learn latent causal representations and the latent causal graph underlying them. In this framework, by uncovering novel connections between score functions and CRL, we have established identifiability results for linear and general transformations without restricting the latent causal models, and designed LSCALE-I and GSCALE-I algorithms that achieve these identifiability guarantees. There are several exciting directions for future work.

As discussed in 0.0.3.4 after defining the identifiability objectives, and emphasized when presenting the corresponding results through the paper, our results are tight given a complete set of atomic interventions. A natural direction for future work is relaxing the atomic intervention requirement. In this aspect, the paper [57] extends the score-based framework for linear transformations into the setting of unknown multi-node interventions by using combinations of multi-node score differences. Given a sufficiently diverse multi-node intervention set, similar identifiability results to those in 0.0.5 are shown. That being said, establishing the necessary conditions for multi-node interventions and the case of general transformations remains an open problem. Next, we note that our LSCALE-I algorithm for linear transformations is agnostic to the intervention type and latent causal model, meaning that it can be used for both soft and hard interventions and different causal models. Designing similar universal CRL algorithms that can handle general transformations under different sizes and types of interventional environments is the ultimate goal for studying the identifiability of CRL from interventions. A missing component of existing CRL literature is the finite-sample analysis. Probabilistic identifiability results for a given number of interventional data samples can be useful, especially in applications where performing interventions is costly. In this direction, the recent paper [58] establishes finite sample guarantees for CRL under linear transformations. Extending such analysis to general transformations can provide insights for bridging the gap between theory and practice in CRL.

1 Appendix↩︎

1.0.1 Proofs of Score Function Properties and Transformations↩︎

We start by providing the following facts that will be used repeatedly in the proofs.

If there exists \(\bz \in \R^n\) such that \(f(\bz)\neq g(\bz)\), then \(f(\bz)-g(\bz)\) is nonzero over a nonzero-measure set due to continuity. Then, \(\E\big[|f(\bZ)-g(\bZ)|^{\alpha}\big] \neq 0\) since \(p\) (pdf of \(\bZ\)) has full support. On the other direction, if \(f(\bz)=g(\bz)\) for all \(\bz \in \R^n\), then \(\E\big[|f(\bZ)-g(\bZ)|^{\alpha}\big]=0\). This means that \(\E\big[|f(\bZ)-g(\bZ)|^{\alpha}\big]\neq 0\) implies that there exists \(\bz \in \R^n\) such that \(f(\bz) \neq g(\bz)\).

1.0.1.1 Proof of Lemma 1

Our score-based methodology builds on the changes in score functions under interventions. For proving 1, we start by showing that \(\E\Big[\big|\bss(\bZ) - \bss^{m}(\bZ)\big|_i\Big] \neq 0 \implies i \in \Paplus(I^m)\), which holds true regardless of the causal model and the intervention type.

1.0.1.1.0.1 Proof of \(\E\Big[\big|\bss(\bZ) - \bss^{m}(\bZ)\big|_i \Big] \neq 0 \implies i \in \Paplus(I^m)\):

Let \(\ell\) denote the node intervened in \(\mcE^m\), i.e., \(I^m = \ell\). Recalling 16 and 17 , the latent scores \(\bss(\bz)\) and \(\bss^{m}(\bz)\) are decomposed as \[\begin{align} \bss(\bz) &= \nabla_{\bz} \log p_{\ell}(z_{\ell} \med \bz_{\Pa(\ell)}) + \sum_{i \neq \ell} \nabla_{\bz} \log p_i(z_i \med \bz_{\Pa(i)}) \;, \tag{70} \\ and \quad \bss^{m}(\bz) &= \nabla_{\bz} \log q_{\ell}(z_{\ell} \med \bz_{\Pa(\ell)}) + \sum_{i \neq \ell} \nabla_{\bz} \log p_i(z_i \med \bz_{\Pa(i)}) \;. \tag{71} \end{align}\] Hence, \(\bss(\bz)\) and \(\bss^{m}(\bz)\) differ in only the causal mechanism of node \(\ell\). Next, we check the derivatives of \(p_{\ell}(z_{\ell} \med \bz_{\Pa(\ell)})\) and \(q_{\ell}(z_{\ell} \med \bz_{\Pa(\ell)})\) in their \(i\)-th coordinates. Note that these two depend on \(\bZ\) only through \(\{Z_j : \; j \in \Paplus(\ell) \}\). Therefore, if \(i \notin \Paplus(\ell)\), \[\require{physics} \begin{align} \pdv{z_i} \log p_{\ell}(z_{\ell} \med \bz_{\Pa(\ell)}) = \pdv{z_i} \log q_{\ell}(z_{\ell} \med \bz_{\Pa(\ell)}) = 0 \;, \label{eq:s95z95zero95indices} \end{align}\tag{72}\] which indicates that if \(i \notin \Paplus(\ell)\), then \([\bss(\bz)]_i = [\bss^{m}(\bz)]_i\) for all \(\bz\). This, equivalently, means that if \(\E\big[|\bss(\bZ) - \bss^{m}(\bZ)|_i \big] \neq 0\), then \(i \in \Paplus({\ell})\).

For the reverse direction, we will use the following intermediate result which formalizes the weakest possible requirement for a meaningful intervention and shows that it is a property of (i) hard interventions under any causal model, and (ii) additive noise model under either soft or hard interventions.

See 1.0.1.3.

1.0.1.1.0.2 Case (i) and Case (ii).

\(\E\Big[\big|\bss(\bZ) - \bss^{m}(\bZ)\big|_i\Big] \neq 0 \implies i \in \Paplus(I^m)\) is already shown above. For the reverse direction, we give the proof for soft interventions on additive noise models, Case (ii). We will use interventional regularity since 12 shows that it is satisfied for additive noise models. The proof for hard interventions, Case (i), follows from similar arguments since interventional regularity is also satisfied for hard interventions by 12.

1.0.1.1.0.3 Proof of \(i \in \Paplus(I^m) \implies \E\Big[\big|\bss(\bZ)-\bss^{m}(\bZ)\big|_i\Big] \neq 0\):

Note that the two score functions \(s\) and \(s^m\) are equal in their coordinate \(i\in \Paplus(\ell)\) only if \[\begin{align} 0 &= \frac{\partial \log q_{\ell}(z_{\ell} \mid \bz_{\Pa(\ell)})}{\partial z_i} - \frac{\partial \log p_{\ell}(z_{\ell} \mid \bz_{\Pa(\ell)})}{\partial z_i} = \frac{\partial}{\partial z_i} \log \frac{q_{\ell}(z_{\ell} \mid \bz_{\Pa(i)})}{p_{\ell}(z_{\ell} \mid \bz_{\Pa(\ell)})} \;. \label{eq:derivate-pq-parent-step} \end{align}\tag{74}\] However, 74 contradicts with interventional regularity. Therefore, if \(i \in \Paplus(\ell)\), \([\bss(\bz)_i]\) and \([\bss^{m}(\bz)]_i\) are not identical and by [prop:continuity-argument], \(\E\big[|\bss(\bZ)-\bss^{m}(\bZ)|_i \big] \neq 0\).

1.0.1.1.0.4 Case (iii) Coupled environments.

Suppose that \(I^m=\tilde{I}^m= \ell\). Following 14 , we have \[\begin{align} \bss^{m}(\bz) &= \nabla_{\bz} \log q_{\ell}(z_{\ell}) + \sum_{i \neq \ell} \nabla_{\bz} \log p_i(z_i \med \bz_{\Pa(i)}) \;, \tag{75} \\ and \quad \tilde{\bss}^{m}(\bz) &= \nabla_{\bz} \log \tilde{q}_{\ell}(z_{\ell}) + \sum_{i \neq \ell} \nabla_{\bz} \log p_i(z_i \med \bz_{\Pa(i)})\;. \tag{76} \end{align}\] Then, subtracting 76 from 75 and looking at \(i\)-th coordinate, we have \[\begin{align} \label{eq:coupled-hard-sz-sz-tilde} \big[\bss^{m}(\bz) - \tilde{\bss}^{m}(\bz)\big]_i &= \dfrac{\partial \log q_{\ell}(z_{\ell})}{\partial z_i} - \dfrac{\partial \log \tilde{q}_{\ell}(z_{\ell})}{\partial z_i} \;. \end{align}\tag{77}\] If \(i\neq \ell\), the right-hand side is zero and we have \(\big[\bss^{m}(\bz) - \tilde{\bss}^{m}(\bz)\big]_i=0\) for all \(\bz\). On the other hand, if \(i=\ell\), since \(q_{\ell}(z_{\ell})\) and \(\tilde{q}_{\ell}(z_{\ell})\) are distinct, there exists \(\bz \in \R^n\) such that \(q_{\ell}(z_{\ell})\neq \tilde{q}_{\ell}(z_{\ell})\). Subsequently, by [prop:continuity-argument], we have \(\E\big[|\bss^{m}(\bZ)-\tilde{\bss}^{m}(\bZ)|_i\big] \neq 0\).

1.0.1.1.0.5 Case (iv) Uncoupled environments.

Suppose that \(I^m = \ell\) and \(\tilde{I}^m= j\), and \(\ell \neq j\). Following 14 , we have \[\begin{align} \bss^{m}(\bz) &= \nabla_{\bz} \log q_{\ell}(z_{\ell}) + \nabla_{\bz} \log p_j(z_j \med \bz_{\Pa(j)}) + \sum_{k \in [n] \setminus \{\ell,j\}} \nabla_{\bz} \log p_k(z_k \med \bz_{\Pa(k}) \;, \tag{78} \\ and \quad \tilde{\bss}^{m}(\bz) &= \nabla_{\bz} \log q_j(z_j) + \nabla_{\bz} \log p_{\ell}(z_{\ell} \med \bz_{\Pa(\ell)}) + \sum_{k \in [n] \setminus \{\ell,j\}} \nabla_{\bz} \log p_k(z_k \med \bz_{\Pa(k)})\;. \tag{79} \end{align}\] Then, subtracting 79 from 78 we have \[\begin{align} \bss^{m}(\bz) - \tilde{\bss}^{m}(\bz) &= \; \nabla_{\bz} \log q_{\ell}(z_{\ell}) + \nabla_{\bz} \log p_j(z_j \med \bz_{\Pa(j)}) - \nabla_{\bz} \log q_j(z_j) - \nabla_{\bz} \log p_{\ell}(z_{\ell} \med \bz_{\Pa(\ell)}) \;. \label{eq:score-difference-proof-1} \end{align}\tag{80}\] Scrutinizing the \(i\)-th coordinate, we have \[\begin{align} \big[\bss^{m}(\bz) - \tilde{\bss}^{m}(\bz)\big]_i &= \dfrac{\partial \log q_{\ell}(z_{\ell})}{\partial z_i} + \dfrac{\partial \log p_j(z_j \med \bz_{\Pa(j)}) }{\partial z_i} - \dfrac{\partial \log q_j(z_j)}{\partial z_i} - \dfrac{\partial \log p_{\ell}(z_{\ell} \med \bz_{\Pa(\ell)})}{\partial z_i} \;. \label{eq:score-difference-proof-k} \end{align}\tag{81}\]

1.0.1.1.0.6 Proof of \(\E\Big[\big|\bss^{m}(\bZ) - \tilde{\bss}^{m}(\bZ)\big|_i\Big] \neq 0\;\implies i\in \Paplus(\ell,j)\):

Suppose that \(i \notin \Paplus({\ell},j)\). Then, none of the terms in the RHS of 81 is a function of \(z_i\). Therefore, all the terms in the RHS of 81 are zero, and we have \(\big[\bss^{m}(\bz) - \tilde{\bss}^{m}(\bz)\big]_i = 0\) for all \(\bz\). By [prop:continuity-argument], \(\E\big[|\bss^{m}(\bZ)-\tilde{\bss}^{m}(\bZ)|_i\big]=0\). This, equivalently, means that if \(\E\big[|\bss^{m}(\bZ)-\tilde{\bss}^{m}(\bZ)|_i\big] \neq 0\), then \(i \in \Paplus({\ell},j)\).

1.0.1.1.0.7 Proof of \(\E\Big[\big|\bss^{m}(\bZ)-\tilde{\bss}^{m}(\bZ)\big|_i\Big] \neq 0 \impliedby i\in \Paplus(\ell,j)\):

We prove it by contradiction. Assume that \(\big[\bss^{m}(\bz) - \tilde{\bss}^{m}(\bz)\big]_i = 0\) for all \(\bz\). Without loss of generality, let \(\ell \notin \Paplus(j)\).

1. If \(i=\ell\). In this case, 81 is simplified to \[\begin{align} 0 = \big[\bss^{m}(\bz) - \tilde{\bss}^{m}(\bz)\big]_{\ell} = \dfrac{\partial \log q_{\ell}(z_{\ell})}{\partial z_{\ell}} - \dfrac{\partial \log p_{\ell}(z_{\ell} \med \bz_{\Pa(\ell)})}{\partial z_{\ell}} \;. \label{eq:proof-k-equals-i} \end{align}\tag{82}\] If \(\ell\) is a root node, i.e., \(\Pa(\ell)=\emptyset\), 82 implies that \((\log q_{\ell})' (z_{\ell}) = (\log p_{\ell})'(z_{\ell})\) for all \(z_{\ell}\). Integrating, we get \(p_{\ell}(z_{\ell}) = \alpha q_{\ell}(z_{\ell})\) for some constant \(\alpha\). Since both \(p_{\ell}\) and \(q_{\ell}\) are pdfs, they both integrate to one, implying \(\alpha = 1\) and \(p_{\ell}(z_{\ell}) = q_{\ell}(z_{\ell})\), which contradicts the premise that observational and interventional mechanisms are distinct. If \(\ell\) is not a root node, consider some \(k \in \Pa(\ell)\). Then, taking the derivative of 82 with respect to \(z_k\), we have \[\begin{align} \label{eq:proof-pi-by-zi-zl} 0 = \dfrac{\partial^2 \log p_{\ell}(z_{\ell} \med \bz_{\Pa(\ell)})}{\partial z_{\ell} \partial z_k} \;. \end{align}\tag{83}\] Recall the equation \(Z_{\ell} = f_{\ell}(\bZ_{\Pa(\ell)}) + N_{\ell}\) for additive noise models specified in 7 . Denote the pdf of the noise term \(N_{\ell}\) by \(p_N\). Then, the conditional pdf \(p_{\ell}(z_{\ell} \med \bz_{\Pa(\ell)})\) is given by \(p_{\ell}(z_{\ell} \med \bz_{\Pa(\ell)}) = p_N(z_{\ell} - f_{\ell}(\bz_{\Pa(\ell)}))\). Denoting the score function of \(p_N\) by \(r_p\), \[\begin{align} \label{eq:score-of-additive-noise} r_p(u) \triangleq \dfrac{\rm d}{{\rm d} u} \log p_N(u) \;, \end{align}\tag{84}\] we have \[\begin{align} \label{eq:additive-component-score} \dfrac{\partial \log p_{\ell}(z_{\ell} \med \bz_{\Pa(\ell)})}{\partial z_{\ell}} = \dfrac{\partial \log p_{N}(z_{\ell} - f_{\ell}(\bz_{\Pa(\ell)}))}{\partial z_{\ell}} = r_p(z_{\ell} - f_{\ell}(\bz_{\Pa(\ell)})) \;. \end{align}\tag{85}\] Substituting this into 83 , we obtain \[\begin{align} 0 &= \dfrac{\partial r_p\big(z_{\ell} - f_{\ell}(\bz_{\Pa(\ell)})\big)}{\partial z_k} = - \dfrac{\partial f_{\ell}(\bz_{\Pa(\ell)})}{\partial z_k} \cdot r_p'\big(z_{\ell}-f_{\ell}(\bz_{\Pa(\ell)})\big) \;, \quad \forall \bz \in \R^n \;. \label{eq:proof-pi-by-zizl-additive} \end{align}\tag{86}\] Since \(k\) is a parent of \(\ell\), there exists a fixed \(\bZ_{\Pa(\ell)} = \bz_{\Pa(\ell)}^*\) realization for which \(\partial f_{\ell}(\bz_{\Pa(\ell)}^*)/\partial z_k\) is nonzero. Otherwise, \(f_{\ell}(\bz_{\Pa(\ell)})\) would not be sensitive to \(z_k\) which is contradictory to \(k\) being a parent of \(\ell\). Note that \(Z_{\ell}\) can vary freely after fixing \(\bZ_{\Pa(\ell)}\). Therefore, for 86 to hold, the derivative of \(r_p\) must always be zero. However, the score function of a valid pdf with full support cannot be constant. Therefore, \(\big[\bss^{m}(\bz)_i - \tilde{\bss}^{m}(\bz)\big]_i\) is not always zero, and we have \(\E\big[|\bss^{m}(\bZ) - \tilde{\bss}^{m}(\bZ)|_i\big] \neq 0\).

2. If \(i\neq \ell\). In this case, 81 is simplified to \[\begin{align} 0 = \big[\bss^{m}(\bz) - \tilde{\bss}^{m}(\bz)\big]_i = \dfrac{\partial \log p_j(z_j \med \bz_{\Pa(j)}) }{\partial z_i} - \dfrac{\partial \log q_j(z_j)}{\partial z_i} - \dfrac{\partial \log p_{\ell}(z_{\ell} \med \bz_{\Pa(\ell)})}{\partial z_i} \;. \label{eq:proof-k-notequals-i-1} \end{align}\tag{87}\] We investigate case by case and reach a contradiction for each case. First, suppose that \(i \notin \Pa({\ell})\). Then, we have \(i \in \Paplus(j)\), and 87 becomes \[\begin{align} 0 = \big[\bss^{m}(\bz) - \tilde{\bss}^{m}(\bz)\big]_i = \dfrac{\partial \log p_j(z_j \med \bz_{\Pa(j)}) }{\partial z_i} - \dfrac{\partial \log q_j(z_j)}{\partial z_i} \;. \label{eq:proof-k-notequals-i-2} \end{align}\tag{88}\] If \(i = j\), the impossibility of 88 directly follows from the impossibility of 82 . The remaining case is \(i \in \Pa(j)\). In this case, taking the derivative of the right-hand side of 88 with respect to \(z_j\), we obtain \[\begin{align} 0 = \dfrac{\partial^2 \log p_j(z_j \med \bz_{\Pa(j)})}{\partial z_i \partial z_j} \;, \end{align}\] which is a realization of 83 for \(i \in \Pa(j)\) and \(j\) in place of \(k \in \Pa(\ell)\) and \(\ell\), which we proved to be impossible in \(i=\ell\) case. Therefore, \(i \notin \Pa(\ell)\) is not viable. Finally, suppose that \(i \in \Pa({\ell})\). Then, taking the derivative of the right-hand side of 87 with respect to \(z_{\ell}\), we obtain \[\begin{align} \label{eq:proof-pi-by-zi-zk} 0 = \dfrac{\partial^2 \log p_{\ell}(z_{\ell} \med \bz_{\Pa(\ell)})}{\partial z_i \partial z_{\ell}} \;, \end{align}\tag{89}\] which is again a realization of 83 for \(k = i\), which we proved to be impossible.

Hence, we showed that \(\big[\bss^{m}(\bz) - \tilde{\bss}^{m}(\bz)\big]_i\) cannot be zero for all \(\bz\) values. Then, by [prop:continuity-argument] we have \(\E\big[|\bss^{m}(\bZ) - \tilde{\bss}^{m}(\bZ)|_i\big] \neq 0\), and the proof is concluded.

1.0.1.2 Proof of Lemma 2

Let us recall the setting. Consider random vectors \(\bY_1,\bY_2\in\R^r\) and \(\bW_1\), \(\bW_2 \in \R^s\) that are related through \(\bY_1=f(\bW_1)\) and \(\bY_2=f(\bW_2)\) such that \(r \geq s\), probability measures of \(\bW_1,\bW_2\) are absolutely continuous with respect to the \(s\)-dimensional Lebesgue measure and \(f: \R^s \to \R^r\) is an injective and continuously differentiable function.

In this setting, the realizations of \(\bW_1\) and \(\bY_1\), and that of \(\bW_2\) and \(\bY_2\), are related through \(\by = f(\bw)\). Since \(f\) is injective and continuously differentiable, volume element \({\rm d}\bw\) in \(\R^{s}\) gets mapped to \(\left|\det([J_{f}(\bw)]^{\top} \cdot J_{f}(\bw))\right|^{1/2}\;{\rm d}\bw\) on \(\image(f)\). Since \(\bW_1\) has density \(p_{\bW_1}\) absolutely continuous with respect to the \(s\)-dimensional Lebesgue measure, using the area formula [59], we can define a density for \(\bY_1\), denoted by \(p_{\bY_1}\), supported only on manifold \(\mcM \triangleq \image(f)\) which is absolutely continuous with respect to the \(s\)-dimensional Hausdorff measure: \[p_{\bY_1}(\by) = p_{\bW_1}(\bw) \cdot \left|\det([J_{f}(\bw)]^{\top} \cdot J_{f}(\bw))\right|^{-1/2} \;, \quad where \quad \by = f(\bw) \;. \label{eq:px-py-via-g-jacobian}\tag{90}\] Densities \(p_{\bY_2}\) and \(p_{\bW_2}\) of \(\bY_2\) and \(\bW_2\) are related similarly. Subsequently, score functions of \(\{\bW_1,\bW_2\}\) and \(\{\bY_1,\bY_2\}\) are specified similarly to 9 and 12 , respectively. Denote the Jacobian matrix of \(f\) at point \(\bw \in \R^{s}\) by \(J_{f}(\bw)\), which is an \(r \times s\) matrix with entries given by \[\require{physics} \big[J_{f}(\bw)\big]_{i,j} = \pdv{\big[f(\bw)\big]_i}{w_j}(x) = \pdv{y_i}{w_j} \;, \quad \forall i \in [r] \, , j \in [s] \;. \label{eq:g-jacobian-defn}\tag{91}\] Next, consider a function \(\phi \colon \mcM \to \R\). Since the domain of \(\phi\) is a manifold, its differential, denoted by \(D\phi\), is defined according to 11 . By noting \(\by = f(\bw)\), we can also differentiate \(\phi\) with respect to \(\bw \in \R^{s}\) as [49] \[\nabla_{w} \phi(\by) = \nabla_{w} (\phi \circ f)(\bw) = \big[J_{f}(\bw)\big]^{\top} \cdot D \phi(\by) \;. \label{eq:nabla-x-fy}\tag{92}\] Next, given the identities in 90 and 92 , we find the relationship between score functions of \(\bW_1\) and \(\bY_1\) as follows. \[\begin{align} \bss_{\bW_1}(\bw) &= \nabla_{\bw} \log p_{\bW_1}(\bw) \tag{93} \\ \overset{\eqref{eq:px-py-via-g-jacobian}}&{=} \nabla_{\bw} \log p_{\bY_1}(\by) + \nabla_{w} \log \left|\det([J_{f}(\bw)]^{\top} \cdot J_{f}(\bw))\right|^{1/2} \\ \overset{\eqref{eq:nabla-x-fy}}&{=} \big[J_{f}(\bw)\big]^{\top} \cdot D \log p_{\bY_1}(\by) + \nabla_{\bw} \log \left|\det([J_{f}(\bw)]^{\top} \cdot J_{f}(\bw))\right|^{1/2} \\ &= \big[J_{f}(\bw)\big]^{\top} \cdot \bss_{\bY_1}(\by) + \nabla_{w} \log \left|\det([J_{f}(\bw)]^{\top} \cdot J_{f}(\bw))\right|^{1/2} \;. \tag{94} \end{align}\] Following the similar steps that led to 94 for \(\bW_2\) and \(\bY_2\), we obtain \[\begin{align} \bss_{\bW_2}(\bw) &= \big[J_{f}(\bw)\big]^{\top} \cdot \bss_{\bY_2}(\by) + \nabla_{\bw} \log \left|\det([J_{f}(\bw)]^{\top} \cdot J_{f}(\bw))\right|^{1/2} \;. \label{eq:sxm-szm-for-g} \end{align}\tag{95}\] Subtracting 95 from 94 , we obtain the desired result \[\begin{align} \bss_{\bW_1}(\bw)-\bss_{\bW_2}(\bw) = \big[J_f(\bw)\big]^{\top} \cdot \big[\bss_{\bY_1}(\by) - \bss_{\bY_2}(\by)\big] \;. \label{eq:sw-from-sy-in-proof} \end{align}\tag{96}\]

1.0.1.2.0.1 Proof of the reverse direction.

Multiplying 96 from left with \(\big[[J_f(\bw)]^{\dagger}\big]^{\top}\), we obtain \[\Big[\big[J_f(\bw)\big]^{\dagger}\Big]^{\top} \cdot \big[\bss_{\bW_1}(\bw)-\bss_{\bW_2}(\bw)\big] = \Big[\big[J_f(\bw)\big]^{\dagger}\Big]^{\top} \cdot \big[J_f(\bw)\big]^{\top} \cdot \big[\bss_{\bY_1}(\by) - \bss_{\bY_2}(\by)\big] \;. \label{eq:jac-rel-left-multiplied}\tag{97}\] Note that \[\Big[\big[J_f(\bw)\big]^{\dagger}\Big]^{\top} \cdot \big[J_f(\bw)\big]^{\top} = J_f(\bw) \cdot \big[J_f(\bw)\big]^{\dagger} \;. \label{eq:jjw-col-proj-eq}\tag{98}\] By properties of the Moore-Penrose inverse, for any matrix \(\bA\), we have \(\bA \cdot \bA^{\dagger} \cdot \bA = \bA\). This means that \(\bA \cdot \bA^{\dagger}\) acts as a left identity for vectors in the column space of \(\bA\). By definition, \(\bss_{\bY_1}\) and \(\bss_{\bY_2}\) have values in \(T_{\bw} \image(f)\), the tangent space of the image manifold \(f\) at point \(\bw\). This space is equal to the column space of the matrix \(J_f(\bw)\). Therefore, \(J_f(\bw) \cdot [J_f(\bw)]^{\dagger}\) acts as a left identity for \(\bss_{\bY_1}(\by)\) and \(\bss_{\bY_2}(\by)\), and we have \[J_f(\bw) \cdot \big[J_f(\bw)\big]^{\dagger} \cdot \big[\bss_{\bY_1}(\by) - \bss_{\bY_2}(\by)\big] = \bss_{\bY_1}(\by) - \bss_{\bY_2}(\by) \;. \label{eq:sy-in-col-jfw}\tag{99}\] Substituting 98 and 99 into 97 completes the proof.

1.0.1.2.0.2 Proof of 1

For a given linear transform \(\bF\), we have \(J_f(\bw)=\bF\), which is independent of \(\bw\). Then, in the proof of 2, 94 reduces to \(\bss_{\bW}(\bw) = \bF^{\top} \cdot \bss_{\bY}(\by)\). Finally, we note that the score difference of \(\bY\) can be similarly written in terms of the score difference of \(\bW\).

1.0.1.3 Proof of Lemma 12

We will prove that causal mechanisms \(p_i\) and \(q_i\) satisfy interventional regularity for (i) hard interventions and (ii) additive noise models. To this end, we first define \[\begin{align} \psi(z_i, \bz_{\Pa(i)}) \triangleq \frac{q_i(z_i\mid \bz_{\Pa(i)})}{p_i(z_i\mid \bz_{\Pa(i)})} \;. \label{eq:h-pq-parent-dependence} \end{align}\tag{100}\] We start by showing that \(\psi(z_i, \bz_{\Pa(i)})\) varies with \(z_i\). We prove it by contradiction. Assume the contrary, i.e., let \(\psi(z_i, \bz_{\Pa(i)}) = \psi(\bz_{\Pa(i)})\). By rearranging 100 we have \[\begin{align} \label{eq:misc1} q_i(z_i \mid \bz_{\Pa(i)}) = \psi(\bz_{\Pa(i)}) \, p_i(z_i \mid \bz_{\Pa(i)}) \;. \end{align}\tag{101}\] Fix a realization of \(\bz_{\Pa(i)}^*\) and integrate both sides of 101 with respect to \(z_i\). Since both \(p_i\) and \(q_i\) are pdfs, we have \[\require{physics} \begin{align} 1 &= \int_{\R} q_i(z_i \mid \bz_{\Pa(i)}^*) \, \dd z_i = \int_{\R} \psi(\bz_{\Pa(i)}^*) \, p_i(z_i \mid \bz_{\Pa(i)}^*) \mathrm{d}{z_i} \dd z_i \\ &= \psi(\bz_{\Pa(i)}^*) \int_{\R} p_i(z_i \mid \bz_{\Pa(i)}^*) \dd z_i \\ &= \psi(\bz_{\Pa(i)}^*)\;. \end{align}\] This identity implies that \(p_i(z_i \mid \bz_{\Pa(i)}^*) = q_i(z_i \mid \bz_{\Pa(i)}^*)\) for any arbitrary realization \(\bz_{\Pa(i)}^*\). This contradicts the premise that observational and interventional distributions are distinct. As a result, to check if a model satisfies interventional regularity for node \(i\), it suffices to investigate whether the function \(\psi\) is not invariant with respect to \(z_k\) for \(k\in\Pa(i)\). To this end, from 100 we know that \(\psi(z_i, \bz_{\Pa(i)})\) varies with \(z_k\) if and only if \[\begin{align} \label{eq:int-regularity-suff} \frac{\partial}{\partial z_i} \log \psi(z_i, \bz_{\Pa(i)}) = \frac{\partial q_i(z_i \mid \bz_{\Pa(i)})}{\partial z_i} \cdot \frac{1}{q_i(z_i \mid \bz_{\Pa(i)})} - \frac{\partial p_i(z_i \mid \bz_{\Pa(i)})}{\partial z_i} \cdot \frac{1}{p_i(z_i \mid \bz_{\Pa(i)})} \neq 0 \;. \end{align}\tag{102}\] Next, we investigate the sufficient conditions listed in 12.

1.0.1.3.1 Hard Interventions

Under hard interventions, note that for any \(k\in\Pa(i)\), \[\begin{align} \frac{\partial}{\partial z_k} p_i(z_i \mid \bz_{\Pa(i)}) \neq 0 \;, \quad and \quad \frac{\partial}{\partial z_k} q_i(z_i) = 0\;. \end{align}\] Then, it follows directly from 102 that \[\begin{align} \frac{\partial}{\partial z_k} \log \psi(z_i, \bz_{\Pa(i)}) &= \underset{= \, 0}{\underbrace{\frac{\partial q_i(z_i)}{\partial z_k}}} \cdot \frac{1}{q_i(z_i)} - \underset{\neq \, 0}{\underbrace{\frac{\partial p_i(z_i \mid \bz_{\Pa(i)})}{\partial z_k}}} \cdot \frac{1}{p_i(z_i \mid \bz_{\Pa(i)})} \neq 0\;. \end{align}\] Thus, hard interventions on any latent causal model satisfy interventional regularity.

1.0.1.3.2 Additive Noise Models

The additive noise model for node \(i\) is given by \[Z_i = f_i(\bZ_{\Pa(i)}) + N_i\;, \label{eq:additive-model-intreg-proof-obs}\tag{103}\] as specified in 7 . When node \(i\) is soft intervened, \(Z_i\) is generated according to \[Z_i = \bar f_i(\bZ_{\Pa(i)}) + \bar N_i\;, \label{eq:additive-model-intreg-proof-int}\tag{104}\] in which \(\bar f_i\) and \(\bar N_i\) specify the interventional mechanism for node \(i\). Then, denoting the pdfs of \(N_i\) and \(\bar N_i\) by \(p_N\) and \(q_N\), respectively, 103 and 104 imply that \[\begin{align} p_i(z_i \mid \bz_{\Pa(i)}) = p_N\big(z_i - f_i(\bz_{\Pa(i)})\big)\;, \quad and \quad q_i(z_i \mid \bz_{\Pa(i)}) = q_N\big(z_i - \bar f_i(\bz_{\Pa(i)})\big)\;. \label{eq:additive-model-pq} \end{align}\tag{105}\] Denote the score functions associated with \(p_N\) and \(q_N\) by \[\begin{align} \label{eq:def-r-additive} r_p(u) \triangleq \frac{{\rm d}}{{\rm d} u} \log p_N(u) = \frac{p_N'(u)}{p_N(u)}\;, \quad and \quad r_q(u) \triangleq \frac{{\rm d}}{{\rm d} u} \log q_N(u) = \frac{q_N'(u)}{q_N(u)}\;. \end{align}\tag{106}\] We will prove that 102 holds by contradiction. Assume the contrary and let \[\begin{align} \label{eq:int-regularity-additive-contrary} \frac{\partial q_i(z_i \mid \bz_{\Pa(i)})}{\partial z_k} \cdot \frac{1}{q_i(z_i \mid \bz_{\Pa(i)})} = \frac{\partial p_i(z_i \mid \bz_{\Pa(i)})}{\partial z_k} \cdot \frac{1}{p_i(z_i \mid \bz_{\Pa(i)})}\;. \end{align}\tag{107}\] From 105 and 106 , for the numerators in 107 we have, \[\begin{align} \frac{\partial p_i(z_i \mid \bz_{\Pa(i)})}{\partial z_k} &= - \frac{\partial f_i(\bz_{\Pa(i)})}{\partial z_k} \cdot p_N'\big(z_i - f_i(\bz_{\Pa(i)})\big) \;, \\ and \qquad \frac{\partial q_i(z_i \mid \bz_{\Pa(i)})}{\partial z_k} &= - \frac{\partial \bar f_i(\bz_{\Pa(i)})}{\partial z_k} \cdot q_N'\big(z_i - \bar f_i(\bz_{\Pa(i)})\big) \;. \end{align}\] Hence, the identity in 107 can be written as \[\frac{\partial f_i(\bz_{\Pa(i)})}{\partial z_k} \cdot r_p\big(z_i - f_i(\bz_{\Pa(i)})\big) = \frac{\partial \bar f_i(\bz_{\Pa(i)})}{\partial z_k} \cdot r_q\big(z_i - \bar f_i(\bz_{\Pa(i)})\big) \;. \label{eq:pre-rearrange-additive-proof}\tag{108}\] Define \(n_i\) and \(\bar n_i\) as the realizations of \(N_i\) and \(\bar N_i\) when \(Z_i=z_i\) and \(\bZ_{\Pa(i)}=\bz_{\Pa(i)}\). By defining \(\delta(\bz_{\Pa(i)}) \triangleq f_i(\bz_{\Pa(i)}) - \bar f_i(\bz_{\Pa(i)})\), we have \(\bar n_i = n_i + \delta(\bz_{\Pa(i)})\). Then, 108 can rewritten as \[\label{eq:pre-rearrange-additive-proof2} \frac{\partial f_i(\bz_{\Pa(i)})}{\partial z_k} \cdot r_p(n_i) = \frac{\partial \bar f_i(\bz_{\Pa(i)})}{\partial z_k} \cdot r_q(n_i + \delta(\bz_{\Pa(i)})) \;.\tag{109}\] Note that \(\frac{\partial f_i(\bz_{\Pa(i)})}{\partial z_k}\) is a nonzero continuous function. Hence, there exists an interval \(\Phi\subseteq\R^{|\Pa(i)|}\) over which \(\frac{\partial f_i(\bz_{\Pa(i)})}{\partial z_k} \neq 0\). Likewise, \(r_p(n_i)\) cannot be constantly zero over all possible intervals \(\Omega \subseteq \R\). This is because otherwise, it would have to necessarily be a constant zero function (since it is analytic), which is an invalid score function. Hence, there exists an open interval \(\Omega \subseteq \mathbb{R}\) over which \(r_p(n_i)\) is nonzero for all \(n_i \in \Omega\). Then, we can rearrange 109 as \[\begin{align} \label{eq:rp-rq-equals-d-fq-d-fp} \frac{r_p(n_i)}{r_q(n_i+\delta(\bz_{\Pa(i)}))} &= \frac{\frac{\partial \bar f_i(\bz_{\Pa(i)})}{\partial z_k}}{\frac{\partial f_i(\bz_{\Pa(i)})}{\partial z_k}} \;, \quad \forall (n_i,\bz_{\Pa(i)})\in \Omega \times \Phi \;. \end{align}\tag{110}\] Note that the RHS of 110 is not a function of \(n_i\). Then, taking the derivative of both sides with respect to \(n_i\), we get \[\label{eq:d-log-rp-d-log-rq-with-delta} \frac{r_p'(n_i)}{r_p(n_i)} = \frac{r_q'(n_i + \delta(\bz_{\Pa(i)}))}{r_q(n_i + \delta(\bz_{\Pa(i)}))} \;.\tag{111}\] In the next step, we show that \(\delta\) is not a constant function. We prove this by contradiction. Suppose that \(\delta(\bz_{\Pa(i)}) = \delta^*\) is a constant function. Then, the gradients of \(f_i\) and \(\bar f_i\) are equal. From 110 , this implies that \[r_p(n_i) = r_q(n_i + \delta^*) \;, \quad \forall n_i \in \Omega \;.\] Since \(r_p(n_i)\) and \(r_q(n_i + \delta^*)\) are analytic functions that agree on an open interval of \(\R\), they are equal for all \(n_i \in \R\). This implies that \(p_N(n_i) = \eta \cdot q_N(n_i+\delta^*)\) for some constant \(\eta \in \R\). Since \(p_N\) and \(q_N\) are pdfs, \(\eta=1\) is the only choice that maintains \(p_N\) and \(q_N\) are pdfs. Therefore, \(p_N(n_i)=q_N(n_i+\delta^*)\). However, using 105 , this implies that \(p_i(z_i \mid \bz_{\Pa(i)}) = q_i(z_i \mid \bz_{\Pa(i)})\), which contradicts the premise that an intervention changes the causal mechanism of target node \(i\). Therefore, \(\delta\) is a continuous, non-constant function, and its image over \(\bz_{\Pa(i)}\in\Phi\) includes an open interval \(\Theta\subseteq\R\). With this result in mind, we return to 111 . Consider a fixed realization \(n_i = n_i^*\) and denote the value of the left-hand side (LHS) for \(n_i^*\) by \(\alpha\). By defining \(u \triangleq \delta(\bz_{\Pa(i)})\), we get \[\begin{align} \alpha &= \frac{r_q'(n_i^* + u)}{r_q(n_i^* + u)} \;, \qquad \forall u\in\Theta \;. \end{align}\] This is only possible if \(r_q\) is an exponential function, i.e., \(r_q(u) = k_1 \exp(\alpha u)\) over interval \(u \in \Theta\). Since \(r_q\) is an analytic function, it is, therefore, exponential over entire \(\R\). Then, the associated pdf must have the form \(q_N(u) = k_2 \exp((k_1/\alpha) \exp(\alpha u))\). However, the integral of this function over the entire domain \(\R\) diverges. Hence, it is not a valid pdf, rendering a contradiction. Hence, the additive noise model satisfies interventional regularity.

1.0.2 Proofs of the Results for Linear Transformations↩︎

1.0.2.1 Auxiliary Results

First, we provide a linear algebraic property, which will be used in the proofs.

First, note that \(\bL_{\sur} \preccurlyeq \bL_{\Pa} \preccurlyeq \bL_{\An}\). Hence, we start with considering a generic lower triangular matrix \(\bL\) such that \(\bI_{n \times n} \preccurlyeq \mathbb{1}\{\bL\} \preccurlyeq \bL_{\Pa}\). Matrix \(\bL\) can be decomposed as \[\bL = \bE \cdot (\bI_{n \times n} + \bLambda) \;,\] in which \(\bE\) is a diagonal matrix and \(\bLambda\) is a strictly lower triangular matrix that satisfies \(\mathbb{1}\{\bLambda\}\preccurlyeq \bL_{\Pa}\). Subsequently, since \(\bE\) is a diagonal matrix, we have \[\bL^{-1} = (\bI_{n \times n} + \bLambda)^{-1} \cdot \bE^{-1} \;, \quad and \quad \mathbb{1}\{\bL^{-1}\} = \mathbb{1}\big\{(\bI_{n \times n}+\bLambda)^{-1}\big\} \;.\] Note that \(\bLambda\) is a strictly lower triangular \(n \times n\) matrix, which implies that \(\bLambda^{n}\) is a zero matrix. Therefore, the inverse of \((\bI_{n \times n} + \bLambda)\) can be expanded as \[(\bI_{n \times n}+\bLambda)^{-1} = \bI_{n \times n} - \bLambda + \bLambda^{2} - \dots + (-1)^{n-1}\bLambda^{n-1} \;. \label{eq:inverse-expansion}\tag{112}\] If \(\big[\bI_{n \times n}+\bLambda\big]^{-1}_{i,j} \neq 0\) for some \(i, j \in [n]\) with \(j < i\), by 112 we have \([\bLambda^{k}]_{i,j} \neq 0\) for some \(k \in [n-1]\). Expanding matrix \(\bLambda^{k}\) yields that \([\bLambda^{k}]_{i,j}\) is equal to the sum of the products with \(k\) terms, i.e., \[[\bLambda^{k}]_{i,j} = \sum_{i > a_1 > \dots > a_{k-1} > j} [\bLambda]_{i,a_1} [\bLambda]_{a_1,a_2} \dots [\bLambda]_{a_{k-1},j} \;.\] Therefore, if \([\bLambda^{k}]_{i,j} \neq 0\), there exists a sequence of entries \(\big([\bLambda]_{i,a_1},[\bLambda]_{a_1,a_2},\dots,[\bLambda]_{a_{k-1},j}\big)\) in which all terms are nonzero. Next, we prove the two cases as follows.

1.0.2.1.0.1 Case 1:

Since \(\mathbb{1}\{\bLambda\}\preccurlyeq \bL_{\Pa}\), we have \([\bL_{\Pa}]_{i,a_1}=[\bL_{\Pa}]_{a_1,a_2}=\dots=[\bL_{\Pa}]_{a_{k-1},j}=1\). By the definition of \(\bL_{\Pa}\), this means that there is a path \(j \rightarrow a_{k-1} \rightarrow \dots \rightarrow a_1 \rightarrow i\) in \(\mcG\), which implies that \(j \in \An(i)\) and \([\bL_{\An}]_{i,j}=1\). Therefore, we conclude that \[\big[(\bI_{n \times n}+\bLambda)^{-1}\big]_{i,j}\neq 0 \quad \implies \quad [\bL_{\An}]_{i,j}=1 \;.\] Hence, \(\bI_{n \times n} \preccurlyeq \mathbb{1}\{\bL^{-1}\} \preccurlyeq \bL_{\An}\).

1.0.2.1.0.2 Case 2:

Since \(\mathbb{1}\{\bLambda\}\preccurlyeq \bL_{\sur}\), we have \([\bL_{\sur}]_{i,a_1}=[\bL_{\sur}]_{a_1,a_2}=\dots=[\bL_{\sur}]_{a_{k-1},j}=1\). By the definition of \(\bL_{\sur}\), this means that \[\Chplus(i) \subseteq \Ch(a_1) \subset \dots \subset \Ch(a_{k-1}) \subset \Ch(j) \;,\] which implies \([\bL_{\sur}]_{i,j}=1\) since \(\Chplus(i) \subseteq \Chplus(j)\). Therefore, we conclude that if \[\big[(\bI_{n \times n}+\bLambda)^{-1}\big]_{i,j}\neq 0 \quad \implies \quad [\bL_{\sur}]_{i,j}=1 \;.\] Hence, \(\bI_{n \times n} \preccurlyeq \mathbb{1}\{\bL^{-1}\} \preccurlyeq \bL_{\sur}\).

1.0.2.2 Proof of Theorem 3

1.0.2.2.0.1 Proof of the scaling consistency.

First, by 2, encoder estimate \(\hat{\bH}\) at the end of Stage L1 of 3 satisfies \[\begin{align} \label{eq:before-unmixing} \hat{\bH} \cdot \bG = \bP_{\mcI} \cdot \bL \;, \end{align}\tag{113}\] where \(\bP_{\mcI}\) is the permutation matrix of the intervention order \((I^1,\dots,I^n)\), \(\bL\) has nonzero diagonal entries and satisfies \(\bL_{i,j} = 0\) for all \(j \notin \Paplus(i)\). We will show that the output of Stage L3 of 3 for hard interventions satisfies that \[\begin{align} \label{eq:unmixing95goal95node} \mathbb{1}\big([\hat{\bH} \cdot \bG]_m \big) = \be_{I^m}^{\top} \;, \quad \forall m \in [n] \;, \end{align}\tag{114}\] which will imply scaling consistency as \[\hat{\bZ}(\bX; \hat{\bH}) = \hat{\bH} \cdot \bX = \hat{\bH} \cdot \bG \cdot \bZ = \bP_{\mcI} \cdot \bC_{\rm s} \cdot \bZ \;,\] where \(\bC_{\rm s}\) is a constant diagonal matrix with nonzero diagonal entries. We prove this as follows.

First, consider the topological order \(\pi\) of \(\hat{\mcG}\). By 2, transitive closures of \(\hat{\mcG}\) and \(\mcG\) are the same under relabeling of the nodes with permutation \((1,\dots,n) \rightarrow (I^1,\dots,I^n)\). Then, since \((1,\dots,n)\) is assumed to be topological order of \(\mcG\), \(\rho \triangleq (I^{\pi_1},\dots,I^{\pi_n})\) is also a topological order of \(\mcG\). According to this notation, in the \(\pi_t\)-th environment, node \(\rho_t\) is intervened.

Next, since \(\rho_1\) has no parent in \(\mcG\), by 2 we already know that \(\hat{Z}_{\pi_1} = b_1 \cdot Z_{\rho_1}\) for some constant \(b_1 \in \R_+\). Consider \(m = \pi_k\) step of the algorithm, and assume that for all \(t \in [k-1]\), the updated encoder rows satisfy \[\label{eq:unmixing-induction-step} \hat{Z}_{\pi_t} = \hat{\bH}_{\pi_t} \cdot \bG \cdot \bZ = b_t \cdot Z_{\rho_t} \;,\tag{115}\] for some constant \(b_t \in \R_+\). We will show that the update in step \(m = \pi_k\) will ensure that \(\hat{Z}_{\pi_k}\) also satisfies 115 . Before the update, by 2 we know that \(\hat{Z}_{\pi_k}\) is a linear function of \(Z_{\rho_k}\) and \(\bZ_{\Pa(\rho_k)}\). Then, using 115 , we have \[\hat{Z}_{\pi_k} = c_0 \cdot Z_{\rho_k} + \bc \cdot \hat{\bZ}_{\hat{\Pa}(\pi_k)} \;,\] for \(|\hat{\Pa}(\pi_k)|\)-dimensional vector \(\bc\) and nonzero constant \(c_0 \in \R_+\). We will use [fact:hard-indep] in environment \(\mcE^{\pi_k}\). For brevity, using \(m = \pi_k\) and letting \(\bV \triangleq \hat{\bZ}^{m}_{\hat{\Pa}(\pi_k)}\), \[\hat{Z}^{m}_{m} = c_0 \cdot Z^{m}_{\rho_k} + \bc \cdot \bV \;.\] Since node \(\rho_k\) is intervened in environment \(\mcE^m\), by [fact:hard-indep], we know that \(Z^{m}_{\rho_k} \ci \bV\), which implies that \[\Cov(Z^{m}_{\rho_k}, \bV) = \boldsymbol{0}_{1 \times |\bV|} \;.\] Consider a \(|\bV|\)-dimensional row vector \(\bu\) and let \[Y = \big(\hat{\bH}_{m} - \bu \cdot \hat{\bH}_{\hat{\Pa}(m)}\big) \cdot \bX^{m} = \hat{Z}^{m}_{m} - \bu \cdot \bV \;.\] Note that \(\bu = \bc\) would yield that \(Y = c_0 \cdot Z^{m}_{\rho_k}\), and subsequently \(Y \ci \bV\). On the other hand, \(Y \ci \bV\) implies that \[\boldsymbol{0}_{1 \times |\bV|} = \Cov(Y, \bV) = \Cov(\hat{Z}^{m}_{m} - \bu \cdot \bV, \bV) \;,\] which has a unique solution \[\bu = \Cov(\hat{Z}^{m}_{m}, \bV) \cdot \big(\Cov(\bV)\big)^{-1} \;.\] i.e., linear minimum mean square error (LMMSE) estimator. Finally, \(\Cov(\bV)\) is invertible since the causal relationships among the entries in \(\hat{Z}^{m}_{\hat{\Pa}(\pi_k)}\) are not deterministic. Using this row vector \(\bu\), we update \[\hat{\bH}_m \gets \bH_m - \bu \cdot \bH_{\hat{\Pa}(\pi_k)} \;,\] and achieve \(\hat{Z}_m = c_0 \cdot Z_{\rho_k}\). Therefore, 115 holds for \(t= k\) as well, and by induction we obtain \[\label{eq:z-hat-z-unmixed} \hat{\bZ} = \hat{\bH} \cdot \bG \cdot \bZ = \bP_{\mcI} \cdot \bC_{\rm s} \cdot \bZ \;,\tag{116}\] for a diagonal matrix \(\bC_{\rm s}\) with nonzero diagonal entries.

1.0.2.2.0.2 Proof of the perfect DAG recovery.

We show that the graph construction process in 40 achieves perfect DAG recovery as follows. First, using 2 and 116 , we have \[\begin{align} \bsd_{\hat{\bZ}}^{m}(\hat{\bZ}; \hat{\bH}) = (\bP_{\mcI} \cdot \bC_{\rm s})^{-\top} \cdot \bsd_{\bZ}^{m}(\bZ) = \bP_{\mcI} \cdot \bC_{\rm s}^{-\top} \cdot \bsd_{\bZ}^{m}(\bZ) \;. \end{align}\] Subsequently, using the fact that \(\bC_{\rm s}\) is a diagonal matrix, for all \(i, m \in [n]\) we obtain \[\begin{align} \big[\bsd_{\hat{\bZ}}^{m}(\hat{\bZ}; \hat{\bH})\big]_i \neq 0 \; \; \iff \;\; \big[\bsd_{\bZ}^{m}(\bZ)\big]_{I^i} \neq 0 \;. \end{align}\] Then, 1 and 40 , we have \[\begin{align} I^i \in \Paplus(I^m) \; &\iff \; \E\Big[\big|\bsd_{\bZ}^{m}(\bZ)\big|\Big]_{I^i} \neq 0 \; \iff \; \E\Big[\big|\bsd_{\hat{\bZ}}^{m}(\hat{\bZ}; \hat{\bH})\big|_i \Big] \; \iff \; i \in \hat{\Paplus}(m) \;. \end{align}\] This concludes the proof that \(\hat{\mcG}\) and \(\mcG\) are related through a graph isomorphism by permutation \(\mcI\), which denotes the intervention order.

1.0.2.3 Proof of Theorem 4

1.0.2.3.0.1 Latent graph estimation.

We continue from the intuitions provided in 0.0.5.4 and prove the perfect recovery of the latent graph as follows. First, recall that output \(\hat{\mcG}\) of 3 for soft interventions satisfies that transitive closures of \(\hat{\mcG}\) and \(\mcG\) are the same under permutation \((I^1,\dots,I^n)\). Since \((1,\dots,n)\) is assumed to be a topological order of \(\mcG\), \(\rho \triangleq (I^{\pi_1},\dots,I^{\pi_n})\) is also a topological order of \(\mcG\). Then, for any \(k = \pi_u\) and \(t = \pi_v \in \{\pi_1,\dots,\pi_{u-1}\}\), 42 becomes \[\label{eq:full-rank-parent-dim} \dim\big( \colspace(\bR_{\bX}^{t}) \cap \colspace(\bR_{\bX}^{k})\big) = |\Pa(\rho_u) \cap \Pa(\rho_v)| + \mathbb{1}\big\{\rho_v \in \Pa(\rho_u) \big\} \;.\tag{117}\] By induction, we will prove that the parent sets generated by 4 satisfy \[\label{eq:full-rank-parent-condition} \pi_u \in \hat{\Pa}(\pi_v) \; \iff \; \rho_u \in \Pa(\rho_v) \;.\tag{118}\] At the base case, consider \(k=\pi_2\). The only possible parent is \(t=\pi_1\). Since \(\rho_1\) is a root node in \(\mcG\), 117 implies \[\dim\big( \colspace(\bR_{\bX}^{\pi_1}) \cap \colspace(\bR_{\bX}^{\pi_2})\big) = \mathbb{1}\{ \rho_1 \in \Pa(\rho_2) \} \;.\] Since \(\hat{\Pa}(\pi_1) = \hat{\Pa}(\pi_2) = \emptyset\) at this step, the algorithm adds \(\pi_1\) to \(\hat{\Pa}(\pi_2)\) if and only if \(\rho_1 \in \Pa(\rho_2)\), and 118 is satisfied for \(k=\pi_2\). Next, as the induction step, assume that the algorithm output \(\hat{\Pa}(k)\) satisfies 118 for \(k \in \{\pi_1,\dots,\pi_{u-1}\}\). We prove that \(k = \pi_u\) also satisfies the condition by induction again. At the base case, consider \(t = \pi_1\). Similar to the previous base case, 117 implies \[\dim\big( \colspace(\bR_{\bX}^{\pi_1}) \cap \colspace(\bR_{\bX}^{\pi_u})\big) = \mathbb{1}\{ \rho_1 \in \Pa(\rho_u) \} \;,\] and since \(\hat{\Pa}(\pi_u)\) does not contain \(\pi_1\) yet, the algorithm adds \(\pi_1\) to \(\hat{\Pa}(\pi_u)\) if and only if \(\rho_1 \in \Pa(\rho_u)\). Next, as the induction step, assume that for \(\pi_j \in \{\pi_1,\dots,\pi_{v-1}\}\) where \(v < u\), the algorithm correctly identified the parents, i.e., \(\pi_j \in \hat{\Pa}(\pi_u)\) if and only if \(\rho_j \in \Pa(\rho_u)\). Consider \(k = \pi_u\) and \(t = \pi_v\), for which 117 implies \[\dim\big( \colspace(\bR_{\bX}^{\pi_v}) \cap \colspace(\bR_{\bX}^{\pi_u})\big) = |\Pa(\rho_u) \cap \Pa(\rho_v)| + \mathbb{1}\{ \rho_v \in \Pa(\rho_u)\} \;.\] By the induction hypothesis, we know that \(|\hat{\Pa}(\pi_u) \cap \hat{\Pa}(\pi_v)| = |\Pa(\rho_u) \cap \Pa(\rho_v)|\). Therefore, the algorithm adds \(\pi_v\) to \(\hat{\Pa}(\pi_v)\) if and only if \(\rho_v \in \Pa(\rho_u)\). Therefore, by the inner induction, \(\hat{\Pa}(k)\) satisfies 118 . Then, by the outer induction, 118 is satisfied for all parent set estimates, which concludes the proof of perfect graph recovery.

1.0.2.3.0.2 Markov property of consistency up to mixing with surrounding parents.

Let \(\hat{\bZ}\) denote \(\hat{\bZ}(\bX; \hat{\bH})\). We further investigate the properties of consistency up to mixing with surrounding parents to prove that \(\hat{\bZ}\) is Markov with respect to \(\hat{\mcG}\). Using the results of the encoder estimation step, we have \[\label{eq:hat-z-surrounding} \hat{\bZ} = \bP_{\mcI} \cdot \bC_{\rm sur} \cdot \bZ \;,\tag{119}\] where \(\bC_{\rm sur}\) has nonzero diagonal entries and satisfies \([\bC_{\rm sur}]_{i,j} = 0\) for all \(j \notin \overline{\sur}(i)\). Recall the following SCM specified in 6 \[Z_i = f_i(\bZ_{\Pa(i)},N_i) \;.\] Let \(\tau\) be the permutation that maps \(\{1,\dots,n\}\) to \(\mcI\), i.e., \(I^{\tau_i}=i\) for all \(i\in[n]\). Since \(\sur(i) \subseteq \Pa(i)\), 119 implies \[\begin{align} \hat{Z}_{\tau_i} &= [\bC_{\sur}]_{i,i} \cdot Z_i + \sum_{j \in \sur(i)} [\bC_{\sur}]_{i,j} \cdot Z_j \tag{120} \\ &= [\bC_{\sur}]_{i,i} \cdot f_i(\bZ_{\Pa(i)},N_i) + \sum_{j \in \sur(i)} [\bC_{\sur}]_{i,j} \cdot Z_j \;, \tag{121} \end{align}\] which is a function of \(\bZ_{\Pa(i)}\) and \(N_i\). To prove Markov property, we need to specify \(Z_{\tau_i}\) in terms of \(\hat{\bZ}_{\hat{\Pa}(\tau_i)}\). To this end, it suffices to show that for any \(k \in \Pa(i)\), \(Z_k\) is a function of \(\hat{\bZ}_{\hat{\Pa}}(\tau_i)\). Using 119 , we have \[\bZ = \bC_{\sur}^{-1} \cdot \bP_{\tau} \cdot \hat{\bZ} \;. \label{eq:z-from-hat-z-pre}\tag{122}\] By 13, we know that \([\bC_{\sur}^{-1}]_{i,j}=0\) if \(j \notin \sur(i)\) for distinct \(i\) and \(j\). Subsequently, as counterpart of 120 , we have \[\label{eq:rho-next-4} Z_k = [\bC_{\sur}^{-1}]_{k,k} \cdot \hat{Z}_{\tau_k} + \sum_{j \in \sur(k)} [\bC_{\sur}^{-1}]_{k,j} \cdot \hat{Z}_{\tau_j} \;.\tag{123}\] Since \(k \in \Pa(i)\), we have \(\tau_k \in \Pa(\tau_i)\) and \(\hat{Z}_{\tau_k}\) is in \(\hat{\bZ}_{\hat{\Pa}(\tau_i)}\). Note that if \(j \in \sur(k)\), \(j\) is also in \(\Pa(i)\). Therefore, every term in the RHS of 123 belongs to \(\hat{\bZ}_{\hat{\Pa}(\tau_i)}\) and \(Z_k\) is a function of \(\hat{\bZ}_{\hat{\Pa}(\tau_i)}\). Then, using 121 , we know that \(\hat{Z}_{\tau_i}\) is a function of only \(\hat{\bZ}_{\hat{\Pa}(\tau_i)}\) and \(N_i\), which concludes the proof that \(\hat{\bZ}\) is Markov with respect to \(\hat{\mcG}\).

1.0.3 Proofs of the Results for General Transformations↩︎

The proof of 6 is already given in 0.0.6.2. In this section, we first prove 7, that is, identifiability for coupled interventions without observational environment. Then, we prove the results in 0.0.6.3 related to the uncoupled interventions.

1.0.3.1 Proof of Theorem 7

We will show that if \(p(\bZ)\) is adjacency-faithful to \(\mcG\) and the latent causal model is an additive noise model, then we can recover \(\mcG\) without having access to observational environment \(\mcE^0\). Let \(\tau\) be the permutation that maps \(\{1,\dots,n\}\) to \(\mcI\), i.e., \(I^{\tau_i}=i\) for all \(i\in[n]\) and \(\bP_{\tau}\) to denote the permutation matrix that corresponds to \(\tau\), i.e., \[\begin{align} \label{eq:intervention-order} [\bP_{\tau}]_{i,m} = \begin{cases} 1 \;, &m = \rho_i \;, \\ 0 \;, & \text{else} \;. \end{cases} \end{align}\tag{124}\] By 1(iv), true latent score changes across \(\{\mcE^{\tau_i},\tilde{\mcE}^{\tau_j}\}\) gives us \(\Paplus(i,j)\) for \(i\neq j\). First, we use the perfect latent recovery result to show that 1(iv) also applies to estimated latent score changes. Denote \(f = \hat{h} \circ g\). Using score transform between \(\bZ\) and \(\hat{\bZ}\) in \(\eqref{eq:score-difference-z-zhat}\), and recalling \(\mathbb{1}\{J_{f}^{-1}\}=\bP_{\tau}\), we have \[\begin{align} \big[\bss_{\hat{\bZ}}^{\tau_i}(\hat{\bz};\hat{h}) - \tilde{\bss}_{\hat{\bZ}}^{\tau_j}(\hat{\bz}; \hat{h})\big]_{\tau_k} &= \big[J_{f}^{-\top}(\bz)\big]_{\tau_k} \cdot \big[\bss^{\tau_i}(\bz) - \tilde{\bss}^{\tau_j}(\bz)\big] \\ &= \big[J_{f}^{-\top}(\bz)\big]_{\tau_k,k} \cdot \big[\bss^{\tau_i}(\bz) - \tilde{\bss}^{\tau_j}(\bz)\big]_k \;. \end{align}\] Using \(\big[J_{f}^{-\top}(\bz)\big]_{\tau_k,k} \neq 0\) for all \(\bz \in \R^n\), we have \[\begin{align} \label{eq:graph-iso-step-1} \E\Big[\big|\bss_{\hat{\bZ}}^{\tau_i}(\hat{\bZ};\hat{h}) - \tilde{\bss}_{\hat{\bZ}}^{\tau_j}(\hat{\bZ}; \hat{h})\big|_{\tau_k}\Big] \neq 0 \quad \iff \quad \E\Big[\big|\bss^{\tau_i}(\bZ) - \tilde{\bss}^{\tau_j}(\bZ)\big|_{\tau_k}\Big] \neq 0 \;. \end{align}\tag{125}\] Hence, by 1(iv), \[\begin{align} \label{eq:graph-iso-step-2} \E\Big[\big|\bss_{\hat{\bZ}}^{\tau_i}(\hat{\bZ};\hat{h}) - \tilde{\bss}_{\hat{\bZ}}^{\tau_j}(\hat{\bZ}; \hat{h})\big|_{\tau_k}\Big] \neq 0 \quad \iff \quad k \in \Paplus(i,j)\;. \end{align}\tag{126}\] Let us denote the graph \(\mcG_{\tau}\) that is related to \(\mcG\) by permutation \(\tau\), i.e., \(i\in\Pa(j)\) if and only if \(\tau_i \in \Pa_{\tau}(\tau_j)\) for which \(\Pa_{\tau}(\tau_j)\) denotes the parents of node \(\tau_j\) in \(\mcG_{\tau}\). Using 126 , we have \[\begin{align} \E\Big[\big|\bss_{\hat{\bZ}}^{\tau_i}(\hat{\bZ};\hat{h}) - \tilde{\bss}_{\hat{\bZ}}^{\tau_j}(\hat{\bZ}; \hat{h})\big|_{\tau_k}\Big] \neq 0 \quad \iff \quad \tau_k \in \Paplus_{\tau}(\tau_i,\tau_j)\;. \end{align}\] In the rest of the proof, we will show how to obtain \(\{\Pa_{\tau}(i) : i \in [n]\}\) using \(\{\Paplus_{\tau}(i,j) : i,j \in [n], \; i \neq j\}\). Since \(\mcG_{\tau}\) is a graph isomorphism of \(\mcG\), this problem is equivalent to obtaining \(\{\Pa(i) : i \in [n]\}\) using \(\{\Paplus(i,j) : i,j \in [n], \; i \neq j\}\). Note that \(\hat{Z}_i\) (which corresponds to node \(i\) in \(\mcG_{\tau}\)) is intervened in environments \(\mcE^i\) and \(\tilde{\mcE}^i\). We denote the set of root nodes by \[\mcK \triangleq \{i\in[n] : \Pa(i) = \emptyset \} \;,\] and also define \[\mcB_i \triangleq \, \bigcap_{j \neq i} \Paplus(i,j) \;, \quad \forall i \in [n] \;, \quad and \quad \mcB \triangleq \{i : |\mcB_i|=1\} \;.\] Note that \(\Paplus(i) \subseteq \mcB_i\). Hence, \(|\mcB_i|=1\) implies that \(i\) is a root node. We investigate the graph recovery in three cases.

1. \(|\mcB|\geq 3\): For any node \(i\in[n]\), we have \[\begin{align} \Paplus(i) \subseteq \mcB_i \subseteq \bigcap_{j \in \mcK \setminus \{i\}} \Paplus(i,j) = \Paplus(i) \cup \Big\{\bigcap_{j \in \mcK \setminus \{i\}} \{ j \} \Big\} = \Paplus(i) \;. \end{align}\] Note that the last equality is due to \(\textstyle\bigcap_{j \in \mcK \setminus \{i\}} \{ j \} = \emptyset\) since there are at least two root nodes excluding \(i\). Then, \(\mcB_i = \Paplus(i)\) for all \(i \in [n]\) and we are done.

2. \(|\mcB|=2\): The two nodes in \(\mcB\) are root nodes. If there were at least three root nodes, we would have at least three nodes in \(\mcB\). Hence, the two nodes in \(\mcB\) are the only root nodes. Subsequently, every \(i \notin \mcB\) is also not in \(\mcK\) and we have \[\begin{align} \Paplus(i) \subseteq \mcB_i \subseteq \bigcap_{j \in \mcK} \Paplus(i,j) = \Paplus(i) \cup \Big\{\bigcap_{j \in \mcK}\{j\}\Big\} = \Paplus(i) \;. \end{align}\] Hence, \(\mcB_i = \Paplus(i)\) for every non-root node \(i\) and we already have the two root nodes in \(\mcB\), which completes the graph recovery.

3. \(|\mcB|\leq 1\): First, consider all \((i,j)\) pairs such that \(|\Paplus(i,j)|=2\). For such an \((i,j)\) pair, at least one of the nodes is a root node; otherwise \(\Paplus(i,j)\) would contain a third node. Using these pairs, we identify all root nodes as follows. Note that a hard intervention on node \(i\) makes \(Z_i\) independent of all of its non-descendants, and all conditional independence relations are preserved under componentwise diffeomorphisms such as \(f\). Then, using the adjacency-faithfulness assumption, we infer that

   • if \(\hat{Z}_i \ci \hat{Z}_j\) in \(\mcE^i\) and \(\hat{Z}_i \ci \hat{Z}_j\) in \(\tilde{\mcE}^j\), then both \(i\) and \(j\) are root nodes.

   • if \(\hat{Z}_i \notci \hat{Z}_j\) in \(\mcE^i\), then \(i \rightarrow j\) and \(i\) is a root node.

   • if \(\hat{Z}_i \notci \hat{Z}_j\) in \(\tilde{\mcE}^j\), then \(j \rightarrow i\) and \(j\) is a root node.

   This implies that we can determine whether \(i\) and \(j\) nodes are root nodes by using at most two independence tests. Hence, we identify all root nodes by using at most \(n\) independence tests. We also know that there are at most two root nodes. If we have two root nodes, then \(\mcB_i = \Paplus(i)\) for all non-root nodes, and the graph is recovered. If we have only one root node \(i\), then for any \(j \neq i\) we have \[\begin{align} \Paplus(j) \subseteq \mcB_j \subseteq \Paplus(i,j) = \Paplus(j) \cup \{i\} \;. \end{align}\] Finally, if \(\hat{Z}_j \ci \hat{Z}_i \mid \{\hat{Z}_\ell : \ell \in \mcB_j \setminus \{i\}\}\) in \(\tilde{\mcE}^j\), we have \(i \notin \Paplus(j)\) due to adjacency-faithfulness. Otherwise, we conclude that \(i \in \Paplus(j)\). Hence, an additional \((n-1)\) conditional independence tests ensure the recovery of all \(\Paplus(j)\) sets, and the graph recovery is complete.

1.0.3.2 Proof of Lemma 10

Let us start by scrutinizing the constraints. For the true encoder \(g^{-1}\), 1 gives us \[\begin{align} [\bD(g^{-1})]_{i,m} \neq 0 \;\; &\iff \;\; i \in \Paplus(I^m) \;, \\ and \quad [\tilde{\bD}(g^{-1})]_{i,m} \neq 0 \;\; &\iff \;\; i \in \Paplus(\tilde{I}^m) \;. \end{align}\] Therefore, using \(\tilde{\mcI} = \sigma \circ \mcI\), we have \(\mathbb{1}\{\tilde{\bD}(g^{-1})\} = \mathbb{1}\{\bD(g^{-1})\} \cdot \bP_{\sigma}\). Note that acyclicity of graph \(\mcG\) also implies that \(\mathbb{1}\{\bD(g^{-1})\} \odot \mathbb{1}\{\bD(g^{-1})^{\top}\} = \bI_{n \times n}\).

Next, note that \(\bD_{\rm t}(g^{-1};\sigma)\) is equal to the true score change matrix \(\bD_{\rm t}(g^{-1})\) defined for the coupled interventions, which satisfies \(\mathbb{1}\{\bD_{\rm t}(g^{-1};\sigma)\} = \bP_{\mcI}^{\top}\). Let \(h^* = \big[\bD_{\rm t}(g^{-1};\sigma)\big]^{\top} \cdot g^{-1}\). Since \(\bD_{\rm t}(g^{-1};\sigma)\) is a scaled permutation matrix, using 48 we obtain \[\begin{align} \bD_{\rm t}(h^*;\sigma) &= [\bD_{\rm t}(g^{-1};\sigma)]^{-1} \cdot \bD_{\rm t}(g^{-1};\sigma) = \bI_{n \times n} \;, \\ \bD(h^*) &= [\bD_{\rm t}(g^{-1};\sigma)]^{-1} \cdot \bD(g^{-1}) \\ \tilde{\bD}(h^*) &= [\bD_{\rm t}(g^{-1};\sigma)]^{-1} \cdot \tilde{\bD}(g^{-1}) \;. \end{align}\] Thus, using \(\mathbb{1}\{\tilde{\bD}(g^{-1})\} = \mathbb{1}\{\bD(g^{-1})\} \cdot \bP_{\sigma}\), we have \[\begin{align} \mathbb{1}\{\tilde{\bD}(h^*)\} = \mathbb{1}\{\bD(h^*)\} \cdot \bP_{\sigma} \;, \end{align}\] which satisfies the first constraint. Also, using \(\mathbb{1}\{\bD(h)\}= \bP_{\sigma} \cdot \mathbb{1}\{\tilde{\bD}(h)\}\), \[\mathbb{1}\{\bD(h^*)\} = \bP_{\mcI} \cdot \mathbb{1}\{\bD(g^{-1})\} \;.\] Therefore, we have \(\mathbb{1}\{\bD(h^*)\} \odot \mathbb{1}\{\bD(h^*)^{\top}\} = \bI_{n \times n}\). Finally, \((h^*)^{-1} \triangleq g \cdot \big[\bD_{\rm t}(g^{-1};\sigma)\big]^{-\top}\) satisfies \((h^*)^{-1}\circ h(\bX)=\bX\). Hence, \((h^*, \sigma)\) minimizes the objective in 57 at value zero.

1.0.3.3 Proof of Lemma 11

We will prove the desired result by contradiction. Suppose that \((h^*, \pi)\) is a solution to the optimization problem specified in 57 and let \(\pi \neq \sigma\). Denote \(f = h^* \circ g\), so we have \(\hat{\bZ} = f(\bZ)\). If \(\pi_m = \sigma_m\) for some \(m \in [n]\), this means \(I^m = \tilde{I}^{\pi_m} = \ell\) for some node \(\ell \in [n]\). We follow the same approach in the proof of 5. Specifically, using score difference transformation property in 30 and one-sparse property of \([\bss^{m}(\bZ) - \tilde{\bss}^{\pi_m}(\bZ)]\) via 1(iii), we have \[\begin{align} [\bD_{\rm t}(h^*)]_{i,m} &= \E \Big[\big| \bss_{\hat{\bZ}}^{m}(\hat{\bZ}) - \tilde{\bss}_{\hat{\bZ}}^{\pi_m}(\hat{\bZ}) \big|_i \Big] \\ &= \E \bigg[ \Big| \big[J_{f}^{-\top}(\bZ)\big]_i \cdot \big[\bss^{m}(\bZ) - \tilde{\bss}^{\pi_m}(\bZ)\big] \Big| \bigg] \\ &= \E \bigg[ \Big[J_{f}^{-1}(\bZ)\big]_{\ell,i} \cdot \big[\bss^{m}(\bZ) - \tilde{\bss}^{\pi_m}(\bZ)\big]_{\ell} \Big] \;. \end{align}\] By definition of \(h^*\), \([\bD_{\rm t}(h^*)]_{i,m} = 0\) for all \(i \neq m\). Also, interventional discrepancy between \(q_\ell(z_\ell)\) and \(\tilde{q}_\ell(z_\ell)\) implies that \(\big[\bss^{m}(\bz) - \tilde{\bss}^{\pi_m}(\bz)\big]_{\ell} \neq 0\) except for a null set. Then, \([\bD_{\rm t}(h^*)]_{i,m} = 0\) implies that \(\big[J_{f}^{-1}(\bz)\big]_{\ell,i} = 0\) except for a null set. Since \(J_{f}^{-1}\) is a continuous function, this implies that \(\big[J_{f}^{-1}(\bz)\big]_{\ell,i} = 0\) for all \(\bz \in \R^n\). Furthermore, since \(J_{f}^{-1}\) is invertible for all \(\bz\), none of its columns can be a zero vector. Hence, for all \(\bz \in \R^n\), \(\big[J_{f}^{-1}(\bz)\big]_{\ell,m} \neq 0\). To summarize, if \(\pi_m = \sigma_m\), then \[\begin{align} \label{eq:matching-node-column} \forall \bz\in \R^n \;\; \big[J_f^{-1}(\bz)\big]_{I^m,i} \neq 0 \quad \iff \quad i = m \;. \end{align}\tag{127}\] Next, consider the set of mismatched nodes \[\begin{align} \mcA \triangleq \{I^m : I^m \neq \tilde{I}^{\pi_m}\} \;. \end{align}\] Let \(I^a \in \mcA\) be a non-descendant of all the other nodes in \(\mcA\). There exist nodes \(I^b, I^c \in \mcA\), not necessarily distinct, such that \[\label{eq:rho-abc} I^a = \tilde{I}^{\pi_b} \;, \quad and \quad I^c = \tilde{I}^{\pi_a} \;.\tag{128}\] In four steps, we will show that \([\bD(h^*)]_{a,b}\neq 0\) and \([\bD(h^*)]_{b,a}\neq 0\), which violates the constraint \(\mathbb{1}\{\bD(h^*)\} \odot \mathbb{1}\{\bD^{\top}(h^*)\} = \bI_{n \times n}\) and will conclude the proof by contradiction. Before giving the steps, we provide the following argument which we repeatedly use in the rest of the proof. For any continuous function \(f:\R^n \to \R\), we have \[\begin{align} \E\Big[\big|f(\bZ)\big|\Big] \neq 0 \;\; &\iff \;\; \E\Big[\big|f(\bZ) \cdot \big[\bss(\bZ) - \bss^{a}(\bZ)\big]_{I^a} \big|\Big] \neq 0 \;, \tag{129} \\ and \quad \E\Big[\big|f(\bZ)\big|\Big] \neq 0 \;\; &\iff \;\; \E\Big[\big|f(\bZ) \cdot \big[\bss(\bZ) - \tilde{\bss}^{\pi_b}(\bZ)\big]_{I^a} \big|\Big] \neq 0 \;. \tag{130} \end{align}\] First, suppose that \(\E\big[|f(\bZ)|\big] \neq 0\). Then, there exists an open set \(\Psi \subseteq \R^n\) for which \(f(\bz)\neq 0\) for all \(z\in \Psi\). Due to the interventional discrepancy condition, there exists an open set within \(\Psi\) for which \([\bss(\bZ) - \bss^{a}(\bZ)]_{I^a} \neq 0\). This implies that \[\E\Big[\big|f(\bZ) \cdot \big[\bss(\bZ) - \bss^{a}(\bZ)\big]_{I^a} \big|\Big] \neq 0\] For the other direction, suppose that \(\E\Big[\big|f(\bZ) \cdot \big[\bss(\bZ) - \bss^{a}(\bZ)\big]_{I^a} \big|\Big] \neq 0\), which implies that there exists an open set \(\Psi\) for which both \(f(\bz)\) and \([\bss(\bz) - \bss^{a}(\bz)]_{I^a}\) are nonzero. Then, \(\E\big[|f(\bZ)|\big] \neq 0\), and we have 129 . Similarly, due to \(\tilde{I}^{\pi_b} = I^a\) and interventional discrepancy, we obtain 130 .

1.0.3.3.0.1 Step 1: Show that \(\E\Big[\big|[J_f^{-1}(\bZ)]_{I^a,a}\big|\Big] \neq 0\).

First, using 30 and 1(i), we have \[\begin{align} \big[\bD(h^*)\big]_{a,a} &= \E\bigg[\Big|\big[J_f^{-\top}(\bZ)\big]_{a} \cdot \big[\bss(\bZ)-\bss^{a}(\bZ)\big]\Big|\bigg] \\ &= \E\bigg[\Big|\sum_{I^j \in \Paplus(I^a)} \big[J_f^{-1}(\bZ)\big]_{I^j,a} \cdot \big[\bss(\bZ)-\bss^{a}(\bZ)\big]_{I^j}\Big|\bigg] \;. \label{eq:no-solution-step1} \end{align}\tag{131}\] Note that \(\Paplus(I^a) \cap \mcA = \{I^a\}\) since \(I^a\) is non-descendant of the other nodes in \(\mcA\). Consider \(I^j\in\Pa(I^a)\), which implies that \(I^j \notin \mcA\) and \(I^j=\tilde{I}^{\pi_j}\). By 127 , we have \([J_f^{-1}(\bZ)]_{I^j,a} = 0\). Then, 131 becomes \[\big[\bD(h^*)\big]_{a,a} = \E\bigg[\Big|\big[J_f^{-1}(\bZ)\big]_{I^a,a} \cdot \big[\bss(\bZ)-\bss^{a}(\bZ)\big]_{I^a}\Big|\bigg] \neq 0 \;,\] since diagonal entries of \(\bD(h^*)\) are nonzero due to the last constraint in 57 . Then, 129 implies that \(\E\Big[\big|[J_f^{-1}(\bZ)]_{I^a,a}\big|\Big] \neq 0\).

1.0.3.3.0.2 Step 2: Show that \(\big[\tilde{\bD}(h^*)\big]_{a,b}\neq 0\).

Next, we use \(I^a=\tilde{I}^{\pi_b}\) and 1(i) to obtain \[\begin{align} \big[\tilde{\bD}(h^*)\big]_{a,b} &= \E\bigg[\Big|\big[J_f^{-\top}(\bZ)\big]_{a} \cdot \big[\bss(\bZ)-\tilde{\bss}^{\pi_b}(\bZ)\big]\Big|\bigg]\\ &= \E\bigg[\Big|\sum_{I^j \in \Paplus(I^a)} \big[J_f^{-1}(\bZ)\big]_{I^j,a} \cdot \big[\bss(\bZ) - \tilde{\bss}^{\pi_b}(\bZ)\big]_{I^j} \Big|\bigg] \\ &= \E\bigg[\Big|\big[J_f^{-1}(\bZ)\big]_{I^a,a} \cdot \big[\bss(\bZ)-\tilde{\bss}^{\pi_b}(\bZ)\big]_{I^a} \Big|\bigg] \;. \end{align}\] Using 130 and Step 1 result, we have \(\big[\tilde{\bD}(h^*)\big]_{a,b}\neq 0\).

1.0.3.3.0.3 Step 3: Show that \(\E\Big[\big|[J_f^{-1}(\bZ)]_{I^a,b}\big|\Big] \neq 0\).

Using 30 and 1(i), we have \[\begin{align} \big[\tilde{\bD}(h^*)\big]_{b,b} &= \E\bigg[\Big|\big[J_f^{-1}(\bZ)\big]_{b} \cdot \big[\bss(\bZ)-\tilde{\bss}^{\pi_b}(\bZ)\big]\Big|\bigg] \\ &= \E\bigg[\Big|\sum_{I^j \in \Paplus(I^a)} \big[J_f^{-1}(\bZ)\big]_{I^j,b} \cdot \big[\bss(\bZ)-\tilde{s}^{\pi_b}(\bZ)\big]_{I^j}\Big|\bigg] \\ &= \E\bigg[\Big|\big[J_f^{-1}(\bZ)]_{I^a,b} \cdot \big[\bss(\bZ)-\tilde{s}^{\rho_c}(\bZ)\big]_{I^a} \Big|\bigg] \;. \end{align}\] Since \(\mathbb{1}\{\bD(h^*)\}=\mathbb{1}\{\tilde{\bD}(h^*)\}\), the diagonal entry \(\big[\tilde{\bD}(h^*)\big]_{b,b}\) is nonzero. Then, using 130 we have \(\E\Big[\big|[J_f^{-1}(\bZ)]_{I^a,b}\big|\Big] \neq 0\).

1.0.3.3.0.4 Step 4: Show that \(\big[\bD(h^*)\big]_{b,a}\neq 0\).

Next, we use \(I^a=\tilde{I}^{\pi_b}\) and 1(i) to obtain \[\begin{align} \big[\bD(h^*)\big]_{b,a} &= \E\bigg[\Big|\big[J_f^{-\top}(\bZ)\big]_{b} \cdot \big[\bss(\bZ)-\bss^{a}(\bZ)\big]\Big|\bigg]\\ &= \E\bigg[\Big|\sum_{I^j \in \Paplus(I^a)} \big[J_f^{-1}(\bZ)\big]_{I^j,b} \cdot \big[\bss(\bZ) - \bss^{a}(\bZ)\big]_{I^j} \Big|\bigg] \\ &= \E\bigg[\Big|\big[J_f^{-1}(\bZ)\big]_{I^a,b} \cdot \big[\bss(\bZ)- \bss^{a}(\bZ)\big]_{I^a} \Big|\bigg] \;. \end{align}\] Using 129 and Step 3 result, we have \([\bD(h^*)]_{b,a}\neq 0\).

Finally, using the constraint \(\mathbb{1}\{\bD(h^*)\}=\mathbb{1}\{\tilde{\bD}(h^*)\}\), Step 2 implies that \([\bD(h^*)]_{a,b}\neq 0\). Combining with Step 4 result, we have \([\bD(h^*) \odot \bD^{\top}(h^*)]_{a,b} \neq 0\), which violates the last constraint in 57 . Therefore, if \((h^*,\pi)\) is a minimizer of 57 , then \(\pi\) must be the correct coupling \(\sigma\).

1.0.4 Analysis of the Assumptions↩︎

1.0.4.1 Analysis of Assumption 1

In this subsection, we prove 6 statement, i.e., 1 is satisfied for additive noise models under hard interventions, as follows. Consider a hard interventional environment \(\mcE^m\) and let \(I^m = i\) be the intervened node. Recall 23 , which becomes \[[\bsd_{\bZ}^{m}(\bz)]_k = \big[\bss(\bz) - \bss^{m}(\bz)\big]_k = \frac{\partial}{\partial z_k} \log p_{i}(z_i \mid \bz_{\Pa(i)}) \;,\] for the hard intervention on node \(i\) and parent node \(k \in \Pa(i)\). Then, the ratio in 1 is given by \[\label{eq:ratio-SN-2} \frac{[\nabla_{\bz} \log p(z_i \mid \bz_{\Pa(i)})]_k }{[\nabla_{\bz} \log p(z_i \mid \bz_{\Pa(i)})]_i - [\nabla_{\bz} \log q(z_i)]_i } = \frac{[\bsd_{\bZ}^{m}(\bz)]_k}{[\bsd_{\bZ}^{m}(\bz)]_i} \;.\tag{132}\] We will prove the desired result by contradiction. Assume that there exists a nonzero constant \(c \in \R\) such that \([\bsd_{\bZ}^{m}(\bz)]_k = c \cdot [\bsd_{\bZ}^{m}(\bz)]_i\) for all \(\bz \in \R^n\).

1.0.4.1.0.1 Additive noise model.

The additive noise model for node \(i\) is given by \(Z_i = f_i(\bZ_{\Pa(i)}) + N_i\) as specified in 7 . When node \(i\) is hard intervened, \(Z_i\) is generated according to \(Z_i = \bar N_i\) in which \(\bar N_i\) denotes the exogenous noise term for the interventional model. Then, denoting the pdfs of \(N_i\) and \(\bar N_i\) by \(p_N\) and \(q_N\), respectively, we have \[\begin{align} \label{eq:additive-pdf-connection} p_i(z_i \mid \bz_{\Pa(i)}) = p_N(z_i - f_i(\bz_{\Pa(i)}))\;, \quad and \quad q_i(z_i) = q_N(z_i)\;. \end{align}\tag{133}\] Denote the score functions associated with \(p_N\) and \(q_N\) by \[\begin{align} \label{eq:def-noise-scores} r_p(u) \triangleq \frac{{\rm d}}{{\rm d} u} \log p_N(u) = \frac{p_N'(u)}{p_N(u)}\;, \quad and \quad r_q(u) \triangleq \frac{{\rm d}}{{\rm d} u} \log q_N(u) = \frac{q_N'(u)}{q_N(u)}\;. \end{align}\tag{134}\] Define \(n_i\) and \(\bar n_i\) as the realizations of \(N_i\) and \(\bar N_i\) when \(Z_i=z_i\) and \(\bZ_{\Pa(i)}=\bz_{\Pa(i)}\). Then, we have \(\bar n_i = n_i + f_i(\bz_{\Pa(i)})\). Using 133 and 134 , we can express 132 as \[\require{physics} \begin{align} \label{eq:rank2-proof-additive-score-diff} \frac{[\bsd_{\bZ}^{m}(\bz)]_k}{[\bsd_{\bZ}^{m}(\bz)]_i} = \frac{-\pdv{f_i(\bz_{\Pa(i)})}{z_k} \cdot r_p(n_i) }{r_p(n_i) - r_q(n_i + f_i(\bz_{\Pa(i)}))} \;. \end{align}\tag{135}\] Then, \([\bsd_{\bZ}^{m}(\bz)]_k = c \cdot [\bsd_{\bZ}^{m}(\bz)]_i\) for all \(\bz \in \R^n\) becomes \[\label{eq:additive-rank2-step2} c \cdot \frac{\partial f_i(\bz_{\Pa(i)})}{\partial z_k} + 1 = \frac{r_q(n_i + f_i(\bz_{\Pa(i)}))}{r_p(n_i)} \;.\tag{136}\] We scrutinize 136 in two cases.

1.0.4.1.0.2 Case 1: \(\frac{\partial f_i(\bz_{\Pa(i)})}{\partial z_k}\) is constant.

In this case, let \(\eta = \frac{\partial f_i(\bz_{\Pa(i)})}{\partial z_k}\) for some \(\eta \in \R\) for all \(\bz_{\Pa(i)} \in \R^{|\Pa(i)|}\). Then, the RHS of 136 is also constant for all \(\bz_{\Pa(i)} \in \R^{|\Pa(i)|}\) and \(n_i \in \R\). Fix a realization \(n_i = n_i^*\) and note that \(f_i(\bz_{\Pa(i)})\) needs to be constant, denoted by \(\delta^*\), for all \(\bz_{\Pa(i)} \in \R^{|\Pa(i)|}\). However, this implies that \(\frac{\partial f_i(\bz_{\Pa(i)})}{\partial z_k} = 0\), and we have \(r_q(n_i + \delta^*) = r_p(n_i)\) for all \(n_i \in \R\). This implies that \(p_N(n_i) = \eta^* q_N(n_i+\delta^*)\) for some constant \(\eta^*\). Since \(p_N\) and \(q_N\) are pdfs, the only choice is \(\eta^*=1\) and \(p_N(n_i)=q_N(n_i+\delta^*)\), then \(p_i(z_i \mid \bz_{\Pa(i)}) = q_i(z_i)\), which contradicts the premise that an intervention changes the causal mechanism of the target node \(i\).

1.0.4.1.0.3 Case 2: \(\frac{\partial f_i(\bz_{\Pa(i)})}{\partial z_k}\) is not constant.

In this case, note that LHS of 136 is not a function of \(n_i\). Then, taking the derivative of both sides with respect to \(n_i\) and rearranging, we obtain \[\label{eq:additive-rank2-ratio} \frac{r_p'(n_i)}{r_p(n_i)} = \frac{r_q'(n_i + f_i(\bz_{\Pa(i)}))}{r_q(n_i + f_i(\bz_{\Pa(i)}))} \;, \quad \forall (n_i,\bz_{\Pa(i)})\in \R \times \R^{|\Pa(i)|} \;.\tag{137}\] Next, consider a fixed realization \(n_i=n_i^*\), and denote the value of LHS by \(\alpha\). Since \(f_i(\bz_{\Pa(i)})\) is continuous and not constant, its image contains an open interval \(\Theta \subseteq \mathbb{R}\). Denoting \(u\triangleq f_i(\bz_{\Pa(i)})\), we have \[\alpha = \frac{r_q'(n_i^* + u)}{r_q(n_i^* + u)} \;, \quad \forall u \in \Theta \;.\] The only solution to this equality is that \(r_q(n_i^* + u)\) is an exponential function, \(r_q(u) = k_1\exp(\alpha u)\) over interval \(u \in \Theta\). Since \(r_q\) is an analytic function that equals to an exponential function over an interval, it is exponential over entire \(\mathbb{R}\). This implies that the pdf \(q_N(u)\) is of the form \(q_N(u) = k_2 \exp((k_1/\alpha) \exp(\alpha u))\). However, this cannot be a valid pdf since its integral over \(\R\) diverges. Then, \([\bsd_{\bZ}^{m}(\bz)]_k = c \cdot [\bsd_{\bZ}^{m}(\bz)]_i\) cannot be true, which concludes the proof.

1.0.4.2 Analysis of Assumption 2

In this section, we establish the necessary and sufficient conditions under which 2 holds for additive models. Furthermore, we show that a large class of nonlinear models in the latent space satisfy these conditions, including the two-layer neural networks. We have focused on such NNs since they effectively approximate continuous functions [60]. Readily, the necessary and sufficient conditions can be investigated for other choices of nonlinear functions.

1.0.4.2.1 Interpreting Assumption 2

The implication of 2 is that the effect of an intervention is not lost in any linear combination of the varying coordinates of the scores. To formalize this, for each node \(i\in[n]\), we define \[\begin{align} \mcC_i \triangleq \{\bc\in\R^n\;:\; \exists j \in \Paplus(i) \;\; such that \; c_j\neq 0 \}\;, \label{eq:C} \end{align}\tag{138}\] Let \(I^m = i\). Note that \(\bR_{\bZ}^{m}\) is positive semi-definite and \([\bR_{\bZ}^{m}]_{j,k} = 0\) for all \((j,k) \notin \Paplus(i) \times \Paplus(i)\). Therefore, for any \(\bc \in \mcC_{i}\), \[\bc^{\top} \cdot \bR_{\bZ}^{m} \cdot \bc = \bE\big[ (\bc^{\top} \cdot \bsd_{\bZ}^m) \cdot (\bc^{\top} \cdot \bsd_{\bZ}^m)^{\top} \big] \neq 0 \;,\] since \(\bR_{\bZ}^{m}\) has rank \(|\Paplus(i)|\). The reverse direction also holds true, i.e., if \(\bR_{\bZ}^{m}\) has rank less than \(|\Paplus(i)|\), then there exists \(\bc \mcC_{i}\) which makes \(\bc^{\top} \cdot \bR_{\bZ}^{m} \cdot \bc = 0\). Therefore, 2 is equivalent to the following statement: \[\label{eq:assumption-full-rank-original} \bE\big[|\bc^{\top} \cdot \bsd_{\bZ}^{m}|\big] \neq 0 \;, \quad \forall \bc \in \mcC_{i} \;.\tag{139}\]

1.0.4.2.1.1 Necessary and sufficient conditions.

Consider the additive noise model for node \(i\) \[\begin{align} Z_i &= f_i(\bZ_{\Pa(i)}) + N_i\;, \label{eq:additive-model-discussion-proof-obs} \end{align}\tag{140}\] as specified in 7 . When node \(i\) is soft intervened, \(Z_i\) is generated according to \[\begin{align} Z_i &= \bar f_i(\bZ_{\Pa(i)}) + \bar N_i\;, \label{eq:additive-model-discussion-proof-int} \end{align}\tag{141}\] in which \(\bar f_i\) and \(\bar N_i\) specify the interventional mechanism for node \(i\). The following lemma characterizes the necessary and sufficient conditions under which 139 (equivalently, 2) is satisfied. In this subsection, we use \(\varphi\) as the shorthand for \(\bz_{\Pa(i)}\).

See 1.0.4.2.2.
To provide some intuition about the conditions in 14, we consider a node \(i \in [n]\) and discuss the conditions in the context of a few examples. Note that by sweeping \(\varphi\in \R^{|\Pa(i)|}\) we generate a continuum of linear equations of the form: \[\begin{align} \label{eq:NS95equations2} c_i - \bc^{\top} \cdot \nabla_{\bz} f_i(\varphi) = 0\;, \qquad and \qquad c_i - \bc^{\top} \cdot \nabla_{\bz} \bar f_i(\varphi) = 0 \;. \end{align}\tag{142}\] Note that for all \(j\not\in \Pa(i)\) we have \([\nabla_{\bz} f_i(\varphi)]_j=[\nabla_{\bz} \bar f_i(\varphi)]_j=0\). Hence, in finding the solutions to 142 only the coordinates \(\{j\in \Paplus(i)\}\) of \(\bc\) are relevant. Let us define \[\begin{align} \label{eq:varphi} u(\varphi)\triangleq \nabla_{\varphi} f_i(\varphi) \;, \quad and \quad \bar u(\varphi)\triangleq \nabla_{\varphi} \bar f_i(\varphi) \;, \end{align}\tag{143}\] which are the gradients of \(f_i\) and \(f_i\) by considering only the coordinates of \(\bz\) in \(\{j\in {\Pa}(i)\}\). Accordingly, we also define \(\bb\) by concatenating only the coordinates of \(\bc\) with their indices in \(\{j\in {\Pa}(i)\}\). Next, consider \(w\) distinct choices of \(\varphi\) and denote them by \(\{\varphi^t\in\R^{|\Pa(i)|}:t\in[w]\}\). By concatenating the two equations in 143 specialized to these realizations, we get the following linear system with \(2w\) equations and \(|\Pa(i)|+1\) unknown variables. \[\begin{align} \label{eq:discussion-full-rank} \underset{\triangleq \, \bV \, \in \, \R^{2w\times (|\Pa(i)|+1)}}{\underbrace{ \begin{bmatrix} [u(\varphi^1)]^{\top} & -1\\ [\bar u(\varphi^1)]^{\top} & -1\\ \vdots & \vdots\\ [u(\varphi^w)]^{\top} & -1\\ [\bar u(\varphi^w)]^{\top} & -1\\ \end{bmatrix}}} \begin{bmatrix} \bb \\ c_i \end{bmatrix} =\boldsymbol{0}_{2w}\;. \end{align}\tag{144}\] When \(\bV\) is full-rank, i.e., \(\rank(\bV)=|\Pa(i)|+1\), the system has only the trivial solutions \(c_i=0\) and \(\bb=\boldsymbol{0}\). Then, we make the following observations.

1. If \(f_i\) and \(\bar f_i\) are linear functions, the vector spaces generated by \(u\) and \(\bar u\) have dimensions 1. Subsequently, we always have \(\rank(\bV)\leq 2\), rendering an underdetermined system when \(|\Pa(i)|\geq 2\). Hence, when the maximum degree of \(\mcG\) is at least 2, a linear causal model does not satisfy 2.

2. If \(f_i\) and \(\bar f_i\) are quadratic with full-rank matrices, i.e., \(f_i(\varphi) = \varphi^{\top} \bA \varphi\) and \(\bar f_i(\varphi) = \varphi^{\top} \bar \bA \varphi\) where \(\rank(\bA)=\rank(\bar \bA)=|\Pa(i)|\), there is a choice of \(w \geq |\Pa(i)| + 1\) and realizations \(\{\varphi^t \in \R^{|\Pa(i)|} : t \in [w] \}\) for which \(\rank(\bV)=|\Pa(i)|+1\) and the system in 142 admits only the trivial solutions \(a_i=0\) and \(\bb=\boldsymbol{0}\). Hence, quadratic causal models satisfy 2.

3. If \(f_i\) and \(\bar f_i\) are two-layer NNs with a sufficiently large number of hidden neurons, they also render a fully determined system, and as a result, they satisfy 2.

We investigate the last example in detail as follows. Assume that \(f_i\) and \(\bar f_i\) are two-layer NNs with \(|\Pa(i)|\) inputs, \(w_i\) and \(\bar w_i\) hidden nodes, respectively, and with sigmoid activation functions. Denote the weight matrices between input and hidden layers in \(f_i\) and \(\bar f_i\) by \(\bW^i \in \R^{w_i \times |\Pa(i)|}\) and \(\bar \bW^i \in \R^{\bar w_i \times |\Pa(i)|}\), respectively. Furthermore, define \(\bnu \in \R^{w_i}\) and \(\bar \bnu \in \R^{w_i}\) as the weights between the hidden layer and output in \(f_i\) and \(\bar f_i\), respectively. Finally, define \(\nu_0\) and \(\bar \nu_0\) as the bias terms. Hence, we have \[\begin{align} f_i(\varphi) & = \bnu^{\top} \cdot \sigma(\bW^i \cdot \varphi) + \nu_0 = \sum_{j=1}^{w_i} \nu_j \cdot \sigma(\bW^i_j \cdot \varphi) + \nu_0 \;, \tag{145} \\ and \quad \bar f_i(\varphi) & = \bar \bnu^{\top} \cdot \sigma(\bar \bW^i \cdot \varphi) + \nu_0 = \sum_{j=1}^{\bar w_i} \bar \nu_j \cdot \sigma(\bar \bW^i_j \cdot \varphi) + \bar \nu_0 \;, \tag{146} \end{align}\] in which activation function \(\sigma\) is applied element-wise.

1.0.4.2.1.2 7

Consider NNs \(f_i\) and \(\bar f_i\) specified in 145 and 146 . If for all \(i \in[n]\) we have \[\begin{align} \max\{\rank(\bW^i)\; , \; \rank(\bar \bW^i)\} = |\Pa(i)|\;, \end{align}\] then 2 holds.
See 1.0.4.2.3.

1.0.4.2.2 Proof of Lemma 14

We show that for node \(I^m = i\), the condition in 139 \[\begin{align} \E\Big[\big| \bc^{\top} \cdot \bsd_{\bZ}^{m} \big|\Big] \neq 0 \;, \;\; \forall c \in \mcC_i\;, \label{eq:assumption95sc:node95level} \end{align}\tag{147}\] holds if and only if the following continuum of equations admit their solutions \(\bc\) in \(\R^n \setminus \mcC_i\). \[\begin{align} \left\{ \begin{array}{l} c_i - \bc^{\top} \cdot \nabla_{\bz} f_i(\varphi) = 0 \\ \\ c_i - \bc^{\top} \cdot \nabla_{\bz} \bar f_i(\varphi) = 0 \end{array} \right.\;, \qquad \qquad \forall \varphi \in\R^{|\Pa(i)|}\;, \label{eq:assumption95sc:condition95node95level} \end{align}\tag{148}\] in which shorthand \(\varphi\) is used for \(\bz_{\Pa(i)}\). We first note that the condition in 147 can be equivalently stated using 16 and 17 as \[\begin{align} \label{eq:equivalent95statement2-pre} \E\bigg[\Big| \bc^{\top}\cdot \big[\nabla_{\bz} \log p_i(Z_i \mid \bZ_{\Pa(i)}) - \nabla_{\bz} \log q_i(Z_i \mid \bZ_{\Pa(i)}) \big]\Big|\bigg] \neq 0 \;, \quad \forall \bc \in \mcC_{i} \;. \end{align}\tag{149}\] Also, using [prop:continuity-argument], 149 is equivalent to \[\begin{align} \label{eq:equivalent95statement2} \forall \bc \in \mcC_i, \;\; \exists \bz : \quad \bc^{\top}\cdot \nabla_{\bz} \log p_i(z_i \mid \bz_{\Pa(i)}) \neq \bc^{\top}\cdot \nabla_{\bz} \log q_i(z_i \mid \bz_{\Pa(i)}) \;. \end{align}\tag{150}\] The additive noise model for node \(i\) is given by \[\begin{align} Z_i &= f_i(\bZ_{\Pa(i)}) + N_i\;, \label{eq:additive-model-aux-proof-obs} \end{align}\tag{151}\] as specified in 7 . When node \(i\) is soft intervened, \(Z_i\) is generated according to \[\begin{align} Z_i &= \bar f_i(\bZ_{\Pa(i)}) + \bar N_i\;, \label{eq:additive-model-aux-proof-int} \end{align}\tag{152}\] in which \(\bar f_i\) and \(\bar N_i\) specify the interventional mechanism for node \(i\). Then, denoting the pdfs of \(N_i\) and \(\bar N_i\) by \(p_N\) and \(q_N\), respectively, 151 and 152 imply that \[\begin{align} \label{eq:p-h-connection} p_i(z_i \mid \bz_{\Pa(i)}) = p_N\big(z_i - f_i(\bz_{\Pa(i)})\big)\;, \quad and \quad q_i(z_i \mid \bz_{\Pa(i)}) = q_N\big(z_i - \bar f_i(\bz_{\Pa(i)})\big)\;. \end{align}\tag{153}\] Denote the score functions associated with \(p_N\) and \(q_N\) by \[\begin{align} \label{eq:def-r} r_p(u) \triangleq \frac{{\rm d}}{{\rm d} u} \log p_N(u) = \frac{p_N'(u)}{p_N(u)}\;, \quad and \quad r_q(u) \triangleq \frac{{\rm d}}{{\rm d} u} \log q_N(u) = \frac{q_N'(u)}{q_N(u)}\;. \end{align}\tag{154}\] Define \(n_i\) and \(\bar n_i\) as the realizations of \(N_i\) and \(\bar N_i\) when \(Z_i=z_i\) and \(\bZ_{\Pa(i)}=\bz_{\Pa(i)}\). By defining \(\delta(\bz_{\Pa(i)}) \triangleq f_i(\bz_{\Pa(i)}) - \bar f_i(\bz_{\Pa(i)})\), we have \(\bar n_i = n_i + \delta(\bz_{\Pa(i)})\). Using 154 and 153 , we can express the relevant entries of \(\nabla_{\bz} \log p_i(z_i \mid \bz_{\Pa(i)})\) and \(\nabla_{\bz} \log q_i(z_i \mid \bz_{\Pa(i)})\) as \[\require{physics} \begin{align} \big[\nabla_{\bz} \log p(z_i \mid \bz_{\Pa(i)})\big]_j &= \begin{cases} r_p(n_i) \;, &j = i \;, \\ -\pdv{f_i(\bz_{\Pa(i)})}{z_j} r_p(n_i) \;, &j \in \Pa(i) \;, \\ 0 \;, &j \notin \Paplus(i) \;, \end{cases} \tag{155} \\ and \quad \big[\nabla_{\bz} \log q_i(z_i \mid \bz_{\Pa(i)})\big]_j &= \begin{cases} r_q(\bar n_i) \;, &j = i \;, \\ -\pdv{\bar f_i(\bz_{\Pa(i)})}{z_j} r_q(\bar n_i) \;, &j \in \Pa(i) \; , \\ 0 \;, &j \notin \Paplus(i) \;. \end{cases} \tag{156} \end{align}\] By substituting 155 –156 in 150 and rearranging the terms, the statement in 150 becomes equivalent to following statement. For all \(\bc \in \mcC_i\), there exist \(n_i \in \R\) and \(\bz_{\Pa(i)}\in\R^{|\Pa(i)|}\) such that \[\require{physics} \begin{align} \label{eq:equivalent95statement3} r_p(n_i) \cdot \bigg(c_i- \sum_{j \in {\Pa}(i)} c_j \cdot \pdv{f_i(\bz_{\Pa(i)})}{z_j}\bigg) \neq r_q(n_i+\delta(\bz_{\Pa(i)})) \cdot \bigg(c_i- \sum_{j \in {\Pa}(i)} c_j\cdot \pdv{\bar f_i(\bz_{\Pa(i)})}{z_j}\bigg)\;, \end{align}\tag{157}\] which by using the shorthand \(\varphi\) for \(\bz_{\Pa(i)}\) can be compactly presented as follows. For all \(\bc\in \mcC_i\), there exists \(n_i \in \R\) and \(\varphi \in \R^{|\Pa(i)|}\) such that \[\begin{align} \label{eq:equivalent95statement4} r_p(n_i) \cdot \left[c_i - \bc^{\top} \cdot \nabla_{\bz} f_i(\varphi)\right] \neq r_q\big(n_i+\delta(\varphi)\big) \cdot \left[c_i - \bc^{\top} \cdot \nabla_{\bz} \bar f_i(\varphi)\right] \;. \end{align}\tag{158}\] Hence, 2 is equivalent to the statement in 158 , which we use for the rest of the proof.

1.0.4.2.2.1 Sufficient condition.

We show that if 148 admit solutions \(\bc\) only in \(\R^n\setminus \mcC_i\), then the statement in 158 holds. By contradiction, assume that there exists \(\bc^*\in\mcC_i\) such that for all \(n_i \in\R\) and \(\varphi \in \R^{|\Pa(i)|}\), \[\begin{align} \label{eq:equivalent95statement5} r_p(n_i) \cdot \left[c^*_i-(\bc^*)^{\top} \cdot \nabla_{\bz} f_i(\varphi)\right] = r_q\big(n_i+\delta(\varphi)\big) \cdot \left[c^*_i-(\bc^*)^{\top} \cdot \nabla_{\bz} \bar f_i(\varphi)\right] \; . \end{align}\tag{159}\] We show that \(\bc^*\in\mcC_i\) is also a solution to 148 , contradicting the premise. In order to show that \(\bc^*\in\mcC_i\) is also a solution to 148 , suppose, by contradiction, that 159 holds, and there exists \(\varphi^*\in\R^{|\Pa(i)|}\) corresponding to which \[\begin{align} \label{eq:contradiction95premise} (\bc^*)^{\top}\cdot \nabla_{\bz} \bar f_i(\varphi^*)- c^*_i \neq 0\;. \end{align}\tag{160}\] Note that \(f_i\) is a continuously differentiable function and also a function of \(z_j\) for all \(j \in \Pa(i)\). Hence, there exists an open set \(\Phi \subseteq \mathbb{R}^{|\Pa(i)|}\) for \(\varphi\) for which \((\bc^*)^{\top}\cdot \nabla_{\bz} f_i(\varphi^*)- c^*_i\) is nonzero everywhere in \(\Phi\). Likewise, \(r_p\) cannot constantly be zero over all possible intervals \(\Omega\). This is because otherwise, it would have to necessarily be a constant zero function (since it is analytic), which is an invalid score function. Hence, there exists an open interval \(\Omega \subseteq \mathbb{R}\) over which \(r_p(n_i)\) is nonzero for all \(n_i \in\Omega\). Subsequently, the left-hand side of 159 is nonzero over the Cartesian product \((n_i, \varphi) \in \Omega \times \Phi\). This means that if 159 is true, then both functions on its right-hand side must also be nonzero over \(\Omega \times \Phi\). Hence, by rearranging the terms in 159 we have \[\begin{align} \frac{r_q\big(n_i+\delta(\varphi)\big)}{r_p(n_i)} = \frac{c^*_i - (\bc^*)^{\top}\cdot \nabla_{\bz} f_i(\varphi)}{ c^*_i - (\bc^*)^{\top}\cdot \nabla_{\bz} \bar f_i(\varphi)}\;, \qquad \forall (n_i, \varphi) \in \Omega \times \Phi \;. \label{eq:ratio} \end{align}\tag{161}\] In two steps, we show that 161 cannot be valid.

Step 1. First, we show that function \(\delta\) cannot be a constant over \(\Phi\) (interval specified above). Suppose the contrary and assume that \(\delta(\varphi)= \delta^*\) for all \(\varphi \in \Phi\). Hence, the gradient of \(\delta\) is zero. Using the definition of \(\delta\), this implies that \[\begin{align} \nabla_{\bz} f_i(\varphi) = \nabla_{\bz} \bar f_i(\varphi) \;, \quad \forall \varphi \in \Phi \;. \end{align}\] Then, by leveraging 161 , we conclude that \(r_p(n_i) = r_q(n_i + \delta^*)\) for all \(n_i \in \Omega\). Since \(r_p(n_i)\) and \(r_q(n_i + \delta^*)\) are analytic functions that agree on an open interval of \(\mathbb{R}\), they are equal for all \(n_i \in \mathbb{R}\) as well. This implies that \(p_N(n_i) = \eta \cdot q_N(n_i+\delta^*)\) for some constant \(\eta\). Since \(p_N\) and \(q_N\) are pdfs, the only choice is \(\eta=1\) and \(p_N(n_i)=q_N(n_i+\delta^*)\), then \(p_i(z_i \mid \bz_{\Pa(i)}) = q_i(z_i \mid \bz_{\Pa(i)})\), which contradicts the premise that observational and interventional mechanisms are distinct.

Step 2. Finally, we will show that 161 cannot be true when \(\delta\) is not a constant function over \(\Phi\). Note that the right-hand side of 161 is not a function of \(n_i\). Then, taking the derivative of both sides with respect to \(n_i\) and rearranging, we obtain \[\begin{align} \label{eq:ratio-derivatives} \frac{r_p'(n_i)}{r_p(n_i)} &= \frac{r_q'(n_i + \delta(\varphi))}{r_q(n_i + \delta(\varphi))} \;, \quad \forall (n_i,\varphi)\in \Omega \times \Phi \;. \end{align}\tag{162}\] Next, consider a fixed realization \(n_i=n_i^*\), and denote the value of LHS by \(\alpha\). Since \(\delta\) is continuous and not constant over \(\Phi\), its image contains an open interval \(\Theta \subseteq \mathbb{R}\). Denoting \(u\triangleq \delta(\varphi)\), we get \[\begin{align} \alpha = \frac{r_q'(n_i^* + u)}{r_q(n_i^* + u)} \;, \quad \forall u \in \Theta \;. \end{align}\] The only solution to this equality is that \(r_q(n_i^* + u)\) is an exponential function, \(r_q(u) = k_1\exp(\alpha u)\) over interval \(u \in \Theta\). Since \(r_q\) is an analytic function that equals to an exponential function over an interval, it is exponential over entire \(\mathbb{R}\). This implies that the pdf \(q_N(u)\) is of the form \(q_N(u) = k_2 \exp((k_1/\alpha) \exp(\alpha u))\). However, this cannot be a valid pdf since its integral over \(\R\) diverges. Hence, 161 is not true, and the premise that there exists \(\bc^* \in \mcC_i\) and \(\varphi\in\R^{|\Pa(i)|}\) corresponding to which 160 holds is invalid, concluding that for all \(\varphi\in\R^{|\Pa(i)|}\) we have \(\bc^{\top}\cdot \nabla_{\bz} \bar f_i(\varphi) = c_i\). Proving the counterpart identity \(\bc^{\top}\cdot \nabla_{\bz} f_i(\varphi) = c_i\) follows similarly. Therefore, 159 implies \[\begin{align} \left\{ \begin{array}{l} c^*_i - (\bc^*)^{\top} \cdot \nabla_{\bz} f_i(\varphi) = 0 \\ \\ c^*_i - (\bc^*)^{\top} \cdot \nabla_{\bz} \bar f_i(\varphi) = 0 \end{array} \right.\;, \qquad \qquad \forall \varphi \in\R^{|\Pa(i)|}\;, \end{align}\] which means that we have found a solution to 148 that is not in \(\R^n\setminus \mcC_i\), contradicting the premise.

1.0.4.2.2.2 Necessary condition.

Assume that 2, and equivalently, the statement in 158 holds but 148 has a solution \(\bc^*\) in \(\mcC_i\). Then, \[\begin{align} c^*_i - (\bc^*)^{\top} \cdot \nabla_{\bz} f_i(\varphi) = c^*_i - (\bc^*)^{\top} \cdot \nabla_{\bz} \bar f_i(\varphi)= 0 \;, \quad \forall \varphi \in \R^{|\Pa(i)|} \;. \end{align}\] Multiplying the left side by \(r_p(n_i)\) and the right side by \(r_q(n_i + \delta(\varphi))\), we obtain that for all \(n_i \in \R\) and \(\varphi\in\R^{|\Pa(i)|}\) \[\begin{align} r_p(n_i) \cdot \left[c^*_i-(\bc^*)^{\top} \cdot \nabla_{\bz} f_i(\varphi)\right] = r_q(n_i+\delta(\varphi)) \cdot \left[c^*_i-(\bc^*)^{\top} \cdot \nabla_{\bz} \bar f_i(\varphi)\right] = 0 \;, \end{align}\] which implies that for all \(\bc \in \mcC_i\), \(n_i \in \R\) and \(\varphi\in\R^{|\Pa(i)|}\) \[\begin{align} r_p(n_i) \cdot \left[c_i-\bc^{\top} \cdot \nabla_{\bz} f_i(\varphi)\right] = r_q(n_i+\delta(\varphi)) \cdot \left[c_i-\bc^{\top} \cdot \nabla_{\bz} \bar f_i(\varphi)\right] \;. \end{align}\] This contradicts 158 , and equivalently 2. Hence, the proof is complete.

1.0.4.2.3 Proof of Lemma 7

1.0.4.2.3.1 Approach.

We will use the same argument as at the beginning of the proof of 14. Specifically, we will show that for any node \(i\) and two-layer NNs \(f_i\) and \(\bar f_i\) with weight matrices \(\bW^i\) and \(\bar \bW^i\), the following continuum of equations admit their solutions \(\bc\) in \(\R^n\setminus \mcC_i\). \[\begin{align} \left\{ \begin{array}{l} c_i -\bc^{\top} \cdot \nabla_{\bz} \bar f_i(\varphi) = 0 \\ \\ c_i -\bc^{\top} \cdot \nabla_{\bz} f_i(\varphi) = 0 \end{array} \right.\;, \qquad \qquad \forall \varphi \in\R^{|\Pa(i)|}\;, \end{align}\] in which shorthand \(\varphi\) is used for \(\bz_{\Pa(i)}\). Hence, by invoking 14, 2 holds.

1.0.4.2.3.2 Definitions.

Define \(d_i \triangleq |\Pa(i)|\). Since \(\max\{\rank(\bW^i), \rank(\bar \bW^i)\}=d_i\), without loss of generality, suppose that \(\bW^i \in \R^{w_i \times d_i}\) has rank \(d_i\). The rest of the proof follows similarly for the case of \(\rank(\bar \bW^i)=d_i\). We use shorthand \(\{\bW, w\}\) for \(\{\bW^i, w_i\}\) when it is obvious from context.

1.0.4.2.3.3 Parameterization.

Note that \(f_i\) can be represented by different parameterizations, some containing more hidden nodes than others. Without loss of generality, let \(w\) be the fewest number of nodes that can represent \(f_i\). This implies that the entries of \(\bnu\) are nonzero. Otherwise, if \(\nu_i=0\), we can remove \(i\)-th hidden node and still have the same \(f_i\). Similarly, the rows of \(\bW\) are distinct. Otherwise, if there exist rows \(\bW_i = \bW_j\) for distinct \(i,j \in [w]\), removing \(j\)-th hidden node and using \((\nu_i + \nu_j)\) in place of \(\nu_i\) results the same function as \(f_i\) with \((w-1)\) hidden nodes. Similarly, we have \(\bW_i \neq \boldsymbol{0}\) for all \(i \in [w]\). Otherwise, we have \[\nu_i \cdot \sigma(\bW_i \cdot \varphi) + \nu_0 = (\nu_0 + \frac{\nu_i}{2}) \;,\] and by removing \(i\)-th hidden node and using \((\nu_0 + \frac{\nu_i}{2})\) instead of \(\nu_0\), we reach the same function as \(f_i\) with \((w-1)\) hidden nodes. Finally, we have \(\bW_i + \bW_j \neq \boldsymbol{0}\). Otherwise, we have \[\nu_i \cdot \sigma(\bW_i \varphi) + \nu_j \sigma(\bW_j \cdot \varphi) + \nu_0 = (\nu_i - \nu_j) \cdot \sigma(\bW_i \cdot \varphi) + (\nu_0+\nu_j) \;,\] and by removing \(j\)-th hidden node and using \((\nu_i-\nu_j)\) instead of \(\nu_i\) and \((\nu_0+\nu_j)\) instead of \(\nu_0\), we reach the same function as \(f_i\) with \((w-1)\) hidden nodes. In summary, we have \[\begin{align} \bW_i \neq \boldsymbol{0} \;, \;\; \nu_i \neq 0 \;, \quad \forall i \in [w] \;, \quad and \quad \bW_i \pm \bW_j \neq \boldsymbol{0} \;, \quad \forall i, j \in [w] : i \neq j \;. \label{eq:W95construction} \end{align}\tag{163}\] We will show that there does not exist \(\bc \in \mcC_i\) such that \(\bc^{\top} \cdot \nabla_{\bz} f_i(\varphi) = c_i\) for all \(\varphi \in \R^{d_i}\). Assume the contrary, and assume there exist \(\bc^* \in \mcC_i\) such that \[\begin{align} (\bc^*)^{\top} \cdot \nabla_{\bz} f_i(\varphi) = c^*_i \;, \quad \forall \varphi \in \R^{d_i} \;. \label{eq:nn-contrary} \end{align}\tag{164}\] This is equivalent to showing that there exists nonzero \(\bb^* \in \R^{d_i}\) such that \[\begin{align} (\bb^*)^{\top} \cdot \nabla_{\varphi} f_i(\phi) = c_i^{*} \;, \quad \forall \varphi \in \R^{d_i} \;. \end{align}\] Based on 145 , the gradient of \(f_i(\varphi)\) is \[\begin{align} \nabla_{\varphi} f_i(\varphi) = \bW^{\top} \cdot \diag(\bnu) \cdot \sdot(W \cdot \varphi) \;, \label{eq:u-expression} \end{align}\tag{165}\] where \(\diag(\bnu)\) is the diagonal matrix with \(\bnu\) as its diagonal elements, and \(\sdot\) is the derivative of the sigmoid function, applied element-wise to its argument. Hence, based on 164 , the contradiction premise is equivalent to having a nonzero \(\bb^* \in \R^{d_i}\) such that \[\begin{align} \big[\sdot(\bW \cdot \varphi)\big]^{\top} \cdot \diag(\bnu) \cdot \bW \cdot \bb^* = c_i \;, \quad \forall \varphi \in \R^{d_i} \;. \label{eq:nn-constant-requirement} \end{align}\tag{166}\] We note that since \(\bW\) is full-rank and \(\nu_i\) is nonzero for all \(i \in [w]\), the matrix \(\diag(\bnu) \cdot \bW\) is full-rank as well and it has a trivial null space. Subsequently, \(\bb^* \in \R^{d_i}\) is nonzero if and only if \((\diag(\bnu) \cdot \bW \cdot \bb^*) \in \R^{w}\) is nonzero. We will use the following lemma to show that 166 cannot be true. This establishes that the contradiction premise is not correct, and completes the proof.

See 1.0.4.2.4.

Let us define \(\xi \triangleq \bW \cdot \varphi\). We will show that there exists \(\varphi^* \in \R^{d_i}\) such that \(\xi = \bW \cdot \varphi^*\) satisfies the conditions in 15. Then, using 15 with the choice of \(\bu=\bW \cdot \varphi^*\), \(d=\diag(\bnu) \cdot \bW \cdot \bb^*\), and \(a=c_i\), we find that there exists \(\alpha \in \R\) such that \[\begin{align} \label{eq:nn-constant-requirement-contrary} \big[\sdot(\alpha \cdot \bW \cdot \varphi^*)\big]^{\top} \cdot \diag(\bnu) \cdot \bW \cdot \bb^* \neq c_i \;. \end{align}\tag{167}\] Hence, 166 is false since it is violated for \(\varphi=\alpha \varphi^*\) and the proof is completed. We show the existence of such \(\varphi^*\) as follows. We first construct the set of \(\varphi\) values for which conditions of 15 on \(w\) are not satisfied. The set in question is the union of the following cases: (i) \(\xi_i = \bW_i \cdot \varphi = 0\) for some \(i \in [w]\), (ii) \(|\xi_i| = |\xi_j|\) for some distinct \(i,j \in [w]\), or equivalently, \((\bW_i \pm \bW_j) \varphi = 0\). Note that \(\bW_i \neq \boldsymbol{0}\) and \(\bW_i \pm \bW_j \neq \boldsymbol{0}\) by 163 . For a nonzero vector \(\by \in \R^{d_i}\), the set \(\{\varphi \in \R^{d_i} \;: \;\by^{\top} \varphi = 0\}\) is a \((d_i-1)\)-dimensional subspace of \(\R^{d_i}\). Then, there is \(w\) number of \((d_i-1)\)-dimensional subspaces that fall under case (i), and \(w(w-1)\) number of \((d_i-1)\)-dimensional subspaces that fall under case (ii). Therefore, there are \(w^2\) lower-dimensional subspaces for which the conditions of 15 do not hold. However, \(\R^{d_i}\) cannot be covered by a finite number of lower-dimensional subspaces of itself. Therefore, there exists \(\varphi^* \in \R^{d_i}\) such that \(\xi = \bW \cdot \varphi^*\) satisfies the conditions of 15, and the proof is completed.

1.0.4.2.4 Proof of Lemma 15

Assume the contrary and suppose that there exists a nonzero \(\bc\) and \(\alpha\) for a given \(\bu\). Define the function \(g_u(\alpha) \triangleq\bc^{\top} \cdot \sdot(\alpha \bu)\). Note that \[\begin{align} \sdot(x) = \frac{1}{1+e^{-x}+e^{x}} \end{align}\] is an even analytic function, and its Taylor series expansion at \(0\) has the domain of convergence \(\{x \in \R : |x| < \frac{\pi}{2}\}\). Thus, for all \(\alpha \in (-\frac{\pi}{2 \max_{i} |u_{i}|}, \frac{\pi}{2 \max_{i} |u_{i}|})\), \[\begin{align} a = g_u(\alpha) &= \sum_{j=1}^{p} c_j \cdot \sdot(\alpha u_j) = \sum_{i=0}^{\infty} \gamma_i \sum_{j=1}^{p} c_j (\alpha u_j)^{2i} = \sum_{i=0}^{\infty} (\gamma_i \sum_{j=1}^{p} c_j u_j^{2i}) \alpha^{2i} \;. \end{align}\] Note that \(g_u(\alpha)\) is a constant function of \(\alpha \in (-\frac{\pi}{2 \max_{i} |u_{i}|}, \frac{\pi}{2 \max_{i} |u_{i}|})\). Thus, its Taylor coefficients, i.e., \((\gamma_i \sum_{j=1}^{p} c_j u_j^{2i})\), are zero for all \(i \in \N^+\). However, the coefficients of even powers in Taylor expansion of \(\sdot\), i.e., \(\gamma_i\)’s, are nonzero [61]. Therefore, we have \(\sum_{i=1}^{p}(c_j u_j^{2i}) = 0\) for all \(i \in \N^+\). Next, construct the following system of linear equations, \[\begin{align} \begin{bmatrix} u_1^2 & u_1^4 & \dots & u_1^{2w} \\ \vdots & \vdots & \dots & \vdots \\ u_p^2 & u_p^4 & \dots & u_p^{2w} \end{bmatrix} \begin{bmatrix} c_1 \\ \vdots \\ c_p \end{bmatrix} = \boldsymbol{0} \;. \end{align}\] This is equivalent to \[\begin{align} \diag\big([u_1^2,\dots,u_p^2]^{\top}\big) \cdot \underset{{\rm Vandermonde}}{\underbrace{ \begin{bmatrix} 1 & u_1^2 & \dots & u_1^{2(p-1)} \\ \vdots & \vdots & \dots & \vdots \\ 1 & u_p^2 & \dots & u_p^{2(p-1)} \end{bmatrix}}} \cdot \begin{bmatrix} c_1 \\ \vdots \\ c_p \end{bmatrix} = \boldsymbol{0} \;. \label{eq:vandermonde} \end{align}\tag{168}\] Note that the Vandermonde matrix in 168 has determinant \(\prod_{1 \leq i \leq j < p}(u_i^2 - u_j^2)\), which is nonzero since \(|u_i|\neq |u_j|\) for \(i\neq j\). Multiplying an invertible matrix with a diagonal invertible matrix generates another invertible matrix, and \(\bc\) must be a zero vector which is a contradiction. Hence, there does not exist such \(\bc\) for which \(\bc^{\top} \cdot \sdot(\alpha \bu)\) is constant for every \(\alpha \in \R\).

1.0.4.3 Assumptions of [10]

We clarify that Assumption 2 of [10], which is referred to as “linear interventional faithfulness”, requires nonlinear SCMs. Subsequently, both of their main results (Theorem 1 and Theorem 2) require this assumption. First, note that [10] implicitly explain that their results are for nonlinear SCMs: “In general, we show in Appendix B that a large class of nonlinear SCMs and soft interventions satisfy this assumption.” Next, we show that this assumption is violated for linear additive noise models as follows. First, we quote their assumption:

1.0.4.3.0.1 Linear interventional faithfulness

 [10]: Intervention \(I\) on node \(i\) satisfies linear interventional faithfulness if for every \(j \in \Ch(i) \cup \{i\}\) such that \(\Pa(j) \cap \De(i) = \emptyset\), it holds that \(\mathbb{P}(Z_j + \bc^{\top} \bZ_S) \neq \mathbb{P}^I(Z_j + \bc^{\top} \bZ_S)\) for all constant vectors \(\bc \in \R^{|S|}\) where \(S = [n] \setminus \De(i)\).

1.0.4.3.0.2 Linear SCMs violate linear interventional faithfulness.

Next, we show that linear SCMs violate the above linear interventional faithfulness assumption. Consider an intervention on node \(i\), and let \(j \in \Ch(i)\) such that \(\Pa(i) \cap \De(i) = \emptyset\). Let \[Z_j = \bw^\top \bZ_{\Pa(j)} + N_j \;,\] where \(\bw \in \R^{|\Pa(j)|}\) denotes the weight vector. Note that we have \(\Pa(j) \subseteq S = [n] \setminus \De(i)\) since \(\Pa(i) \cap \De(i) = \emptyset\). Consider vector \(\bc \in \R^{|S|}\) such that \(\bc_{\Pa(j)} = - \bW\) and \(\bc_{S \setminus \Pa(j)} = \boldsymbol{0}\) where \(\bc_{\Pa(j)}\) denotes the entries corresponding to nodes in \(\Pa(j)\). Then, we have \[Z_j + \bc^{\top}\bZ_S = N_j \;,\] which remains invariant after an intervention on node \(i\) and yields \(\mathbb{P}(Z_j + \bc^{\top} \bZ_S) = \mathbb{P}^I(Z_j + \bc^{\top} \bZ_S)\). Therefore, linear SCMs violate linear interventional faithfulness.

1.0.5 Empirical Evaluations: Details and Additional Results↩︎

1.0.5.1 Implementation Details of LSCALE-I Algorithm

1.0.5.1.0.1 Preprocessing: Dimensionality reduction.

Since \(\bX = \bG \cdot \bZ\), we can compute \(\image(\bG)\) using \(n\) random samples of \(\bX\) almost surely. Specifically, \(\image(\bG)\) equals the column space of the sample covariance matrix of \(n\) samples of \(\bX\). In LSCALE-I, as a preprocessing step, we compute this subspace and express samples of \(\bX\) and \(\bs_\bX\) in this basis. This procedure effectively reduces the dimension of \(\bX\) and \(\bs_\bX\) from \(d\) to \(n\). Then, we perform steps of LSCALE-I using this \(n\) dimensional observed data.

1.0.5.1.0.2 Implementation details of 3.

In Stage L1 of 3, choosing any \(\bh \in \colspace(\bR_\bX^m)\) suffices. To make our algorithm deterministic, we choose the top eigenvector of \(\colspace(\bR_\bX^m)\) that corresponds to the largest eigenvalue. In Stage L2, we use a nonzero threshold \(\lambda_{\mcG}\), i.e., \[\label{eq:linear-hard-parent-set-empirical} \hat{\Pa}(m) \triangleq \Big\{ i \neq m : \E\Big[\big| \bsd_{\hat{\bZ}}^{m}(\hat{\bZ}; \hat{\bH})\big|_i \Big] \geq \lambda_{\mcG} \Big\} \;.\tag{169}\] In 15, we list all \(\lambda_{\mcG}\) values used in the experiments.

1.0.5.1.0.3 Implementation details of 4.

The latent graph estimation part of 4 involves rank tests, i.e., determining the dimension of the intersection of column spaces. Therefore, we also need a threshold, denoted by \(\lambda_{\rm eigv}\) to determine the numerical ranks in practice. In experiments reported in 8, we set \(\lambda_{\rm eigv}=0.01\).

1.0.5.1.0.4 Using the algorithms of related work.

For the comparisons in 0.0.7.3, we used the shared sources codes of [8]¹¹ and [10]¹².

1.0.5.2 Implementation Details of GSCALE-I Algorithm

We perform experiments for the coupled interventions setting, and solve the following optimization problem in Step G2 of GSCALE-I, \[\label{eq:gscale-i-l1} \min_{h} \big\|\bD_{\rm t}(h) - \bI_{n \times n}\big\|_{1,1} + \lambda \E\Big[ \big\| h^{-1}(h(\bX)) - \bX \big\|_2^2\Big] \;.\tag{170}\] where \(\lambda > 0\) is a regularization parameter to ensure injectivity. In the following, first, we show why this problem is equivalent to solving 51 . Then, we describe the computation of the ground truth score differences for \(\bX\) and discuss other implementation details.

1.0.5.2.0.1 Implementation steps.

We use \(n_{\rm s}\) samples from the observational environment to compute empirical expectations. Since encoder \(h\) is parameterized by \(\bH\), we use gradient descent to learn this matrix. To do so, we minimize the loss described above in 170 . We denote the final parameter estimate by \(\bH^*\) and the encoder parameterized by \(\bH^*\) by \(h^*\).

In the simulation results reported in 0.0.7.2, we set \(\lambda = 1\) and solve 170 using RMSprop optimizer with learning rate \(10^{-3}\) for \(3 \times 10^{4}\) steps for \(n=5\) and \(4 \times 10^{4}\) steps for \(n=8\). We also use early stopping when the training converges before the maximum number of steps. Recall that the latent graph estimate \(\hat{\mcG}\) is constructed using \(\mathbb{1}\{\bD(h^*)\}\). Similarly to 169 , We use a threshold \(\lambda_{\mcG}\) to obtain the graph from the upper triangular part of \(\bD(h^*)\).

1.0.5.2.0.2 Architecture details and hyperparameter selection for image experiments.

In 16 and 17, we list the architecture details and training hyperparameters used in the image experiments in 0.0.7.5. Specifically, 16 specifies the autoencoder architecture we used to estimate the latent variables, and 17 details the classification-based score difference network we trained. Finally, 18 lists the first level CNN-based autoencoder we used for image experiments with 2 balls.

1.0.5.3 Experimental Methodology Comparison with DiscrepancyVAE

One major difference between the LSCALE-I algorithm and the publicly available DiscrepancyVAE implementation is that DiscrepancyVAE tries to learn the intervention effects as an encoding of the intervention labels. Specifically, for each intervention target, the VAE model learns a vector of latent intervention targets (a neural network with softmax activation) to denote which latent node is intervened on. In practice, learning an intervention encoding runs the risk of learning encodings where different interventions affect the same latent node, consequently, having nodes that are not affected by any intervention encoding. We demonstrate this phenomenon for DiscrepancyVAE in 11, plotting the learned intervention embeddings for each single-node intervention. We observe that despite having exhaustive single-node interventions, the intervention embedding network fails to learn exhaustive intervention targets/effects. Therefore, the learned graph is not immediately interpretable as a DAG on \(n\) causal latent variables since the learned intervention targets (i.e., node labels) are not surjective.

In contrast, in LSCALE-I, the effect of an intervention is estimated using score function differences. While the distinction may seem small, this approach ensures that the estimated graph’s nodes are always interpretable since score differences between observational and interventional environments directly correspond to the intervention targets.

1.0.5.4 Visualizations for Intervention Extrapolation on Biological Data

In Figure 12, we plot certain double-node intervention samples, both actual and sampled using our score-difference Langevin method in a UMAP [62] plot. We plot the double-node interventions that [10] plots in their double-node intervention extrapolation results that have more than \(1000\) samples for both single and double-node environments.sec:discrvae-invert-enc-problem

We see across intervention pairs, score-based intervention extrapolation method leads to a sampling pattern that is relatively compliant with the actual double node interventional data. However, we observe that in some pairs, intervention effects are overestimated (particularly for CEBPE+SPI and ETS2+MAPK1).

1.0.5.5 Details of Score Functions Used in the Experiments

1.0.5.5.0.1 Score function of the linear Gaussian model.

For a Gaussian random vector \(\bY \sim \mcN(\mu, \Sigma)\), the score function of \(\bY\) is given by \[\bs_\bY(\by) = - \Sigma^{-1} \cdot (\by - \mu) \;.\] Using the linear model specified in 61 , in environment \(\mcE^0\), \[\bZ = (\bI_{n \times n} - \bA)^{-1} \cdot N_i \;.\] Since \(N_i\) is a zero-mean Gaussian random vector and \((\bI_{n \times n} - \bA)^{-1}\) is a full rank matrix, \(\bZ\) is a zero-mean Gaussian random vector. Hence, \[\bs(\bz) = - \big[\Cov(\bZ)\big]^{-1} \cdot \bz \;.\] Note that the covariance of \(\bZ\) has the following form \[\Cov(\bZ) = (\bI_{n \times n} - \bA)^{-1} \cdot \diag\big([\sigma_{1}^2,\dots,\sigma_{n}^2]^{\top}\big) \cdot \big((\bI_{n \times n} - \bA)^{-1}\big)^{\top} \;.\] Similarly, the score function of \(\bZ\) in \(\mcE^m\) is given by \[\bs^m(z) = - \big[\Cov(\bZ^m)\big]^{-1} \cdot \bz \;, \qquad \forall m \in [n] \;,\] in which \(\bZ^m\) denotes the latent variables \(\bZ\) in environment \(\mcE^m\). Finally, by setting \(f=\bG\), 1 specifies the score differences of \(\bX\) in terms of \(\bZ\) under different environment pairs as \[\bs_\bX(\bx) - \bs_\bX^m(\bx) = \big(\bG^{\dagger}\big)^{\top} \cdot \big[\bs(\bz) - \bs^m(\bz)\big] \;.\]

1.0.5.5.0.2 Score estimation for a linear Gaussian model.

For all environments \(\mcE^m \in \mcE\), given \(n_{\rm s}\) i.i.d.samples of \(\bX\), we first compute the sample covariance matrix denoted by \(\hat{\Sigma}^m\). Then, we compute the sample precision matrix as \[\hat{\Theta}^m \triangleq (\hat{\Sigma}^m)^{\dagger} \;,\] which leads to the score function estimate given by \[\hat{\bs}_\bX^m(\bx) \triangleq - \hat{\Theta}^m \cdot \bx \;.\]

1.0.5.5.0.3 Score function of the quadratic model.

Following 14 , score functions \(s^m\) and \(\tilde{s}^m\) are decomposed as follows. \[\begin{align} \bs^{m}(\bz) &= \nabla_\bz \log q_{\ell}(z_{\ell}) + \sum_{i \neq \ell} \nabla_\bz \log p_i(z_i \med \bz_{\Pa(i)}) \;, \tag{171} \\ and \quad \tilde{\bs}^{m}(\bz) &= \nabla_\bz \log \tilde{q}_{\ell}(z_{\ell}) + \sum_{i \neq \ell} \nabla_\bz \log p_i(z_i \med \bz_{\Pa(i)})\;. \tag{172} \end{align}\] For additive noise models, the terms in 171 and 172 have closed-form expressions. Specifically, using 85 and denoting the score functions of the noise terms \(\{N_i : i \in [n]\}\) by \(\{r_i : i \in [n]\}\), we have \[[\bs(\bz)]_i = r_i(n_i)-\sum_{j \in \Ch(i)} \frac{\partial f_i(\bz_{\Pa(j)})}{\partial z_i} \cdot r_j(n_j) \;.\] Recall that we consider a quadratic latent model with \[f_i(\bz_{\Pa(i)}) = \sqrt{\bz_{\Pa(i)}^{\top} \cdot \bQ_i \cdot \bz_{\Pa(i)}} \;, \quad and \quad N_i \sim \mcN(0, \sigma_i^{2}) \;,\] which implies \[\frac{\partial f_j(\bz_{\Pa(j)})}{\partial z_i} = \frac{[\bQ_j]_i \cdot \bz_{\Pa(j)}}{\sqrt{\bz_{\Pa(j)}^{\top} \cdot \bQ_j \cdot \bz_{\Pa(j)}}} \;, \quad and \quad r_i(n_i) = -\frac{n_i}{\sigma_i^2} \;, \quad \forall i \in [n] \;.\] Components of the score functions \(\bs^m\) and \(\tilde{\bs}^m\) can be computed similarly. Subsequently, using 1 of 2, we can compute the score differences of observed variables as follows. \[\begin{align} \bs_\bX(\bx) - \bs_\bX^m(\bx) &= \Big[\big[J_g(\bz)\big]^{\dagger}\Big]^{\top} \cdot \big[\bs(\bz) - \bs^m(\bz)\big] \;, \\ \bs_\bX(\bx) - \tilde{\bs}_\bX^m(\bx) &= \Big[\big[J_g(\bz)\big]^{\dagger}\Big]^{\top} \cdot \big[\bs(\bz) - \tilde{\bs}^m(\bz)\big] \;, \\ \bs_\bX^m(\bx) - \tilde{\bs}_\bX^m(\bx) &= \Big[\big[J_g(\bz)\big]^{\dagger}\Big]^{\top} \cdot \big[\bs^m(\bz) - \tilde{\bs}^m(\bz)\big] \;. \end{align}\]

References↩︎

---

1. Equal contribution↩︎

2. The work was completed when B.V. was at Rensselaer Polytechnic Institute.↩︎

3. Contributions of A.K. were done when the author was at Google Brain.↩︎

4. For details of where this has been used, see 1.0.1.2.↩︎

5. Transitive closure of a DAG \(\mcG\), denoted by \(\mcG_{\rm tr}\), is a DAG with parents denoted by \(\Pa_{\rm tr}(i)=\An(i)\) for each node \(i\). The transitive reduction of \(\mcG\) is the DAG with the fewest edges that preserves the same reachability relation as \(\mcG\).↩︎

6. It is also called topological ordering* or topological sort in the literature.*↩︎

7. Surrounded node concept is first defined by [48], and later adopted by [13] when considering soft interventions.↩︎

8. Some examples include generic interventions in [8], no pure shift interventions condition in [9], and the genericity condition in [15].↩︎

9. The “linear interventional faithfulness” [10] implies nonlinearity, which we elaborate in 1.0.4.3.↩︎

10. The codebase for the algorithms and simulations is available at:
    https://github.com/acarturk-e/score-based-crl.↩︎

11. https://github.com/csquires/linear-causal-disentanglement-via-interventions↩︎

12. https://github.com/uhlerlab/discrepancy_vae↩︎