Learning Invariant Representations under General Interventions on the Response

Kang Du and Yu Xiang K. Du and Y. Xiang are with the Department of Electrical and Computer Engineering, University of Utah, Salt Lake City, UT, 84112, USA. This paper has supplementary downloadable material available at http://ieeexplore.ieee.org., provided by the authors. The material includes all the proofs. Contact {kang.du, yu.xiang}@utah.edu for further questions about this work.

Abstract

It has become increasingly common nowadays to collect observations of feature and response pairs from different environments. As a consequence, one has to apply learned predictors to data with a different distribution due to distribution shifts. One principled approach is to adopt the structural causal models to describe training and test models, following the invariance principle which says that the conditional distribution of the response given its predictors remains the same across environments. However, this principle might be violated in practical settings when the response is intervened. A natural question is whether it is still possible to identify other forms of invariance to facilitate prediction in unseen environments. To shed light on this challenging scenario, we focus on linear structural causal models (SCMs) and introduce invariant matching property (IMP), an explicit relation to capture interventions through an additional feature, leading to an alternative form of invariance that enables a unified treatment of general interventions on the response as well as the predictors. We analyze the asymptotic generalization errors of our method under both the discrete and continuous environment settings, where the continuous case is handled by relating it to the semiparametric varying coefficient models. We present algorithms that show competitive performance compared to existing methods over various experimental settings including a COVID dataset.

Index Terms:

Multi-environment domain adaptation, invariance, structural causal models, semiparametric varying coefficient model.

1 Introduction

How to make reliable predictions in unseen environments that are different from training environments is a challenging problem, which is fundamentally different from the classical machine learning settings [1, 2, 3]. Modeling these distribution shifts in a principled way is of great importance in many fields including robotics, medical imaging, and environmental science. Apparently, this problem is ill-posed without any constraints on the relationship between training and test distributions, as the test distribution may be arbitrary. Consider the problem of predicting the response $Y$ given its predictors $X=(X_{1},...,X_{d})^{\top}$ in unseen environments. To model distribution changes across different environments (or training and test distributions), we follow the approach of using structural causal models (SCMs) [4, 5] to model different data-generating mechanisms. The common assumption is that the assignment for $Y$ does not change across environments (or $Y$ is not intervened), which allows for natural formulations of the invariant conditional distribution of $Y$ given a subset of $X$ [6, 7, 8, 5, 9, 10, 11, 12, 13]. The underlying principle is known as invariance, autonomy or modularity [14, 15, 4, 16].

Following this principle, the invariance-based causal prediction initiated in [8] (also see [17] and [11] and references therein) assumes that the conditional distribution of $Y$ given a set of predictors $X_{S}\subseteq\{X_{1},...,X_{d}\}$ is invariant in all environments, i.e., $\mathcal{P}_{e}(Y|X_{S})=\mathcal{P}_{h}(Y|X_{S})$ for environments $e$ and $h$ , where $(X,Y)$ is generated according to the joint distribution $\mathcal{P}_{e}:=\mathcal{P}_{e}^{X,Y}$ . Focusing on linear SCMs, it assumes the existence of a linear model that is invariant across environments, with an unknown noise distribution and arbitrary dependence among predictors (see extensions to nonlinear [18] and time series [19] settings). Following this framework, theoretical guarantees for domain adaptation have been developed in [9, 20]. More recently, a multi-environment regression method for domain adaptation called the stabilized regression [21] explicitly enforces stability (based on a weaker version of invariance $\mathsf{E}_{\mathcal{P}_{e}}[Y|X_{S}=x_{s}]=\mathsf{E}_{\mathcal{P}_{h}}[Y|X_{S}=x_{s}]$ ) by introducing the stable blanket, which is a refined version of the Markov blanket to promote generalization. The tradeoff between predictive performance on training and test data has been studied via regularization under shift interventions [22]. Motivated by [8], the invariant risk minimization (IRM) [23] only uses data from the training environments (i.e., the out-of-distribution generalization setting), and imposes $\mathcal{P}_{e}(Y|\phi(X))=\mathcal{P}_{h}(Y|\phi(X))$ , where $\phi$ is invariant across environments, leading to a bi-leveled optimization problem that is not practical. Several relaxed versions of IRM have been proposed in [23], but they behave very differently from the original IRM (see, e.g., [24, 25]). For a framework of the out-of-distribution setting from a causal perspective with a focus on minimizing the worst-case risk, see [26] and references therein. In this line of invariance-based work, the fundamental assumption is that interventions on the target variable $Y$ are not allowed.

In practical settings, however, the structural assignment of $Y$ might change across environments, namely, $Y$ might be intervened. How to relax this assumption in a principled way is one of the main motivations in our work. We propose to explore alternative forms of invariance and make an attempt in this direction by focusing on linear SCMs. Concretely, the assignment for $Y$ allows general interventions

\displaystyle Y^{e}=(\beta^{e})^{\top}X^{e}+\varepsilon_{Y}^{e},

where $Y$ can be intervened through coefficient $\beta^{e}$ and/or the noise $\varepsilon_{Y}^{e}$ , to capture the dependence of structural assignment across different environments (preliminary results have been reported by the same authors in [27]). We consider a multi-environment regression setting for domain adaption: There are multiple training data $(X^{e},Y^{e})$ for $e\in\mathcal{E}^{\text{train}}$ that are generated from a training model and one test data (indexed by $e^{\text{test}}$ ) from a test model; we assume the training model and test model follow SCMs with the same but unknown graph structure, but we allow $\beta^{e}$ and the mean and variance of $\varepsilon_{Y}^{e}$ to be arbitrarily different under the two models. To avoid the setting to be ill-posed, a key necessary condition is that $Y^{e}$ needs to have at least one child in the SCMs, as prediction is not possible otherwise given that $Y^{e}$ may change arbitrarily over environments. The main challenge lies in whether it is still possible to identify other forms of invariance to facilitate prediction in the test environment. We propose an alternative form of invariance $\mathcal{P}_{e}(Y|\phi_{e}(X))=\mathcal{P}_{h}(Y|\phi_{h}(X))$ that is enabled by a family of conditional invariant transforms $\Phi\ni\phi_{e},\phi_{h}$ . Under general interventions on $Y$ , we provide explicit constructions of such transforms by developing invariant matching property (IMP), a deterministic relation between an estimator of $Y$ and $X$ along with an additional predictor constructed from $X$ . To enable a systematic way of constructing the IMP, we provide a natural decomposition of it and demonstrate this when only $Y$ is intervened or both $X$ and $Y$ are intervened. The IMP comes with several appealing features: (1) it does not vary over environments, making it applicable in unseen environments, and (2) the identification of the IMP follows directly from the fact that the training data contains multiple environments. We study the asymptotic generalization error for both the discrete environment setting and continuous counterpart, which is the more challenging setting. Interestingly, we reveal a connection between the continuous environment setting with the semiparametric varying models, which makes the asymptotic generalization analysis possible. We believe that our results open up new possibilities for multi-environment regression methods for domain adaptation under the structure causal models.

1-A Further Related Work

In [28], the authors have provided a systematic treatment of domain adaptation using the SCMs to enable analysis and comparisons of domain adaption methods, which leads to the conditional invariant residual matching (CIRM) method. The CIRM and its variants combine the domain invariant projection (DIP)-type methods (see [29, 30] and the generalized label shift to handle target label perturbation [31, 32]) with the idea of conditional invariance penalty (appeared in [33, 34] under slightly different settings) that assumes the existence of conditional invariant components (CICs) in the anticausal setting where $Y$ causes $X$ . Theoretical guarantees have been provided for the prediction performance under shift interventions on $Y$ , while numerical studies are provided for interventions on the noise variance of $Y$ [28]. It has also been pointed out that the general mixed-causal-anticausal domain adaptation problem remains open. We aim to shed light on this challenging setting by constructing explicit conditional invariant transforms.

The role of causality in facilitating domain adaptation problems is first articulated in [6], focusing on causal and anticausal predictions. Reweighting methods have been extensively studied for covariate shift [1, 35, 36], which assumes that only the feature distribution changes over environments while the conditionals remain the same. The label shift, which aligned with the anticausal setting, has attracted much attention recently [37, 38, 39]. Many other interesting domain adaption methods have been developed but they are less related to this work. The performance bounds using Vapnik-Chervonenkis (VC) theory have been initiated in [40]. There are fundamental works from the robust statistics perspective including distributional robust learning [41, 42, 43, 44, 45, 46] and adversarial machine learning [47, 48].

1-B Contribution and Structure

There are four main contributions in this work. (I) We formulate a general invariance property and the corresponding conditional invariant transforms for analyzing general interventions on linear SCMs (Section 2). (II) We tackle this problem by introducing the invariant matching property (IMP) and providing a systematic approach for establishing explicit characterization of it (Section 3 and Section 4). (III) To handle the continuous environment setting, we bridge our framework with the profile likelihood estimators developed in the semiparametric literature, leading to asymptotic performance guarantees under this challenging setting (Section 6). (IV) Motivated by our theoretical results, we develop efficient algorithms that show competitive empirical performance over a variety of simulation settings including a COVID dataset (Section 5 and Section 7). All the technical details are deferred to Appendix.

2 Background and Problem Formulation

Consider a response $Y\in\mathbb{R}$ and a vector of predictors $X=(X_{1},\ldots,X_{d})^{\top}\in\mathcal{X}\subseteq\mathbb{R}^{d}$ following an acyclic linear SCM,

	$\displaystyle\mathcal{M}:$	$\displaystyle X=\gamma Y+BX+\varepsilon_{X}$
	$\displaystyle\mathcal{M}:$	$\displaystyle Y=\beta^{\top}X+\varepsilon_{Y},$

where $\beta,\gamma\in\mathbb{R}^{d}$ , $B\in\mathbb{R}^{d\times d}$ , $\varepsilon_{X}=(\varepsilon_{X_{1}},\ldots,\varepsilon_{X_{d}})^{\top}$ , and the noise variables $\{\varepsilon_{X_{1}},\ldots,\varepsilon_{X_{d}}\}$ and $\varepsilon_{Y}$ are jointly independent. We use $\mathcal{G}(\mathcal{M})$ to denote the directed acyclic graph induced by $\mathcal{M}$ , with edges determined by the non-zero coefficients in $\mathcal{M}$ . We denote the parents, children, descendants, and Markov blanket of a variable $Z\in\{X_{1},\ldots,X_{d},Y\}$ as $PA(Z)$ , $CH(Z)$ , $DE(Z)$ , and $MB(Z)$ , respectively. When $(X,Y)$ is observed in different environments (e.g., different experiment settings for data collection), the parameters $\{\beta,\gamma,B\}$ and the distributions of $\{\varepsilon_{X},\varepsilon_{Y}\}$ may change. In the following, we use interventions on the SCM $\mathcal{M}$ to model such changes.

Let $\mathcal{E}^{\text{all}}$ denote the set of all possible environments ¹¹1We use training (or test) environments and observable (or unseen) environments interchangeably., which are partitioned into multiple training environments $\mathcal{E}^{\text{train}}$ and one test environment $\{e^{\text{test}}\}$ such that $\mathcal{E}^{\text{all}}=\mathcal{E}^{\text{train}}\cup\{e^{\text{test}}\}$ . In each environment $e\in\mathcal{E}^{\text{all}}$ , $(X,Y)$ is generated according to the joint distribution $\mathcal{P}_{e}:=\mathcal{P}_{e}^{X,Y}$ , and to simplify notation, we write $\mathcal{P}_{\text{test}}:=\mathcal{P}_{e^{\text{test}}}^{X,Y}$ . A variable from $\{X_{1},\ldots,X_{d},Y\}$ is intervened if the parameters or noise distribution in its assignment changes over different $e\in\mathcal{E}^{\text{all}}$ . For instance, the changes of $\beta$ and/or the distribution of $\varepsilon_{Y}$ correspond to the intervention on $Y$ (see different intervention cases below). Importantly, we allow both $Y$ and any subset of $\{X_{1},\ldots,X_{d}\}$ to be intervened in $\mathcal{E}^{\text{all}}$ .

Now, we introduce the linear SCM $\mathcal{M}$ with parameters that change with environments. For each $e\in\mathcal{E}^{\text{all}}$ , the linear SCM $\mathcal{M}$ is modified to be

	$\displaystyle\mathcal{M}^{e}:$	$\displaystyle X^{e}=\gamma^{e}Y^{e}+B^{e}X^{e}+\varepsilon_{X}^{e}$			(1)
	$\displaystyle\mathcal{M}^{e}:$	$\displaystyle Y^{e}=(\beta^{e})^{\top}X^{e}+\varepsilon_{Y}^{e}.$			(2)

This formulation is fairly general. From the structural perspective, this consists of causal, anticausal, and mixed-causal-anticausal settings [6]. It should be noted that we only adopt the linear SCM rather than the fully specified SCMs as in [49], since learning the functional forms can be more complicated than the prediction problem we aim to solve. Regarding the intervention types, we discuss several special cases to put them into perspective.

1.

Shift interventions on $X$ or $Y$ : A variable $X_{j}$ is intervened through a shift if the mean of the noise variable $\varepsilon_{X,j}^{e}$ changes with $e\in\mathcal{E}^{\text{all}}$ . For the shift intervention on $Y$ , the mean of $\varepsilon_{Y}^{e}$ changes.
2.

Interventions on the coefficients of $X$ or $Y$ : A variable $X_{j}$ is intervened through coefficients if the coefficients $\{\gamma_{j}^{e},B_{j\cdot}^{e}\}$ change with $e\in\mathcal{E}^{\text{all}}$ . For $Y$ , the change is on the coefficient vector $\beta^{e}$ .
3.

Interventions on the noise variance of $X$ or $Y$ : Similar to shift interventions, a variable $X_{j}$ or $Y$ is intervened if its noise variance changes.

We observe $n^{e}$ i.i.d. samples $\{(x_{1},y_{1}),...,(x_{n^{e}},y_{n^{e}})\}$ from each training environment distribution $\mathcal{P}_{e}$ for $e\in\mathcal{E}^{\text{train}}$ , but in the test environment $e^{\text{test}}$ we only observe $m$ i.i.d. samples $\{x_{1},...,x_{m}\}$ from $\mathcal{P}_{\text{test}}^{X}$ . The goal is to learn a function $f:\mathcal{X}\to\hat{\mathcal{Y}}$ that works well on $e^{\text{test}}$ in the sense that it minimizes the test population loss

\mathcal{L}_{\text{test}}(f):=\mathsf{E}_{(X,Y)\sim\mathcal{P}_{\text{test}}}[l(Y,f(X))].

(3)

where $l$ is the square loss function $l(y,\hat{y})=(y-\hat{y})^{2}$ . The optimal function is $f(x)=\mathsf{E}_{\mathcal{P}_{\text{test}}}[Y|X=x]$ , which cannot be learned from the observed data in general when $X$ and/or $Y$ is intervened. Without any constraints on the relationship between $\mathcal{P}_{\text{test}}$ and $\mathcal{P}_{e}$ for $e\in\mathcal{E}^{\text{train}}$ , the test population loss can be arbitrarily large. To make this problem tractable, we assume that $(X,Y)$ under $\mathcal{P}_{\text{test}}$ and $\mathcal{P}_{e}$ are generated according to the SCMs described in (2) but we do not assume that the causal graph is known and we allow for general types of interventions.

It is well-known that if $Y$ is not intervened, a general form of invariance principle applies, assuming the existence of some subset $S\subseteq\{1,...,d\}$ such that

\mathcal{P}_{e}(Y|X_{S})=\mathcal{P}_{h}(Y|X_{S})

(4)

holds for any $e,h\in\mathcal{E}^{\text{all}}$ . Under this assumption, the causal function $\mathsf{E}_{\mathcal{P}_{e}}[Y|X_{PA(Y)}]$ is invariant across different environments and minimax under the class of all possible interventions on $X$ [9]. If not every predictor is intervened arbitrarily, previous works that are motivated by (4) (as mentioned in Section 1) aim to improve upon the causal function. The main challenge in our setting comes from the general interventions on $Y$ , making the traditional invariance principle not applicable. Importantly, the causal function can change arbitrarily with environments. In this work, we propose to exploit an alternative form of invariance to tackle this problem.

Definition 1.

A function $\phi:(\mathcal{E}^{\text{all}},\mathbb{R}^{d})\to\mathbb{R}^{q}$ is called a conditional invariant transform if the following invariance property holds for any $e,h\in\mathcal{E}^{\text{all}}$

\mathcal{P}_{e}(Y|\phi_{e}(X))=\mathcal{P}_{h}(Y|\phi_{h}(X)).

(5)

Under general intervention settings, we denote this class of conditional invariant transforms as $\Phi$ , and we provide explicit characterizations of it via the invariant matching property (IMP) (see Definition 2). For each $\phi\in\Phi$ , the invariance property (5) enables us to compute

f_{\phi_{e}}(x)=g\circ\phi_{e}(x)=\mathsf{E}_{\mathcal{P}_{\text{e}}}[Y|\phi_{e}(x)],

(6)

for any $e\in\mathcal{E}^{\text{all}}$ , where the function $g:\mathbb{R}^{q}\to\hat{\mathcal{Y}}$ is invariant across environments and is nonlinear in general. Equivalently, this solves a relaxed version of (3) by minimizing $\mathcal{L}_{\text{test}}(f_{\phi})$ over $\{\phi\in\Phi\}$ . This formulation allows us to treat the general mixed-causal-anticausal problem under general interventions on $Y$ in a unified manner.

3 Invariant Matching Property and Theoretical Guarantees

Among the intervention settings below (2), the interventions on only $Y$ are important but rarely studied, which concerns the changes of $\beta^{e}$ and the distribution of $\varepsilon^{e}_{Y}$ . Our method can be motivated in this setting by the following observation: If $X$ includes any descendants of $Y$ , then $\beta$ and $\varepsilon_{Y}$ (that may change arbitrarily in the unseen environments) will be passed on to the descendants. Thus, the changes might be revealed by the changes of certain statistical properties of $X$ , leading to our proposed invariant matching property detailed in this section. We start with a few toy examples.

3-A Motivating Examples

Example:

Consider $(Y^{e},X^{e}),e\in\mathcal{E}^{\text{toy}}=\{1,2\}$ , with $X^{e}:=(X_{1}^{e},X_{2}^{e},X_{3}^{e})^{\top}$ satisfying the following linear acyclic SCM (illustrated in Fig. 1),

Figure 1: Directed acyclic graph

\mathcal{G}(\mathcal{M}_{\text{toy}}^{e})

\displaystyle\mathcal{M}_{\text{toy}}^{e}:\begin{cases}Y^{e}=a^{e}X_{1}^{e}+X_{2}^{e}+N_{Y}^{e}\\ X_{3}^{e}=Y^{e}+X_{1}^{e}+N_{3}^{e},\end{cases}

(7)

where $X_{1}^{e},X_{2}^{e},N_{3}^{e}$ , and $N_{Y}^{e}$ are independent and $\mathcal{N}(0,1)$ -distributed for every $e\in\mathcal{E}^{\text{toy}}$ . Since $(Y^{e},X^{e})$ is multivariate Gaussian, the MMSE estimator of $Y^{e}$ given $X^{e}$ is

	$\displaystyle\mathsf{E}_{\mathcal{P}_{e}}[Y\|X]=$	$\displaystyle X^{\top}\left(\mathsf{E}_{\mathcal{P}_{e}}[XX^{\top}]\right)^{-1}\mathsf{E}_{\mathcal{P}_{e}}[XY]$
	$\displaystyle=$	$\displaystyle\frac{1}{2}(a^{e}-1)X_{1}^{e}+\frac{1}{2}X_{2}^{e}+\frac{1}{2}X_{3}^{e},$

Refer to caption — Figure 2: Illustrations of invariant relations in Example Example. In (a), we illustrate the linear invariant relation (8) by visualizing estimates³³3For $e\in\{1,2\}$ , $X_{4}^{e}$ and $\mathsf{E}_{\mathcal{P}_{e}}[Y|X]$ are estimated using OLS. Then, $\eta$ (denoted as $\eta_{\text{ⓐ}}$ in Figure 2.(a)) is estimated using OLS by regressing $\mathsf{E}_{\mathcal{P}_{e}}[Y|X]$ on $(X_{4}^{e},X^{e})$ using the pooled data. of the tuple $(X_{4}^{e}$ , $\eta^{\top}X^{e},\mathsf{E}_{\mathcal{P}_{e}}[Y|X])$ (corresponding to the $\{x,y,z\}$ axes). Similarly, for (b) and (c), we verify such a relation for $X_{5}^{e}$ and $X_{6}^{e}$ . The overlaps between the red dots ( $e=1$ ) and the blue dots ( $e=2$ ) in (a) and (b) indicate the invariant linear relations. However, the red and blue dots in (c) are not aligned, implying that no invariant relations as in (8) and (9) hold for $(\mathsf{E}_{\mathcal{P}_{e}}[Y|X],X_{6}^{e})$ .

which is not directly applicable for predicting $Y^{e^{\text{test}}}$ as $a^{e}$ can change arbitrarily. Similarly, one can compute $\mathsf{E}_{\mathcal{P}_{e}}[X_{3}|X_{1},X_{2}]=(1+a^{e})X_{1}^{e}+X_{2}^{e}$ . It is noteworthy that $\mathsf{E}_{\mathcal{P}_{e}}[Y|X]$ and $\mathsf{E}_{\mathcal{P}_{e}}[X_{3}|X_{1},X_{2}]$ can each be generated by a linear combination of $\mathsf{E}_{\mathcal{P}_{e}}[Y|X_{PA(Y)}]=a^{e}X_{1}^{e}+X_{2}^{e}$ and $\{X_{1}^{e},X_{2}^{e},X_{3}^{e}\}$ , and we will highlight this observation by introducing two invariant relations (see Def. 4 and Def. 5). As a result, there exists a deterministic linear relation, which we refer to as matching,

\displaystyle\mathsf{E}_{\mathcal{P}_{e}}[Y|X]

\displaystyle=\lambda\mathsf{E}_{\mathcal{P}_{e}}[X_{3}|X_{1},X_{2}]+\eta^{\top}X^{e},

(8)

with coefficients $\lambda=1/2$ and $\eta=(-1,0,1/2)^{\top}$ that are invariant with respect to the environment (illustrated in Fig. 2.(a)-(b)). For simplicity of notation, we denote $X_{4}^{e}:=\mathsf{E}_{\mathcal{P}_{e}}[X_{3}|X_{1},X_{2}]$ in Fig. 2. Moreover, one can verify that $\mathcal{P}_{e}(Y|X,X_{4})$ is invariant since $\mathsf{E}_{\mathcal{P}_{e}}[Y|X,X_{4}]$ and $\mathop{\rm Var}\nolimits_{\mathcal{P}_{e}}(Y|X,X_{4})$ are invariant. Thus $\phi_{e}(X)=(X^{\top},X_{4})^{\top}$ satisfies the invariance property (5). A prediction model in (8) with invariant coefficients is often not unique when it exists. One can show that

\displaystyle\mathsf{E}_{\mathcal{P}_{e}}[Y|X]=-X_{1}^{e}+\frac{1}{2}X_{2}^{e}+X_{3}^{e}-\frac{3}{2}X_{5}^{e},

(9)

with $X_{5}^{e}:=\mathsf{E}_{\mathcal{P}_{e}}[X_{2}|X_{1},X_{3}]$ . However, invariant relations in (8) and (9) do not hold for $X_{6}^{e}:=\mathsf{E}_{\mathcal{P}_{e}}[X_{1}|X_{2},X_{3}]$ , since

\displaystyle X_{6}^{e}=-\frac{a^{e}+1}{(a^{e}+1)^{2}+2}X_{2}^{e}+\frac{a^{e}+1}{(a^{e}+1)^{2}+2}X_{3}^{e}

depends on $e$ in a more complicated form so that there is no linear invariant relation between $(X_{6}^{e},X^{e})$ and $\mathsf{E}_{\mathcal{P}_{e}}[Y|X]$ , as illustrated in Fig. 2.(c). $\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\blacksquare$

In the next example, we extend Example Example to allow for interventions on $X_{1}$ , $X_{2}$ , and $Y$ through the means and/or variances of the noise variables (see Fig. 3).

Example (additional settings):

Consider model $\mathcal{M}_{\text{toy}}^{e}$ in Example Example with additional shift interventions and/or interventions on the noise variances. The results are summarized in the following.

1.

Under shift interventions on $Y$ and $X_{1}$ and an intervention on the variance of $X_{1}$ as illustrated in Fig. 3(a), the two invariant relations (8) and (9) hold.
2.

When $Y$ is intervened through the variance of $N_{Y}$ as illustrated in Fig. 1, the relations (8) and (9) will not hold. In this case, new relations can be established if $\mathsf{E}_{\mathcal{P}_{e}}[Y|X]$ in (8) and (9) is replaced by $\mathsf{E}_{\mathcal{P}_{e}}[Y|X_{PA(Y)}]=\mathsf{E}_{\mathcal{P}_{e}}[Y|X_{1},X_{2}]$ .
3.

When $X_{2}$ is intervened through either the mean or variance as illustrated in Fig. 3(b), a relation as in (9) that is based on $X_{5}^{e}$ will not hold.
4.

Combining the interventions above along with an intervention on $X_{3}$ through the noise variance as illustrated in Fig. 3 (c), there will be one invariant relation left,

$\mathsf{E}_{\mathcal{P}_{e}}[Y|X_{1},X_{2}]=X_{4}^{e}-X_{1}^{e}+b,$ (10)

for some intercept $b\in\mathbb{R}$ . It is noteworthy that (10) will fail to hold if $X_{3}$ is intervened through the coefficients or noise mean. However, due to the intervention on the noise variance of $Y$ , $\mathcal{P}_{e}(Y|X_{S},X_{4})$ is not invariant for any $S\subseteq\{1,2,3\}$ , since $\mathop{\rm Var}\nolimits_{\mathcal{P}_{e}}(Y|X_{S},X_{4})$ changes with $e$ .

Remark 1.

In [50], the authors have shown that a form of varying filter connecting the features and response (as a special case of the varying coefficients captured by $\beta^{e}$ in (2)) is effective for causal inference tasks, by adopting estimators from [51].

$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\blacksquare$

Example (cyclic graph setting):

Consider model $\mathcal{M}_{\text{toy}}^{e}$ in Example Example with $X_{2}^{e}$ generated by $X^{e}_{2}=aX^{e}_{3}+N^{e}_{2}$ such that $a\neq 1$ , which creates a cyclic model as illustrated in Fig. 4(a). One can compute that the relation in (8) still holds. Moreover, the same relation holds even when $X_{2}$ is intervened through the coefficient $a$ (i.e., when it changes over $e$ , denoted by $a^{e}$ ) and/or the noise $N^{e}_{2}$ , illustrated in Fig. 4(b); this example can be further generalized to more complicated cyclic models. This cyclic setting highlights the point that the matching described in (8) can hold without assuming acyclic causal models, which will be left for future work. Focusing on acyclic causal models, in the subsequent sections, we formalize the matching conditions and provide a systematic analysis of sufficient conditions through a natural decomposition. $\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\blacksquare$

Remark 2.

The condition on $a$ avoids the cyclic model to be ill-defined, since combining the structural equations under $a=1$ yields a contradiction where the independent noise variables satisfy a linear constraint.

3-B Invariant Matching Property

In this section, we generalize the invariant relations observed in Example Example and Example Example to a class of such relations for $\mathcal{M}^{e}$ , $e\in\mathcal{E}^{\text{all}}$ without assuming joint Gaussian distributions, and connect this with the invariance property in (5). For the identification of such relations, we show that even two training environments suffice (see Proposition 1).

To handle non-Gaussian cases (beyond Example Example), we choose to adopt the linear MMSE (or LMMSE) estimators for constructing linear invariant relations. For a target variable $Y\in\mathbb{R}$ given a vector of predictors $X\in\mathbb{R}^{p}$ , the LMMSE estimator is defined as

\mathsf{E}_{l}[Y|X]:=(\theta^{\text{ols}})^{\top}(X-\mathsf{E}[X])+\mathsf{E}[Y],

where $\theta^{\text{ols}}:=\mathop{\rm Cov}\nolimits(X,X)^{-1}\mathop{\rm Cov}\nolimits(X,Y)$ is called the population ordinary least squares (OLS) estimator. With a slight abuse of notation, we write $\mathsf{E}_{l,\mathcal{P}_{e}}[Y\,|\,X]$ to denote the LMMSE of $Y$ given $X$ with respect to $(X,Y)\sim\mathcal{P}_{e}$ . To simplify presentation, we focus on $(X^{e},Y^{e})$ with zero means for each $e\in\mathcal{E}^{\text{all}}$ (or equivalently all the noise variables have zero means), while the non-zero mean settings can be handled by introducing the constant one as an additional predictor.

Definition 2.

For $k\in\{1,\ldots,d\}$ , $R\subseteq\{1,\ldots,d\}\setminus k$ , and $S\subseteq\{1,\ldots,d\}$ , we say that the tuple $(k,R,S)$ satisfies the invariant matching property (IMP) if, for every $e\in\mathcal{E}^{\text{all}}$ ,

\displaystyle\mathsf{E}_{l,\mathcal{P}_{e}}[Y\,|\,X_{S}]=\lambda\mathsf{E}_{l,\mathcal{P}_{e}}[X_{k}\,|\,X_{R}]+\eta^{\top}X^{e},

(11)

for some $\lambda\in\mathbb{R}$ and $\eta\in\mathbb{R}^{d}$ that do not depend on $e$ . We denote $\mathcal{I}_{\mathcal{M}}:=\{(k,R,S):~{}\eqref{invar_coeff}\text{ holds}\}$ for model $\mathcal{M}$ , and we call $(\eta^{\top},\lambda)^{\top}$ the matching parameters.

For a tuple $(K,R,S)$ such that the IMP does not hold, we simply call it a non-IMP. Observe that $\mathsf{E}_{l,\mathcal{P}_{e}}[Y|X_{S}]$ is not directly applicable to the test environment due to its components depending on $e$ , but those components are fully captured by $\mathsf{E}_{l,\mathcal{P}_{e}}[X_{k}|X_{R}]$ . If the matching parameters are identified from the training environments, the IMP is applicable to the test environment since $\mathsf{E}_{l,\mathcal{P}_{\text{test}}}[X_{k}|X_{R}]$ is completely determined by the distribution of $X^{e^{\text{test}}}$ without the need of the target $Y^{e^{\text{test}}}$ . Since computing the additional feature $\mathsf{E}_{l,\mathcal{P}_{e}}[X_{k}|X_{R}]$ is simply the prediction of $X_{k}$ (as if $X_{k}$ is not observed), the IMP indicates that the prediction of $Y$ can benefit from the predictions of certain predictors. We formally define this class of additional features as follows.

Definition 3.

For any $k\in\{1,\ldots,d\}$ and $R\subseteq\{1,\ldots,d\}\setminus k$ , we call $\mathsf{E}_{l,\mathcal{P}_{e}}[X_{k}\,|\,X_{R}]$ a prediction module. If a prediction module satisfies an IMP for some $S\subseteq\{1,\ldots,d\}$ , we call it a matched prediction module for $S$ .

Now we discuss the relationship between the IMP and the invariance property $\mathcal{P}_{e}(Y|\phi_{e}(X))=\mathcal{P}_{h}(Y|\phi_{h}(X))$ in (5), we rewrite (11) in a compact form as

\mathsf{E}_{l,\mathcal{P}_{e}}[Y|X_{S}]=\theta^{\top}\widetilde{X}^{e},

(12)

where

\widetilde{X}^{e}:=(X_{1}^{e},\ldots,X_{d}^{e},\mathsf{E}_{l,\mathcal{P}_{e}}[X_{k}|X_{R}])^{\top},

(13)

and $\theta=(\eta^{\top},\lambda)^{\top}$ denotes the matching parameter. Define

\phi_{e}^{(k,R,S)}(X^{e}):=(X_{S^{\prime}}^{e},\mathsf{E}_{l,\mathcal{P}_{e}}[X_{k}|X_{R}])^{\top},

(14)

where $X_{S^{\prime}}^{e}$ is a row vector for some $S^{\prime}\subseteq\{1,\ldots,d\}$ . For notational convenience, we introduce the shorthand $\phi_{e}(X^{e}):=\phi_{e}^{(k,R,S)}(X^{e})$ . Note that $\phi_{e}(X^{e})$ is a linear transform of $X^{e}$ .

In general, however, the invariance of the matching parameter $\theta$ does not imply the invariance property (5). In Section 4, we will characterize a class of IMPs that each satisfies (5). When the invariance property holds, one can apply the general conditional expectation $f_{\phi_{e}}(x)=\mathsf{E}_{\mathcal{P}_{\text{e}}}[Y|\phi_{e}(x)]$ as in (6), since the linear estimator from the IMP is in general sub-optimal for the non-Gaussian cases. We will focus on linear estimators in this work as the extension can be handled via nonlinear regression methods in a straightforward manner.

It is noteworthy that since $\mathsf{E}_{l,\mathcal{P}_{e}}[X_{k}|X_{R}]$ is a linear function of $X_{R}^{e}$ , the matching parameter $\theta$ is not unique given a single environment $e\in\mathcal{E}^{\text{train}}$ . This causes issues when one aims to identify the possible IMPs given the distribution of $(X^{e},Y^{e})$ . However, we show that two environments in the training data are sufficient to identify $\theta$ under a mild assumption on the matched prediction module.

Proposition 1.

For a tuple $(k,R,S)$ that satisfies an IMP, the matching parameter $\theta$ can be uniquely identified in $\mathcal{E}^{\text{train}}$ if $|\mathcal{E}^{\text{train}}|\geq 2$ and

\mathsf{E}_{l,\mathcal{P}_{e}}[X_{k}|X_{R}=x]\neq\mathsf{E}_{l,\mathcal{P}_{h}}[X_{k}|X_{R}=x]

(15)

for some $e,h\in\mathcal{E}^{\text{train}}$ and $x\in\mathcal{X}_{R}$ .

This proposition shows how the heterogeneity of the data generating process can be helpful for identifying important invariant relations.

3-C A Decomposition of the IMP

In our toy examples, recall that the IMPs are derived by first computing $\mathsf{E}_{\mathcal{P}_{e}}[Y|X_{S}]$ and $\mathsf{E}_{\mathcal{P}_{e}}[X_{k}|X_{R}]$ separately and then fitting a linear relation from $(\mathsf{E}_{\mathcal{P}_{e}}[X_{k}|X_{R}],X_{S})$ to $\mathsf{E}_{\mathcal{P}_{e}}[Y|X_{S}]$ . These two steps reveal a natural decomposition of the IMP, which we term as the first and second matching properties below.

Definition 4.

We say that $S\subseteq\{1,\ldots,d\}$ satisfies the first matching property if, for every $e\in\mathcal{E}^{\text{all}}$ ,

\mathsf{E}_{l,\mathcal{P}_{e}}[Y|X_{S}]=\lambda_{Y}\mathsf{E}_{\mathcal{P}_{e}}[Y|X_{PA(Y)}]+\eta_{Y}^{\top}X^{e},

(16)

for some $\lambda_{Y}\in\mathbb{R}$ and $\eta_{Y}\in\mathbb{R}^{d}$ that do not depend on $e$ .

First, observe that the first matching property holds for $S=PA(Y)$ since

\displaystyle\mathsf{E}_{l,\mathcal{P}_{e}}[Y|X_{PA(Y)}]=\mathsf{E}_{\mathcal{P}_{e}}[Y|X_{PA(Y)}]=(\beta^{e})^{\top}X^{e}.

The first matching property concerns the set $S$ such that the components in $\mathsf{E}_{l,\mathcal{P}_{e}}[Y|X_{S}]$ that depend on $e$ are fully captured by the causal function $\mathsf{E}_{\mathcal{P}_{e}}[Y|X_{PA(Y)}]$ . However, this invariant relation is not directly useful for the prediction of $Y^{e^{\text{test}}}$ , since the causal function can change arbitrarily with $e$ . To this end, we identify another invariant relation from $\mathcal{M}^{e}$ which is called the second matching property.

Definition 5.

For $k\in\{1,\ldots,d\}$ and $R\subseteq\{1,\ldots,d\}\setminus k$ , we say that a tuple $(k,R)$ satisfies the second matching property if, for every $e\in\mathcal{E}^{\text{all}}$ ,

\mathsf{E}_{l,\mathcal{P}_{e}}[X_{k}|X_{R}]=\lambda_{X}\mathsf{E}_{\mathcal{P}_{e}}[Y|X_{PA(Y)}]+\eta_{X}^{\top}X^{e},

(17)

for some $\lambda_{X}\in\mathbb{R}$ and $\eta_{X}\in\mathbb{R}^{d}$ that do not depend on $e$ .

It is straightforward to see that, if $\lambda_{X}\neq 0$ in the second matching property, the first and second matching properties imply the IMP as follows,

	$\displaystyle\mathsf{E}_{l,\mathcal{P}_{e}}[Y\|X_{S}]$	$\displaystyle=\frac{\lambda_{Y}}{\lambda_{X}}\mathsf{E}_{l,\mathcal{P}_{e}}[X_{k}\|X_{R}]+\left(\eta_{Y}-\frac{\lambda_{Y}}{\lambda_{X}}\eta_{X}\right)^{\top}X^{e}$
		$\displaystyle:=\lambda\mathsf{E}_{l,\mathcal{P}_{e}}[X_{k}\|X_{R}]+\eta^{\top}X^{e}.$

For prediction tasks under SCMs, the causal function often plays a central role. Our first and second matching properties show how the LMMSE estimator $\mathsf{E}_{l,\mathcal{P}_{e}}[Y|X_{S}]$ and the matched prediction module $\mathsf{E}_{l,\mathcal{P}_{e}}[X_{k}|X_{R}]$ are connected with the causal function, respectively. Together, the two individual connections make up the IMP (illustrated in Fig. 5).

4 Characterization of Invariant Matching Properties

4-A Interventions on the Response

First, we consider model $\mathcal{M}^{e}$ with interventions only on $Y$ through the coefficients, i.e.,

	$\displaystyle\mathcal{M}^{e,1}:$	$\displaystyle X^{e}=\gamma Y^{e}+BX^{e}+\varepsilon_{X}^{e}$			(18)
	$\displaystyle\mathcal{M}^{e,1}:$	$\displaystyle Y^{e}=(\alpha^{e}+\beta)^{\top}X^{e}+\varepsilon_{Y}^{e}.$			(19)

To distinguish the parents of $Y$ with varying and invariant coefficients, we decompose $\beta^{e}$ in $\mathcal{M}^{e}$ into two parts $\alpha^{e}$ and $\beta$ . Without loss of generality, we assume that $\alpha_{j}^{e}\neq 0$ if and only if $\alpha_{j}^{e}$ is a non-constant function of $e$ , and we define the following subset of parents of $Y$ ,

PE=\{j\in\{1,\ldots,d\}:\alpha_{j}^{e}\neq 0\}.

Remark 3.

Note that, in the non-zero mean settings, this model covers the shift intervention on $Y$ through the varying coefficient of the predictor which is a constant one.

Recall that prediction modules do not rely on the response $Y$ but on the relations between the predictors for each environment. When $Y$ is unobserved (or equivalently, substituting $Y$ in (19) into (18)), the relations between the predictors are as follows,

X^{e}=\left(\gamma\left(\alpha^{e}+\beta\right)^{\top}+B\right)X^{e}+\gamma\varepsilon_{Y}+\varepsilon_{X},

(20)

where $\gamma\varepsilon_{Y}+\varepsilon_{X}$ a vector of dependent random variables when $\gamma$ is not a zero vector. If $\alpha^{e}$ vanishes from (20), the distribution of $X^{e}$ becomes invariant with respect to environments. As a consequence, the condition (15) in Proposition 1 will not be satisfied. Observe that $\alpha^{e}$ is non-vanishing in (20) only if $\gamma$ is not a zero vector, which brings up the following key assumption.

Assumption 1.

When $Y$ is intervened, we assume that $Y$ has at least one child.

Note that if $Y$ has no children, by the Markov property of SCMs [4], the test data sampled from $\mathcal{P}_{\text{test}}^{X}$ provides no information about $\mathcal{P}_{\text{test}}^{Y|X}$ and thus the observed $X^{e^{\text{test}}}$ may correspond to two arbitrarily different $Y^{e^{\text{test}}}$ ’s, which makes the problem ill-posed.

The first and second matching properties enable us to characterize the tuples $(k,R,S)$ ’s that satisfy IMPs through the characterizations of $S$ (for the first matching property) and $(k,R)$ (for the second matching property) separately. In the following theorem, we show that a class of IMPs implied by the first and second invariant matching properties satisfy the invariance property (5).

Theorem 1.

For model $\mathcal{M}^{e,1}$ , the first and second matching properties hold in the following cases.

1.

On the first MP: For each $S\subseteq\{1,\ldots,d\}$ such that $PE\subseteq S$ , the first matching property holds.
2.

On the second MP: For each $k\in\{1,\ldots,d\}\setminus PE$ and $R\subseteq\{1,\ldots,d\}\setminus k$ such that $PE\subseteq R$ , the second matching property holds.

For any tuple $(k,R,S)$ above such that $R\subseteq S$ , if $\lambda_{X}\neq 0$ in the second matching property, then $\phi_{e}(X^{e})=(X_{S}^{e},\,\mathsf{E}_{l,\mathcal{P}_{e}}[X_{k}|X_{R}])^{\top}$ satisfies (5). Furthermore, $\mathcal{L}_{\text{test}}(f_{\phi})$ is minimized by any $\phi$ with $S=\{1,...,d\}$ .

As a concrete example, recall (8) from our motivating example, where $R\subseteq S$ and $\lambda_{X}\neq 0$ are satisfied, $\mathsf{E}_{\mathcal{P}_{e}}[Y|X]$ can be represented using invariant coefficients. It is noteworthy that Assumption 1 is a necessary condition for $\lambda_{X}\neq 0$ , and we provide a sufficient condition for $\lambda_{X}\neq 0$ in a concrete setting with $S=\{1,\ldots,d\}$ below.

Corollary 1.

For model $\mathcal{M}^{e,1}$ , the first and second matching properties hold in the following cases.

1.

On the first MP: The first matching property holds for $S=\{1,\ldots,d\}$ .
2.

On the second MP: For each $k\in\{j\in MB(Y):\alpha_{j}^{e}=0\}$ and $R=-k:=\{1,\ldots,d\}\setminus k$ , the second matching property holds.

Proposition 2.

Under Assumption 1, for each $(k,R)$ in Corollary 1, we have $\lambda_{X}\neq 0$ in the second matching property if $B_{-k,k}$ is not in the following hyperplane,

w^{\top}x+b=0,

(21)

where $w\in\mathbb{R}^{d-1}$ and $b\in\mathbb{R}$ are determined by the parameters in $\mathcal{M}^{e,1}$ other than $B_{-k,k}$ .

The explicit expressions of $w$ and $b$ using the parameters in $\mathcal{M}^{e,1}$ are provided in the proof of Proposition 2. For generic choices of the parameters, the second matching property holds with $\lambda_{X}\neq 0$ since $B_{k,-k}$ is not necessarily on the hyperplane described in (21).

When $Y$ is additionally intervened through the noise variance, the proof of Theorem 1 will break down in general (see Remark 8 in Appendix C). However, recall that the first matching property holds for $S=PA(Y)$ by definition. In this case, we provide an example for the second matching property in the following corollary.

Corollary 2.

Under Assumption 1, if $Y$ is intervened through the noise variance in model $\mathcal{M}^{e,1}$ , the second matching property holds for $k\in CH(Y)$ such that $k\not\in DE(i)$ for any $i\in CH(Y)\setminus k$ , and $R=\{1,\ldots,d\}\setminus DE(Y)$ .

The resulting IMPs no longer satisfy the invariance property (5), but we can use the IMP directly for the prediction of $Y^{e^{\text{test}}}$ .

Remark 4.

To sum up, the class of IMPs constructed under interventions on $Y$ only through the coefficients and shifts will in general imply the invariance property (5), but it is not the case under interventions on $Y$ through the noise variance. For the characterizations of IMP, we focus on sufficient conditions that are relatively simple to evaluate, while they are not necessary conditions in general. One way to find necessary and sufficient conditions for the first/second IMP is through explicit (but tedious) calculations as in the proof of Proposition 2. However, such conditions are hard to evaluate and provide limited insights into our methodology. We thus do not pursue them in this work. When $Y$ is intervened through the coefficients of all its parents, i.e., $PE=PA(Y)$ , a recent work [52] shows that $PA(Y)\subseteq S$ and $PA(Y)\subseteq R$ are also necessary for the IMP under mild assumptions, leading to a novel algorithm for identifying direct causes of $Y$ in the multi-environment setting.

4-B Interventions on both Predictors and Response

To generalize the setting when only $Y$ is intervened to the general setting when $X$ and $Y$ are both intervened, an idea is to merge the setting when only $Y$ is intervened with the one when only $X$ is intervened. The latter setting has been studied in the stabilized regression framework [21]. The following set of predictors is identified (see Definition 3.4 therein),

\displaystyle X^{\text{int}}(Y)=CH^{I}(Y)\cup\bigl{\{}j\in\{1,\ldots,d\}\,|\,\exists i\in CH^{I}(Y)\text{ such that }j\in DE(X_{i})\bigr{\}},

which contains the intervened children of $Y$ (denoted by $CH^{I}(Y)$ ) and the descendants of such children. This useful notion can be defined for each $X_{j}\in\{X_{1},\ldots,X_{d}\}$ , denoted by $X^{\text{int}}(X_{j})$ for each $X_{j}$ . When only $X$ is intervened, the invariance principle (4) holds for $S^{*}=\{1,\ldots,d\}\setminus X^{\text{int}}(Y)$ ; the Markov blanket of $Y$ defined with respect to $X_{S^{*}}$ is called the stable blanket of $Y$ in [21]. In other words, by excluding the predictors in $X^{\text{int}}(Y)$ , the target $Y$ is blocked from the interventions on $X$ when conditioning on $X_{S^{*}}$ . This holds when $Y$ is additionally intervened, as if only $Y$ is intervened given $X_{S^{*}}$ . In order for $S^{*}$ to include as least one child of $Y$ as in Assumption 1, we need the following assumption.

Assumption 2.

When $Y$ is intervened, we assume that $Y$ has at least one child that is not intervened and that child is not a descendant of some intervened child of $Y$ .

Based on the observation above, we identify an important class of IMPs for the general setting in the following Theorem.

Theorem 2.

For the training model $\mathcal{M}^{e}$ without the intervention on the noise variance of $Y$ , the first and second matching properties hold in the following cases.

1.

On the first MP: For $S\subseteq\{1,\ldots,d\}\setminus X^{\text{int}}(Y)$ such that $PA(Y)\subseteq S$ , the first matching property holds.
2.

On the second MP: For each $k\in\{1,\ldots,d\}\setminus\{PE\cup X^{\text{int}}(Y)\}$ , and $R\subseteq\{1,\ldots,d\}\setminus\{k,X^{\text{int}}(X_{k})\cup X^{\text{int}}(Y)\}$ such that $PA(X_{k})\setminus Y\subseteq R$ , the second matching property holds.

Furthermore, if $\lambda_{X}\neq 0$ in the second matching property, then $\phi_{e}(X^{e})=(X_{S}^{e},\,\mathsf{E}_{l,\mathcal{P}_{e}}[X_{k}|X_{R}])^{\top}$ satisfies (5).

Remark 5.

In general, there exist subsets of $S$ and $R$ from Theorem 2 such that the first and second matching properties hold, respectively. In particular, $PA(Y)\subseteq S$ and $PA(X_{k})\subseteq R$ are not necessarily satisfied. Due to the Markov property of SCMs, $Y$ is independent of its ancestors when conditioning on $PA(Y)$ , thus we focus on the IMPs that are more predictive by including $PA(Y)$ in $S$ .

When $Y$ is intervened through its noise variance, recall again that the first matching property always holds for $S=PA(Y)$ by definition. The following proposition follows from Remark 9 in Appendix E.

Proposition 3.

When $Y$ is additionally intervened through the noise variance in Theorem 2, the second matching property in Theorem 2 still holds.

Similar to the argument in the proof of Theorem 1, the class of $\phi$ ’s from Theorem 2 will lead to the same test population loss, as they depend on the same $S$ that is fixed in this setting. Assumption 2 is necessary for $\lambda_{X}\neq 0$ , while sufficient conditions for $\lambda_{X}\neq 0$ can be found similarly as in Proposition 2.

5 Algorithms

For each $e\in\mathcal{E}^{\text{train}}$ , we are given the i.i.d. training data $\boldsymbol{X}_{e}\in\mathbb{R}^{n_{e}\times d},\boldsymbol{Y}_{e}\in\mathbb{R}^{n_{e}}$ , and we observe the i.i.d. test data $\boldsymbol{X}^{\tau}\in\mathbb{R}^{m\times d}$ and aim to predict $\boldsymbol{Y}^{\tau}\in\mathbb{R}^{m}$ . Let $\boldsymbol{X}\in\mathbb{R}^{n\times d}$ with $n:=\sum_{i=1}^{|\mathcal{E}^{\text{train}}|}n_{e}$ denote the pooled data of $\boldsymbol{X}_{e}$ ’s. In this section, we present the implementation of our method starting with the case when $e$ is sampled from a discrete distribution with a finite support. In this setting, we expect to have $n_{e}\gg 1$ for every $e\in\mathcal{E}^{\text{train}}$ in general, thus it is possible to do estimation based on the data from each single environment. The challenging setting of continuous environments will be handled afterward.

5-A Discrete Environments

To implement our method, the main task is to identify the set of IMPs $\mathcal{I}_{\mathcal{M}}$ from the training data. For each tuple $(k,R,S)$ in Algorithm 1, we test the following null hypothesis

\mathcal{H}_{0}\text{ : There exists }\theta\in\mathbb{R}^{d+1}\text{ such that}~{}\eqref{imp_theta}\text{ holds}.

We propose two test procedures.

Test of the Deterministic Relation: Since the IMPs are linear and deterministic (i.e., noiseless), we test whether the residual vector $\boldsymbol{R}\in\mathbb{R}^{n}$ of fitting an IMP on $(k,R,S)$ is a zero vector or not using the test statistics,

T=\frac{1}{n}\boldsymbol{R}^{\top}\boldsymbol{R},

where $\boldsymbol{R}$ is a pooled data vector of $\boldsymbol{R}_{e}$ ’s defined below. To fit an IMP, we first estimate the two LMMSE estimators in (11) using OLS for each environment,

	$\displaystyle\boldsymbol{\hat{L}}_{e,1}$	$\displaystyle:=(\boldsymbol{X}_{e,S}^{\top}\boldsymbol{X}_{e,S})^{-1}\boldsymbol{X}_{e,S}^{\top}\boldsymbol{Y}_{e}\,,$
	$\displaystyle\boldsymbol{\hat{L}}_{e,2}$	$\displaystyle:=(\boldsymbol{X}_{e,R}^{\top}\boldsymbol{X}_{e,R})^{-1}\boldsymbol{X}_{e,R}^{\top}\boldsymbol{X}_{e,k}\,.$

Let $\boldsymbol{\hat{L}}_{1}\in\mathbb{R}^{n}$ and $\boldsymbol{\hat{L}}_{2}\in\mathbb{R}^{n}$ denote the pooled data of $\boldsymbol{\hat{L}}_{e,1}$ ’s and $\boldsymbol{\hat{L}}_{e,2}$ ’s, respectively. It is noteworthy that $\boldsymbol{\hat{L}}_{2}$ only depends on $\boldsymbol{X}$ , thus $\boldsymbol{\hat{L}}_{2}^{\tau}$ for the test data can be computed similarly using $\boldsymbol{X}^{\tau}$ . Next, we estimate the matching parameter using OLS on the pooled data (recall that the matching parameter cannot be identified using the data from a single environment). The OLS estimator of the matching parameter is

\hat{\theta}:=(\hat{\eta}^{\top},\hat{\lambda})^{\top}=([\boldsymbol{X}_{S},\boldsymbol{\hat{L}}_{2}]^{\top}[\boldsymbol{X}_{S},\boldsymbol{\hat{L}}_{2}])^{-1}[\boldsymbol{X}_{S},\boldsymbol{\hat{L}}_{2}]^{\top}\boldsymbol{\hat{L}}_{1}\,.

For each $e\in\mathcal{E}^{\text{train}}$ , we obtain the residual vector of fitting an IMP

\displaystyle\boldsymbol{R}_{e}=\boldsymbol{\hat{L}}_{e,1}-\hat{\lambda}\boldsymbol{\hat{L}}_{e,2}-\boldsymbol{X}_{e,S}\hat{\eta}\,.

2.

Approximate Test of Invariant Residual Distributions: According to the invariance property (6), we test whether the residual when regressing $\boldsymbol{Y}$ on $[\boldsymbol{X}_{S},\boldsymbol{L}_{2}]$ has constant mean and variance. Specifically, we use the t-test and F-test with corrections for multiple hypothesis testing from [8] (see Section 2.1 Method II). The test yields a p-value.

The test statistic from the first procedure and the p-value from the second procedure quantify how likely an IMP holds (i.e., the smaller the more likely), and thus we will refer to either one of them as an IMP score denoted by $s_{\text{IMP}}$ . Let $\hat{\mathcal{I}}=\{(k,R,S):s_{\text{IMP}}(k,R,S)<c_{\text{IMP}}\}$ denote the set of IMPs identified from the training data, where $c_{\text{IMP}}$ is some cutoff parameter. Then, since IMPs are not equally predictive in general, we focus on the most predictive ones by introducing the mean squared prediction error as a prediction score $s_{\text{pred}}$ , and we select the set of IMPs that are more predictive $\widehat{\mathcal{I}_{\text{pred}}}=\{(k,R,S)\in\hat{\mathcal{I}}:s_{\text{pred}}(k,R,S)<c_{\text{pred}}\}$ with some cutoff parameter $c_{\text{pred}}$ . For the second IMP score that is a p-value, the cutoff parameter $c_{\text{IMP}}$ is simply a significance level that is fixed to $0.05$ in this work. For choosing the rest of the cutoff parameters, we follow a bootstrap procedure from [19] with one subtle difference: We sample the same amount of bootstrap samples from each environment rather than sampling over the pool data as in [19] since our procedure involves estimations using the data from each environment.

Algorithm 1 Invariant Prediction using the IMP (discrete)

procedure Identify IMPs from the Training data

for

k\in\{1,\ldots,d\}

S\subseteq\{1,\ldots,d\}

R\subseteq S\setminus k

Compute the IMP score

s_{\text{imp}}

and the prediction

score

s_{\text{pred}}

for

i=(k,R,S)

Regress

\boldsymbol{Y}

[\boldsymbol{X}_{S},\boldsymbol{\hat{L}}_{2}]

to obtain

f_{i}

Identify

\hat{\mathcal{I}}

and

\widehat{\mathcal{I}_{\text{pred}}}

procedure Prediction on the Testing data

\hat{\boldsymbol{Y}}^{\tau}=\frac{1}{|\widehat{\mathcal{I}_{\text{pred}}}|}\sum_{i\in\widehat{\mathcal{I}_{\text{pred}}}}f_{i}(\boldsymbol{X}_{S}^{\tau},\boldsymbol{\hat{L}}_{2}^{\tau})

In practice, there can be spurious IMPs that have extremely small IMP scores but have large prediction scores, e.g., when $Y$ is independent of $X_{S}$ and $X_{k}$ is independent of $X_{R}$ . To this end, we will pre-select $(k,R,S)$ ’s with prediction scores smaller than the median of all the computed prediction scores before identifying $\hat{\mathcal{I}}$ .

If the regression function $f_{i}$ in Algorithm 1 is chosen to be linear, one can use the IMP directly, i.e.,

\boldsymbol{\hat{Y}}_{\text{IMP}}(k,R,S)=[\boldsymbol{X}_{S},\boldsymbol{\hat{L}}_{2}]\hat{\theta},

(22)

which we call the discrete IMP estimator denoted by $\text{IMP}_{\text{d}}$ . To make use of all the IMPs selected in $\widehat{\mathcal{I}_{\text{pred}}}$ , we use an averaging step for the prediction of $\boldsymbol{Y}^{\tau}$ in Algorithm 1.

Now we discuss the computation complexity of our method. For a given graph of size $d\geq 3$ , the number of $(k,R,S)$ ’s with nonempty $(R,S)$ ’s is given by

\sum_{j=1}^{d}\binom{d}{j}\biggl{(}j\cdot(2^{j-1}-1)+(d-j)\cdot(2^{j}-1)\biggr{)}=d\cdot(2\cdot 3^{d-1}-2^{d}),

which follows by first choosing a set $S$ with $j$ elements, and then considering two settings $k\in S$ and $k\not\in S$ since $R$ has to satisfy $R\subseteq S\setminus k$ . This calculation implies that the exhaustive search step in our method is not applicable for relatively large graphs (e.g., $d\geq 15$ ) due to the exponential time complexity. It is noteworthy that in the ICP framework [8], a similar issue occurs due to an exhaustive search over $S\subseteq\{1,\ldots,d\}$ . We discuss several options to alleviate this issue. (1) One can adopt a preprocessing step for feature selection to reduce the dimension $d$ , using Lasso [53] or Boosting [54], as proposed in ICP [8]. (2) When prior information about the graph structure (e.g., the maximum number of parents of the nodes) is available or the graph structure can be estimated, the search space can be reduced. For instance, one can first adopt existing causal discovery methods (e.g., [55]) for graph structure estimation and then apply our methods according to the sufficient conditions in Theorem 1 and 2. (If no IMP is found, one can still perform an exhaustive search over other $(k,R,S)$ ’s). (3) One can also exploit some intrinsic sparsity regarding IMPs, observe that there is only one matched prediction module (recall Definition 3) in the IMP, while the number of prediction modules grows as $d\cdot 3^{d}$ . This sparsity has been studied when only $Y$ is intervened [56], leveraging a variant of Lasso.

5-B Continuous Environments

To model continuous environments, we introduce an environmental variable $U$ that is a continuous random variable with support $\mathcal{U}$ . Apparently, this is a much more challenging setting compared with the discrete environment case, as we only have one training data sample for each $u\in\mathcal{U}$ , making the OLS a poor estimate of $\mathsf{E}_{l,\mathcal{P}_{u}}[Y|X_{S}]$ . Fortunately, it turns out that we can leverage the semi-parametric varying coefficient (SVC) models [57] (see Appendix G) to remedy this issue. In particular, we estimate $\mathsf{E}_{l,\mathcal{P}_{u}}[Y|X_{S}]$ by fitting,

Y=M+\beta^{\top}Z+N\;\text{ with }\;M=\alpha^{\top}(U)W,

(23)

where $N$ is independent of $U$ and the two vectors of predictors $W\in\mathbb{R}^{p}$ (for the varying coefficient) and $Z\in\mathbb{R}^{q}$ (for the invariant coefficients) with $p+q=|S|$ . Since we assume $N\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}U$ , we focus on the settings when $Y$ is not intervened through the noise variance.

Remark 6.

Our estimation procedure for the discrete environments can also be formulated under the SVC model with a discrete random variable $U$ , where we treat all the coefficients as varying coefficients (i.e., $\beta=\boldsymbol{0}$ ), and $\text{IMP}_{\text{d}}$ becomes an estimate of $M$ .

An SVC model over $(Y^{\tau},W^{\tau},Z^{\tau},N^{\tau})$ for the test data can be defined similarly, where $\sigma^{2}=\mathsf{E}[(N^{\tau})^{2}]$ is the population generalization error of the IMP estimator. Observe that the linear SCM $\mathcal{M}^{u}$ in (2) can be viewed as a collection of SVC models parameterized by $U=u$ . Thus the estimation tasks for the linear SCMs from continuous environments can greatly benefit from the existing theories developed for SVC models. More precisely, we employ the following estimate

\widehat{\mathsf{E}_{l,\mathcal{P}_{u}}[Y|X_{S}]}=\hat{M}|_{U=u}+\hat{\beta}^{\top}Z,

where the profile least-squares estimation of $\beta$ and $M$ proposed in [57] can be found in Appendix G. Similarly, $\mathsf{E}_{l,\mathcal{P}_{u}}[X_{k}|X_{R}]$ can be estimated by fitting another semi-parametric varying coefficient model

\displaystyle V=M_{V}+\beta_{V}^{\top}Z_{V}+N_{V}\;\text{ with }\;M_{V}=\alpha_{V}^{\top}(U)W,

(24)

where $V$ denotes any $X_{k}$ , and $X_{R}$ is divided into $Z_{V}\in\mathbb{R}^{r}$ and $W$ .

It is noteworthy that the two SVC models share the same set of predictors with varying coefficients, which we explain below. A challenge for fitting such models is that the vector of predictors with varying coefficients, namely $W$ , needs to be known. For continuous environments, we focus on discovering IMPs that can be decomposed into the first and second matching properties. Thus, since the causal function captures the predictors with varying coefficients, the first and second matching properties imply that the vector $W$ is simply $X_{PE}$ , i.e., the parents of $Y$ with varying coefficients in $\mathcal{M}^{u}$ , for both models.

Based on this observation and Theorem 1, we replace the exhaustive search over $(k,R,S)$ in Algorithm 1 by a search over $(P,k,R,S)$ according to the conditions in Theorem 1 with $PE=P$ . That is, we choose $(P,k,R,S)$ from

	$\displaystyle P\subseteq\{1,\ldots,d\},\quad k\in\{1,\ldots,d\}\setminus P,$
	$\displaystyle P\subseteq S\subseteq\{1,\ldots d\},\quad P\subseteq R\subseteq S\setminus k,$

such that $W=X_{P}$ , $Z=X_{S\setminus P}$ , $V=X_{k}$ , and $Z_{V}=X_{R\setminus P}$ .

Unlike $\text{IMP}_{\text{d}}$ in (22), we make use of the fact that $\beta$ is invariant and propose the continuous IMP estimator denoted by $\text{IMP}_{\text{c}}$ as follows

\boldsymbol{\hat{Y}}_{\text{IMP}}(P,k,R,S)=[\boldsymbol{W},\boldsymbol{\hat{M}}_{V}]\hat{w}+\boldsymbol{Z}\hat{\beta},

(25)

with the matching parameter $w\in\mathbb{R}^{p+1}$ estimated by

\hat{w}=([\boldsymbol{W},\boldsymbol{\hat{M}}_{V}]^{\top}[\boldsymbol{W},\boldsymbol{\hat{M}}_{V}])^{-1}[\boldsymbol{W},\boldsymbol{\hat{M}}_{V}]^{\top}\boldsymbol{\hat{M}}.

(26)

The data matrices for the two models (e.g., $\boldsymbol{W}\in\mathbb{R}^{n\times p}$ ) can be defined accordingly and we provide the details in Appendix I. In this case, the residual vector for the first IMP score is given by

\boldsymbol{R}=\boldsymbol{M}-[\boldsymbol{W},\boldsymbol{\hat{M}}_{V}]\hat{w}.

Note the second IMP score is not applicable for continuous environments due to the small sample size in each environment, thus we focus on the first IMP score.

Remark 7.

The $\text{IMP}_{\text{d}}$ and $\text{IMP}_{\text{c}}$ are similar in spirit, the $\text{IMP}_{\text{c}}$ relies on the first and second matching properties for identifying $W$ (so we can reuse $Z$ in (25)), whereas $\text{IMP}_{\text{d}}$ directly tests the IMP since the estimation process treats all the coefficients as varying coefficients (also see the proof of Corollary 4).

6 Asymptotic Generalization Error

In this section, we provide the asymptotic generalization errors (as $n,m\to\infty$ ) of the $\text{IMP}_{\text{c}}$ and $\text{IMP}_{\text{d}}$ estimators for $(k,R,S)$ ’s that satisfy IMPs, i.e., $(k,R,S)\in\mathcal{I}_{\mathcal{M}}$ . Recall that $\sigma^{2}=\mathsf{E}[(N^{\tau})^{2}]$ is the population generalization error of the IMP estimators. Due to the estimations on both training and test data, the asymptotic generalization error can be decomposed to the error terms depending on the training data size $n$ and the test data size $m$ as follows. Let $c_{n}=\left\{\frac{\log(1/h)}{nh}\right\}^{1/2}+h^{2}$ , where $h$ is the kernel bandwidth (see Appendix G for details).

Theorem 3.

For any $(k,R,S)\in\mathcal{I}_{\mathcal{M}}$ , under the technical assumptions in Appendix H, the asymptotic generalization error of the $\text{IMP}_{\text{c}}$ estimator is given by

\frac{1}{m}\sum_{i=1}^{m}(\hat{Y}^{\tau}_{i}-Y_{i}^{\tau})^{2}=\sigma^{2}+O_{p}(c_{n}\vee n^{-1/2})+O_{p}(c_{m}\vee m^{-1/2}).

The following corollary considers the setting when the amount of unlabeled training and test data grows in a higher order than that of labels in the training data. The generalization error due to the estimation on the test data disappears.

Corollary 3.

Given i.i.d. training data of size $n$ with $0<l_{n}<n$ labels and test data of size $m$ , if $\max(\frac{l_{n}}{n},\frac{l_{n}}{m})\to 0$ as $\min(m,n)\to\infty$ , under the technical assumptions in Appendix H, the asymptotic generalization error of the $\text{IMP}_{\text{c}}$ estimator is given by

\frac{1}{m}\sum_{i=1}^{m}(\hat{Y}^{\tau}_{i}-Y_{i}^{\tau})^{2}=\sigma^{2}+O_{p}(c_{l_{n}}\vee l_{n}^{-1/2}).

The setting of discrete environments can be viewed as a special case of continuous environments, where the error term $c_{n}\vee n^{-1/2}$ due to the kernel estimation procedure is replaced by an error term from multiple OLS estimations.

Corollary 4.

For any $(k,R,S)\in\mathcal{I}_{\mathcal{M}}$ , under the technical assumptions in Appendix H, the asymptotic generalization error of the $\text{IMP}_{\text{d}}$ estimator is given by

\frac{1}{m}\sum_{i=1}^{m}(\hat{Y}^{\tau}_{i}-Y_{i}^{\tau})^{2}=\sigma^{2}+O_{p}(a_{n})+O_{p}(a_{m}),

where $a_{n}=(\min_{e\in\mathcal{E}^{\text{train}}}n_{e})^{-1/2}$ and $a_{m}=m^{-1/2}$ .

This asymptotic generalization error heavily depends on the environment with the smallest sample size, which also supports the fact that $\text{IMP}_{\text{d}}$ should not be employed for continuous environment settings.

7 Experiments

The prediction performance is measured by the mean residual sum of squares (RSS) on the test environments. We compare our method with several baseline methods: Ordinary Least Squares (OLS), stabilized regression (SR) [21], anchor regression (AR) [22]. We have also compared with domain invariant projection (DIP) [29], conditional invariance penalty (CIP) [34], conditional invariant residual matching (CIRM) [28], and invariant risk minimization (IRM) [23]; it turns out that the empirical performance of these methods is not as competitive as the other baselines in our experimental settings, thus we do not report them below.

The two IMP scores lead to two versions of our algorithm, and we refer to the first one as IMP and the second one as $\text{IMP}_{\text{inv}}$ (since it tests the invariance of the noise mean and variance). We focus on linear functions $f_{i}$ ’s in Algorithm 1, namely, we use the IMP estimators. For the profile likelihood estimation, we adopt the Epanechnikov kernel $k(u)=0.75\max(1-u^{2},0)$ with the bandwidth fixed to be $0.1$ . We test DIP, CIP, CIRM, and their variants provided in [28] with the default parameters. For the anchor regression, we use a $5$ -fold cross-validation procedure to select the hyper-parameter $\gamma$ from $\{0,0.05,0.1,\ldots,0.5\}$ . The significance levels are fixed to be $0.05$ for all methods. We randomly simulate $500$ data sets for each experiment, if not mentioned otherwise.

7-A Discrete environments

First, we generate linear SCMs $\mathcal{M}^{e}$ ’s without interventions. For each $e\in\mathcal{E}^{\text{train}}=\{1,\ldots,5\}$ or $e\in\mathcal{E}^{\text{test}}=\{6,\ldots,10\}$ , we randomly generate a linear SCM with 9 variables as follows. The graph $\mathcal{G}(\mathcal{M}^{e})$ is specified by a lower triangular matrix of i.i.d. $\mathrm{Bernoulli}(1/2)$ random variables. The response $Y$ is randomly selected from the $9$ variables and we require that $Y$ has a least one parent and one child in $\mathcal{G}(\mathcal{M}^{e})$ . When $X$ and $Y$ are both intervened, we randomly choose a child of $Y$ to not be intervened. For each linear SCM, the non-zero coefficients are sampled from $\mathrm{Unif}[-1.5,-0.5]\cup[0.5,1.5]$ and the noise variables are standard normal. For each training or test environment, we simulate i.i.d. data of sample size $300$ .

7-A1 Interventions on $X$

Since the baseline methods have been examined extensively under shift interventions, we focus on shift interventions on $X$ for comparison. The general interventions on $X$ will be considered in Section 7-C. Specifically, for each training environment, we randomly selected $4$ predictors to be intervened through shifts sampled from $\mathrm{Unif}[-2,2]$ . For each test environment, the shifts are sampled from $\mathrm{Unif}[-10,10]$ . In Fig. 6, $\text{IMP}_{\text{inv}}$ performs similarly to SR since they share a similar idea when only $X$ is intervened. Due to the averaging steps of IMP, $\text{IMP}_{\text{inv}}$ , and SR, they have smaller variances compared with OLS and AR (a similar result has been reported in [19]).

7-A2 Interventions on $Y$

We consider the response $Y$ to be intervened through both the coefficients and shifts. We randomly select $n_{p}\sim\mathrm{Unif}\{1,\ldots,|PA(Y)|\}$ of parents of $Y$ to have varying coefficients. For each training environment, we add perturbation terms sampled from $\mathrm{Unif}[-2,2]$ to the original coefficients. For each test environment, the perturbations are sampled from $\mathrm{Unif}[-10,10]$ . The shift intervention on $Y$ is the same as the shift interventions on $X$ in Section 7-A1. In this setting, since none of the baseline methods allow interventions on $Y$ through the coefficients, they cannot even improve upon OLS. In Fig. 7, IMP performs slightly better than $\text{IMP}_{\text{inv}}$ , which may be due to the fact that the IMP method aims to find all possible IMPs, but the $\text{IMP}_{\text{inv}}$ only looks for IMPs that imply invariance. To further examine our test procedures for identifying IMPs, we check if the sufficient conditions in Theorem 1 are satisfied for the estimated IMPs $(\hat{k},\hat{R},\hat{S})$ ’s, summarized in Table 1. We observe that the conditions on $(k,S)$ are satisfied for the majority of cases, while the condition on $R$ holds with a noticeably lower empirical probability. Note that our sufficient conditions might be conservative. Moreover, due to the randomly selected coefficients, the intervention strength can be quite weak in some cases, making it challenging to distinguish true IMPs and non-IMPs. Fortunately, the averaging step helps to mitigate inaccurate predictions made by non-IMPs.

	$\hat{P}(\hat{k}\not\in PE)$	$\hat{P}(PE\subseteq\hat{R})$	$\hat{P}(PE\subseteq\hat{S})$
IMP	0.9610 (0.0796)	0.7503 (0.2511)	0.9376 (0.1079)
$\text{IMP}_{\text{inv}}$	0.9603 (0.0850)	0.7060 (0.2617)	0.9373 (0.1090)

TABLE 1: Empirical results on the estimated IMP vs. conditions in Theorem 1 in terms of empirical probability and standard deviation (in parenthesis): Only

Y

is intervened.

7-A3 Interventions on both $X$ and $Y$

The setting of interventions on both $X$ and $Y$ is simply a combination of the two settings above. In this challenging setting, our method outperforms the baselines by a large margin. Since the sufficient conditions in Theorem 2 are not necessary for some settings as discussed in Remark 10, we verify slightly relaxed conditions according to Remark 10. The estimation task is more challenging in this setting compared with the setting with interventions only on $X$ or $Y$ . The results in Table 2 indicate that a large percentage of estimated IMPs conform to the conditions. In Fig. 8, we present an illustration of the difference between the predictions made by the estimated IMPs and other $(k,R,S)$ ’s that are rejected by our test procedures (which we call the estimated non-IMPs). For the estimated IMP, the gap between the prediction error of the training environments and that of the test environments is relatively small. However, the gap can be huge for the majority of the estimated non-IMPs.

	$\hat{P}(\hat{k}\not\in PE\cup X^{\text{int}})$	$\hat{P}(A)$	$\hat{P}(B)$
IMP	0.8767 (0.2749)	0.7170 (0.3515)	0.7003 (0.4412)
$\text{IMP}_{\text{inv}}$	0.8558 (0.3144)	0.6571 (0.4085)	0.6853 (0.4460)

TABLE 2: Empirical results on the estimated IMP vs. conditions in Theorem 2 in terms of empirical probability and standard deviation (in parenthesis): Both

X

and

Y

are intervened. Events

A

and

B

are defined by

A=\{PE\subseteq\hat{R}\}\cap\{(X^{\text{int}}(Y)\cup X^{\text{int}}(X_{\hat{k}}))\cap\hat{R}=\varnothing\}

and

B=\{PE\subseteq\hat{S}\}\cap\{X^{\text{int}}(Y)\cap\hat{S}=\varnothing\}

, respectively.

7-B Interventions on both $X$ and $Y$ (continuous)

In this more challenging setting, we only compare with OLS, since AR only considers shift interventions while the other baselines are not developed for continuous environments, and also recall that $\text{IMP}_{\text{inv}}$ is proposed for discrete environments. First, we define $\{U_{1},\ldots,U_{800}\}$ sampled from $\mathrm{Unif}[0,1]$ and $\{U_{1}^{\tau},\ldots,U_{800}^{\tau}\}$ sampled from $\mathrm{Unif}[1,2]$ as the training and test environments, respectively. Similar to the setting of discrete environments, we randomly generate the linear SCMs without interventions first and then add interventions to the model. Due to the high computational complexity, we focus on graphs with $5$ nodes where $2$ predictors are intervened. The interventions on the coefficients and shift interventions are defined by adding a perturbation term $a\sin(2\pi wU_{i})$ , where $w$ is sampled from $\mathrm{Unif}[0.5,2]$ . The parameter $a$ is fixed to be $2$ for the training environments (i.e., the same range as for the discrete case) and $5$ for the test environments. Since the parameter space is much smaller than in the previous experiments, we only generated $100$ data sets. Overall, the performance of our IMP algorithm is similar to that in the discrete setting.

7-C Robustness

In the previous experiments, we only consider shift interventions on $X$ and interventions on $Y$ other than the noise variance. In this experiment, we consider an extreme case when only a child of $Y$ with the highest causal ordering is not intervened (i.e., Assumption 2 is satisfied), all other variables are intervened through every parameter. Specifically, a shift intervention or an intervention on the coefficient is defined by adding a perturbation term to the original parameter. The perturbation term is sampled from $\mathrm{Unif}[-2,2]$ for the training data, and from $\mathrm{Unif}[-5,5]$ for the test data. The intervened noise variances are sampled from $\mathrm{Unif}[0.75,1.25]$ for each training environment, and from $\mathrm{Unif}[0.5,1.5]$ for each test environment. To test how sensitive our method is with respect to Assumption 2 in this challenging setting, we gradually add interventions to the child of $Y$ that is not intervened, where the shifts and coefficient interventions are sampled from $\mathrm{Unif}[-2\lambda,2\lambda]$ and $\mathrm{Unif}[-5\lambda,5\lambda]$ for the training and test environments, respectively. The intervened noise variances are sample from $\mathrm{Unif}[1-0.25\lambda,1+0.25\lambda]$ for training, and from $\mathrm{Unif}[1-0.5\lambda,1+0.5\lambda]$ for testing. The parameter $\lambda\in[0,1]$ controls the intervention strength. Due to the results in Section 7-A3, it would not be informative to compare with the baseline methods, so we focus on our IMP method. Note that our $\text{IMP}_{\text{inv}}$ is also not included since IMPs will not imply invariance in this case (see Remark 4). Overall, as shown in Fig. 11, the median of the mean RSS is not too sensitive with respect to mild interventions on the child that is not intervened, but the variance increases rapidly.

7-D COVID Data Set

For observational data, there is often no ground truth regarding the set of variables that are intervened. If the data is not collected under a carefully designed experimental setting (e.g., the gene perturbation experiments from [8]), it is reasonable to assume that every variable is more or less intervened, especially the response. To examine our methods under real-world distribution shifts, we consider a COVID data set [58] collected at $3142$ US counties from January 22, 2020, to June 10, 2021. The data set consists of 46 predictive features that are relevant to the number of COVID cases. Some of the features are treated as fixed or time-invariant, e.g., population density, age distribution, and education level, and others are temporal features, e.g., temperature, social distancing grade, and the number of COVID tests performed.

There are plenty of prediction tasks that can be potentially designed for this data set. In our experiment, we focus on the prediction of the number of COVID cases ⁴⁴4The number of COVID cases is measured in thousands. using the $12$ temporal features for the time interval from March 1, 2020, to September 30, 2020 ( $214$ days). For different regions of the US, features such as population density can differ greatly, making it challenging for a prediction model to perform well on a diverse set of regions. Specifically, we take $5$ major cities from the Western U.S. (Los Angeles, San Francisco, San Diego, Seattle, and Phoenix) as $5$ training environments, and we examine our methods and the baselines regarding $8$ cities/counties that are mostly from the Eastern U.S. (Baltimore, Boston, Philadelphia, New York County, Queens County ⁵⁵5Only the data from Manhattan (New York County) and Queens are available for New York City in this data set., Houston, Miami, and Chicago). According to the total number of COVID cases in the considered time interval, we partition (using $50$ k as the threshold) the cities/counties for testing into two groups in Table 3, which makes the results more informative as most of the methods perform worse on cities/counties with higher total cases. Since no invariance in terms of mean and variance can be identified by $\text{IMP}_{\text{inv}}$ , we focus on the IMP method. Overall, IMP and CIP show similar performance for the lower total cases, while the performance of CIP degrades more over the higher total cases compared with IMP. Since real-world interventions are often beyond shift interventions, AR does not improve upon OLS for this data set. It is noteworthy that no stable sets are found by SR at the significance level of $0.05$ , implying that the invariance principle is likely to be violated. IRM shows less competitive performance compared with other methods, thus its results are not presented in Table 3.

Cities/Counties (total cases)	IMP	OLS	AR	SR	CIP	DIP	CIRM
Baltimore (18k)	0.1292	0.0937	0.0379	0.0978	0.0184	0.9199	0.3503
Boston (26k)	0.0238	0.0503	0.0351	0.0543	0.0192	0.5820	0.4605
Philadelphia (32k)	0.1532	0.1191	0.0701	0.0970	0.0353	1.0407	0.4217
NY County (34k)	0.1611	0.4088	0.5606	0.3449	0.3537	2.4055	2.6372
average mean RSS	0.1168	0.1680	0.1759	0.1485	0.1066	1.2370	0.9674
Queens County (73k)	0.1126	0.1020	0.0882	0.1109	0.0673	0.9389	0.9647
Houston (143k)	1.0884	1.2042	1.2200	1.1673	1.3705	1.7802	1.4982
Chicago (146k)	0.3687	0.4489	0.4026	0.4099	0.4562	0.6328	0.3499
Miami (171k)	0.4534	0.4621	0.5170	0.4025	0.6011	1.4465	0.8004
average mean RSS	0.5058	0.5543	0.5548	0.5226	0.6238	1.1996	0.9033

TABLE 3: Experiment results of the COVID data set.

8 Discussion

To deal with general interventions on the response, we introduce the IMP that allows for an alternative form of invariance when the traditional invariance principle fails. We provide explicit characterizations of the IMP under different intervention settings and provide asymptotic generalization error analysis for both the discrete environment and the more challenging continuous environment settings, supported by our empirical studies.

Our work is motivated by allowing general interventions on the response, and it can be extended into several directions. First, we observe a connection between the IMP and self-supervised learning schemes (e.g., [59] from natural language processing). Specifically, the prediction of $X_{k}$ using $X_{R}$ corresponds to the pretext (or pretrained) task, and the linear relation in the IMP solves the downstream prediction task. We believe that methods developed through causal modeling including our IMP method will bring new opportunities for self-supervised learning. We have discussed several heuristic approaches to alleviate the computation complexity issue; it would be worthwhile to analyze such methods and provide theoretical guarantees (e.g., exploiting sparsity beyond the setting in [56]). Our asymptotic analysis of the generalization error is developed given that IMPs are correctly identified. It would be interesting to further analyze the consistency of the set of IMPs $\hat{\mathcal{I}}$ identified by the proposed algorithms or new efficient algorithms.

References

[1] J. Quinonero-Candela, M. Sugiyama, A. Schwaighofer, and N. D. Lawrence, Dataset shift in machine learning. MIT Press, 2008.
[2] K. Weiss, T. M. Khoshgoftaar, and D. Wang, “A survey of transfer learning,” Journal of Big Data, vol. 3, no. 1, pp. 1–40, 2016.
[3] G. Csurka, “Domain adaptation for visual applications: A comprehensive survey,” arXiv preprint arXiv:1702.05374, 2017.
[4] J. Pearl, Causality. Cambridge University Press, 2009.
[5] J. Peters, D. Janzing, and B. Schölkopf, Elements of causal inference: foundations and learning algorithms. The MIT Press, 2017.
[6] B. Schölkopf, D. Janzing, J. Peters, E. Sgouritsa, K. Zhang, and J. Mooij, “On causal and anticausal learning,” in 29th International Conference on Machine Learning (ICML 2012). International Machine Learning Society, 2012, pp. 1255–1262.
[7] K. Zhang, B. Schölkopf, K. Muandet, and Z. Wang, “Domain adaptation under target and conditional shift,” in International Conference on Machine Learning. PMLR, 2013, pp. 819–827.
[8] J. Peters, P. Bühlmann, and N. Meinshausen, “Causal inference by using invariant prediction: identification and confidence intervals,” Journal of the Royal Statistical Society. Series B (Statistical Methodology), pp. 947–1012, 2016.
[9] M. Rojas-Carulla, B. Schölkopf, R. Turner, and J. Peters, “Invariant models for causal transfer learning,” The Journal of Machine Learning Research, vol. 19, no. 1, pp. 1309–1342, 2018.
[10] C. Heinze-Deml and N. Meinshausen, “Conditional variance penalties and domain shift robustness,” arXiv preprint arXiv:1710.11469, 2017.
[11] P. Bühlmann, “Invariance, causality and robustness,” Statistical Science, vol. 35, no. 3, pp. 404–426, 2020.
[12] A. Ghassami, N. Kiyavash, B. Huang, and K. Zhang, “Multi-domain causal structure learning in linear systems,” Advances in neural information processing systems, vol. 31, 2018.
[13] A. Ghassami, S. Salehkaleybar, N. Kiyavash, and K. Zhang, “Learning causal structures using regression invariance,” Advances in Neural Information Processing Systems, vol. 30, 2017.
[14] T. Haavelmo, “The probability approach in econometrics,” Econometrica: Journal of the Econometric Society, pp. iii–115, 1944.
[15] J. Aldrich, “Autonomy,” Oxford Economic Papers, vol. 41, no. 1, pp. 15–34, 1989.
[16] G. W. Imbens and D. B. Rubin, Causal inference in statistics, social, and biomedical sciences. Cambridge University Press, 2015.
[17] N. Meinshausen, “Causality from a distributional robustness point of view,” in 2018 IEEE Data Science Workshop (DSW). IEEE, 2018, pp. 6–10.
[18] C. Heinze-Deml, J. Peters, and N. Meinshausen, “Invariant causal prediction for nonlinear models,” Journal of Causal Inference, vol. 6, no. 2, 2018.
[19] N. Pfister, P. Bühlmann, and J. Peters, “Invariant causal prediction for sequential data,” Journal of the American Statistical Association, vol. 114, no. 527, pp. 1264–1276, 2019.
[20] S. Magliacane, T. Van Ommen, T. Claassen, S. Bongers, P. Versteeg, and J. M. Mooij, “Domain adaptation by using causal inference to predict invariant conditional distributions,” Advances in Neural Information Processing Systems, vol. 31, 2018.
[21] N. Pfister, E. G. Williams, J. Peters, R. Aebersold, and P. Bühlmann, “Stabilizing variable selection and regression,” The Annals of Applied Statistics, vol. 15, no. 3, pp. 1220–1246, 2021.
[22] D. Rothenhäusler, N. Meinshausen, P. Bühlmann, and J. Peters, “Anchor regression: Heterogeneous data meet causality,” Journal of the Royal Statistical Society: Series B (Statistical Methodology), vol. 83, no. 2, pp. 215–246, 2021.
[23] M. Arjovsky, L. Bottou, I. Gulrajani, and D. Lopez-Paz, “Invariant risk minimization,” arXiv preprint arXiv:1907.02893, 2019.
[24] E. Rosenfeld, P. Ravikumar, and A. Risteski, “The risks of invariant risk minimization,” in International Conference on Learning Representations, vol. 9, 2021.
[25] P. Kamath, A. Tangella, D. Sutherland, and N. Srebro, “Does invariant risk minimization capture invariance?” in International Conference on Artificial Intelligence and Statistics. PMLR, 2021, pp. 4069–4077.
[26] R. Christiansen, N. Pfister, M. E. Jakobsen, N. Gnecco, and J. Peters, “A causal framework for distribution generalization,” IEEE Transactions on Pattern Analysis and Machine Intelligence, 2021.
[27] K. Du and Y. Xiang, “An invariant matching property for distribution generalization under intervened response,” in European Conference on Signal Processing (EUSIPCO) 2022 arXiv:2205.09162, 2022.
[28] Y. Chen and P. Bühlmann, “Domain adaptation under structural causal models,” Journal of Machine Learning Research, vol. 22, pp. 1–80, 2021.
[29] M. Baktashmotlagh, M. T. Harandi, B. C. Lovell, and M. Salzmann, “Unsupervised domain adaptation by domain invariant projection,” in Proceedings of the IEEE International Conference on Computer Vision, 2013, pp. 769–776.
[30] Y. Ganin, E. Ustinova, H. Ajakan, P. Germain, H. Larochelle, F. Laviolette, M. Marchand, and V. Lempitsky, “Domain-adversarial training of neural networks,” The Journal of Machine Learning Research, vol. 17, no. 1, pp. 2096–2030, 2016.
[31] Y. Li, M. Murias, S. Major, G. Dawson, and D. Carlson, “On target shift in adversarial domain adaptation,” in The 22nd International Conference on Artificial Intelligence and Statistics. PMLR, 2019, pp. 616–625.
[32] R. Tachet des Combes, H. Zhao, Y.-X. Wang, and G. J. Gordon, “Domain adaptation with conditional distribution matching and generalized label shift,” Advances in Neural Information Processing Systems, vol. 33, pp. 19 276–19 289, 2020.
[33] M. Gong, K. Zhang, T. Liu, D. Tao, C. Glymour, and B. Schölkopf, “Domain adaptation with conditional transferable components,” in International Conference on Machine Learning. PMLR, 2016, pp. 2839–2848.
[34] C. Heinze-Deml and N. Meinshausen, “Conditional variance penalties and domain shift robustness,” Machine Learning, vol. 110, no. 2, pp. 303–348, 2021.
[35] A. Storkey, “When training and test sets are different: characterizing learning transfer,” Dataset shift in machine learning, vol. 30, pp. 3–28, 2009.
[36] M. Sugiyama and M. Kawanabe, Machine learning in non-stationary environments: Introduction to covariate shift adaptation, 2012.
[37] Z. Lipton, Y.-X. Wang, and A. Smola, “Detecting and correcting for label shift with black box predictors,” in International Conference on Machine Learning. PMLR, 2018, pp. 3122–3130.
[38] K. Azizzadenesheli, A. Liu, F. Yang, and A. Anandkumar, “Regularized learning for domain adaptation under label shifts,” in International Conference on Learning Representations, 2018.
[39] S. Garg, Y. Wu, S. Balakrishnan, and Z. Lipton, “A unified view of label shift estimation,” Advances in Neural Information Processing Systems, vol. 33, pp. 3290–3300, 2020.
[40] S. Ben-David, J. Blitzer, K. Crammer, and F. Pereira, “Analysis of representations for domain adaptation,” Advances in Neural Information Processing Systems, vol. 19, 2006.
[41] J. A. Bagnell, “Robust supervised learning,” in AAAI, 2005, pp. 714–719.
[42] Z. Hu and L. J. Hong, “Kullback-Leibler divergence constrained distributionally robust optimization,” Available at Optimization Online, pp. 1695–1724, 2013.
[43] A. Sinha, H. Namkoong, and J. Duchi, “Certifying some distributional robustness with principled adversarial training,” in International Conference on Learning Representations, 2018.
[44] R. Gao, X. Chen, and A. J. Kleywegt, “Distributional robustness and regularization in statistical learning,” arXiv preprint arXiv:1712.06050, 2017.
[45] J. Lee and M. Raginsky, “Minimax statistical learning with Wasserstein distances,” Advances in Neural Information Processing Systems, vol. 31, 2018.
[46] J. C. Duchi and H. Namkoong, “Learning models with uniform performance via distributionally robust optimization,” The Annals of Statistics, vol. 49, no. 3, pp. 1378–1406, 2021.
[47] I. Goodfellow, P. McDaniel, and N. Papernot, “Making machine learning robust against adversarial inputs,” Communications of the ACM, vol. 61, no. 7, pp. 56–66, 2018.
[48] A. Raghunathan, J. Steinhardt, and P. S. Liang, “Semidefinite relaxations for certifying robustness to adversarial examples,” Advances in Neural Information Processing Systems, vol. 31, 2018.
[49] J. Pearl and E. Bareinboim, “External validity: From do-calculus to transportability across populations,” in Probabilistic and Causal Inference: The Works of Judea Pearl, 2022, pp. 451–482.
[50] K. Du and Y. Xiang, “Causal inference from slowly varying nonstationary processes,” arXiv preprint arXiv:2012.13025, 2020.
[51] Y. Xiang, J. Ding, and V. Tarokh, “Estimation of the evolutionary spectra with application to stationarity test,” IEEE Transactions on Signal Processing, vol. 67, no. 5, pp. 1353–1365, 2019.
[52] K. Du, Y. Xiang, and I. Soloveychik, “Identifying direct causes using intervened target variable,” arXiv preprint arXiv:2307.07736, 2023.
[53] R. Tibshirani, “Regression shrinkage and selection via the lasso,” Journal of the Royal Statistical Society: Series B (Methodological), vol. 58, no. 1, pp. 267–288, 1996.
[54] P. Bartlett, Y. Freund, W. S. Lee, and R. E. Schapire, “Boosting the margin: A new explanation for the effectiveness of voting methods,” The annals of statistics, vol. 26, no. 5, pp. 1651–1686, 1998.
[55] S. Shimizu, P. O. Hoyer, A. Hyvärinen, A. Kerminen, and M. Jordan, “A linear non-gaussian acyclic model for causal discovery.” Journal of Machine Learning Research, vol. 7, no. 10, 2006.
[56] K. Du and Y. Xiang, “Generalized invariant matching property via lasso,” in ICASSP 2023-2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2023, pp. 1–5.
[57] J. Fan and T. Huang, “Profile likelihood inferences on semiparametric varying-coefficient partially linear models,” Bernoulli, vol. 11, no. 6, pp. 1031–1057, 2005.
[58] A. Haratian, H. Fazelinia, Z. Maleki, P. Ramazi, H. Wang, M. A. Lewis, R. Greiner, and D. Wishart, “Dataset of COVID-19 outbreak and potential predictive features in the USA,” Data in Brief, vol. 38, p. 107360, 2021.
[59] J. Devlin, M.-W. Chang, K. Lee, and K. Toutanova, “BERT: Pre-training of deep bidirectional transformers for language understanding,” arXiv preprint arXiv:1810.04805, 2018.
[60] M. S. Bartlett, “An inverse matrix adjustment arising in discriminant analysis,” The Annals of Mathematical Statistics, vol. 22, no. 1, pp. 107–111, 1951.
[61] Y.-p. Mack and B. W. Silverman, “Weak and strong uniform consistency of kernel regression estimates,” Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete, vol. 61, no. 3, pp. 405–415, 1982.

Appendix A Proof of Proposition 1

Let $\mathcal{E}^{\text{train}}=\{e_{1},\ldots,e_{n}\}$ . For a tuple $(k,S,R)$ that satisfies the IMP, let

\displaystyle\hat{\boldsymbol{Y}}=(\mathsf{E}[Y^{e_{1}}|X_{S}^{e_{1}}],\ldots,\mathsf{E}[Y^{e_{n}}|X_{S}^{e_{n}}])^{\top},

\boldsymbol{\tilde{X}}=\begin{bmatrix}(X^{e_{1}})^{\top}&\mathsf{E}_{l}[X_{k}^{e_{1}}|X_{R}^{e_{1}}]\\ \vdots&\vdots\\ (X^{e_{n}})^{\top}&\mathsf{E}_{l}[X_{k}^{e_{n}}|X_{R}^{e_{n}}]\end{bmatrix}:=\begin{bmatrix}\boldsymbol{X}&v\end{bmatrix},

where the rows of $\boldsymbol{\tilde{X}}$ are independent. According to (12), we have $\hat{\boldsymbol{Y}}=\boldsymbol{\tilde{X}}\theta$ . Then, if $\mathsf{E}[\boldsymbol{\tilde{X}}^{\top}\boldsymbol{\tilde{X}}]$ is invertible, we have

\displaystyle(\mathsf{E}[\boldsymbol{\tilde{X}}^{\top}\boldsymbol{\tilde{X}}])^{-1}\mathsf{E}[\boldsymbol{\tilde{X}}^{\top}\hat{\boldsymbol{Y}}]=(\mathsf{E}[\boldsymbol{\tilde{X}}^{\top}\boldsymbol{\tilde{X}}])^{-1}\mathsf{E}[\boldsymbol{\tilde{X}}^{\top}\boldsymbol{\tilde{X}}]\theta=\theta.

(27)

Now, we prove the invertibility. Observe that

\mathsf{E}[\boldsymbol{\tilde{X}}^{\top}\boldsymbol{\tilde{X}}]=\begin{bmatrix}\mathsf{E}[\boldsymbol{X}^{\top}\boldsymbol{X}]&\mathsf{E}[\boldsymbol{X}^{\top}v]\\ \mathsf{E}[v^{\top}\boldsymbol{X}]&\mathsf{E}[v^{\top}v]\end{bmatrix},

where $\mathsf{E}[\boldsymbol{X}^{\top}\boldsymbol{X}]=\sum_{i=1}^{n}\mathsf{E}[X^{e_{i}}(X^{e_{i}})^{\top}]$ is inveritble since it is a sum of positive-definite matrices. Then, $\mathsf{E}[\boldsymbol{\tilde{X}}^{\top}\boldsymbol{\tilde{X}}]$ is invertible if and only if

\mathsf{E}[v^{\top}v]-\mathsf{E}[v^{\top}\boldsymbol{X}]\mathsf{E}^{-1}[\boldsymbol{X}^{\top}\boldsymbol{X}]\mathsf{E}[\boldsymbol{X}^{\top}v]\neq 0.

This is equivalent to

\mathsf{E}[(v-\boldsymbol{X}\beta)^{\top}(v-\boldsymbol{X}\beta)]\neq 0,

where $\beta:=\mathsf{E}^{-1}[\boldsymbol{X}^{\top}\boldsymbol{X}]\mathsf{E}[\boldsymbol{X}^{\top}v]$ . This is true since there is no $b\in\mathbb{R}^{d}$ such that $v=\boldsymbol{X}b$ almost surely by our assumption in (15). Therefore, $\theta$ is uniquely determined by (27).

Appendix B Problem Formulation Using an environmental random variable

By introducing an environmental random variable $E\in\mathcal{E}^{\text{all}}$ , we define a mixture of $\mathcal{M}^{e}$ ’s, $e\in\mathcal{E}^{\text{all}}$ , as follows,

	$\displaystyle\mathcal{M}:$	$\displaystyle X=\gamma(E)Y+B(E)X+\varepsilon_{X}(E)$
	$\displaystyle\mathcal{M}:$	$\displaystyle Y=\beta^{\top}(E)X+\varepsilon_{Y}(E),$

where $E$ in a root node in $\mathcal{G}(\mathcal{M})$ and the noise variables are jointly independent given $E$ . We do not specify the distribution of $E$ but assume that $E$ does not have a degenerate distribution. Under this formulation, the invariance property (5) is equivalent to

Y\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}E\,|\,\phi(E,X),

(28)

and the invariant, first, and second matching properties can be equivalently written as

$\displaystyle\mathsf{E}_{l}[Y\|X_{S},E=e]$	$\displaystyle=\lambda\mathsf{E}_{l}[X_{k}\|X_{R},E=e]+\eta^{\top}X,$	(29)
$\displaystyle\mathsf{E}_{l}[Y\|X_{S},E=e]$	$\displaystyle=\lambda_{Y}\mathsf{E}[Y\|X_{PA(Y)},E=e]+\eta_{Y}^{\top}X,$	(30)
$\displaystyle\mathsf{E}_{l}[X_{k}\|X_{R},E=e]$	$\displaystyle=\lambda_{X}\mathsf{E}[Y\|X_{PA(Y)},E=e]+\eta_{X}^{\top}X,$	(31)

for $e\in\mathcal{E}^{\text{all}}$ . As a special case of $\mathcal{M}$ , we define a mixture of $\mathcal{M}^{e,1}$ ’s as

	$\displaystyle\mathcal{M}^{1}:$	$\displaystyle X=\gamma Y+BX+\varepsilon_{X}$			(32)
	$\displaystyle\mathcal{M}^{1}:$	$\displaystyle Y=(\alpha(E)+\beta)^{\top}X+\varepsilon_{Y},$			(33)

where only $Y$ is intervened through the coefficients.

The proofs of Theorems 1, 2 and Corollaries 1, 2 will be presented under this formulation.

Appendix C Proof of Theorem 1

The problem formulation and necessary notation for the proof are introduced in Appendix B. For the first part, let $Z(E)=\alpha^{\top}(E)X$ denote an additional node in the acyclic graph $\mathcal{G}(\mathcal{M}^{1})$ , then the assignment of $Y$ in (33) becomes

Y=Z(E)+\beta^{\top}X+\varepsilon_{Y},

(34)

where $E$ is no longer a parent of $Y$ . Note that $Z(E)=\mathsf{E}[Y|X_{PA(Y)},E]-\beta^{\top}X$ as we assume both $X$ and $Y$ have zero means. Since $E$ is a root node, observe that $E$ and $Y$ can be d-connected through only two types of paths as follows,

1.

$E\to Z(E)\to Y$ ,
2.

$E\to Z(E)\leftarrow X_{i}\to\cdots\to X_{l}\leftarrow\cdots\leftarrow Y$ ,

where $i\in PE$ , and $X_{i}\to\cdots\to X_{l}\leftarrow\cdots\leftarrow Y$ is a V-structure for some $l\in DE(i)\cap DE(Y)$ . Note that the second type of path does not exist if Assumption 1 is not satisfied or $CH(i)\setminus Y$ is empty.

We start by showing that the d-separation $Y\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}_{\mathcal{G}}E\,|\,\{Z(E),X_{S}\}$ holds given $PE\subseteq S$ . First, the first path is immediately blocked by $Z(E)$ . Second, for any $s\in CH(i)\setminus Y$ , the path $E\to Z(E)\leftarrow X_{i}\to X_{s}$ is blocked by $\{Z(E),X_{i}\}$ . Thus, the second path is blocked given $Z(E)$ and $X_{S}$ . According to the Markov property of SCMs [4], the d-separation $Y\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}_{\mathcal{G}}E\,|\,\{Z(E),X_{S}\}$ implies

Y\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}E\,|\,\{Z(E),X_{S}\}.

(35)

Now, we prove that the above conditional independence implies the first matching property (30).

By definition, the LMMSE estimators only rely on the (finite) first two moments of the variables. Thus we start with the case when $(X,Y)|_{E=e}$ is jointly Gaussian, for each $e\in\mathcal{E}^{\text{all}}$ . First we have

	$\displaystyle\mathsf{E}_{l}[Y\|X_{S},E=e]$	$\displaystyle\overset{(a)}{=}\mathsf{E}[Y\|X_{S},E=e]$
		$\displaystyle\overset{(b)}{=}\mathsf{E}[Y\|Z(e),X_{S},E=e],$

where $(a)$ follows from the Gaussian assumption on $(X,Y)|_{E=e}$ and $(b)$ from the fact that $Z(e)$ is a function of $\{X_{S},E=e\}$ given our assumption that $PE\subseteq S$ . This implies that

\displaystyle\mathsf{E}_{l}[Y|X_{S},E]=\mathsf{E}[Y|Z(E),X_{S},E]\overset{(a)}{=}\mathsf{E}[Y|Z(E),X_{S}],

where $(a)$ follows from the conditional independence relation (35). Using the Gaussian assumption again, we have $\mathsf{E}[Y|Z(e),X_{S}]=\mathsf{E}_{l}[Y|Z(e),X_{S}]$ . Thus putting all the pieces together, we obtain

	$\displaystyle\mathsf{E}_{l}[Y\|X_{S},E=e]$	$\displaystyle=\mathsf{E}_{l}[Y\|Z(e),X_{S}]$		(36)
		$\displaystyle=aZ(e)+b^{\top}X_{S},$		(37)

where $a\in\mathbb{R}$ and $b\in\mathbb{R}^{|S|}$ that are not functions of $E$ . When $(X,Y)|_{E=e}$ is non-Gaussian, one can replace it with Gaussian random variables with the matching first and second moments. Then the same argument leading to (37) still holds. Then, the first matching property follows from the fact that $Z(E)=\mathsf{E}[Y|X_{PA(Y)},E]-\beta^{\top}X$ .

Similarly, for the second part, we first show the following d-separation $X_{k}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}_{\mathcal{G}}\,E\,|\,\{Z(E),X_{R}\}$ . Observe that $X_{k}$ and $E$ can be d-connected through two types of paths as follows,

1.

$E\to Z(E)\to Y\cdots X_{k}$ ,
2.

$E\to Z(E)\leftarrow X_{i}\cdots X_{k}$ ,

with $i\in\{j\in\{1,\ldots,d\}:\alpha_{j}(E)\neq 0\}$ , where $Y\cdots X_{k}$ denotes any directed path between $Y$ and $X_{k}$ , and similarly for $X_{i}\cdots X_{k}$ . The two types of paths are immediately blocked by $\{Z(E),X_{R}\}$ under our assumption that $PE\subseteq R$ . We thus have the following when $X|_{E=e}$ is jointly Gaussian,

$\displaystyle\mathsf{E}_{l}[X_{k}\|X_{R},E=e]$	$\displaystyle=\mathsf{E}_{l}[X_{k}\|Z(e),X_{R},E=e]$
	$\displaystyle=\mathsf{E}_{l}[X_{k}\|Z(e),X_{R}]$
	$\displaystyle=cZ(e)+d^{\top}X_{R}$	(38)

for some $c\in\mathbb{R}$ and $d\in\mathbb{R}^{|R|}$ do not depend on $E$ . The non-Gaussian cases can be handled in the same way as before, and thus the second property follows again from $Z(E)=\mathsf{E}[Y|X_{PA(Y)},E]-\beta^{\top}X$ . Observe that we have $\lambda_{X}=c$ in the second matching property and Assumption 1 is a necessary condition for $\lambda_{X}\neq 0$ .

Finally, given $R\subseteq S$ , we have that (38) with $c\neq 0$ (i.e., $\lambda_{X}\neq 0$ ) provides a one-to-one mapping between $\{Z(E),X_{S}\}$ and $\{\mathsf{E}_{l}[X_{k}|X_{R},E],X_{S}\}$ . Therefore, the conditional independence (35) is equivalent to

Y\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}E\,\big{|}\,\{\mathsf{E}_{l}[X_{k}|X_{R},E],\,X_{S}\}.

(39)

This implies that, using our previous notation, $\phi_{e}(X^{e})=(X_{S}^{e},\,\mathsf{E}_{l,\mathcal{P}_{e}}[X_{k}|X_{R}])^{\top}$ satisfies the invariance property (5), which implies⁶⁶6We use the shorthand $\phi_{e}(X)$ for $\phi_{e}(X^{e})$ since the expectation is with respect to $\mathcal{P}_{e}$ . $\mathsf{E}_{\mathcal{P}_{e}}[Y|X_{S}]=\mathsf{E}_{\mathcal{P}_{e}}[Y|\phi_{e}(X)]$ . Effectively, $\mathsf{E}_{\mathcal{P}_{e}}[Y|\phi_{e}(X)]$ serves as a representation of $\mathsf{E}_{\mathcal{P}_{e}}[Y|X_{S}]$ that is invariant. Note that $\Phi$ is not empty when $\lambda_{X}\neq 0$ holds with $R\subseteq S$ . This implies that any $\phi\in\Phi$ with $S=\{1,\ldots,d\}$ minimizes $\mathcal{L}_{\text{test}}(f_{\phi})$ , since the optimality of $\phi\in\Phi$ only relies on the corresponding $S$ .

Remark 8.

When $\varepsilon_{Y}$ is replace by $\varepsilon_{Y}(E)$ , we consider the following two cases.

1.

The mean of $\varepsilon_{Y}^{e}$ is a function of $e$ , and its variance is a constant.
2.

The variance of $\varepsilon_{Y}^{e}$ is a function of $E$ , and its mean can be either a function of $e$ or a constant.

For the first case, we can introduce $X_{d+1}:=1$ that is a parent of $Y$ with $\alpha_{d+1}:=\mathsf{E}[\varepsilon_{Y}(E)|E]$ . Additionally, $X_{d+1}$ is a parent of every $X_{j}$ such that $\varepsilon_{X,j}$ has a non-zero mean. Specifically, the coefficient of $X_{d+1}$ in the assignment of $X_{j}$ will be $\mathsf{E}[\varepsilon_{X,j}]$ . Thus, the problem reduces to the setting when $\varepsilon_{X}^{e}$ and $\varepsilon_{Y}^{e}$ have zero means, which has been proved in Theorem 1. For the second case, however, the varying variance of $\varepsilon_{Y}(E)$ cannot be separated as the mean, thus $E$ is always a parent of $Y$ and the path $E\to Y$ cannot be blocked by $Z$ or any $X_{S}\subseteq\{X_{1},\ldots,X_{d}\}$ , i.e., the proof of Theorem 1 breaks down.

Appendix D Proof of Proposition 2

The following lemma is a slight extension of Lemma 3.6 from [21], where we consider linear models with dependent noise variables rather than linear SCMs considered in [21].

Lemma 1.

Consider $V\in\mathbb{R}$ and $X=(X_{1},\ldots,X_{p})^{\top}\in\mathbb{R}^{p}$ satisfying a linear model,

\begin{bmatrix}V\\ X\end{bmatrix}=\begin{bmatrix}0&h^{\top}\\ g&A\end{bmatrix}\begin{bmatrix}V\\ X\end{bmatrix}+\begin{bmatrix}e_{V}\\ e_{X}\end{bmatrix}:=\tilde{A}\begin{bmatrix}V\\ X\end{bmatrix}+e,

(40)

where $e_{V}\in\mathbb{R}$ , $g,h,e_{X}\in\mathbb{R}^{p}$ , and $A\in\mathbb{R}^{p\times p}$ . Assume that $I-\tilde{A}$ is invertible, the population OLS estimator when regressing $V$ on $X$ is given by

\displaystyle\theta^{\text{ols}}

\displaystyle=\left[c-g^{\top}\tilde{\Sigma}^{-1}\left(I-c\sigma^{2}gg^{\top}\tilde{\Sigma}^{-1}\right)v\right]h+\left(I-A^{\top}\right)\left[c\sigma^{2}\tilde{\Sigma}^{-1}g+\tilde{\Sigma}^{-1}\left(I-c\sigma^{2}gg^{\top}\tilde{\Sigma}^{-1}\right)v\right],

where $c:=\left(1+\sigma^{2}g^{\top}\tilde{\Sigma}^{-1}g\right)^{-1}$ , $\tilde{\Sigma}:=\Sigma+vg^{\top}+gv^{\top}$ , and

\mathop{\rm Cov}\nolimits(e,e):=\begin{bmatrix}\sigma^{2}&v^{\top}\\ v&\Sigma\end{bmatrix},

with $\sigma^{2}:=\mathop{\rm Var}\nolimits(e_{V})$ , $\Sigma:=\mathop{\rm Cov}\nolimits(e_{X},e_{X})$ , and $v:=\mathop{\rm Cov}\nolimits(e_{X},e_{V})$ .

Proof:

Let $P=[0_{1\times d},I_{d\times d}]$ , we continue from Equation A.7 in [21] as follows,

	$\displaystyle\theta^{\text{ols}}$	$\displaystyle=\mathop{\rm Cov}\nolimits^{-1}(X,X)\mathop{\rm Cov}\nolimits(X,V)$
		$\displaystyle=\left\{P(I-\tilde{A})^{-1}\begin{bmatrix}\sigma^{2}&v^{\top}\\ v&\Sigma\end{bmatrix}\left[(I-\tilde{A})^{-1}\right]^{\top}P^{\top}\right\}^{-1}P(I-\tilde{A})^{-1}\begin{bmatrix}\sigma^{2}\\ v\end{bmatrix}+h,$

and it has been shown that

P(I-\tilde{A})^{-1}=[w,M],

for $M=(I-A-gh^{\top})^{-1}$ and $w=Mg$ . Then,

	$\displaystyle\theta^{\text{ols}}$	$\displaystyle=\left[\sigma^{2}ww^{\top}+wv^{\top}M^{\top}+Mvw^{\top}+M\Sigma M^{\top}\right]^{-1}\left(\sigma^{2}w+Mv\right)+h$
		$\displaystyle:=\left[\sigma^{2}ww^{\top}+M\tilde{\Sigma}M^{\top}\right]^{-1}\left(\sigma^{2}w+Mv\right)+h,$

where $\tilde{\Sigma}:=\Sigma+vg^{\top}+gv^{\top}$ and it was computed in [21] using the Sherman-Morrison formula [60] that

		$\displaystyle\left[\sigma^{2}ww^{\top}+M\tilde{\Sigma}M^{\top}\right]^{-1}$
	$\displaystyle=$	$\displaystyle(I-A-gh^{\top})^{\top}\tilde{\Sigma}^{-1}\left(I-c\sigma^{2}gg^{\top}\tilde{\Sigma}^{-1}\right)M^{-1},$

with $c:=\left(1+\sigma^{2}g^{\top}\tilde{\Sigma}^{-1}g\right)^{-1}$ . Then, some simple algebra leads to

	$\displaystyle\theta^{\text{ols}}$	$\displaystyle=h+c\sigma^{2}(I-A-gh^{\top})^{\top}\tilde{\Sigma}^{-1}g+(I-A-gh^{\top})^{\top}\tilde{\Sigma}^{-1}\left(I-c\sigma^{2}gg^{\top}\tilde{\Sigma}^{-1}\right)v$
		$\displaystyle=\left[c-g^{\top}\tilde{\Sigma}^{-1}\left(I-c\sigma^{2}gg^{\top}\tilde{\Sigma}^{-1}\right)v\right]h+\left(I-A^{\top}\right)\left[c\sigma^{2}\tilde{\Sigma}^{-1}g+\tilde{\Sigma}^{-1}\left(I-c\sigma^{2}gg^{\top}\tilde{\Sigma}^{-1}\right)v\right].$

$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\blacksquare$

D-A Proof of Proposition 2

For each tuple $(k,R,S),k\in\{j\in MB(Y):\alpha_{j}^{e}=0\}$ . $R=-k:=\{1,\ldots,d\}\setminus k$ , $S=\{1,\ldots,d\}$ , we prove that $\lambda_{X}\neq 0$ . Equivalently, we prove that $\mathsf{E}_{l}[X_{k}|X_{-k},E=e]$ is a non-constant function of $e$ .

First, recall that, when the target variables $Y$ is unobserved, the relations between the predictors in $\mathcal{M}^{e,1}$ , are described by

X=\left(\gamma(\beta+\alpha^{e})^{\top}+B\right)X+\gamma\varepsilon_{Y}+\varepsilon_{X}.

Now, we rewrite the above equation in the same form as the linear model (40) as follows,

\displaystyle\begin{bmatrix}X_{k}\\ X_{-k}\end{bmatrix}=\begin{bmatrix}0&\gamma_{k}(\beta+\alpha^{e})_{-k}^{\top}+B_{k,-k}\\ \gamma_{-k}\beta_{k}+B_{-k,k}&\gamma_{-k}(\beta+\alpha^{e})_{-k}^{\top}+B_{-k,-k}\end{bmatrix}\begin{bmatrix}X_{k}\\ X_{-k}\end{bmatrix}+\varepsilon_{X}+\gamma\varepsilon_{Y},

where the top-left element of the coefficient matrix is zero, i.e., $\gamma_{k}(\beta_{k}+\alpha_{k}^{e})+B_{k,k}=0$ since $\alpha_{k}^{e}=0$ (by assumption), $B_{k,k}=0$ (due to acyclicity), and $\gamma_{k}\beta_{k}=0$ (since $X_{k}$ cannot be both a child and a parent of $Y$ ). Now, by Lemma 1, the population OLS estimator when regressing $X_{k}$ on $X_{-k}$ given $E=e$ is

	$\displaystyle\theta^{\text{ols},k}(e)$	$\displaystyle=\left[c-g^{\top}\tilde{\Sigma}^{-1}\left(I-c\sigma^{2}gg^{\top}\tilde{\Sigma}^{-1}\right)v\right]h+\left(I-A^{\top}\right)\left[c\sigma^{2}\tilde{\Sigma}^{-1}g+\tilde{\Sigma}^{-1}\left(I-c\sigma^{2}gg^{\top}\tilde{\Sigma}^{-1}\right)v\right]$
		$\displaystyle:=ah+(I-A^{\top})b,$		(41)

where $h=\gamma_{k}(\beta+\alpha^{e})_{-k}+B_{k,-k}^{\top}$ , $v=\gamma_{k}\sigma_{Y}^{2}\gamma_{-k}$ , $g=\beta_{k}\gamma_{-k}+B_{-k,k}$ , $A=\gamma_{-k}(\beta+\alpha^{e})_{-k}^{\top}+B_{-k,-k}$ , $\Sigma=\mathop{\rm Cov}\nolimits(N_{X_{-k}},N_{X_{-k}})$ , $\tilde{\Sigma}=\Sigma+vg^{\top}+gv^{\top}$ , and $c=\left(1+\sigma^{2}g^{\top}\tilde{\Sigma}^{-1}g\right)^{-1}$ . Note that $a\in\mathbb{R}$ and $b\in\mathbb{R}^{(d-1)}$ are not functions of $e$ and

a=g^{\top}b+c(1+\sigma^{2}g^{\top}\tilde{\Sigma}^{-1}g)=g^{\top}b+1.

(42)

$k\in CH(Y)$ : First, observe that $\beta_{k}=0$ implies $g=B_{-k,k}$ . By plugging $h$ and $A$ into (41), we obtain

\displaystyle\theta^{\text{ols},k}(e)

\displaystyle=(a-\gamma^{\top}_{-k}b)\alpha_{-k}^{e}+a\gamma_{k}\beta_{-k}+(I-\beta_{-k}\gamma_{-k}^{\top}-B_{-k,-k}^{\top})b,

where $\alpha_{-k}^{e}$ in the first term is non-vanishing only if

a-\gamma^{\top}_{-k}b=1+(B_{-k,k}-\gamma_{-k})^{\top}b\neq 0,

where we use (42).

$k\in PA(Y)$ : Observe that $v=0_{(d-1)}$ and $\gamma_{k}=0$ , then

	$\displaystyle\theta^{\text{ols},k}(e)$	$\displaystyle=c\alpha^{e}+c\beta+c\sigma^{2}\left(I-A^{\top}\right)\Sigma^{-1}g$
		$\displaystyle=c(1+\sigma^{2}\beta_{k}\gamma_{-k}^{\top}\Sigma^{-1}\gamma_{-k}+\sigma^{2}\gamma_{-k}^{\top}\Sigma^{-1}B_{-k,k})\alpha_{-k}^{e}+c\beta+c\sigma^{2}(I-\beta_{-k}\gamma_{-k}^{\top}-B_{-k,-k}^{\top})\Sigma^{-1}g,$

where the first term is not vanishing only if

1+\sigma^{2}\beta_{k}\gamma_{-k}^{\top}\Sigma^{-1}\gamma_{-k}+\sigma^{2}\gamma_{-k}^{\top}\Sigma^{-1}B_{-k,k}\neq 0.

3.

$k\in\{j:\exists i\in CH(X_{j})\text{ such that }i\in CH(Y)\}$ : Again, we have $v=0_{(d-1)\times 1}$ and $\gamma_{k}=0$ . Additionally, we have $\beta_{k}=0$ , then $\alpha_{-k}^{e}$ in $\theta^{\text{ols},k}(e)$ is non-vanishing only if

$1+\sigma^{2}\gamma_{-k}^{\top}\Sigma^{-1}B_{-k,k}\neq 0.$

Appendix E Proof of Corollary 2

First, we extract the assignments of $(X_{k},X_{R})$ ’s from (20) with $E=e$ , where $k\in CH(Y)$ and $k\not\in DE(X_{i})$ for any $i\in CH(Y)\setminus k$ and $R=\{1,\ldots,d\}\setminus DE(Y)$ as follows,

\displaystyle\begin{bmatrix}X_{k}\\ X_{R}\end{bmatrix}

\displaystyle=\begin{bmatrix}0&\gamma_{k}(\beta+\alpha^{e})_{R}^{\top}+B_{k,R}\\ 0_{|R|\times 1}&B_{R,R}\end{bmatrix}\begin{bmatrix}X_{k}\\ X_{R}\end{bmatrix}+\begin{bmatrix}\varepsilon_{X,k}+\gamma_{k}\varepsilon_{Y}^{e}\\ \varepsilon_{X,R}\end{bmatrix},

where the zeros in the coefficient matrix are due to the fact that all the descendants of $X_{k}$ are excluded from $X_{R}$ since $DE(X_{k})\subseteq DE(Y)$ and $R=\{1,\ldots,d\}\setminus DE(Y)$ . Note that $X_{k}$ is a child of $Y$ that is not a descendant of any other child of $Y$ , thus any removed node $j\in\{1,\ldots,d\}\setminus\{k,R\}$ can not be the parent of any remaining node $i\in\{k,R\}$ .

Now, using Lemma 1, the population OLS estimator when regressing $X_{k}$ on $X_{R}$ given $E=e$ is

\theta^{\text{ols}}(e)=\gamma_{k}(\beta+\alpha^{e})_{R}+B_{k,R}^{\top}.

This implies

\mathsf{E}_{l}[X_{k}|X_{R},E=e]=\gamma_{k}(\beta+\alpha^{e})^{\top}X+B_{k,\cdot}X,

where we use the fact that $\beta_{j}=\alpha^{e}_{j}=B_{k,j}=0$ for $j\not\in R$ (recall that $\{1,\ldots,d\}\setminus R$ does not contain either the parents of $X_{k}$ or parents of $Y$ ). Therefore,

\mathsf{E}_{l}[X_{k}|X_{R},E=e]=\gamma_{k}\mathsf{E}[Y|X_{PA(Y)},E=e]+B_{k,\cdot}X,

which is the second matching property.

Remark 9.

Observe that $\theta^{\text{ols}}(e)$ does not depend on $B_{R,R}$ or the covariance of $\varepsilon_{X,R}$ , thus the second matching property holds even if every $X_{j}\in X_{R}$ is intervened.

Appendix F Proof of Theorem 2

For the first part, we only need to prove that $Y$ and $E$ are d-connected only through the arrow $E\to Y$ when conditioning on $X_{S}$ , $S\subseteq\{1,\ldots,d\}\setminus X^{\text{int}}(Y)$ such that $PA(Y)\subseteq S$ . The rest is simply the proof of the first part of Theorem 1. Since $E$ is a root node, $Y$ and $E$ can only be d-connected through the two types of paths,

1.

$E\to\cdots\to Y$ ,
2.

$Y\to\cdots\to X_{j}\leftarrow\cdots\leftarrow E$ ,

where $j\in DE(Y)\cap DE(E)$ , and $E\to\cdots\to Y$ denotes any directed path from $E$ to $Y$ . For the first type of paths, since $PA(Y)\subseteq S$ , the only unblocked path when conditioning on $X_{S}$ is $E\to Y$ . For the second type of paths that are V-structures, we have $j\not\in S$ since $X^{\text{int}}(Y)$ contains all the intervened children of $Y$ and the descendants of such children, implying that the V-structures are always blocked when conditioning on $X_{S}$ . Thus, the only unblocked path between $E$ and $Y$ is $E\to Y$ when conditioning on $X_{S}$ .

For the second part, similarly, for $k\in\{1,\ldots,d\}\setminus\{PE\cup X^{\text{int}}(Y)\}$ and $R\subseteq\{1,\ldots,d\}\setminus\{k,X^{\text{int}}(X_{k})\cup X^{\text{int}}(Y)\}$ such that $PA(X_{k})\setminus Y\subseteq R$ , we prove that $X_{k}$ is d-connected with $E$ only through paths that contain the subpath $E\to Y$ when conditioning on $X_{R}$ (i.e., the second type of path below). Following the same idea as in the first part, $Y$ and $E$ are d-connected only through the arrow $E\to Y$ when conditioning on $X_{R}$ since $X^{\text{int}}(Y)\cap R=\varnothing$ and $PA(Y)\subseteq R$ . Then, since $X_{k}$ is not intervened (i.e., $E\not\in PA(X_{k})$ ), the variables $E$ and $X_{k}$ can only be d-connected through three types of paths as follows,

1.

$E\to\dots\to X_{k}$ that does not contain $Y$ ,
2.

$E\to Y\to\cdots\to X_{k}$ ,
3.

$E\to\cdots\to X_{i}\leftarrow\cdots\leftarrow X_{k}$ ,

where $i\in R$ and we use the fact that $k\not\in X^{\text{int}}(Y)$ (i.e., the node $k$ is not an intervened child of $Y$ or a descendant of an intervened child of $Y$ ). The first type of path is blocked immediately due to $PA(X_{k})\setminus Y\subseteq R$ .

To handle the third type of paths, we will need two technical results. (I) If $i\in R$ , then $j\in R$ for any $j\in PA(i)$ such that $j\neq k$ . This can be proved by contradiction. If $j\not\in R$ (i.e., $j\in X^{\text{int}}(Y)\cup X^{\text{int}}(X_{k})$ ), then we have $i\in X^{\text{int}}(Y)\cup X^{\text{int}}(X_{k})$ by the definition of $X^{\text{int}}(Y)$ and $X^{\text{int}}(Y)$ , i.e., $i\not\in R$ . (II) For $i\in R$ , observe that $X_{i}$ can only be a child of $X_{k}$ , otherwise the path $X_{i}\leftarrow X_{j}\leftarrow\cdots\leftarrow X_{k}$ is blocked by $X_{j}\in PA(X_{i})$ when conditioning on $X_{R}$ , since $i\in R$ implies $j\in R$ .

Now, we proved that the third type of paths are always blocked when conditioning on $X_{R}$ , by focusing on $E\to\cdots\to X_{i}$ . For any $X_{i}\in CH(X_{k})$ , the subpath $E\to\dots\to X_{i}$ cannot be $E\to X_{i}$ , since $i\in R$ implies that $X_{i}$ is not an intervened child of $X_{k}$ . Finally, the path $E\to\cdots\to X_{l}\to X_{i}$ is blocked by $X_{l}\in PA(X_{i})$ when conditioning on $X_{R}$ , since $i\in R$ implies that $l\in R$ .

Remark 10.

We require $PA(Y)$ to be included $S$ in order to block every path from intervened ancestors to $Y$ , but this is not necessary when the ancestors are not intervened. While it is important to include $PE$ in $S$ to handle interventions on $Y$ . Also, including $PE$ will block paths from intervened ancestors that pass through variables in $PE$ . Similar arguments apply to $R$ as well.

Appendix G Semi-parametric Varying Coefficient Models and Profile Least-Squares Estimation

First, we introduce the semi-parametric varying coefficient model following most of the notation in [57]. Consider $Y\in\mathbb{R}$ , $U\in\mathcal{U}$ , and two vectors of predictors $W=(W_{1},\ldots,W_{p})^{\top}$ and $Z=(Z_{1},\ldots,Z_{q})^{\top}$ such that $Z_{j}$ ’s have invariant coefficients, a semi-parametric varying coefficient model over $(U,Y,W,Z)$ is defined by

Y=M+\beta^{\top}Z+N,\quad M=\alpha^{\top}(U)W,

(43)

where $N$ is independent of $(U,W,Z)$ .

We briefly introduce the profile least-squares estimator of $\beta$ proposed in [57]. Denote $n$ i.i.d. samples of $(U,Y,Z,N)$ as $\boldsymbol{U}=(U_{1},\ldots,U_{n})^{\top}$ , $\boldsymbol{Y}=(Y_{1},\ldots,Y_{n})^{\top}$ , $\boldsymbol{W}=(\boldsymbol{W_{1}},\ldots,\boldsymbol{W_{n}})^{\top}$ , $\boldsymbol{W}_{i}=(W_{i1},\ldots,W_{ip})^{\top}$ , $\boldsymbol{Z}=(\boldsymbol{Z}_{1},\ldots,\boldsymbol{Z}_{n})^{\top}$ , $\boldsymbol{Z}_{i}=(Z_{i1},\ldots,Z_{iq})^{\top}$ , and $\boldsymbol{N}=(N_{1},\ldots,N_{n})^{\top}$ . We thus have $\boldsymbol{M}=(M_{1},\ldots,M_{n})^{\top}$ with $M_{i}=\alpha^{\top}(U_{i})\boldsymbol{W}_{i}$ . Let $\boldsymbol{K}_{u}=\text{diag}(K_{h}(U_{1}-u),\ldots,K_{h}(U_{n}-u))$ for some kernel function $K_{h}(\cdot)=K(\cdot/h)/h$ with bandwidth $h$ , and

\boldsymbol{\tilde{W}}_{u}=\begin{bmatrix}\boldsymbol{W_{1}}^{\top}&\frac{U_{1}-u}{h}\boldsymbol{W_{1}}^{\top}\\ \vdots&\vdots\\ \boldsymbol{W_{n}}^{\top}&\frac{U_{n}-u}{h}\boldsymbol{W_{n}}^{\top}\end{bmatrix}\in\mathbb{R}^{n\times 2p}.

For $u$ in a neighborhood of $u_{0}\in\mathcal{U}$ , assume that each function $\alpha_{i}(u)$ in (43) can be approximated locally by the first-order Taylor expansion $\alpha_{i}(u)\approx\alpha_{i}+b_{i}(u-u_{0})$ . Then the varying coefficients $\alpha(U_{1}),\ldots,\alpha(U_{n})$ can be estimated by solving the weighted local least squares problem

\displaystyle\min_{\{a_{i},b_{i}\}}\sum_{k=1}^{n}\left[\boldsymbol{Y}_{k}-\sum_{j=1}^{q}\beta_{j}Z_{kj}-\sum_{i=1}^{p}(a_{i}+b_{i}(U_{k}-u))W_{ki}\right]^{2}\cdot K_{h}(U_{k}-u),

(44)

which has the solution

\displaystyle[\hat{a}_{1}(u),\ldots,\hat{a}_{p}(u),h\hat{b}_{1}(u),\ldots,h\hat{b}_{p}(u)]=(\boldsymbol{\tilde{X}}_{u}^{\top}\boldsymbol{K}_{u}\boldsymbol{\tilde{X}}_{u})^{-1}\boldsymbol{\tilde{X}}_{u}^{\top}\boldsymbol{K}_{u}(\boldsymbol{Y}-\boldsymbol{Z}\beta),

for $u\in\{U_{1},\ldots,U_{n}\}$ . Then, each variable $\boldsymbol{M}_{i}$ can be estimated by $\boldsymbol{\hat{M}}_{i}=\sum_{j=1}^{p}\hat{a}_{j}(u)\boldsymbol{X}_{ij}$ , and thus the vector $\boldsymbol{M}$ can be estimated by

	$\displaystyle\boldsymbol{\hat{M}}(\beta)=$
	$\displaystyle\begin{bmatrix}[\boldsymbol{W}_{1}^{\top}\quad 0_{1\times p}]\{\boldsymbol{\tilde{W}}_{u_{1}}^{\top}\boldsymbol{K}_{u_{1}}\boldsymbol{\tilde{W}}_{u_{1}}\}^{-1}\boldsymbol{\tilde{W}}_{u_{1}}^{\top}\boldsymbol{K}_{u_{1}}\\ \vdots\\ [\boldsymbol{W}_{n}^{\top}\quad 0_{1\times p}]\{\boldsymbol{\tilde{W}}_{u_{n}}^{\top}\boldsymbol{K}_{u_{n}}\boldsymbol{\tilde{W}}_{u_{n}}\}^{-1}\boldsymbol{\tilde{W}}_{u_{n}}^{\top}\boldsymbol{K}_{u_{n}}\\ \end{bmatrix}(\boldsymbol{Y}-\boldsymbol{Z}\beta)$
	$\displaystyle:=A(\boldsymbol{Y}-\boldsymbol{Z}\beta),$		(45)

which depends on the unknown parameter $\beta$ that will be estimated below. Substituting $\boldsymbol{\hat{M}}(\beta)$ into the vector form of (43), we obtain $(I-A)\boldsymbol{Y}=(I-A)\boldsymbol{Z}\beta+\boldsymbol{N}$ . The profile least-squares estimator of $\beta$ is given by

\displaystyle\hat{\beta}

\displaystyle=\{\boldsymbol{Z}^{\top}(I-A)^{\top}(I-A)\boldsymbol{Z}\}^{-1}\cdot\boldsymbol{Z}^{\top}(I-A)^{\top}(I-A)\boldsymbol{Y}.

(46)

Finally, by replacing $\beta$ in (45) with $\hat{\beta}$ , the final form of the estimator for $\boldsymbol{M}$ is given by

\displaystyle\boldsymbol{\hat{M}}

\displaystyle=A(\boldsymbol{Y}-\boldsymbol{Z}\hat{\beta}).

(47)

Appendix H Technical Lemmas for the Proof of Theorem 3

The semi-parametric varying coefficient model and profile least-squares estimation are introduced in Appendix G. First, we present some technical assumptions and two technical lemmas from [57]. Let $c_{n}^{\prime}=\left\{\frac{\log(1/h)}{nh}\right\}^{1/2}$ and $c_{n}=c_{n}^{\prime}+h^{2}$ .

1.

$U$ has a bounded support $\mathcal{U}$ and has density function $f(\cdot)$ that is Lipschitz continuous and bounded away from $0$ .
2.

For each $U_{i}$ , the matrix $\mathsf{E}[W^{\top}W|U_{i}]$ is non-singular, and the matrices $\mathsf{E}[W^{\top}W|U_{i}]$ , $\left(\mathsf{E}[W^{\top}W|U_{i}]\right)^{-1}$ , and $\mathsf{E}[WZ^{\top}|U_{i}]$ are Lipschitz continuous.
3.

$\alpha_{1}(u),\ldots,\alpha_{p}(u)$ have continuous second derivatives.
4.

$K(\cdot)$ is a symmetric density function.
5.

There exists an $s>2$ such that $\mathsf{E}[||X||^{2s}]<\infty$ , $\mathsf{E}[||Z_{Y}||^{2s}]<\infty$ , and $\mathsf{E}[||Z_{V}||^{2s}]<\infty$ and there exists $\varepsilon<2-s^{-1}$ such that $n^{2\varepsilon-1}h\to\infty$ .
6.

The bandwidth $h$ satisfies $nh^{8}\to 0$ and $nh^{2}/(\log(n))^{2}\to\infty$ .

Lemma 2 ([61]).

Let $(U_{1},Y_{1}),\ldots,(U_{n},Y_{n})$ be i.i.d. random vectors in $\mathbb{R}^{2}$ . Assume that $\mathsf{E}[|Y|^{s}]<\infty$ and $\sup_{x}\int|y|^{s}f(u,y)dy<\infty$ , where $f(u,y)$ is the density of $(U,Y)$ . Let $K$ be a bounded positive function with a bounded support that satisfies a Lipschitz condition. Given that $n^{2\varepsilon-1}h\to\infty$ for some $\varepsilon<1-s^{-1}$ , then

\displaystyle\sup_{u}\left|\frac{1}{n}\sum_{i=1}^{n}K_{h}(U_{i}-u)Y_{i}-\mathsf{E}[K_{h}(U_{i}-u)Y_{i}]\right|=O_{p}(c_{n}^{\prime}).

Lemma 3 ([57]).

Under assumptions $1)\sim 6)$ , $\sqrt{n}(\hat{\beta}-\beta)\xrightarrow{d}\mathcal{N}(0,\Sigma)$ as $n\to\infty$ , where $\Sigma=\mathop{\rm Var}\nolimits(N)\cdot C^{-1}$ with

\displaystyle C=\mathsf{E}[ZZ^{\top}]-\mathsf{E}\big{[}\mathsf{E}[ZW^{\top}|U]\mathsf{E}[WW^{\top}|U]\mathsf{E}[WZ^{\top}|U]\big{]}.

In the following lemma, we provide the rates of several quantities needed for the proof of Theorem 3.

Lemma 4.

Under the same assumptions as in Lemma 3,

	$\displaystyle R_{1}$	$\displaystyle=\frac{1}{n}(\hat{\boldsymbol{M}}-\boldsymbol{M})^{\top}(\hat{\boldsymbol{M}}-\boldsymbol{M})=O_{p}(c_{n}^{2}\vee n^{-1}),$
	$\displaystyle R_{2}$	$\displaystyle=\frac{1}{n}(\hat{\boldsymbol{M}}-\boldsymbol{M})^{\top}[\boldsymbol{W},\boldsymbol{M}]=O_{p}(c_{n}\vee n^{-1/2}),$
	$\displaystyle R_{3}$	$\displaystyle=\frac{1}{n}(\hat{\boldsymbol{M}}-\boldsymbol{M})^{\top}\boldsymbol{Z}=O_{p}(c_{n}\vee n^{-1/2}),$
	$\displaystyle R_{4}$	$\displaystyle=\frac{1}{n}(\hat{\boldsymbol{M}}-\boldsymbol{M})^{\top}\boldsymbol{N}=O_{p}(c_{n}).$

Proof:

First, (47) and the vector form of (43) give

	$\displaystyle\boldsymbol{\hat{M}}-\boldsymbol{M}$	$\displaystyle=A(\boldsymbol{Y}-\boldsymbol{Z}\hat{\beta})-\boldsymbol{M}$
		$\displaystyle=(A-I)\boldsymbol{M}+A\boldsymbol{Z}(\beta-\hat{\beta})+A\boldsymbol{N}.$

Observe that $R_{2}\sim R_{4}$ can be defined through

	$\displaystyle I_{1}$	$\displaystyle=\frac{1}{n}\boldsymbol{M}^{\top}(A^{\top}-I)\boldsymbol{P},$
	$\displaystyle I_{2}$	$\displaystyle=\frac{1}{n}(\beta-\hat{\beta})^{\top}\boldsymbol{Z}^{\top}A^{\top}\boldsymbol{P},$
	$\displaystyle I_{3}$	$\displaystyle=\frac{1}{n}\boldsymbol{N}^{\top}A^{\top}\boldsymbol{P},$

where $\boldsymbol{P}$ can be replaced by $\boldsymbol{W}$ , $\boldsymbol{Z}$ , $\boldsymbol{M}$ , or $\boldsymbol{N}$ (note that the rows of $\boldsymbol{P}$ are i.i.d. ). We also have $R_{1}=\sum_{i=4}^{9}I_{i}$ , where

	$\displaystyle I_{4}$	$\displaystyle=\frac{1}{n}\boldsymbol{M}^{\top}(A^{\top}-I)(A-I)\boldsymbol{M},$
	$\displaystyle I_{5}$	$\displaystyle=\frac{1}{n}(\beta-\hat{\beta})^{\top}\boldsymbol{Z}^{\top}A^{\top}A\boldsymbol{Z}(\beta-\hat{\beta}),$
	$\displaystyle I_{6}$	$\displaystyle=\frac{1}{n}\boldsymbol{N}^{\top}A^{\top}A\boldsymbol{N},$
	$\displaystyle I_{7}$	$\displaystyle=\frac{1}{n}\boldsymbol{M}^{\top}(A^{\top}-I)\bigg{\{}A\boldsymbol{Z}(\beta-\hat{\beta})+A\boldsymbol{N}\bigg{\}},$
	$\displaystyle I_{8}$	$\displaystyle=\frac{1}{n}(\beta-\hat{\beta})^{\top}\boldsymbol{Z}^{\top}A^{\top}\bigg{\{}(A-I)\boldsymbol{M}+A\boldsymbol{N}\bigg{\}},$
	$\displaystyle I_{9}$	$\displaystyle=\frac{1}{n}\boldsymbol{N}^{\top}A^{\top}\bigg{\{}(A-I)\boldsymbol{M}+A\boldsymbol{Z}(\beta-\hat{\beta})\bigg{\}}.$

It can be shown that $I_{1},I_{3}=O_{p}(c_{n})$ (we use this as a shorthand for $I_{1}=O_{p}(c_{n}),I_{3}=O_{p}(c_{n})$ , as $I_{1}$ and $I_{3}$ are not equal), $I_{2}=O_{p}(n^{-1/2})$ , $I_{5}=O_{p}(n^{-1})$ , $I_{4},I_{6}=O_{p}(c_{n}^{2})$ , $I_{8}=O_{p}(c_{n}n^{-1/2})$ , $I_{7},I_{9}=O_{p}(c_{n}^{2}\vee c_{n}n^{-1/2})$ , which implies $R_{1}=O_{p}(c_{n}^{2}\vee n^{-1})$ , $R_{2},R_{3}=O_{p}(c_{n}\vee n^{-1/2})$ and $R_{4}=O_{p}(c_{n})$ . The techniques for proving the rates of $I_{1}\sim I_{9}$ are similar; observe that all the components in $I_{4}\sim I_{9}$ are already computed to obtain $I_{1}\sim I_{3}$ , thus we only provide the proof for $I_{1}\sim I_{3}$ for simplicity of presentation.

Using Lemma 2, it has been shown in [57] that $\boldsymbol{\tilde{W}}_{u_{i}}^{\top}\boldsymbol{K}_{u_{i}}\boldsymbol{\tilde{W}}_{u_{i}}$ can be equivalently expressed as

\displaystyle\begin{bmatrix}B_{1}&B_{2}\\ B_{3}&B_{4}\end{bmatrix}=nf(U_{i})\mathsf{E}[WW^{\top}|U_{i}]\otimes\begin{pmatrix}1&0\\ 0&\mu_{2}\end{pmatrix}\{1+O_{p}(c_{n})\},

(48)

where $\mu_{2}=\int_{\mathcal{U}}u^{2}K(u)du$ , and the four block matrices are $B_{1}=S_{0}(U_{i})$ , $B_{2}=B_{3}=S_{1}(U_{i})$ , and $B_{4}=S_{2}(U_{i})$ , with respect to

\displaystyle S_{k}(U_{i})

\displaystyle=\sum_{j=1}^{n}\left(\frac{U_{j}-U_{i}}{h}\right)^{k}\boldsymbol{W}_{j}\boldsymbol{W}_{j}^{\top}K_{h}(U_{j}-U_{i}).

Since the techniques for proving (48), omitted in [57], will be used repeatedly in the rest of this paper, we provide the proof for completeness. Applying Lemma 2, it holds uniformly in $u\in\mathcal{U}$ that $S_{k}(u)$ can be expressed as

	$\displaystyle\mathsf{E}\left[\left(\frac{U-u}{h}\right)^{k}\mathsf{E}[WW^{\top}\|U]K_{h}(U-u)\right]\{1+O_{p}(c_{n}^{\prime})\}$
$\displaystyle=$	$\displaystyle\int_{\mathcal{V}}v^{k}\mathsf{E}[WW^{\top}\|U=hv+u]K(v)f(hv+u)dv\cdot\{1+O_{p}(c_{n}^{\prime})\}$	(49)
$\displaystyle=$	$\displaystyle\int_{\mathcal{V}}v^{k}K(v)\left[\mathsf{E}[WW^{\top}\|U=u]+vO(h)\right]\cdot\left[f(u)+vO(h)\right]dv\{1+O_{p}(c_{n}^{\prime})\}$	(50)
$\displaystyle=$	$\displaystyle\begin{cases}\mathsf{E}[WW^{\top}\|U=u]f(u)\mu_{k}\{1+O_{p}(c_{n})\},&k\text{ even}\\ O(h)+O_{p}(c_{n}^{\prime}),&k\text{ odd}\end{cases}$	(51)

where (49) is due to the change of variable $V=(U-u)/h$ , (50) uses the Lipschitz continuity assumptions on $\mathsf{E}[XX^{T}|U]$ and $f(\cdot)$ , and (51) is by the symmetry of the kernel function $K(\cdot)$ . Similarly, we obtain

	$\displaystyle\boldsymbol{\tilde{W}}_{u_{i}}^{\top}\boldsymbol{K}_{u_{i}}\boldsymbol{M}$
$\displaystyle=$	$\displaystyle\begin{bmatrix}\sum_{j=1}^{n}\boldsymbol{W}_{j}\boldsymbol{W}_{j}^{\top}\alpha(U_{j})K_{h}(U_{j}-U_{i})\\ \sum_{j=1}^{n}\frac{U_{j}-U_{i}}{h}\boldsymbol{W}_{j}\boldsymbol{W}_{j}^{\top}\alpha(U_{j})K_{h}(U_{j}-U_{i})\end{bmatrix}$
$\displaystyle=$	$\displaystyle nf(U_{i})\mathsf{E}[WW^{\top}\|U_{i}]\alpha(U_{i})\otimes\begin{pmatrix}1&0\end{pmatrix}^{\top}\{1+O_{p}(c_{n})\}.$	(52)

Recall the expression of $A$ in (45), we have

		$\displaystyle(A-I)\boldsymbol{M}$
	$\displaystyle=$	$\displaystyle\begin{bmatrix}[\boldsymbol{W}_{1}^{\top}\quad\boldsymbol{0}]\{\boldsymbol{\tilde{W}}_{u_{1}}^{\top}\boldsymbol{K}_{u_{1}}\boldsymbol{\tilde{W}}_{u_{1}}\}^{-1}\boldsymbol{\tilde{W}}_{u_{1}}^{\top}\boldsymbol{K}_{u_{1}}\boldsymbol{M}\\ \vdots\\ [\boldsymbol{W}_{n}^{\top}\quad\boldsymbol{0}]\{\boldsymbol{\tilde{W}}_{u_{n}}^{\top}\boldsymbol{K}_{u_{n}}\boldsymbol{\tilde{W}}_{u_{n}}\}^{-1}\boldsymbol{\tilde{W}}_{u_{n}}^{\top}\boldsymbol{K}_{u_{n}}\boldsymbol{M}\\ \end{bmatrix}-\boldsymbol{M}$

Using (48) and (Proof), we obtain

	$\displaystyle I_{1}$	$\displaystyle=\frac{1}{n}\boldsymbol{M}^{\top}(A^{\top}-I)\boldsymbol{P}$
		$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}(M_{i}\{1+O_{p}(c_{n})\}-M_{i})\boldsymbol{P}_{i}=O_{p}(c_{n}),$

where $\boldsymbol{P}_{i}$ denotes the $i^{\text{th}}$ row of $\boldsymbol{P}$ and the last equality is due to the law of large numbers. For $I_{2}$ , similarly as above, we compute

A\boldsymbol{Z}=\begin{bmatrix}\boldsymbol{W_{1}}^{\top}\left\{\mathsf{E}[WW^{T}|U_{1}]\right\}^{-1}\mathsf{E}[WZ^{\top}|U_{1}]\\ \vdots\\ \boldsymbol{W_{n}}^{\top}\left\{\mathsf{E}[WW^{T}|U_{n}]\right\}^{-1}\mathsf{E}[WZ^{\top}|U_{n}]\end{bmatrix}\{1+O_{p}(c_{n})\}.

Since $\beta-\hat{\beta}=O_{p}(n^{-1/2})$ by Lemma 3, we obtain

	$\displaystyle I_{2}$	$\displaystyle=\frac{1}{n}(\beta-\hat{\beta})^{\top}\boldsymbol{Z}^{\top}A^{\top}\boldsymbol{P}$
		$\displaystyle=(\beta-\hat{\beta})^{\top}\frac{1}{n}\sum_{i=1}^{n}\boldsymbol{W_{i}}^{\top}\left\{\mathsf{E}[WW^{T}\|U_{i}]\right\}^{-1}\cdot\mathsf{E}[WZ^{\top}\|U_{i}]\{1+O_{p}(c_{n})\}\boldsymbol{P}_{i}$
		$\displaystyle=O_{p}(n^{-1/2}),$

where, again, the last equality is due to the law of large numbers. Finally, for

\displaystyle A\boldsymbol{N}

\displaystyle=\begin{bmatrix}[\boldsymbol{W}_{1}^{\top}\quad\boldsymbol{0}]\{\boldsymbol{\tilde{W}}_{u_{1}}^{\top}\boldsymbol{K}_{u_{1}}\boldsymbol{\tilde{W}}_{u_{1}}\}^{-1}\boldsymbol{\tilde{W}}_{u_{1}}^{\top}\boldsymbol{K}_{u_{1}}\boldsymbol{N}\\ \vdots\\ [\boldsymbol{W}_{n}^{\top}\quad\boldsymbol{0}]\{\boldsymbol{\tilde{W}}_{u_{n}}^{\top}\boldsymbol{K}_{u_{n}}\boldsymbol{\tilde{W}}_{u_{n}}\}^{-1}\boldsymbol{\tilde{W}}_{u_{n}}^{\top}\boldsymbol{K}_{u_{n}}\boldsymbol{N}\\ \end{bmatrix},

the same argument for (51) leads to

\displaystyle\boldsymbol{\tilde{W}}_{u_{i}}^{\top}\boldsymbol{K}_{u_{i}}\boldsymbol{N}=nf(U_{i})\mathsf{E}[WN^{\top}|U_{i}]\{1+O_{p}(c_{n})\},

where $\mathsf{E}[WN^{\top}|U_{i}]=0$ since $N$ is independent of $W$ and $U_{i}$ , and $N$ has a zero mean. Thus, by the law of large numbers,

\displaystyle I_{3}

\displaystyle=\frac{1}{n}\boldsymbol{N}^{\top}A^{\top}\boldsymbol{P}=\frac{1}{n}\sum_{i=1}^{n}O_{p}(c_{n})\boldsymbol{P}_{i}=O_{p}(c_{n}).

$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\blacksquare$

Appendix I Proofs of Theorem 3 and Corollaries 3 and 4

To reuse model (43) for the prediction of $Y^{\tau}$ , the main challenge comes from $M$ that changes with $U$ (while $\beta$ remains invariant). First, we introduce some notation for the proofs. We define $(\boldsymbol{Y},\boldsymbol{W},\boldsymbol{W}_{i})$ and $(\boldsymbol{Z},\boldsymbol{Z}_{i},\boldsymbol{M},\boldsymbol{N})$ in the same way as in Appendix G. Similarly, we define $\boldsymbol{Z}_{V},\boldsymbol{M}_{V}$ , $\boldsymbol{N}_{V}$ , and all the corresponding data matrices for the test data (e.g., $\boldsymbol{Y}^{\tau}$ ), and let $\sigma^{2}=\mathsf{E}[(N^{\tau})^{2}]$ . With this notation, the IMP (11) implies that

M=\lambda M_{V}+\zeta^{\top}W:=[W,M_{V}]w,

(53)

for some $\zeta\in\mathbb{R}^{p}$ . Let $\boldsymbol{\hat{W}}_{V}:=[\boldsymbol{W},\boldsymbol{\hat{M}}_{V}]$ and $\boldsymbol{\hat{W}}_{V}^{\tau}:=[\boldsymbol{W}^{\tau},\boldsymbol{\hat{M}}_{V}^{\tau}]$ . Then, the OLS estimator of $w$ according to the above equation is given by

\hat{w}=(\boldsymbol{\hat{W}}_{V}^{\top}\boldsymbol{\hat{W}}_{V})^{-1}\boldsymbol{\hat{W}}_{V}^{\top}\boldsymbol{\hat{M}}.

(54)

We predict $\boldsymbol{Y}^{\tau}$ using the continuous IMP estimator

\boldsymbol{\hat{Y}}^{\tau}=\boldsymbol{\tilde{W}}_{V}^{\tau}\hat{w}+\boldsymbol{Z}^{\tau}\hat{\beta},

(55)

where $\hat{\beta}$ , $\boldsymbol{\hat{M}}$ , and $\boldsymbol{\hat{M}}_{V}$ are provided in Appendix G.

Lemma 5.

Under assumptions 1) $\sim$ 4) in Appendix H, it holds that

\hat{w}-w=O_{p}(c_{n}\vee n^{-1/2}).

Proof:

Using the fact that

\boldsymbol{\hat{W}}_{V}=[\boldsymbol{W},\boldsymbol{M}_{V}]+[\boldsymbol{0},\boldsymbol{\hat{M}}_{V}-\boldsymbol{M}_{V}],

we obtain

	$\displaystyle\frac{1}{n}\boldsymbol{\hat{W}}_{V}^{\top}\boldsymbol{\hat{W}}_{V}$
	$\displaystyle=\frac{1}{n}[\boldsymbol{W},\boldsymbol{M}_{V}]^{\top}[\boldsymbol{W},\boldsymbol{M}_{V}]+\frac{1}{n}[\boldsymbol{W},\boldsymbol{M}_{V}]^{\top}[\boldsymbol{0},\boldsymbol{\hat{M}}_{V}-\boldsymbol{M}_{V}]$
	$\displaystyle\hskip 10.00002pt+\frac{1}{n}[\boldsymbol{0},\boldsymbol{\hat{M}}_{V}-\boldsymbol{M}_{V}]^{\top}[\boldsymbol{W},\boldsymbol{M}_{V}]+\frac{1}{n}[\boldsymbol{0},\boldsymbol{\hat{M}}_{V}-\boldsymbol{M}_{V}]^{\top}\cdot[\boldsymbol{0},\boldsymbol{\hat{M}}_{V}-\boldsymbol{M}_{V}]$
	$\displaystyle=\mathsf{E}\bigg{[}(W^{\top},M_{V})^{\top}(W^{\top},M_{V})\bigg{]}\{1+O_{p}(c_{n}\vee n^{-1/2})\},$

where we use $R_{1}$ and $R_{2}$ from Lemma 4 and the law of large numbers. Similarly,

	$\displaystyle\frac{1}{n}\boldsymbol{\hat{W}}_{V}^{\top}\boldsymbol{\hat{M}}$	$\displaystyle=\frac{1}{n}[\boldsymbol{W},\boldsymbol{M}_{V}]^{\top}\boldsymbol{M}+\frac{1}{n}[\boldsymbol{W},\boldsymbol{M}_{V}]^{\top}(\boldsymbol{\hat{M}}-\boldsymbol{M})$
		$\displaystyle\hskip 50.00008pt+\frac{1}{n}[\boldsymbol{0},\boldsymbol{\hat{M}}_{V}-\boldsymbol{M}_{V}]^{\top}(\boldsymbol{\hat{M}}-\boldsymbol{M})+\frac{1}{n}[\boldsymbol{0},\boldsymbol{\hat{M}}_{V}-\boldsymbol{M}_{V}]^{\top}\boldsymbol{M}$
		$\displaystyle=\mathsf{E}\bigg{[}(W^{\top},M_{V})^{\top}M\bigg{]}\{1+O_{p}(c_{n}\vee n^{-1/2})\}.$

Note that the rate is not directly implied by Lemma 4, but it can be proved using the same techniques demonstrated in the proof Lemma 4. Therefore, using (53),

\displaystyle\hat{w}=(\boldsymbol{\hat{W}}_{V}^{\top}\boldsymbol{\hat{W}}_{V})^{-1}\boldsymbol{\hat{W}}_{V}^{\top}\boldsymbol{\hat{M}}=w\{1+O_{p}(c_{n}\vee n^{-1/2})\}.

$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\blacksquare$

I-A Proof of Theorem 3

Using (55) and $\boldsymbol{Y}^{\tau}=\boldsymbol{M}^{\tau}+\boldsymbol{Z}^{\tau}\beta+\boldsymbol{N}^{\tau}$ , we derive

		$\displaystyle\boldsymbol{\hat{Y}^{\tau}}-\boldsymbol{Y}^{\tau}$
	$\displaystyle=$	$\displaystyle[\boldsymbol{W}^{\tau},\boldsymbol{\hat{M}}^{\tau}_{V}]\hat{w}+\boldsymbol{Z}^{\tau}\hat{\beta}-[\boldsymbol{W}^{\tau},\boldsymbol{M}^{\tau}_{V}]w-\boldsymbol{Z}^{\tau}\beta-\boldsymbol{N}^{\tau}$
	$\displaystyle=$	$\displaystyle[\boldsymbol{0},\boldsymbol{\hat{M}}^{\tau}_{V}-\boldsymbol{M}^{\tau}_{V}](\hat{w}-w)+[\boldsymbol{W}^{\tau},\boldsymbol{M}^{\tau}_{V}](\hat{w}-w)+[\boldsymbol{0},\boldsymbol{\hat{M}}^{\tau}_{V}-\boldsymbol{M}^{\tau}_{V}]w+\boldsymbol{Z}^{\tau}(\hat{\beta}-\beta)-\boldsymbol{N}^{\tau}.$

Then, the generalization error is given by

\displaystyle\frac{1}{m}(\boldsymbol{\hat{Y}^{\tau}}-\boldsymbol{Y}^{\tau})^{\top}(\boldsymbol{\hat{Y}^{\tau}}-\boldsymbol{Y}^{\tau})=\sum_{i=1}^{10}J_{i}

where

	$\displaystyle J_{1}$	$\displaystyle=\frac{1}{m}(\hat{w}-w)^{\top}[\boldsymbol{0},\boldsymbol{\hat{M}}^{\tau}_{V}-\boldsymbol{M}^{\tau}_{V}]^{\top}[\boldsymbol{0},\boldsymbol{\hat{M}}^{\tau}_{V}-\boldsymbol{M}^{\tau}_{V}](\hat{w}-w),$
	$\displaystyle J_{2}$	$\displaystyle=\frac{1}{m}(\hat{w}-w)^{\top}[\boldsymbol{X}^{\tau},\boldsymbol{M}^{\tau}_{V}]^{\top}[\boldsymbol{X}^{\tau},\boldsymbol{M}^{\tau}_{V}](\hat{w}-w),$
	$\displaystyle J_{3}$	$\displaystyle=\frac{1}{m}w^{\top}[\boldsymbol{0},\boldsymbol{\hat{M}}^{\tau}_{V}-\boldsymbol{M}^{\tau}_{V}]^{\top}[\boldsymbol{0},\boldsymbol{\hat{M}}^{\tau}_{V}-\boldsymbol{M}^{\tau}_{V}]w,$
	$\displaystyle J_{4}$	$\displaystyle=\frac{1}{m}(\hat{\beta}-\beta)^{\top}(\boldsymbol{Z}^{\tau})^{\top}\boldsymbol{Z}^{\tau}(\hat{\beta}-\beta),$
	$\displaystyle J_{5}$	$\displaystyle=\frac{1}{m}(\boldsymbol{N}^{\tau})^{\top}\boldsymbol{N}^{\tau},$
	$\displaystyle J_{6}$	$\displaystyle=\frac{1}{m}(\hat{w}-w)^{\top}[\boldsymbol{0},\boldsymbol{\hat{M}}^{\tau}_{V}-\boldsymbol{M}^{\tau}_{V}]^{\top}\bigg{\{}[\boldsymbol{W}^{\tau},\boldsymbol{M}^{\tau}_{V}]$
		$\displaystyle\hskip 50.00008pt\cdot(\hat{w}-w)+[\boldsymbol{0},\boldsymbol{\hat{M}}^{\tau}_{V}-\boldsymbol{M}^{\tau}_{V}]w+\boldsymbol{Z}^{\tau}(\hat{\beta}-\beta)-\boldsymbol{N}^{\tau}\bigg{\}},$
	$\displaystyle J_{7}$	$\displaystyle=\frac{1}{m}(\hat{w}-w)^{\top}[\boldsymbol{W}^{\tau},\boldsymbol{M}^{\tau}_{V}]^{\top}\bigg{\{}[\boldsymbol{0},\boldsymbol{\hat{M}}^{\tau}_{V}-\boldsymbol{M}^{\tau}_{V}]$
		$\displaystyle\hskip 50.00008pt\cdot(\hat{w}-w)+[\boldsymbol{0},\boldsymbol{\hat{M}}^{\tau}_{V}-\boldsymbol{M}^{\tau}_{V}]w+\boldsymbol{Z}^{\tau}(\hat{\beta}-\beta)-\boldsymbol{N}^{\tau}\bigg{\}},$
	$\displaystyle J_{8}$	$\displaystyle=\frac{1}{m}w^{\top}[\boldsymbol{0},\boldsymbol{\hat{M}}^{\tau}_{V}-\boldsymbol{M}^{\tau}_{V}]^{\top}\bigg{\{}[\boldsymbol{0},\boldsymbol{\hat{M}}^{\tau}_{V}-\boldsymbol{M}^{\tau}_{V}]$
		$\displaystyle\hskip 50.00008pt\cdot(\hat{w}-w)+[\boldsymbol{W}^{\tau},\boldsymbol{M}^{\tau}_{V}](\hat{w}-w)+\boldsymbol{Z}^{\tau}(\hat{\beta}-\beta)-\boldsymbol{N}^{\tau}\bigg{\}},$
	$\displaystyle J_{9}$	$\displaystyle=\frac{1}{m}(\hat{\beta}-\beta)^{\top}(\boldsymbol{Z}^{\tau})^{\top}\bigg{\{}[\boldsymbol{0},\boldsymbol{\hat{M}}^{\tau}_{V}-\boldsymbol{M}^{\tau}_{V}](\hat{w}-w)$
		$\displaystyle\hskip 50.00008pt+[\boldsymbol{W}^{\tau},\boldsymbol{M}^{\tau}_{V}](\hat{w}-w)+[\boldsymbol{0},\boldsymbol{\hat{M}}^{\tau}_{V}-\boldsymbol{M}^{\tau}_{V}]w+\boldsymbol{N}^{\tau}\bigg{\}},$
	$\displaystyle J_{10}$	$\displaystyle=-\frac{1}{m}(\boldsymbol{N}^{\tau})^{\top}\bigg{\{}[\boldsymbol{0},\boldsymbol{\hat{M}}^{\tau}_{V}-\boldsymbol{M}^{\tau}_{V}](\hat{w}-w)$
		$\displaystyle\hskip 50.00008pt+[\boldsymbol{W}^{\tau},\boldsymbol{M}^{\tau}_{V}](\hat{w}-w)+[\boldsymbol{0},\boldsymbol{\hat{M}}^{\tau}_{V}-\boldsymbol{M}^{\tau}_{V}]w+\boldsymbol{Z}^{\tau}(\hat{\beta}-\beta)\bigg{\}}.$

We can show the following rates of these terms through simple applications of Lemmas 3, 4, 5, along with the law of large number and the central limit theorem. $J_{2},J_{3}=O_{p}(c_{n}^{2}\vee n^{-1})$ , $J_{4}=O_{p}(n^{-1})$ , $J_{5}=\sigma^{2}+O_{p}(m^{-1/2})$ ,

\displaystyle J_{10}

\displaystyle=O_{p}(c_{m}\vee m^{-1/2})\cdot O_{p}(c_{n}\vee n^{-1/2})+O_{p}(c_{n}\vee n^{-1/2})+O_{p}(c_{m}\vee m^{-1/2})+O_{p}(n^{-1/2}),

and $J_{1},J_{6}\sim J_{9}$ have higher orders compared with either $O_{p}(c_{n}\vee n^{-1/2})$ or $O_{p}(c_{m}\vee m^{-1/2})$ in $J_{10}$ . This completes the proof of Theorem 3.

I-B Proof of Corollary 3

We use $l$ as the shorthand notation for $l_{n}$ in the proof. Since the estimation of $\beta,\boldsymbol{M},$ and $w$ are based on the labeled data, we have $\hat{\beta}-\beta=O_{p}(l^{-1/2})$ in Lemma 3 and $R_{1},R_{2},R_{3}=O_{p}(c_{l}\vee l^{-1/2})$ and $R_{4}=O_{p}(c_{l})$ in Lemma 4. In the proof of Lemma 5, the rate of $\boldsymbol{\hat{W}}_{V}^{\top}\boldsymbol{\hat{W}}_{V}$ remains the same since the estimation of $\boldsymbol{M}_{V}$ only uses the unlabeled data but the rate of $\boldsymbol{\hat{W}}_{V}^{\top}\boldsymbol{\hat{M}}$ now depends on $l$ . This observation implies $\hat{w}-w=O_{p}(c_{l}\vee l^{-1/2})$ . Now, observe that $J_{2}=O_{p}(c_{l}^{2}\vee l^{-1})$ , $J_{3}=O_{p}(c_{m}\vee m^{-1/2})$ , $J_{4}=O_{p}(l^{-1})$ , $J_{5}=\sigma^{2}+O_{p}(m^{-1/2})$ , and $J_{10}=O_{p}(c_{l}\vee l^{-1/2})+O_{p}(m^{-1/2})$ . Since $\max(\frac{l}{n},\frac{l}{m})\to 0$ as $\min(n,m)\to\infty$ , we get $\sum_{i=1}^{10}J_{i}=O_{p}(c_{l}\vee l^{-1/2})$ .

I-C Proof of Corollary 4

According to the definition of the discrete IMP estimator, all the coefficients are treated as varying coefficients (i.e., $\beta=\boldsymbol{0}$ and $\beta_{V}=\boldsymbol{0}$ ), and $\boldsymbol{M}$ is estimated by performing the OLS for each environment and then putting the estimates into one vector, thus $R_{1}=O_{p}(a_{n}^{2})$ and $R_{2},R_{4}=O_{p}(a_{n})$ in Lemma 4, with $a_{n}=(\min_{e\in\mathcal{E}^{\text{train}}}n_{e})^{-1/2}$ by the asymptotic normality of the OLS estimators ( note that $X$ and $Z$ are assumed to have finite fourth moments in Lemma 3). Accordingly, we have $\hat{w}-\hat{w}=O_{p}(a_{n})$ in Lemma 5 using the law of large numbers. Setting $Z^{\tau}=\boldsymbol{0}$ in $J_{1}\sim J_{10}$ , we obtain $J_{2},J_{3}=O_{p}(a_{n}^{2})$ , $J_{4}=0$ , $J_{5}=\sigma^{2}+O_{p}(m^{-1/2})$ , $J_{10}=O_{p}(a_{n})+O_{p}(a_{m})$ , and the other terms have higher orders compared with $O_{p}(a_{n})$ or $O_{p}(a_{m})$ . Similar to Theorem 3, the asymptotic generalization error is dominated by $J_{10}$ .

$\displaystyle\mathsf{E}_{l}[Y\|X_{S},E=e]$	$\displaystyle=\lambda\mathsf{E}_{l}[X_{k}\|X_{R},E=e]+\eta^{\top}X,$	(29)
$\displaystyle\mathsf{E}_{l}[Y\|X_{S},E=e]$	$\displaystyle=\lambda_{Y}\mathsf{E}[Y\|X_{PA(Y)},E=e]+\eta_{Y}^{\top}X,$	(30)
$\displaystyle\mathsf{E}_{l}[X_{k}\|X_{R},E=e]$	$\displaystyle=\lambda_{X}\mathsf{E}[Y\|X_{PA(Y)},E=e]+\eta_{X}^{\top}X,$	(31)

$\displaystyle\mathsf{E}_{l}[X_{k}\|X_{R},E=e]$	$\displaystyle=\mathsf{E}_{l}[X_{k}\|Z(e),X_{R},E=e]$
	$\displaystyle=\mathsf{E}_{l}[X_{k}\|Z(e),X_{R}]$
	$\displaystyle=cZ(e)+d^{\top}X_{R}$	(38)

Learning Invariant Representations under General Interventions on the Response

Abstract

Index Terms:

1 Introduction

1-A Further Related Work

1-B Contribution and Structure

2 Background and Problem Formulation

Definition 1.

3 Invariant Matching Property and Theoretical Guarantees

3-A Motivating Examples

Example:

Example (additional settings):

Remark 1.

Example (cyclic graph setting):

Remark 2.

3-B Invariant Matching Property

Definition 2.

Definition 3.

Proposition 1.

3-C A Decomposition of the IMP

Definition 4.

Definition 5.

4 Characterization of Invariant Matching Properties

4-A Interventions on the Response

Remark 3.

Assumption 1.

Theorem 1.

Corollary 1.

Proposition 2.

Corollary 2.

Remark 4.

4-B Interventions on both Predictors and Response

Assumption 2.

Theorem 2.

Remark 5.

Proposition 3.

5 Algorithms

5-A Discrete Environments

5-B Continuous Environments

Remark 6.

Remark 7.

6 Asymptotic Generalization Error

Theorem 3.

Corollary 3.

Corollary 4.

7 Experiments

7-A Discrete environments

7-A1 Interventions on XX

7-A2 Interventions on YY

7-A3 Interventions on both XX and YY

7-B Interventions on both XX and YY (continuous)

7-C Robustness

7-D COVID Data Set

8 Discussion

References

Appendix A Proof of Proposition 1

Appendix B Problem Formulation Using an environmental random variable

Appendix C Proof of Theorem 1

Remark 8.

Appendix D Proof of Proposition 2

Lemma 1.

Proof:

D-A Proof of Proposition 2

Appendix E Proof of Corollary 2

Remark 9.

Appendix F Proof of Theorem 2

Remark 10.

Appendix G Semi-parametric Varying Coefficient Models and Profile Least-Squares Estimation

Appendix H Technical Lemmas for the Proof of Theorem 3

Lemma 2 (​[61]).

Lemma 3 (​[57]).

Lemma 4.

Proof:

Appendix I Proofs of Theorem 3 and Corollaries 3 and 4

Lemma 5.

Proof:

I-A Proof of Theorem 3

I-B Proof of Corollary 3

I-C Proof of Corollary 4

7-A1 Interventions on $X$

7-A2 Interventions on $Y$

7-A3 Interventions on both $X$ and $Y$

7-B Interventions on both $X$ and $Y$ (continuous)

Lemma 2 ([61]).

Lemma 3 ([57]).