Which Invariance Should We Transfer? A Causal Minimax Learning Approach
Authors: Mingzhou Liu
Xiangyu Zheng
Xinwei Sun
Fang Fang
Yizhou Wang
Venue: Submitted to ICLR 2023 (2023)
Keywords: robustness
minimax
subset selection
causal model
g-equivalence

Abstract: A major barrier to deploy current machine learning models lies in their sensitivity to dataset shifts. To resolve this problem, most existing studies attempted to transfer stable information to unseen environments. Among these, graph-based methods causally decomposed the data generating process into stable and mutable mechanisms. By removing the effect of mutable generation, they identified a set of stable predictors. However, a key question regarding robustness remains: which subset of the whole stable information should the model transfer, in order to achieve optimal generalization ability? To answer this question, we provide a comprehensive minimax analysis that fully characterizes conditions for a subset to be optimal. Particularly in general cases, we propose to maximize over mutable mechanisms (i.e., the source of dataset shifts), which is provable to identify the worst-case risk over all environments. This ensures us to select the optimal subset with the minimal worst-case risk. To reduce computational costs, we propose to search over only equivalent classes in terms of worst-case risk, instead of over all subsets. In cases when the searching space is still large, we turn this subset selection problem into a sparse min-max optimization scheme, which enjoys the simplicity and efficiency of implementation. The utility of our methods is demonstrated on the diagnosis of Alzheimer's Disease and gene function prediction. 

Under review as a conference paper at ICLR 2023
WHICH INVARIANCE SHOULD WE TRANSFER? A CAUSAL MINIMAX LEARNING APPROACH
A Anonymous authors
Paper under double-blind review
ABSTRACT
A major barrier to deploy current machine learning models lies in their sensitivity dataset shifts. To resolve this problem, most existing studies attempted transfer stable information to unseen environments. Among these. graph-based methods causally decomposed the data generating process into stable and mutable mechanisms. By removing the effect of mutable generation, they identified a set of stable predictors. However, key question regarding robustness remains: which subset of the whole stable information should the model transfer, in order to achieve optimal generalization ability? ? To answer this question, we provide a comprehensive minimax analysis that fully characterizes conditions for a subset be optimal. Particularly in general cases, we propose to maximize over mutable mechanisms (i.e., the source of dataset shifts), which is provable identify the worst-case risk over all environments. This ensures us to select the optimal subset with the minimal worst-case risk. To reduce computational costs, we propose search over only equivalent classes in terms of worst-case risk, instead of over all subsets. In cases when the searching space is still large, we turn this subset selection problem into a sparse min-max optimization scheme, which enjoys the simplicity and efficiency of implementation. The utility of our methods is demonstrated on the diagnosis of Alzheimer's Disease and gene function prediction.
INTRODUCTION
Current machine learning systems, which are commonly deployed based on their in-distribution performance, often encounter dataset shifts (Subbasaay al. (2019) such as covariate shift, label shift, etc., due to changes in the data generating process. When such a shift exists in deployment environments, the model may give unreliable prediction results, which can cause severe consequences in safe-critical tasks such as healthcare (Hendrycks et al. 2021). At the heart of this unreliability issue are stability and robustness aspects, which refer to the insensitivity of prediction behavior and generalization errors over shifts, respectively.
For example, consider the system deployed to predict the Functional Activities Questionnaire (FAQ) score, which is commonly adopted Mayo (016)) to measure the severity of Alzheimer's Disease (AD). During prediction, the system can only access biomarkers or volumes of brain regions with anonymous demographic information for privacy consideration. However, the changes in demographics can cause shifts in covariates. To achieve reliability for the deployed model, it is desired for its prediction to be stable against demographic changes, and meanwhile to be constantly accurate over all different populations. To incorporate both aspects, this paper targets at finding the most robust (i.e., min-max Mullrr et al. (2020)) predictor, among the set of stable predictors over all distributions.
MptimalMMüller
To achieve this goal, many studies have proposed to learn invariance to transfer to unseen data, Examples include ICP Peters et al.) (2016) and (Arjovsky et al., 2019 Liu et al.] 2021 Ahuja et al. (02T) that assumed the prediction mechanism given causal features or representations to be invariant; Anchor Regression Rothenhäusler et al. (2021) that explicitly attributed the variation to exogenous variables. Particularly, the Subbaswamy & Saria (2020); (Subbaswamy et al. (2019) causally decomposed the joint distribution into mutable M and stable S sets, with respectively changed and unchanged causal mechanisms. They then proposed to intervene on M obtain a set of stable predictors. Still, a question regarding robustness remains: which subset of stable information should the model transfer, in order to be most robust against dataset shifts? The answer given by

---

Under review as a conference paper at ICLR 2023
30%
26.4%
25.3%
Improvement over vani Vanilla
25%
20%
1:21
15.5%
Max. Test MSE
15%
13.1%
10%
5%6-
2.8%
0%
0%-
-3.6%
Ours Surg(SOTA)
MPM
Out
surgery
(b)
Figure 1: FAQ prediction Alzheimer's Disease. (a) Maximal mean square error (MSE) over test environments; (b) Maximal MSE of predictors that are ranked in ascending order from left to right, respectively according to the estimated worst-case risk of our method and the validation's loss the graph surgery estimator Subbassamy al. (2019). As shown, our method is more reflective of the maximal MSE than the graph surgery method.
Subbaswamy et al. (2019) was to simply search over all subsets S and took the one with the minimal validation loss. However, it lacks theoretical and practical guarantees for the validation's loss to reflect the worst-case risk, as shown by Fig. (b).
To give comprehensive whole stable set to be optimal. This condition can be easily tested via causal discovery. When this condition fails, prove that the worst-case risk be identified by maximizing over the generating mechanism of M, iee, the only source of shift. This conclusion ensures us to select the optimal subset in
first provide graphical condition that is sufficient for the
accurate way. Consider again the example of FAQ prediction in AD diagnosis, Fig. (b) shows that our method is more reflective of the maximal mean squared error (MSE) than Subbaswamy
(2019), which explains our advantage in predicting FAQ across patient groups shown in Fig. (a). Besides, to reduce the searching cost, we propose to search over only equivalent classes in terms of worst-case risk. We however find that in some cases such search can still be expensive. To improve efficiency in these cases, we turn this subset selection task into sparse min-max optimization scheme, which alternates between gradient ascent step on the M's generating function and sparse optimization with Lasso-typp penalty to detect the optimal subset. We demonstrate the utility of our methods on a synthetic dataset and two real-world applications: Alzheimer's Disease diagnosis and gene function prediction.
Contributions. We summarize our contributions as follows:
1. We propose to identify the optimal subset of invariance transfer, guided by a comprehen- sive min-max analysis. To the best of our knowledge, are first to comprehensively study the problem of which part among all sources of invariance should the model transfer, in the literature of robust learning.
2. We introduce the concept of "equivalent relation" in terms of worst-case risk, in order to analyze the computational complexity, and propose a sparse min-max optimization method surrogate scheme to improve efficiency.
3. Our method significantly outperform others in terms of subset selection and generaliza- tion robustness, on Alzheimer's 's Disease diagnosis and gene function prediction.
2 PRELIMINARIES AND BACKGROUND
Problem Setup & Notations. We consider the supervised regression scenario, where the system includes a target variable Y e, a multivariate predictive variable X := [XI..., Xa] E Rd and data collected from heterogeneous environments. In practice, different "environments" can refer to different groups of subjects or different experimental settings. We use (D.le e ErT} to denote training data, with D.: {( vt)) Miad pf (2r, y) being the data from environment with
2

---

Under review as a conference paper at ICLR 2023
sample size ne- The total number of training samples is n Ee Ree We say a predictor f F: y stable it can be learned from the training environments Στy and transferred to a broader family of environments & without any adjustment. We denote the stable predictor set F and the distribution set as P == { (X. F'EE, with P*(X, (X, Y) the distribution over TR x y in environment e. For a causal directed acyelic graph (DAG) G, denote the parents, children, and descendants of the node set V as Pa(V) (), Ch(V) and De( De) respectively. For subset V c V, Gy denotes the sub-graph obtained by deleting edges pointing to any member of V'. We denote conditional independence and d d-separation by 1 and LG, respectively.
Our goal is to find the most robust predictor f among the stable predictor set F using data from En. A commonly used way Peters et al.] (201)); Ahuja et al. (2021) to measure this robustness is to investigate the predictor's worst-case risk, which provides safeguard for deployment in unseen environments. That is, we want f to have the following min-max property:
(x) argmin max Ep- (V,
fEF3 PcCP
Next, we introduce the causal model, Markovian and faithfulness assumptions that our methods are based These assumptions are commonly made in the literature of causal inference and learning Pearl (2009) Spirtes et al. (2000); Arjovsky et al. (2019).
Assumption 2.1 (Causal Model). We assume that YY) is entailed by unknown DAG (V, E) for all E E, where V := X JUD. denotes the node set and E denotes the edge set. Each variable Vi V is associated with structural equation of Vi <-g, (Pa(V)) U2), where U, denotes the exogenous variable. Each edge in E represents a direct causal relationship (Pearl, 1995). Assumption 2.2 (Markovian and Faithfulness). The Markovian means {U, } mutually independent.
Together with faithfulness, it means's V disjoint sets Va, V,, V&: V, 1 V, Vv + V1 La Vvk.
Graph Surgery Estimator. Under assumptions 2.1 2.2 the graph surgery estimator Subbassamy al. 2019) causally decomposed the joint distribution p" (I, y) into disentangled generating factors:
p² v) p(ulpa(v))
II (x:))))(())) EES
II p" (malpa(zx)) ds = |S], das- = IM, +EM
(2)
where 8, M respectively denote stable and mutable sets such that Xx:=(X,IVeeE (x::)e = p(x)|pa(z;))} contains variables with stable mechanisms and XM == {x, e # 62 E E, pd (z; pa();)) * (xx))) contains those with unstable mechanisms. They then re- moved the unstable mechanisms by intervening on XM and obtained a set of stable predictors FS; {fss = Ey|Y|æs__ do(x;))][S__ C S} that independent of In this regard, identifying f" in Eq. Dis is equivalent to selecting the optimal subset S* CS such that f = fs.
To identify S*, they shown that the whole set S is optimal under the degeneration condition: p(ylæs. do(xx)) = p(ylæ=') for some X' C X. In the general cases, they searched over all S_ CS S and selected the best with minimal held-out validation loss. However, this analysis is theoretically incomplete and practically defective: i) it does not provide a procedure to test the degeneration condition, making it inapplicable; ii) the selected subset may not be min-max optimal, as the validation loss does not necessarily reflect the worst-case risk (Fig. (b)); iii) the searching cost exponentially expensive w.r.t. ds, making it hard be applied to large-scale scenarios. In the next section, we will provide a comprehensive min-max analysis to remedy these issues.
3 METHODOLOGY
In this section, we introduce our method to identify S*, Specifically, in Sec. 3.1 we first introduce a comprehensive min-max analysis to identify the S* followed by the learning method in Sec. 3.2 We then analyze the computational complexity in Sec. 3.3 via the lens of g g-equivalence, and show that its searching cost can still be exponentially expensive in some cases. To improve efficiency in these cases, we in Sec. [3.4] introduce a sparse min-max optimization algorithm. which turns the subset selection problem into sparse optimization scheme that enjoys model selection consistency.
3.1 IDENTIFICATION WITH MIN-MAX ANALYSIS
In this section, we introduce our method to identify the min-max optimal subset S* with theoretical guarantees. Our analysis is composed of two main results: Thm. 3.1 and Thm. 3.3 First, Thm.
3

---

Under review as a conference paper at ICLR 2023
3.1 provides testable graphical condition that is sufficient for S* S. When this condition fails, we show that the whole stable set is not necessarily optimal via a counter- counter-example. We then in Thm. 3.3 provide a sufficient and necessary condition for a subset to be optimal in the general cases.
In the following, we first introduce the graphical condition for S = S. which is equivalent to the degeneration condition in (nSSbbaawamyy et al. (2019).
Theorem 3.1 (Graphical Condition for f = fs). Suppose assumptions 2.7 2.2 hold. Denote = (CCY as mutable variables in Y' 'ss children, and K := De(X (xx) M) X "M as descendants
of X ir. Then, p do(xM)) can degenerate to conditional distribution if and only if does not point to any member of K. Further, under either of the two conditions, we have S S.
Example 1. To illustrate, consider the example shown in Fig.2 where the graphical condition holds when the dashed arrow from Y - K does not exist. To see its equivalence to the degeneration condition, we can set K as stable variables. When Y A K, the path from Y to K can be blocked by X2,- We then have Y -xxy K and thus p(ykk, do(2)2)) = p(y) according to the inference rules in Pearl (2009). While when such a dashed arrow of Y + K holds, the K becomes a collider on the path between Y and X making incapable remove "do" in p(y k do(ar))))
Compared to the degeneration condition, our graphical condi- tion is more intuitive and can be easily tested via causal discov- ery, as guaranteed by the following proposition.
Proposition 3.2. Under assumptions 2.1 2.2 we have i) is identifable, and ii) Y + K is testable from joint distribution over training environments.
Y
0
K
Thm.
3. also reminds us that the sufficient condition S may not always hold. Indeed, we provide a counter-example in Sect. [B. in the appendix showing that the whole stable set has a larger worst-case risk than some subsets.
Figure 2: Illustration of Thm.
3.1
To identify S* when the graphical condition fails, we turn to estimate the expected worst-case risk of each subset S_ C S from (De Jeeee By noticing that the variation of unstable mechanisms in Xa is the only source of shifts in (X,XY,,, we propose to parameterize these mechanisms and let them vary arbitrarily, in order to explore the behavior of the worst-case environment. Specifically, we XM and Py P(Y, X s/ld((> X J [(Pa (XPa(aaa as)) ))). By maximizing the population risk of over J for each subset, we can obtain the worst-case risk of this subset, as shown in Thm. 3.3
consider distribution family {P for any J -PaXMr)
Theorem 3.3 (Min-max Property). Denote h hS(s)) == max [Y (x)) as the maximal
expected risk in PJ for subset S_. Then, we have (S-) maxp EP Ep- [( - Js (a:) 2). In this regard, the optimal subset S can be astained via S = argmins cs h" (S_).
This theorem informs us to optimize EP, [[Y (x))²] over J to obtain h- for S_, as it equals the worst-case risk of using subset S_. With this theorem, it is sufficient to compare h² of each subset to identify the optimal one. The following proposition ensures that the optimization is tractable, each component, ie., Pa(XM), P,, fs_ used in optimization is identifiable.
Proposition 3.4. Under assumptions
2.1 2.2 the Pa(XM), PJ, and fs_ are identifiable.
3.2 LEARNING METHOD
According to the last section, we have S* S A K. Otherwise, we can simply search over all subsets of 5 and compare their h² identify the optimal one, as similarly adopted in (Subbaswamy et al. (019). However, this exhaustive search can be redundant, as some subsets are equivalent in terms of prediction. Formally speaking, we introduce g- g-equivalence, i.e., ~G, as follows:
Definition 3.5 (g-equivalence), For two subsets S;, Sy. we say S; ~G Sy if 3S, S S; nS, such that Y -Gxx (Xs. UXs, (Xs., |Xsy XM- We call the elements of the quotient space Pow(S)/ N as g-equivalent classes, and denote NG PPo((S)/ ~G the number of equivalent classes.
Under assumption 2.2, it is easy to see that if S ~g Sj. then we have P do(X P(Y X Ssy do (X.)) and thus S,, Sy have the same efficacy of robustness. In this regard, it suffices
4

---

Under review as a conference paper at ICLR 2023
to search over Pow(S)/ ~G identify the optimal subset, rather than exhaustively searching the Pow(S). To enable this searching, we provide a recovering algorithm that is provable to recover all g g-equivalent classes. For the reason of coherence and space-saving, leave this algorithm and its analysis to Sec. C.1 in the appendix. Equipped with Pow(S)/ ~G. we are now ready to introduce algorithm to identify 5"
Algorithm Identification of the min- min-max optimal subset S and predictor f*
INPUT: The training data (Da e E &n..
Causal discovery to obtain the partially directed acylic graph (PDAG).
2: Detect K and whether Y K.
3: if Y A K then
4: Set 5 = and estimate f fss-
according to Thm. 3.1
5: else
6: Recover Pow(5)/ ~G-
with Alg. Sec. C.1
7: Set Amin = 00, 5" Ø.
8: for Sa E Pow(S)) ~~ do
9: Randomly pick a S_ E Sa, estimate fs and h (S_).
if hS_ min then
11:
Set hmin = h (S_), S_,f" === fs--
according to Thm. 3.3
end if
13:
end for
end if
15:
return S* and f*
As the causal graph is unknown, Alg. ! involves i) causal discovery to detect K and whether Y + K; ii) estimation [s and fii) estimation of h (S").. ). In the following, we roughly introduce the main ideas of our method and left the details to Sec. in the appendix.
Causal discovery to detect K and examine whether Y + K. We first detect a partial directed acyclic graph (PDAG) via the PC algorithm (Spirtes et al. 2000), followed by our method to detect Xas under the assistance of domain index variable E. Specifically, we have X; E XM iff E + X; in the detected PDAG, according Huang et al. (020)). In a similar way, we can identify Pa(X,). Ch(X,) for E M, which is sufficient to detect X := XM n Ch(Y). Applying this method iteratively, we can detect De(X, (e(x)) M)) and Pa(X;) for X; E De(X w) which is sufficient to identify K = De(X)))))) x9, \ X and have Y + K iff EE; E K such that Y E Pa(X,).
M
Estimate fs_. We adopt soft-intervention Eberhardt & Scheines (207 replace (*(XMM|PA(S) Pa(X.)
with P(XMr) and define P(X, Y) == PaY))
Then we have Pe 1" (((a)), with p. n./n. fs_ (x) Ep((Y/ss__.... MM]. Here we set p(#M) Eceen EEeer. Pe
To generate data from P. we first permute XMr in sample-wiss manner to generate data from P(MM). Then, we recursively regenerate data for each variable in De(XM) from its parents in G . by estimating structural equations. This is tractable since De(XM) and the parent nodes for each variable De(XM) are identifiable, mentioned earlier.
Estimate h²" (S-). We first learn h(SS,,)) .J := Ep, [Y fs_ As fs_ can be estimated, we only need to obtain data from Pj. As Pa(X (M)) identifiable, we iteratively regenerate data for X Ms from J (Paa(X.))) (M)) and also for X, E Pa(Xar) from its parents, in order to obtaín samples from p.-
Then we maximize h(S_,J) over J to obtain (S_)).
3.3 COMPLEXITY ANALYSIS
In this section, we discuss the complexity of Alg. I Benefit from the testable condition in Thm. 3.1 algorithm enjoys constant cost when Y does not point to any member of K. This situation happens when the target variable of interest represents the effect of the predictive covariates, egg. the number of bike riding decided by temperature, weather, etc.
When Y does point to K, e.g., represents the disease and K represents its symptoms or biomarkers, Alg. needs search among geequivalent classes and the complexity is O(Na). We will give some
5

---

Under review as a conference paper at ICLR 2023
(a)
(b)
(c)
Figure Examples of Gx = (a) chain (b) skip-chain (c) knot.
examples to show that NG can be polynomial w.r.t. dy in some cases while can also exponentially increase w.r.t. dg, depending on the number of edges in the graph. Before this, we first show (Lemma. D.I in appendix) that NG will not decrease (increase) if we add (delete) edges to (from) G.
Claim 3.6. For chain, we have NG O(ds); for the skip-chain, we have No = (d ) where k is the number of added edges; for the knot graph, we have NG = =O(css for some 1 < <<
We first consider the chain graph, i.e., Y Xs.(1) + - Xs,(ds) (where (1):.., (ds) a permutation of 1,... dy according to the generating order in the chain) in Fig. [ in which we find NG is polynomial as shown in Claim 3.6 This is simply because blocking X $,00 will d-separate Y and Xs.)) for any j > i, making (xs, g g-equivalent to (X).()) )s.(i).
ji S ... j& for any k.. Next, we consider two cases of adding & edges to the chain: i) the skip-chain graph (Fig. (b)) where & does not increasing with ds. In this case, NG is still polynomial since the number of paths between Y and any Xs.()) can be bounded; ii) the knot graph (Fig. 3 (c)) where increases with ds. in this case, NG is shown to exponentially increase w.r.t. ds, because the number of paths between Y and Xs.()) can be exponentially large. Proof of Claim 3.6 and more examples are left to Sec. C.I in the appendix.
3.4 SPARSE MIN-MAX OPTIMIZATION
According to the previous discussion, the overall searching complexity can still be exponentially large. To improve efficiency, we provide alternative method that turns the subset selection problem of Eq. into the following sparse min-max optimization scheme:
min max Ae
+ A|13811.
(3)
0,
where we introduce the coefficient vector 3 and implement lasso-type penalty on 3 with hyperparam- eter À > 0. This penalty regularizes /3 to sparse and its support set, i.e., supp(3) == (i: 3(1) * 0} is used to select the optimal subset. To optimize, we alternatively take a gradient ascent with respect to 8, followed by the minimization over (a, 3). Note that under irrepresentable and restricted convexity conditions Zhao & Yu (2006), we have model selection consistency, i.e., the true support set of 3 can be recovered and fg-consistency properties when ds is fixed, according to Rejchel (2016). showed that this lasso-type
When dg increases with n, under restricted convexity conditions, estimator that belongs to a broader family of M M-estimators Negahban et al. (2012). are Ez-consistent. To further reduce the complexity in the minimization step, we propose to implement Linearized Bregman Iteration (LBI) via differential inclusion, which enjoys the efficiency of generating a whole regularization solution path. In each iteration, the original minimization step can be replaced by gradient descent step followed by soft thresholding step. Details are left to Sec. in the appendix.
4 EXPERIMENT
In this section, we evaluate our method on synthetic data and two real-world applications: diagnosis of Alzheimer's Disease which is one of the most common types of dementia among elder people, and gene function prediction that can potentially help understand the human-disease progress Muñoz Fuentes et al. (2018).
Compared Baselines. We compare our methods with the following baselines: i) ICP (Peters et al. 2016; Bühlmann 2020) that assumed invariance of P( ( iii IC (Rojas-Carulla Carulla et al. 2018) that extended the above assumption to features beyond Pa(Y); ;; iii) Anchor regression (Rothenhäusler
ct al. 202I) that interpolated between ordinary least square (LS) and causal minimax LS, and
constrained the residue in the anchor subspace to be small; iv) IRM (Arjovsky al. 2019) that

---

Under review as a conference paper at ICLR 2023
learned invariant representation to transfer; v) HRM (Liu al. 2021 that extended IRM to the case with unknown domain labels, by exploring the heterogeneity in data via clustering; vi) IB-IRM (Ahuja al. 202T) that leveraged the information bottleneck regularization to supplement the invariance principle in IRM; and vii) Graph Surgery estimator Subbaswamy et al. (2019) that used validation's is loss to identify the optimal subset.
Implementation Details. We leave implementation details to Sec. in the appendix.
4.1 SIMULATION STUDY
Data Generation. We follow Fig- to generate Xs = (X1, 1X,,X3 X2, Xa}
and XM := {X1} via the following structural equations: Z3 f ug
with Ug N(-2, 1); x2 (xx)) + u2 with u2 N(0, 1): y
X3
X
Y
By (x2) + My with My ~ N(0, 1); 24 < (V) + u4 with 8g = e
varied in different domains, u4 N(0, 1); X1
91 (14, 9) + u1 with #1 ~ ((0,0.2) We consider three settings: i) 92 (23) = 0.52g.
X4
(1(-) is changed to (22n ga(y) -0.55 + y. and 91 (x4, y) = (1g + y)- For each setting, we generate 10 environments with 1, 1..., 10, where ne
(y(22) - (y) y, and ø (I4,, -x;; ii)
x4;
In ( x4 2.5g; iii) (2(2x) 10sinc(xs).. (x2)
Figure 4: DAG with
XM == {X4]. Dotted arrow exists in setting-2.3.
200 for each environment. To remove the effect of randomness, we repeat 10 times, and each time we randomly select five domains for training and the rest for testing.
Causal Discovery and Complexity Analysis. We use F1 score, precision, and recall in terms of directed edges, to assess our causal discovery algorithm. We repeat 10 times and the average results are F1 0.99, precision
1.00, and recall 0.98, which suggests the validity of our algorithm. Validations on larger graphs are left in the appendix. As for complexity, in setting-1, the condition in Thm. 3.1 holds and we expect (X1, X3) to be the optimal set. In setting-2.3. the condition is violated and we need to compare h of each equivalent class to find the optimal subset. There are seven equivalent classes, as the only equivalent relation {X2} ~G {X2, X3).
Table 1: Mean Squared Error (MSE) simulation data.
Setting-1
Setting-2
Setting-3
Predictor
h" (S-)
max MSE
h" h*(S-)
max MSE
h" h'(S))
max MSE
Vanilla
1.90
.07
1.72
5.2442888
4.25
5.97
4.12#
4.27 #.01
3.75
H XX
5.20
4.614.u5
.0003 do
.0075+.00
6.84
4.42+1..15
Xa
1.084.03
1.21 .66
1.372.57
1.132.07
1.16
1.10
(Xa)
3.564.06
1.01
5.2642.51
3.51
3.10=2.30
3.11
(X1,2
1.084.03
1.214.06
.064.00
.074.00
1.51
1.20
(X,3)
3.554..6
4.00
.01
.024.00
6.86
4.39
(2.3)
1.064.03
1.194.06
1.23+.26
1.14#.08
1.17
1,10=.05
(X1,23)
1.064.08
1.18#.06
.064.00
.07
1.17
1.21
Results of other baselines are left in the appendix.
Results Analysis. In Tab. [ we report the estimated h" and the maximal mean squared error (MSE) over test sets of each subset S_ and the vanilla regression method. As shown, in setting-1, the whole stable set enjoys the minimal max MSE, which agrees with Thm. 3.T] while setting-2,3, the subset (fX,) } in setting-2, {X2} in setting-3) identified by minimal h- has minimal max MSE. This suggests the effectiveness of our method in finding the optimal subset. Besides, we observe that h estimate the maximal MSE well in most cases, e.g., ht" h*((XX ((X,, X2)) 0.06 0.07 of max MSE in setting-2, h²" ((Xa)) = 10 vs 3.11 of max MSE in setting-3. In addition, equivalent subsets have similar performance, e.g. max MSE of {X2) and {X2,, X33 setting-3 are both 1010, which verifies the searching can be conducted only over equivalent space.
4.2 REAL-WORLD APPLICATIONS
Datasets. We consider the Alzheimer's 's Disease Neuroimaging Initiative (ADNI) et al. (2010) dataset for Alzheimer's Disease diagnosis, and the International Mouse Phenotyping Consortium (IMPC) CRM workshop 2016) dataset for gene function prediction.
7

---

Under review as a conference paper at ICLR 2023
ADNI. The dataset includes 757 patients enrolled in ADNI-GO/1/2 periods. We apply the Automatic Anatomical Labeling atlas [TzuurioMazzoyer et al. 2002) and region index Young et al.) (2018) to partition the whole brain into 9 brain regions: frontal lobe (FL). medial temporal lobe (MTT), parietal lobe (PI.), occipital lobe (OI.), cingulum (CIN), insula (INS), amygdala (AMY), hippocampus (HP), and pallidum (PI.). In addition to brain region volumes,
also include gender (GED) and genetic information (number of ApoE-4 alleles (ApE)). With these covariates, we predict the Functional Activities Questionnaire (FAQ) score Y for each patient. We split the dataset into seven environments according to age (age <60, 60-65, 65-70, 70-75, 75-80. 80-85, >8), which respectively contains ne 27, 59. 90, 240, 182, 117, 42 samples. We repeat 15 times, with each time four domains are randomly selected for training and the rest for testing.
IMPC. The dataset contains the hematology phenotype of both wild-type and mutant mice with 13 kinds of single-gene knockout. To predict the function of the target gene, we knock it out and assess the cell counts of monocyte (MON), with cell counts of neutrophil (NEU), lymphocyte (LYM), eosinophil (EO), basophil (BA), and large unstained cell (LUC) as covariates. Each environment corresponds kind of gene-knockone We repeat 45 times: each time five randomly picked gene knockouts and the wild-type are selected as training sets and the rest eight kinds are left for testing.
GED
CINC
EO
PAL
FL
INS
PL
NEL
AMY
MTL
LYM
OL
HP
MON
BA
Mutable variables
FAO
Stable variables
App
Target
(b)
Figure 5: Learned causal graphs on (a) ADNI and (b) IMPC. ++ to denote undirected edges.
Causal Discovery. We implement our causal discovery algorithm in Sec. 3.2 to learn PDAG respectively shown in Fig. 55 (a.b). For ADNI, Fig. (a) shows that the affection of AD, measured by the FAQ score, firstly shows in the medial temporal lobe (MTL) and the hippocampus (HP), then propagates other brain regions, which echos existing studies that MTL. and HP are early degenerated regions (Barnes et al. 2009 Duara et al. 2008). Besides, we observe that the frontal lobe (FL), pallidum (PAL), and hippocampus (HP) are mutable regions, which agrees with heterogeneity across different age groups found in existing studies Cavedo et al. (2014); Fiford et al. (2018). For IMPC, Fig. (b) shows that the monocyte (MON) affects the number of lymphocytes (LYM) and large unstained cells (LUC), which reflects the activation mechanisms of LYM (Carr et al.) 1994) and LUC (Lee et al. 202T). It also plays a role the activation of basophil (BA) through the causal chain MON LYM + BA, as also found in existing study (Goetzl et al. 1984).
Complexity Analysis. On both ADNI and IMPC, the condition in Thm. 3.1 is violated, as Y (FAQ on ADNI, MON on IMPC) points to K (MTL ADNI, LYM IMPC). So, need to search over g-equivalent classes and compare their ht,, as suggested by Thm. 3.3 The numbers of g-equivalent classes are 98 (out of the 2 256 subsets) on ADNI and 12 (out of the 24 16 subsets) on IMPC.
Results Analysis. Fig. D (a) and Fig. (a) report maximal MSE of our method and baselines. As see, our methods significantly outperform the others (7.8% on ADNI, 9.7% IMPC). Besides, our sparse optimization is comparable to the searching method in Alg. ! These results demonstrate the utility of Thm. 3.3 in identifying the optimal subset, well the effectiveness of
8

---

Under review as a conference paper at ICLR 2023
25%
24.2%
Improvement over vanil
20%
Max. Test MSE
12.7%
10%
8%
5%
1.8%
0.0%
Ours Surg(SOTA)
5%
5.2%
-5.0%
-10%
HAM
Predicters ranked method and surgery estimator
Figure 6: Results on IMPC (a) Maximal MSE over test environments. (b) Maximal MSE of predictors that are ranked in ascending order from left to right, respectively according h" of method and the validation's s loss of the graph surgery estimator.
sparse optimization when the searching cost is expensive. As for advantages over the graph surgery estimator, Fig. (b) and Fig. (b) show that our h" well reflects the worst-case risk, it increases with the worst-case risk; while there no such property for the validation's 's loss in the surgery estimator. Particularly, the top subsets ranked by our h" also have the minimal max MSE; while the one selected fails to identify S* The improvements over ICP, IC, IRM, and their extensions are due our utilization of stable information beyond causal ceaturee/reerreenntationss The advantage over Anchor regression may be due to the relaxation of the linearity assumption.
In addition, Fig. shows that subsets the same equivalent class have similar maximal MSE, which further shows the validity of searching over only g-cquivalent classes in Alg,
1.10
subsett a. (MT)
2.2
subsots - (NEV, LYM, EO) subseess (MEULYYM
ANY)
Max. Test MSE
subsets " [OL, AMY)
Max. Test MSE
2.0
sutsets - (LYM,EO)
subsets (PL)
subsets e (LYY)
1.08
1.8
1.6
1.06
1.4
1.04
1.2
1.0
1.02
i
2
4
Equivalent Classes
Equivalent Classes
(a)
(b)
Figure 7: Max MSE of predictors in the same equivalent class (a) ADNI and (b) IMPC datasets.
5 CONCLUSION
In this paper, we propose a minimax learning approach to identify the optimal subset of invariance to transfer, in order achieve robustness against dataset shift. Among all subsets of stable information, we provide a sufficient and necessary condition for a subset to be min-max optimal. We analyze the searching complexity by introducing the notion of graphical equivalence and propose sparse min-max optimization algorithm when the searching cost is expensive. The subset identified by our method outperforms others in terms of robustness, on Alzheimer's Disease diagnosis and gene function prediction tasks. In the future, we are interested to study the scenarios where the DAG is also allowed vary across domains, which may happen when the number of environments large.

---

Under review as a conference paper at ICLR 2023
REFERENCES
Kartik Ahuja, Ethan Caballero, Dinghuai Zhang, Christopie Gagnon-Audet, Yoshua Bengio, Ioannis Mitliagkas, and Irina Rish. Invariance principle meets information bottleneck for out-of- distribution generalization. Advances in Neural Information Processing Systems, 34, 2021.
Constantin F Aliferis, loannis Tsamardinoss and Alexander Statnikov. Hiton: novel markov blanket algorithm for optimal variable selection. In AMIA annual symposium proceedings. volume 2003, pp. 21. American Medical Informatics Association, 2003.
Martin Arjovsky, Léon Bottou, Ishaan Gulrajani, and David Lopez-Paz. Invariant risk minimization. arXiv preprint arXiv:1907.0288 2019.
John Ashburner. A fast diffeomorphic image registration algorithm. Neuroimage, 38(1):95-113,
2007.
Josephine Barnes, Jonathan W Bartlett, Laura A van de Pol, Clement T Loy, Rachael I Scahill, Chris Frost, Paul Thompson, and Nick C Fox. A meta-analysis of hippocampal atrophy rates in alzheimer's disease. Neurobiology of aging, 30(1 1): 2009.
Peter Bühlmann. Invariance, causality and robustness. Statistical Science, 35(3):404-4260 2020.
M Woldemar Carr, Stephen J Roth, ED Luther, Shayla S Rose, and Timothy A Springer. Monocyte chemoatrractant protein acts as t-lymphocyte chemoattractant. Proceedings of the National Academy of Sciences, 91(9):3552-3656 1994.
Enrica Cavedo, Michela Pievani, Marina Boccardi, Samantha Galluzzi, Martina Bocchetta, Matteo Bonetti, Paul M Thompson, and Giovanni В Frisoni. Medial temporal atrophy early and late-onset alzheimer's s disease. Neurobiology of aging, (509))2004-200011 2014.
Diego Colombo, Marloes H Maathuis, et al. Order-independent constraint- based causal structure learning. J. Mach. Learn. Res., ((1))3341-37822 2014.
CRM workshop. Causal inference challenge. Website, 2016. htttp://ww. crm. umm.tttterrral
ca/201 6/Genet cs 6/ 6ccommptiitionnee.ppll
R Duara, DA Loewenstein, E Potter, J Appel, MT Greig, R Urs, Q Shen, A Raj. B Small, W Barker, et al. Medial temporal lobe atrophy on mri scans and the diagnosis of alzheimer disease. Neurology, 71(24):1986 2008
John C Duchi and Hongseok Namkoong. Learning models with uniform performance via distribu- tionally robust optimization. The Annals of Statistics, 49(3):1378-144) 1406, 2021.
Frederick Eberhardt and Richard Scheines. Interventions and causal inference. Philosophy of science, 74(5):981-995, 2007.
Paul Erdós, Alfréd Rényi, et al. On the evolution of random graphs. Publ. Math. Inst. Hung. Acad. Sci, 5(1): 1960.
Cassidy M Fiford, Gerard R Ridgway, David M Cash, Mare Modat, Jennifer Nicholas, Emily N Manning, lan B Malone, Geert Jan Biessels, Sebastien Ourselin, Owen T Carmichall, et al. Patterns of progressive atrophy vary with age in alzheimer's disease patients. Neurobiology of aging, 63: 22-32, 2018
AmirEmad Ghassami, Negar Kiyavash, Biwei Huang, and Kun Zhang. Multi-domain causal structure learning in linear systems. In Proceedings of the 32nd International Conference on Neural Information Processing Systems, pp. 6269-6279, 2018.
EJ Goetzl, DW Foster, and DG Payan. A Immunology, 53(2):227, 1984.
factor from human lymphocytes.
Dan Hendrycks, Steven Basart, Norman Mu, Saurav Kadavath, Frank Wang, Evan Dorundo, Rahul Desai, Tyler Zhu, Samyak Parajuli, Mike Guo, et al. The many faces of robustness: A critical analysis of distr--disttiuttion generalization. In Proceedings of the IEEE/CVF International Conference Computer Vision, pp. -349, 2021.
10

---

Under review as a conference paper at ICLR 2023
Biwei Huang and Kun Zhang. Specific and shared causal relation modeling and mechanism- based clustering. In Proceedings of the 33rd International Conference on Neural Information Processing Systems, 2019.
Biwei Huang, Kun Zhang, Jiji Zhang, Joseph Ramsey, Ruben Sanchez-Romero, Clark Glymour, and Bernhard Schölkopf. Causal discovery from beccrogenneous/noostttiinnnr data. Journal of Machine Learning Research, 21(89):1-553 2020.
Chendi Huang and Yuan Yao. A unified dynamic approach to sparse model selection. In International conference on artificial intelligence and statistics, pp. 2047-2055. PMLR, 2018.
Maximilian Ilse, Jakub M Tomczak, Christos Louizos, and Max Welling. Diva: Domain invariant variational autoencoders. In Medical Imaging with Deep Learning, pp. 322-348. PMLR, 2020.
Lucy Eunju Lee, Jung Yoon Pyo, Sung Soo Ahn, Jason Jungsik Song, Yong-Beom Park, and Sang- Won Lee. Clinical significance of large unstained cell count in estimating the current activity of antineutrophil cytoplasmic maibodd-associatte vasculitis. International Journal of Clinical Practice, 75(10)=e 14512, 2021.
Jiashuo Liu, Zheyuan Hu, Peng Cui, Bo Li, and Zheyan Shen. Heterogeneous risk minimization, In International Conference on Machine Learning, pp. 6804-6814. PMLR, 2021.
Chaochao Lu, Yuhuai Wu, Jose Miguel Hernández-Looee and Bernhard Schölkopf. Nonlinear invariant risk minimization: A causal approach. arXiv preprint arXiv:2102. 12353, 2021.
Divyat Mahajan, Shruti Tople, and Amit Sharma. Domain generalization using causal matching. arXiv preprint Ttt1:20000075000 2020.
Ann M Mayo. Use of the functional activities questionnaire in older adults with dementia. Hartford Inst Geriatr Nurs, 13:2, 2016
Jovana Mitrovic, Brian McWilliams, Jacob Walker, Lars Buesing, and Charles Blundell. Repre- sentation learning via invariant causal mechanisms. In International Conference on Learning Representations (ICLR), 2021.
Joris M Mooij, Sara Magliacane, and Tom Claassen. Joint causal inference from multiple contexts. Journal of Machine Learning Research, 21(99);1-108, 2020.
Jens Müller, Robert Schmier, Lynton Ardizzone, Carsten Rother, and Ullrich Köthe. Learning robust models using the principle of independent causal mechanisms. arXiv preprint arXiv:2010.07167. 2020.
Violeta Muñoz Muñoz-Fuentes, Pilar Cacheiro, Terrence F Meehan, Juan Antonio Eguilar-Pimente Steve DM Brown, Ann M Flenniken, Paul Flicek, Antonella Galli, Hamed Haseli Mashhadi, Martin Hrabě de Angelis, et al. The international mouse phenotyping consortium (impe): functional catalogue of the mammalian genome that informs conservation. Conservation Genetics, 19(4):995-1005, 2018
Sahand N Negahban, Pradeep Ravikumar, Martin Wainwright, and Bin Yu. unified framework for high- nigh-dimensional analysis of m-estimators with decomposable regularizers. Statistical science, 27(4):538-55) 2012.
Wojeiech Niemiro. Asymptotics for M M-Estimators Defined by Convex Minimization. The Annals of Statistics, 20(3): 1533, 1992. doi: 10. 1214///////1 176348782, BRLhhtts:///oii /rg/
10.1214/aos
Stanley Osher, Martin Burger, Donald Goldfarb, Jinjun Xu, and Wotao Yin. An iterative regularization method for total variation-based image restoration. Multiscale Modeling & Simulation, 4(2):460- 489, 2005. doi: 10. 10.1137/040604122 URL /ttttts:://oi.. .or. http:///011117700000511
Stanley Osher, Feng Ruan, Jiechao Xiong, Yuan Yao, and Wotao Yin. Sparse recovery via differential inclusions. Applied and Computational Harmonic Analysis, 41(2):436 469, 2016. ISSN 1063- 5203. doi: www.i/n//rrriee////////////////e com/ scis
URL www sciencedirect 2000000366600838880 Sparse Representations with Ap- plications in Imaging Science, Data Analysis, and Beyond, Part II.
11

---

Under review as a conference paper at ICLR 2023
Judea Pearl. Causal diagrams for empirical research. Biometrika, 82(4):669 (888,,66--68,, 1995.
Judea Pearl. Causality. Cambridge University Press, 2009.
Jonas Peters, Peter Bühlmann, and Nicolai Meinshausen. Causal inference by using invariant prediction: identification and confidence intervals. Journal of the Royal Statistical Society. Series (Statistical Methodology), pp. 947-1012, 2016.
Ronald Carl Petersen, PS Aisen, Laurel A Beckett, MC Donohue, AC Gamst, Danielle J Harvey, CR Jack, WJ Jagust, LM Shaw, AW Toga, et al. Alzheimer's 's discase neuroimaging initiative (adni): clinical characterization. Neurology, 74(3)2001209, 2010.
Wojciech Rejchel. Lasso with convex loss: Model selection consistency and estimation. Communica- lions in Statistics-Teeory and Methods, (5(():1989-20)) 2016.
Mateo Rojas- Carulla, Bernhard Schölkopf, Richard Turner, and Jonas Peters. Invariant models for causal transfer learning. The Journal of Machine Learning Research, (1t):1309-13422 2018.
Dominik Rothenhäusler, Nicolai Meinshausen, Peter Bühlmann, and Jonas Peters. Anchor regression: Heterogeneous data meet causality, Journal of the Royal Statistical Society: Series (Statistical Methodology), 83(2):215 2222246,,666,-- 2021.
Shiori Sagawa, Pang Wei Koh, Tatsunori Hashimoto, and Percy Liang. Distributionally robust neural networks for group shifts: On the importance of regularization for worst-case generalization. arXiv preprint orXiv: 191.00831, 2019.
Peter Spirtes, Clark N Glymour, Richard Scheines, and David Heckerman. Causation, Prediction, and Search. MIT press, 2000.
Adarsh Subbaswamy and Suchi Saria. I-spec: An end-to-end framework for learning transportable, shift-stable models. arXiv preprint arXiv:2002.08948, 2020.
Adarsh Subbaswamy, Peter Schulam, and Suchi Saria. Preventing failures due to dataset shift: Learning predictive models that transport. In The 22nd International Conference on Artificial Intelligence and Statistics, pp. 3118-3127. PMLR, 2019.
Nathalie Tzourio-Mazoyer Brigitte Landeau, Dimitri Papathanassiou, Fabrice Crivello, Olivier Etard, Nicolas Delcroix, Bernard Mazoyer, and Marc Joliot. Automated anatomical labeling of activations in spm using a macroscopic anatomical parcellation of the MNI MRI single-subject brain. Neuroimage, (5()):733-28, 2002.
Chuanlong Xie, Fei Chen, Yue Liu, and Zhenguo Li. Risk variance penalization: From distributional robustness to causality. arXiv preprint arXiv:2006.077 20040440000 2020.
Alexandra L. Young al. Nature Communications, 9(1):4273, Oct 2018. ISSN 2041-1723, doi: 10.1038/s41467-011 0016677-----009--00
Jiji Zhang. On the completeness of orientation rules for causal discovery in the presence of latent confounders and selection bias. Artificial Intelligence, 172(16-17)1183--188) 2008.
Peng Zhao and Bin Yu. On model selection consistency lasso. The Journal of Machine Learning Research, 7:2541 -2563, 2006.
12

---

Under review as a conference paper at ICLR 2023
REPRODUCIBILITY STATEMENT
Data, code, and instructions to reproduce the main experimental results are provided. Specifically, the ADNI data is available at hlttp://adni.. loni wuw.alllaaii.. the IMPC data set ava- iable at t/http://ww www. crm. umontreal. ctt.caaaac... 2et//////////eiieeeeeee seenetics 6/ competition php the codes are provided in the supplementary materials, the implementation instructions are provided in Sect. in the appendix.
APPENDIX
A Related Works
14
Appendix for Sec. 3.1; Identification with Min-max Property
14
B.1
Proof for Thm. 3.1 Graphical Criterion for f
14
B.2
Proof for Prop. 3.2: Testability of Thm. 3.1
17
B.3 Counter Example of f* + Js
19
Proof Thm. 3.3; Min-max Property
20
B.5
Proof for Prop. 3.4; Identifiability in Thm. 3.3
21
Appendix for Sec. 3.2;
Learning Method
22
C.1
Equivalent Classes and its Recovery Algorithm
22
Details of Causal Discovery to Detect Local Components
................................
27
Details of Estimating fs.
27
D Appendix for Sec. 3.3 Complexity Analysis
28
Appendix for Sec. 3.4; Sparse Min-max Optimization
33
E.
Statistical Consistency
33
Linearized Bregman Iteration
35
Appendix for Sec. 4 Experiment
36
Simulation
36
Alzheimer's s Disease Diagnosis
37
Gene Function Prediction
13

---

Under review as a conference paper at ICLR 2023
A RELATED WORKS
Causal learning for domain generalization. There have been emerging works that consider the
domain generalization problem from causal perspective. One line of work AArjovsy et al. (2019);
Xie et al. (2020); Müller et al. (2020) promoted invariance as a key surrogate feature of causation
where the causal graph more of motivation. Another line of work lise et al. (2020) [ et al. (2021): Mahajan et al. 20:0); Mitrovic et al. (202T) considered domain generalization for unstructured data using specifically designed causal graphs to incorporate priors of the distribution shift, in which the causal features are modeled as latent variables to be inferred for robust prediction. The works most relevant to us pursued robust optimization by making invariance assumptions regarding causal mechanisms (Subbaswamyy et al. (2019); Bühlmann (2020): Peters et al. (2016); [Subbaswamy & Saria (2020), Specifically, the Peters et al. (2016) assumed the generation of Y from its parents is invariant; hence they only utilize Y's s parents for transfer. The The/Subbaswamy et al. (019);; Subbaswamy & Saria (2020) considered a selection diagram framework, in which mutable variables are children of the selection variable. They then remove the unstable mechanism by intervening of XMs and obtain set of stable covariates S. To obtain the optimal subset, they simply search over all subsets in S and took the one with the minimal validation loss. However, their method ís theoretically incomplete and practically defective, as the selected subset may not be min-max optimal and the searching cost expensive in large-scale graphs. In contrast, we provide a comprehensive min-max analysis to guarantee the identification of the optimal subset. For practical employment, we analyze the searching complexity via the lens of g-equivalence. For those graph with expensive searching cost, we provide sparse min-max optimization scheme that can larger improve the efficiency.
Optimization-base domain generalization. There are emerging works that view the domain of
generalization problem as an optimization problem. These methods directly formulate the objective
generalization and conduct optimization for min-max optimum. For example, Distributional Robust Optimization (DRO) Duchi & Namkoong (202T) constrained the distance between test and training distribution with f-divergence or Wasserstein distance and optimized over the min-max objective. One of its popular extensions, GroupDRO Sagawa et al. (2019), provided extra regularization (e.g., weight-decay or early stop) and allowed DRO models to achieve better performances in large neural networks.
However, these methods heavily rely on data- data-n optimization and lack analysis of the source of distributional shifts. For this reason, they have constrain the distributional shifts limited extent, so as to achieve optimization convenience. Such a limitation affects their ability to generalize to broader distribution family, thus limiting their applications in the real-world. In contrast, we consider domain generalization from a causal perspective. Benefiting from the causal framework for distributional shifts, our method can identify the reasons behind distributional shifts and achieve min-max optimum even when the distribution can vary arbitrarily.
Causal discovery in heterogeneous data. Our work also benefits from the recent progress in causal discovery Mooij et al. 2020):1 Huang et al. (2020); Huang & Zhang (2019); Ghassami et al. (2018), a field that focuses on identifying the causal relations from data. Our work shares similar framework with Huang et al. (2020) in formulating the distribution shift. However, they focused on recovering the full causal graph to study relations among variables, we provide a local discovery procedure, which aids the analysis of min-max properties and identification of robust predictors.
B APPENDIX FOR SEC. 3.1 IDENTIFICATION WITH MIN-MAX PROPERTY
B.1 PROOF FOR THM. 3.1 GRAPHICAL CRITERION FOR f" = ffs
Theorem 3.1 (Graphical Criterion for f* = fs). Suppose assumptions 2.2 X Ar :o;:=xx X CC(YY) as mutable variables in Y's children, and K : De( X! of x Then, plylag, do(xxs) can degenerate to conditional distribution if and only Y does not
2.2 hold. Denote M XX M as descendants
point to any member of K. Further, under either of the two conditions, we have f fs.
Proof. Denote the causal DAG as G, the intervened graph that removes all arrowheads into V as Gy. Define X = Xx,, K₂, == (X xx)) (ee)))) )
We firstly prove the equivalence of the following conditions (1), (2), and (3):
14

---

Under review as a conference paper at ICLR 2023
1. Y La x]i
K 22
2. Y and K adjacent in G
3. ply las, ()( can degenerate to conditional distribution.
(1)+(2) Y and K are adjacent, they are also adjacent in G "xy because KnX 2, so Y and
K can not be d-separated by any variable G Gxp. which contradicts with (1).
(2)+(3) Define
I = plykk kk k2, do(x)))) = p(ylpa(y)) exx plz))pa(z))) p(x.(pa(ax))
p Ix p(a, (p((x))) Ixxe p(1; (p(())))
Since PA(Y) [X] K} 2 and AX, E K2 PA(X.) (X M: K} Ø, we have:
plu, k2) II Plzz)lppai (pa(x;))
ply k2) II 11X,eK p(u, pa(r;)d
If and K are not adjacent, then VX, E K. & PA(X)). Therefore, I=
ply.ks) psy,
(k2).
(3)+(1) We will prove by contradiction. Specifically, we will show that i Axxy, K K2 ie., (1) does not hold, then pf (s do(ar)) can not degenerate to any conditional distribution, i.e., (3) does not hold.
We firstly show Y Laxy, K|K2 + pulass do(x))) * of (y k2, do(x))))) then show yf (ylls, do(#r))) + 1f (v)k2, do(x))s) + pe do(xx)) can not degenerate to any condi- tional distribution.
Since Y & (.)() we have:
1d (v, as, wdrldo(
(yllss, who do(x))))
I pe ms w/y do (2 Jyy
(vllpa(y) Il p" (a.ppa(r:)) II, xxxx pf (xv(pa((x))) S (ylpal ()) Ph (x:)(())) [ exx,xx (x)(pa()))) ('l Il ( S (v((vpa())) (1:
Since KUK KUK, Xs U X we have ( = p² k. kg do(z).)) Thus, we can prove: Y LG KKK, = p" (y æs. do( p' k. k2, do(x)))))) # p do(
Next, prove of (yj æs, do( + pr k2) do(x)))) + Pr (y) 2s do( can degenerate to any conditional distribution.
Suppose pf l do(#)s)) of ( |k' R2, do(x)))) ())) We will show if k' * Ø, then the do- operator can not be removed with either Rule 2 (action to observation) or Rule 3 (delction of action). To express do(x)) explicitly, denote X (X and (yyy'' kg, do(x),))
p" kk, k2, do(x) (0...........,. do(2)
Rule 2 can not remove the do-operator of any X X
Note the edge hetween Y and X K only he Y X.
15

---

Under review as a conference paper at ICLR 2023
Recall Rule 2 states that "p(y)do(a) do(z), w) p(y)do(x), Y -Gxz Z W for any
disjoint subsets of variables X, Y, Z, and W
Rule 2 can remove the do-operator of X Md,
E
Xx then
Mr
1a X% (x2,,
M, K2, X {XM,)
(4)
As we have Z (XM,) X XAAM, W K' K2 in the notations of Rule 2.
In the following, we explain why Eq. can not be true. Note XM, E Ch(Y) and the direct edge
Y + XM, M, is reserved in the intervend graph G X M JxS which means that Y and X Mr can not be d-separated by any set of variables in the intervened graph. Thus, Eq. can not be true.
Rule 3 can not remove the do-operator of all XM, E XM..
Recall Rule 3 states that "p(y do (x), do(z), w) = ply do())) w) if Y 1Gx z(w] Z,, w for any disjoint subsets variables X, Y, Z, and W Rule 3 can remove the do-operator X
Mr then:
La (K
X K' U22
because the notations in Rule 3 mean X = 2, Z X
W K'U K2. In the following, we will
M.
show that when K' * Ø, Eq. 5 can not hold. When K'
* Ø, note by definition K' C De (x)r
(r).
X Ar {An(K') U 32} * X Mr. That is X X ( U 22) * Ø. Suppose XM, E X X (K' K2), then the edge + + is in
so An (K') n X M
* Ø. Therefore, X Mr (K' U K2)
the intervened graph G (KUK)) so Y and X can not be d-separated by any variable set. So
Eq. does not hold.
In summary, we have proved that when K' + Ø, the do-operator on X can not be removed entirely by Rule 2 and 3.
Besides, according to Corollary 3.4.2 in Pearl (2009), the inference rules are complete in the sense that if the intervention probability (with do ) can be reduced probability expression (without do ). the "reduction" can be realized by a sequence of transformations. each conforming to one of the Inference Rules 1-3. Note that only Rule and are related to the disappearance of do-operator, is sufficient to prove that Rule 2 and 3 can not remove the do-operator on X M-
We then prove under either of conditions (1), (2), (3), f" fs-
Given any one of the three conditions (1), (2), or (3), (x) Ep [ do(xas)) satisfies the following min-max property:
f (x)) argmin maxpep EplY f(x))
Under any one of the conditions (1)-(3), we have pf (yl as, do( of k2) for Pe EP. For pe P. let of (2:4)= Extevixg of (v) be the marginal distribution of X9 Define Pe as:
(r(w)) = ((n)=p(ylpatff (ylpa((y)
II p" (x.lpa(zz)) X,EK
II D² (z.lpa(z:)) (# X,eK,
by replacing the term II Ixx,ex, (a.1pa(x:)) in (*()) with pa (z1r).
(i) By the definition of P,, P e P and "t (y/z) (3) thy x21) (ylas, do do(x);))
(" k2)
(ii) In the following, we will show 5e (y. k2) = pe (y, k2)
Z () is the set of Z-nodes that ancestors of any W-nde in Gx.
16

---

Under review as a conference paper at ICLR 2023
First, note that K De (Xgr and X Ch(Y), we have K KUXi De(Y). Thus, PA(Y) n
[K U { 2 because otherwise there would be a cycle.
Second, (2( (KUX%,
X } because there exist X; E K K X and also X; E PA())). then K2 n De (XM) A Ø, which contradicts with the definition that K₂ = (x( /ee (X),.
In summary, we have PA(K2 UY) (KUN) 20 which leads to (
pº (k2. p(olpaty))lxx:kkpp)) (silpa(zss) ux (ailpa(x:)) 11 exx, (silpa(ax)) do UxxeK (ailpa(a;)) (xxex, (Iilpa(<;;) de - (x:lpa())) /
(xilpa(zi))
-(ylpo(y)) (zilpa(a:))
and
pº () p(ylpa(y)) (zalpa(xx)) Hxeep (1:1pa(zz)) (a/s) dr
p(ylpa(y)) (aypp((x))) / Uxxekr (r:lpa(zs)) (ale) da dk plulpa(y))llxxxekpp) (a (x)lpa(x;))
Therefore, we have 5r (k2, 8) = (k2,
Note that K2 C X, we have
Varp= (Y(K2) = Ep- [Var (Y|X) + Varp= [Ep= (Y|X))k))
therefore
Epe [arre (2]K2)] Ere (Y|x) Ere [Varre ()((())K.))
and hence Epe [arrp (Y|K))] Ep- (Y
(iii) Because 80 (k2, y) = pd (k2, v) Ep. [Vr pe K2)] = Epe [Varpe (2(K))), we
have Ep. [Varp. (2]
2 Ep- [Var (Y] Besides, since (3((x))) = Ør (ylk2).
Ep. [Var (Y| Ep [Var (Y] 2 Ep- [Varp= (Y [X])..
(iv) In summary, for each pe e P, we may construct Pe such that
Ep. (x 2 Epe pa (Y])))
Denote p { p pe E pP] and p argmaxpep Ep [Varp [X), then P E P. Besides,
note that for any pe E -PEpp.YYa
Eyc (Y vl do (x\r) ) = Ep. (Ylss do(x }} so
(x)x) = P--- Y| xc] (V-lsss do) As Ep- (yp[(V do (x is invariant for all Pe E P.
Therefore, have f*(() Ep- læs, do(xr)))
B.2 PROOF FOR PROP. 3.2: TESTABILITY OF TIM. 3.1
Proposition 3.2. Denote X9, == XM n Ch(Y) as mutable variables in Y's children, and K == De Mr) X as descendants of X M. Under assumptions 2.7 2.2 the K is identifiable; besides, we can determine whether Y + K from the joint distribution over training domains.
Proof. We firstly show K is identifiable. Since all variables in K are descendants of Y. have + X, E K iff X; is adjacent to in the skeleton of DAG (which is identifiable under assumption 2.2). Thus,
determine whether Y K.
17

---

Under review as a conference paper at ICLR 2023
Algorithm 2 Detection of XM and construct the causal skeleton of G
Start with XM = Ø. For E V, test V; T or there exist subset Co.. V such that Vi 1 E v,, If V, y E and there exists such C,,,e, then include Vi, XM =X U V,
2. Start with undirected graph G including edges for any variables in V and the arrows E + Vi for VÚ XM- For each pair of {Vi, Vi If V, 1 V; or there exists a subset Cer.ey CV such that V. VylCon.ey we delete the edge V, V, from G.
Note that K - De n ( } So it suffices to prove the identifiability of X U De(X) M) where Xh, XM Ch(Y). This can be accomplished by three steps: (i) identification of XM, (ii) identification of X and (iii) identification of X U De(X)))
The following algorithm shows step (i), which is the same as Huang et al. (2020).
The following Alg. 3 shows the steps (ii) and (iii), which basically relies on the faithful assumption (conditional independence probability = d-separation in graph).
Algorithm 3 Detection of X and X De(
2.1 Detect X XM Ch(Y):
I: for X; E X and adjacent to Y do
2: If y C, U (X;
3: then X; E XM Ch(Y)
4: end for
2.2 Detect {Ch(Y) n 3 De(C
1: Start with A = B 3- BCCh(Y) n XM and visited(X,) =FALSE
2: while B * Ø do
3:
for X3 B do
4:
for X; E X do
5:
if X, XM and X; 1. F. Ce,r, (X,) Drçe then
6:
A= (X visited (X.) =FALSE then B B U (
8:
end if
10:
end if
11:
if X, XM and X₁ & Adj(Y) and X, Y|C U (X,) \ Dyay then
12:
A= A [X,} visited (X) =FALSE then B = B (X+}
13:
14:
15:
end if end if end for
16:
17:
18:
Let B B \ (X,)
19:
end for
20:
end while
Explanations for 2.2
Line in 2.2 The set A is the final output. The set B only plays part instrumental set that starts with X n Ch(Y) and ends with Ø. During the process, B stores the nodes in Xas O Ch(Y) that has not been searched for the children. Once X, B is searched, it is excluded from the set B (Line and the children of X, are added to it has not been visited (Line S and I4]. which is essentially seadth-first-eaae algorithm.
18

---

Under review as a conference paper at ICLR 2023
Line to 10 (the case when X; & XM): The fact X; XM means E and X; are not adjacent. Besides, note that X, E XM n Ch(Y), there is structure in the form E - .... X, X; where X; and E are not adjacent. In the notation X; - E| r., U {X,} De,,, C..,e denotes separating set such that X; E Cr..e and De.z, denotes the set of variables along the directed clear is determined in the greadth-first-searde process.
Xj. The existence of Cr,,e guaranteed since X; and E are not adjacent, so a separating set has been found when constructing the skeleton. The set De.,, is also
path between F.
Line
to 19 (the case when X; Xas): Firstly, we explain why unnecessary to consider the case when X, E XM and X; E A.)(Y). If X, PA(Y), X; can not be in De(Ch(Y) (.XM) U (Ch(Y) X as it would induce a cycle in this way. If X; E Ch(Y), ) it means X; E Ch(Y) n XM M and has been identified in 2.1 and included in set A in the beginning.
So the remaining case is when X; E XMM and X; & Adj(Y). Note in this case Xj E Ch(Y) or Xj X;, which is the same as E
Xj De(Y), there exists structure Y
X X; in Line to TO
B.3 COUNTER EXAMPLE OF f # fs
Y
Xs
Xm
Figure 8: DAG of the counter example.
Consider the DAG in Fig. 8 in which we denote Y, Xs, Xm are binary variables. We will show that is not min-max
in this scenario, there exists P(Y), P(X, me Y) such that Is optimal. We show this by proving that:
E ElY Ksy do( || E [d]]do(mm))))
Since we have that
E Y Y 1x,, do (o( 2 E E E2 x. do( 2E (E(YYr,, do(xm))l,
and that E [Y - E(Y (do(am)) tion 6 is equivalent to that
)>12 E E| due to that ply)do( = p(y), the Eq. equa-
E
E (P (x)))))]]]
2F E[Y do(xm))) F.2 P..
Besides, we have
Σ p( kmm s)p( v)p(smlu)p(t)
E [ Es do( []]
Σ
( do (5(mm)
(8)
E Y E|Y E. do(xm)))l =
Σ p(xz pm v)p(emlv)p(v)
EE |0,, do(sm))
(9)
Fe.
Since we have plylar, do (m))
p(y)p(us, / (m.))) Σ, ( (t) we have
ply = 1)p(x, Ims y = 1)
a do(xm))
(10)
E, ply))(pp (z Jms v)
Substituting Eq. equation 10 into Eq. equation
equation 9 we have
19

---

Under review as a conference paper at ICLR 2023
R do(
E
Σ plana Im, vpp(am vbl(y)
plv 1p(
Σy p(u))(p) amm
Ne Rm
plu 1)p(zz | y
E|y. E|y do(
Σ
Σ p(ara (am 8) (pp())
Σn p (
2ig Pom
ply y= 1) 11y=1)p(ar | 8=
Σn ((y)p(Ts y)
Denote dy p(u 1), pmam =1y)) == aray p( 1 1/5m 9) aemy. then the left hand side in Eq. equation
has
1 [ v,, do(xm = (2 1m = (Ax110mia (0 + a=100mo
az))
0ya.11
+
4y0a11 + (1 ayJacio
1(x, = 11,mm =0) [a.11(1 Jay a=10 amo) an)]
ay
01y1001 + (1 Bg 0,00
ay a.11)
2
1(x. ((c.=0,am. 1) [(1 d.l1)mm1al mmm kyy]
+
Uy(1 Uoll) (1 Uy [1 - as10)
4y (1 - a.01)
1(,, (1,==0,,,, = 0) (11-a_0)(1-0mm)))) + a100000)(11 (1-m00)(1-mmm) (1-a
By - a001) + (1 0g) (1 a.00)
The right-hand side has
2 0y
2E [YE[Y (z., do(a -| =2 1,
0ydss1 + dy (ss0
dg (1 am1
1(x, 1,
0y0s01 + (1 - (gy (3,00
a (1 0.11 )² )20m1
[(x, = 0, 1)
av as11) + (1 dy Uls0
a aov1 (1 am1)
[(x, = 0, Zm 0)
dy a201) + (1 ay) a.o0) ay.
When ag * 0, the term ot can be removed. Then let ag + 0, the left-hand side approximates to:
22 0ma
20 ama)
(1 - 0.11) 2 amo
(1 - ax01 ( amo)
0x10
a.00
(1 as10)
(1 a,00)
and the right hand side approximates to:
2 aml
2301 (1 am1)
+ (1 am11)
+ (1 - a001)² -am1)
2.
0010
00000
(1 as10)
(1 - as00)
Then the Eq. equation is equivalent to:
2: (amo 22m1)
a (20m1 amo 1] (1 (1-as11) ( 2am1)
(1 (--a00)) )2(mm)
1000m-mmm amo -
1.
a=10
0000
(1-0,10)
(1 (_-a,00)
Let 0,10 0, a,11
1, 0,01
0,10 0.5 and amo 2am1 0, the above inequality holds.
B.4 PROOF FOR THM. MIN-MAX PROPERTY
Theorem 3.3 (Min-max Property), Denote ht hS-) == max Ep, [[Y - (x))² as the maximal J
expected loss over J for S-. Then, we have h² (S-)= = max Ep. [[Y - fs_ (x)) 21, In this regard,
PeEP
the optimal subset S for f -sg- can be attained via S argming 12' (S-).
20

---

Under review as a conference paper at ICLR 2023
Proof. We show the maximum loss is attained when XM is definite function of PA(XM
Let fs_ (as_ xM) = , do(x))) be an invariant predictor. Then
/ (u- fs_ (ass ( p(ylpa(y))
Lp (fs
II RE
((=:ppa(a;)
II 1EM
p
And the maximum loss
S = arg = argmax (p (rappa/ddeettapmaeeddeaaaaaaaaad
Let x X (PA((()))))) and h(xM, pa( fs_ (x)) (xx.
which does not rely on the mutable distribution (p (z:lp((()) E M}. Let m (pa(xr)) = argmax22 h( (a( )))-
Firstly, consider the case of XMM = { Then
(ʃ
/
max Lp (Js max pa(a
h(w, pa(as)))) (#m (pa(as))
II X,EPA(XAr)
p(a: p(d::pa((:))dppp(xxx)
base
f hizw pa(aa (ar (pa((aa))
max
II XxEPA(X)r)
p(=:/pa(x:))dpp(x))
pa(#M) M)
h(m" (pa(æem)) pa(atr))
II X,εPA(X)))
p(x:)pa(x;) (ddpaaw))))))
pa(
When XM is multivariate, cal order { [XM.1, X XM2,, "M,,h where XM3 is a node that is not a parent any other nodes in (Xmali j} in the sub-graph XM- That is, we firstly consider (aM. pa(eMa)) Jam. h(FM.1 pa (FMI1 )de.. and factorize max pe (...) as
consider the maximization sequentially by the topologi-
maxp (ew.l (ew..))()..))) ... (pp, a)) maxp (*M,, (pa(su, .
Note that the sub-graph on XM is always DAG, such topological order always exists.
B.5 PROOF FOR PROP. 3.4 IDENTIFIABILITY THM.
Proposition 3.4. Denote PJ ply. æs do M J (pa(a and fs E[Y d)) Under assumptions 2.1 and 2.2, we have PA(XM). PJ. and fs_ are identifiable.
Proof. To generate data distributed Pj, we need then regenerate De(XM) with structural equations. To estimate fs_, we need to inter- X M. then regenerate De(X,) with structural equations. So, it's suffice to show PA(e(X (A(De(Xar)) XM, De(X (M)) PA(XM),
(M)A(XM)) regenerate XMs
are identifiable.
Identification of X has been shown in Alg.
2
We first show the identification of De(XM).
Line 5 to 10: this case is the same as Algorithm 3 which is based on (i) the structure F .....
X, X; and (ii) X, and are not adjacent.
Line 11 to 19 In this case, we identify the direction between X; and X, by the "Independent Causal Mechanism (ICM) Principle" following Huang et al. (2020), where Ax,-x, and A Ax,+x, are the estimated HSIC (see Eq. 17 Huang et al. (202)) for the detailed formulation of A).
The ICM principle means that "the conditional distribution of each variable given its causes (i.e., its mechanism) does not inform or influence the other mechanisms.". That is, the changes of P(X,(PA(X:)) does not influence the other mechanisms P(X, (PA(X)) for j * i. The ICM principle is implied the definition of "structural causal model" Pearl (2009), where each structural equation represents autonomous physical mechanism.
21

---

Under review as a conference paper at ICLR 2023
Algorithm 4 Detection of De(XM)) UXM
1: Start with A = MM and visited(X, =FALSE while B + 2
3: for X; E B
4:
for X; E Ad)(X;) do
5:
if X; $ XM and X, 1 ClCa,, U1X,) Daye then A= {X+}
if visited X2) =FALSE then
8:
B = BUBSUU (X))
9:
end if end
10:
X. E XM and Ax,-x; < -x then
12:
A = A U ) if visited X2) =FALSE then
13:
14:
B = B (X,)
15:
end if
17:
end for
18:
Let B = B (X,)
19:
end for
20:
end while
We then show the identification of PA) PA(D((X.)))
Algorithm 5 PA((,) for X; XM De(M)))
1:
for X, E Xa U De(Xas) do for X; E Adj(X,) do
2:
3:
if X; g X then
4:
X; E PA X, Y El (Dex ((3))}
5:
clse if X; XM and X; Xm then
6:
X, E PA(X,) when Åx -x, Âx-+Xy. else X; XM and X; Xm then
7:
X; E PA(X;) when / C, U (X,}
9:
end
10:
end for
11:
end for
Line this rule is based the structure E + X, Xx and (X,, E) adjacent.
Line this rule is based the HSIC criterion in Huang (020)).
Line this rule is based on the structure E Xx X, and (E, X3} are not adjacent.
C APPENDIX FOR SEC. 3.2 LEARNING METHOD
C.I EQUIVALENT CLASSES AND ITS RECOVERY ALGORITHM
When the graphical condition Thm. 3.1 fails, Alg. [ needs to search over subsets of the stable set and identify the optimal predictor. However, find exhaustive search of all subsets is redundant, as some subsets are equivalent in the sense predicting Y. Formally speaking.
Definition C.1 (p-equivalence). Two subsets X, and X, (the subscript S is omitted for simplic- ity) of the stable set Xs are probabilistical equivalent, i.e., X; ~p X,, if P(Y X, dx)) P(Y X,, ds)X))))
Remark C.2. It's straight forward to ~p satisfies reflexivity (X, up X.). symmetry (X, ~P
X, Xj ~pp X,), and transitivity (X; ~pp X,,X, ~p X = X; ~p Xx), thus is a legal equivalent relation.
22

---

Under review as a conference paper at ICLR 2023
Under the Markovian assumption, we can infer p-equivalence from structure of the causal graph, espe- cially patterns of d-separations. In the following, we firstly give the notion of graphical equivalence, then show how can we infer p-equivalence from it.
Definition C.3 (g jg-equivalence), Two subsets of vertex X, and X, are graphically equivalent w.r.t the causal graph G, ie., X; NG Xj. if :X C X; n.3 such that Y 1 Cxx X; U UXJ Xzy] XXx, XM-
It's straight forward see that G satisfies reflexivity (X, ~g X,) and symmetry (X G X, + X, ~g X²). To prove it also satisfies transitivity (X; ~G X,, X, ~ X + X; ~ Xt), we need to introduce two properties of the d-separation.
Lemma C.4 (Properties of d-separation), (i) If path p can not be blocked by a vertex set X,, then any of p 's sub-path can not be blocked by X; either. (ii) For two vertex sets X,,X, Xj. and a path p. if can not be blocked by X; and can be blocked by X, K, then, X, must contain a non-collider p.
Proof. The correctness of property-(i) is straightforward, so we focus on proving property-(i). Specifically, there are three possibilities about the path p:
1. all vertices in p are non-colliders. From 'p can not be blocked by X,', we know all vertices in p are not in X, From 'p can be blocked by X; U UX',, we know at least a vertex in p is in X,.
2. all vertices in p are colliders. From 'p can not be blocked by X,', we know VX E p, X E X, or X has a descendant in X, So, AX E p, X E Xx; U X5 or X has a descendant in X; U XX, which means p can not be blocked by X; UX neither.
3. vertices in are colliders and non-colliders. From 'p can not be blocked by X,', we know VX E p. if X is non-collider, X & X,, if X is collider, X or one vertex in Dec(X) is in X,, thus in X, U X, So, X, must contain a non-collider in otherwise, p can not be blocked by X; U U,J O
Equipped with the above properties, we now show ~a also satisfies transitivity, ie,, X, ~ X,.X, Xj NG X = X; ~G X
Proof. Because X, ~ X, X, ~g Xg, by definition, EX, C X, n X, such that Y LG X, U X, Xiy [Xij, and 3X₃ jk X, n nx such that Y La X, X X jk X jkk Different situation of Xij and X5k are discussed below:
1. X, Xij= Xo. Then, we have Y Le X, U UX Xx Xo. So, Y La X, X Xo
and X, C X; Xx. So, X; NG Xx.
2. X 1
9 (a), we have X₁h C CX so, Y Le Xix We also have X.y S X, XX jk. so, Y La Xij (i
X jk 2. As shown by Fig.
In the following, we show any path between Y and Xjv contains at least a collider. We prove this by contradiction, iee., assume there is a path po <<< X2, Xm > between Y and X,k (Xm, E xx)) and every vertex in Po is a non-collider in Po- Because Y IG X g.. so, EX,, i < 1 in Po such that X; E X.z. So, there is path p1 :< X2, ...,.." > between Y and Xiy where every vertex a non-collider. Because Y LG Xiy Xyk= ke so again, EX,, i- in p1 such in
that X, E Xz*- Iterating like this, we have either X1 E X ij or X1 E X5₃k Because X1 E Neig(Y), Y LG X1 given any subset, which contradicts with Y X, XX
jk X5* and Y 1 X, (X,y(X
Because any path between Y and X (X,; can be similarly proved) contains at least a collider, we have Y La X;yl@ and Y La X₁x (0.
In the following, we show any path between Y and X, Xxx contains at least collider. We prove this by contradiction, ie., assume there a path Po :< Y, X1, X2,... Xm > between Y and X, (X. (Xm E X, X..) and every vertex in po is non-collider in po- Because Y LG X, XXyy Xaj, so, i i I m 1 in po such that X, E Xij. So, there is a path pl << X2,XX,..., and X. where every vertex is non-collider, which contradicts with Y -G X,y10.
> between Y
Because any path path between Y and X, X.j contains at least a collider, we have Y Lo XAX,le and similarly Y Ia X X 2 Considering Y Ia X,yl and Y La Xyx (8, we now have Y Lo X: and Y La X Ø0. Because Ø C X, n X and Y Oa X, U Xx Keo0, we have X, NG X
23

---

Under review as a conference paper at ICLR 2023
X
X
X, UX, UX
D
y
Xjk
A
B
Xij
Xjx
Xij U Xyk
Xzj n Xjk
Xij n Xjk
(a)
(b)
(c)
Figure 9
3. X, n nX * Ø. We first define A = Xi Xiy n Xyx and B = X X1 O nXji. as shown by Fig. 9 (b))
In the following, we show any path between Y and A can be blocked by X.y n nx We prove this by contradiction, i.e. assume there a path po Y, X1 X2, Xm > between Y and A such that Po can not be blocked by X,; n jjk- Because A C X, X jk: Po can be blocked by X jk. By Lemma
C.4 X, (X,, n jik (i.e. the subset B in Fig. 9 (b)) must contain a non-collider X,, i I m in
So, there is path pı :<< Y, X1, X2, X; > between Y and Xjx Xiy nx jk: Be Lemma C.4 as Pi a sub-path po. P1 can not be blocked by X.j n XJJ neither. Because X; B X, X P1 can be blocked by X1y, again, X1, XX,, nX3x (i..., the subset A in Fig, 9 (b)) must contain
non-collider in p1. Iterating like this, we have X1 E A or X1 E B. Because X1 is adjacent with Y. we have either Y Lo X,) (yrt/Xa or La Xy (X. Xay not hold.
This contradiction means any path between and A (similarly B) can be blocked by X., O Xjk- Formally speaking, Y La Xis X3 Xi X5 jk X,y O jkk.
Define C := X.y U UX3K X nx fk and == X, U UX U U U Kik,, as shown by Fig, (c). We have already shown that any path between Y and C can be blocked by X,y O Xik- we in the following show any path between Y and D can also be blocked by X.j O Xjjk
Again, prove this by contradiction, L.e. assume there is path po Y,X1, 11, X2, Xm > between Y and D such that po can not be blocked by Xi n Xjk- Because either Xm E X, UXJ XX, or Xzk- Let's assume po can be blocked by Xij (as by j, has a similar analysis). By Lemma C.4. Xij' Xix nnn,X, XXj must contain non-collider
Xm X, U X jk. we have po can be blocked by X iz
X2 in Po- So, we have path Pı :< X,X1, X2, ...X; between Y and X,y XXx n X Because P1 is a sub-path of Po. by Lemma C.4 Pi can not be blocked by X.y n X,s either, which contradict with
Y L. X,j U X S Xi O X jkc
N
This contradiction means any path between Y and D can also be blocked by X,j nXj..
In conclusion, any path between Y and X, U X, X X,y OX5x can be blocked by X, nXi jk Formally speaking, Y Le X,UX, X by n jkl X nx jk. Because X,UX, S X, XXUXX,UX UX and X. 0 X₁ C X, O Xx. we have X; ~g Xx-
To conclude, we have shown that ~y satisfies reflexivity, symmetry, and transitivity. So, it is also a legal equivalent relation.
Recall the Markovian assumption states for any disjoint sets X,, X, Xg. we have X; LG X, X, X, |Xx.. It builds bridge from d-separation in graph conditional independence in probability. As a result, under this assumption, we can infer two subsets are equivalent in predicting Y if they are graphical equivalent the intervened graph Gx Formally speaking,
Proposition C.5. For two subsets X; and X, of the stable set, if X4 ~G Xj, then X, ~p X,-
'The intervened graph means the graph after removing all edges into XMC-
24

---

Under review as a conference paper at ICLR 2023
Remark C.6. Note that the reverse claim X; ~p X, X, NG X, is not true under the faithfulness assumption. Consider the counter example of Y, 1, X2, and the structural equations X1 + Y YN 1), Y (0,, 1). We have (X1} ~o {X2), but do not have {X1} ~ {X2).
Now that we know the definition of two subsets being equivalent and how to infer the equivalence from causal graph, we are ready introduce the notion of g-equivalent class. Denote the power set of the stable variables X,s as Pow (Xx())). s), then elements of the quotient spaces Pow(Xs))o called g g-equivalent classes. Since all predictors in the same g-equivalent class have the same power in predicting Y. the searching for optimal predictor in Alg. should be conducted among g g-equivalent classes.
are
In the following, we introduce an algorithm to recover the Pow(X, //)//// space. The algorithm takes the stable graph H Gg input and recursively explore stable variables in the order of their distance to Y. In each step of exploration, it create sub-graphs to represent conditional independence after including/excluddig
some stable subsets. We use the maximal ancestral graph (MAG) construct these sub-graphs, thanks to its ability to preserve conditional independence when included (selection) or excluded (latent) variables exist, In the following, we omit the subscript in Gg and Xs for brevity.
Algorithm Pg Recover(G))
Input: a causal graph G.
Output: the set of all g-equivalent classes Pg.
1: Let the covariate set of G
2: Find Neig (Y)
3: if Neigo(Y
Ø then
4: return (Pow(X))
5: else
6:
Pg + +0 0
7:
for T in Pow (Neigg (Y)) do
8:
S T.L Neigo (Y) O Neigggy(
9:
G' + MAG S)
Po. Recover(G")
11:
for [X.] in Pay do
for X, XX] do
13:
X, X2 UT
14:
end for end for
15:
16: Pg < Pa UPa
17:
end for
18: return Pa
end if
Algorithm G MAG(G,
Input: a causal graph G over X = OULUS. Output: a causal graph G" over O.
I: for each pair of variables A, B E O, A and B are adjacent in G' if and only if there is an inducing path relative to < L, S > between them in G.
2: for each pair of adjacent vertices A, B in G', orient the edge between them as follows:
3: (a) orient as A + in if A E Auca (B U S) and & Anco AUS);
4: (b) orient it as B + A in G' if A Ance (B U S) and B Anca (AU S);
5: (c) orient as A + in G' if A g Ancg(B U S) and B & Ancg (A (UUS):
6: (d) orient it as A
B in G if A Ance (B S) and B Ance U US;
Proposition C.7. Alg T outputs the correct g-equivalent classes in causal gruph G.
*The stable graph the graph after removing all vertex in X
25

---

Under review as a conference paper at ICLR 2023
Proof. We firstly introduce some notions that will be used in the proof. Denote [X,] (X, X, ~ X,} the equivalent class with representative element X, Denote the set of all equivalent classes as Pow (X)/ ~g. Define length of path the number of edges in it. In causal graph G, we say X; is Y's w-order neighbour if the shortest path between Y and X; has length w. As special case, X; is called O-order neighbour of Y if there is no path between Y and X, Let SQ(G) 0 if does not have any neighbour, let (2(G) = 3, if Y has 1,2, 3, -oddr neighbour, respectively.
Note that if we construct a MAG G' over O by G' MAG(G, O, S,L, then, for any vertices sets V,, V 1 & C O, we have V; La V,J V, S < V; Lo Vş[V.. The proof is available at Sect. 2.3 Zhang 2008
In the following, we prove the proposition by induction on I?G)).
Base. For any causal graph G with (GG)) O, we have Neig(Y) So, for any X, X, XX we have Ø S X, n X2 such that Y La X, n X, 102, which means X; ~G Xj. This means Pow / ~G [[]]], i.e. all subsets of the covariate set are equivalent and there is only one
equivalent class.
Induction Hypotheses. Assume any causal graph Ggw with ((G) S w, Pow(X)/ ~G Recover(G
Step. In the following, we show any causal graph Gw41 with S((G) = w + 1, we have Pow P ~G Recover(Gu-1).
Denote covariates in Neiga (Y) as X2 x Denote the power set 2. Q2 = {X]}
Pow([X) X1 X1 my }) as Q = (Qi, Q2.... Q2n1 J. where Q1
Q2e: = [X] X
n1
Any X; C X can be written as X; = X U Xothee, where x contains all 1-order covariates and Xother contains the others. So, Pow(X] can be partitioned into 2 sets (R1, R2,,R,,... R2,... R2 }. where
R. X Q:).
Now consider an element X; E R, and an element X, E Rj. i + j. Because X * x and X X are connected with Y in the causal graph, we have X; AG Xj. This property means Pow (X)/ ~G.., = - R ~G#+1'
By the aforementioned property of MAG Zhang 2008 if we construct G¹ with ((G < w by G² == MMGGGGw41, O X Neigg.......... (Y, Neigg.., (Y))Q.. Q,), then we have R4/ ~G--1 = RA/ ~ So, we further have Pow(X)/ ~Gu41 - R4/ ~G.t1
U R./ N..
By the induction hypotheses, R²/ ~ Recover(G.). So, the above equation can be further written as Pow(X)/ ~G.11 - R ~G Ra/ ~g = U Recover(")). By de- sign of Alg. (line- have U Recover(G)) (eoovrr(Gw+1)) So, we eventually have Pow / ~G.-- Recover(G..11).
Corollary C.8 (g-equivalent Classes in Sub-graphs). Denote the causal graph G, the covariate set
X. Let Z C X a subset of covariates. Denote covariates in Z as {Xf, X2 Xf}, and the power set Pow( Pow([X;:
= (Q1, Q2..... Q }, with Q1 X²). Construct 2' sub-graphs Gi 1, 2, ..., 2 by G. S Q.
Xf...))))
2, Q2 {Xi) Q2
{X;
ZQQ, O = X(Z). Denote the number of G-equivalent classes in the causal graph G and sub-graph G as No.
Then, we have NG S = No:
Proof. Because we do not restrict the set to Neig(Y), an element X, from R2 and an element X, from R, may be graphical equivalent. So, the equal to mark Prop. C.7 because a greater than or equal to mark.
Indeed, the true causal DAG with complete orientation not identifiable. What we can identify is a partially directed acyclic graph (PDAG), representing all Markovian equivalent graphs of the true
26

---

Under review as a conference paper at ICLR 2023
DAG. The following proposition, which states all the Markovian equivalent graphs have the same 8 g-equivalent classes, assures that our recovery algorithm can be applied to PDAG.
Proposition C.9. Under assumption 2.2 causal graphs in the same Markovian equivalent class have the same g-equivalence.
Proof. By definition, causal graphs in the same Markovian equivalent class have the same probability distribution. Under the Markovian and faithfulness assumptions, this means they have the same set of d-separations. As g-equivalence is defined d-separation, they also have the same g- equivalence.
C.2 DETAILS OF CAUSAL DISCOVERY TO DETECT LOCAL COMPONENTS
In this section, we summarize method detect X De(X XM, De(XM), Blanket(Y), (PA(X.)) xX.UU))X..) PC(Y) == PA(Y) U Ch(Y). The Blanket denotes the Markovian Blanket Y.
The identification of X U De( X M) are in Alg. 3 To distinguish X and De(X)r). ). it suffices identify the direction of X; X; in the case when both X; and X; are in Ar. which can be accomplished by comparing Å xx,+x, and Å Ax,-x, (see Huang et al. 2020) for details). However, it should be noted that distinguishing X and De(X) unnecessary. The identification of XMs and PC(Y) M for the estimation of h(S_, and fs is Alg. 2 where PC(Y) can be obtained from the undirected skeleton. Blanket(Y) can be identified by Aliferis et al. (003)). The identification of XM De(XM) in Alg. and we can distinguish XM from De (e(MM) using the way as in (X Mr De(X)))) The parents {PA(X;) XM U De(X.)) can be identified by Alg.
C.3 DETAILS OF ESTIMATING fs
To estimate fs_. - we adopt soft-intervention replace 1 ("((x[PA(XM)) with P(Xar) and hence define p(x y) = p(y pa(y)) Thes p(=:{pm((z)))p(XX/ Then we have fs Ep [Y| XM). To generate data from P', we first permute X,s in sample-wise manner to generate data from P(XM). We then regenerate data for XM descendants in the intervened graph via estimating structural equations as summarized in Alg. g
Algorithm 8 Estimation of fs
INPUT: training data {a(K), Mo S_ S, XM. De(Xm), and ( (Xu)
OUTPUT: Trained fs-.
Shuffling ((xM) (a) | by randomizing the indices.
2: For X; E Dee (X do
3: Regenerate {(:)) ) as (p)(
4: Train Is the regenerated samples.
Indeed, we only need to regenerate De (XM) n Blanket(Y) since y' (l(y)blankett()) = y (v|=)) To maximally reduce the approximation error in regeneration, consider intervene another variable set X == (De((")))))) CC((Y)) and regenerate variables in Deagy (X))). We prove
Deg (X%)) is the minimum regeneration set in the following proposition.
Proposition C.10. Denote X; == x U (De(X))) \ Ch(Y)). Then: For any admissible set Xáo. we have Deax, (Xao) Blanket(Y) Deazz (Xi,):
x 2. do Dea (Xto), and xPA(X)))x,EDeG (x:.) are identifiable.
Proof. (1) Firstly, we prove that set of variables X is admissible means (d)()) ss (x(x))) (XnnCh(Y)) X and (XsnCh(Y)) n X Ø.
Thhis he achieved because XM, ar), and their parents identifiable, shown in Alg.
27

---

Under review as a conference paper at ICLR 2023
Note that
p (ulp)ay)) IxxxxxnoM(() p(z:)pa(zz))
l plulpaly)) ex,,x p(z; p(a))pa(z)))d pa(x;))
plulpa(u))
Pasly/æ)
p(ubballyy) [x p ddy
can seen (y)a)) = plylas, do(X + X (Ch()) XsnCh(YY, which can be rewritten as
1 XsnCh(Y).
The above equation holds and only if (Cu)C)(Y)) C X and (XsnCh(Y)) n = Ø.
(2) Secondly, prove that X is an admissible set and Dea (Xio) De(XM n Ch(Y))
Xs n Ch(Y). To simplify the notations, let X =X and X₃ x
Dexxx (X
The conditions
[XM (XM(Ch()) CX, and (XsnCh(Y))nXx Ø hold by definition.
(2.1) show x Xo)= C X1
Note Xo C
U M nCh(Y))) we have De(X,) De(X Besides,
since Decx,
xo)=De(X,, (X.) Xo, Then
DeGx, (Xo) -D)xoD(XXo) n =De(X)) (B (on
=De(Xo) (X (X
C
n Ch()) = ( x
(2.2) show X, C CDecx (Xu)
Since X n Ch(Y) Xo. De(X M Ch(Y)) c De(Xg), so X, De(M Ch(Y)) De(X,) and hence X C De(X,))) Xn. Besides, note that X Ø such that X Xo and Decx, (Xo) De(Xo) {X,, we have X, C Deax, (Xo).
(3) given X satisfying the two conditions, we have
XMrCh(YY C X De(XMNCh(Y)) c De(Xoo); X C (XsnCh(Y)) (XsnCh(Y)) cXC
Therefore,
De(X M (Ch(Y)) (XsnCh(Y)) C De(Xat) n X Thus, X, C DoG (Xoo) for any X satisfying X (Ch(Y) X and (xsn(X))))())) a.
(4) The identification of (PA(X)) (x EDe Deoxy x%
(x X and Dea (Xaa) can be readily obtained
in Sec. C2
D APPENDIX FOR SEC. 3.3 COMPLEXITY ANALYSIS
In this section, we given some graphical examples and show the number of g g-equivalent classes (denoted by NG) over them.
28

---

Under review as a conference paper at ICLR 2023
Lemma D.1 (Adding/Deleting Edges). In causal graph G. adding edges does not decrease Ng. deleting edges does not increase Ng.
Proof. 1.
For any causal graph Go. add an edge in it and call the resulted graph as G1. We show No. NG. We show this by proving for any subsets X,, Xj. if X; A.o Xj. then X; Na. Xj.
We prove this by contradiction. Suppose there are X,, X; such that X; *ao X, and X; ~a Xj- By X; ~G1 Xj. we have , So, X, n Y LG, X; U,G X2y Kyy. Because adding an edge does not change the covariate sets, we have Xij CG, X; O X, Because X; AG, X,. we have Y LG. X; X2 XX, In other word, there is path p in Go between Y and X; U Xij such that p can not be blocked by Xij-
This means Xij does not contain any non-collider p. and X,y contains every collider (or its descendants) on p in Gg. Because G1, p still a path between Y and X; UX (X.y. Besides, any collider X. on p in Go still collider on p in G1- Any variable Xa Dec (ee(..), where X, is collider on p in Gn. still a descendant of the collider on p in G1- Any non-collider X,, on p in Go still non-collider on p in G1. We have the path p can be blocked by X.j in G1, neither, which contradicts with the claim that Y Lg. X U X, XX. Xij.
For any causal graph Gj, delete an edge in it and call the resulted graph Go. It is straight forward to prove No, 2 NG. using the conclusion from 1.
Y
Chain
(b) Star
(c) Circle
(d) Knots
Figure 10
Example 2 (Chain). A chain graph is graph whose skeleton a chain, Le, Y -x,,, Xnn- Xn-1-...-K For any chain graph with covariates, NG., 1.
Proof. We prove the claim by induction.
Base. When n NG. 2 n holds.
Induction Hypotheses. Suppose for chain graphs with n covariates, the NG. = + 1.
Step. When there is covariates in the chain graph, i.e., G.+1 has skeleton Y X X,. ., As Neigg (Y) = (Xn+1) we need to discuss the number of g-equivalent classes when including and excluding Xa+1- Xn-1 is collider, then, when including X the X1; when excluding Xn+1, skeleton of the X1- By induction hypotheses, NG = n+1+1 2
X, -x-1
induced MAG G' has skeleton Y X, X,-1
induce MAG becomes Y Xn - X holds. We have similar conclusion if X,+1 is non-collider.
Example 3 (Star). A graph with k-branches is graph whose skeleton is composite of i: disjoint chains, For any star graph with k-branches and n covariates, NG., O
Proof. A star graph is composite of & disjoint chain graphs, each containing a covariates. So, we have NG = a + 1) )
Example (Circle). A circle graph a graph whose skeleton is circle, Le,, Y-X,- -0-xxna X X1, Y - X1- For any eircle graph with covariates, NG... 0(n²).
Proof. If Xn, is collider, then, when including Xn, the induced MAG will have skeleton Y X,-1 X,, Y X1, which is eventually a circle with n covariates; when excluding X
skeleton of the induced MAG becomes chain with covariates Y X1 X2 ... X So.
29

---

Under review as a conference paper at ICLR 2023
have NG., n No,,, which means {NG. hu is arithmetic sequence w.r.t According to summation formula of arithmetic sequence, we have NG.. = =(n²).
X
X
X1
X4
Xs
X
X4
X5
Xs
X
X
X3
X7
X3
Xy
X7
The original DAG
The induced MAG when excluding X,
The incuced MAG when excliding X,.Xa
Figure Recovering g- g-equivalence in knot graph.
Example (Knots). A knots graph (shown in Fig. II (a)) a generalized directed chain graphs, where each knot contains 4 covariates. For knot graph with covariates, have No., Oc")), for some constant
< c <
Proof. We prove the claim by showing the recursion formula of NG w.r.t the knot number k.
For the knot graph shown in Fig. 11 (a), the only neighbour of Y is X1, As X, is non-collider, when including X1, will not have any neighbour in the induced MAG, so, the number G-equivalent classes in this sub-graph will be 1. When Xi excluded, the induced MAG is shown in Fig.
(b), where Y is adjacent to three covariates X2, X3, X4. As a result, we need to consider 2 combinations of covariates including/excluding. Out of the combinations, will induce MAGs where Y has no neighbours while the other 4 combinations (excluding X1) will induce the MAG as shown Fig. [ (c), which is eventually knot graph with k knots.
them (including X,)
So, we have the recursion formula of NG w.r.t the knot number & to be NG. = 1 + (4 + 4 NNG..) This formula means NG increases exponentially w.r.. the knot number k. Because a NG also increases exponentially w.r.t. n.
Y
X1
X2
Xm
/
Xm+1
Xm+z
X
X1
X2
X
Xm+2
Xn
X1
X2
Xm
Xm+z
Xn
(a) The original DAG
(b) The induced MAG when X is blocked
(e) The induced MAG when X open
Figure 12: Recovering g e-euuivalence in lollipop.
Example 6 (lollipop). A lollipop graph, shown in Fig. 12 (a), is constructed by adding edges among Y and the first m covariates in a chain graph. Sugar of the lollipop made up of Y, ".., Xms while the stick is Xm+1 Xm X.. For any lollipop graph with covariates, No.. = O(n).
Proof. Note that covariates X1, X2,.., Xm in the sugar part may play different roles (non-collider or collider) on different paths through them. So, it can be troublesome analyze Y's neighbourhood as we the chain graph. Fortunately, covariates in the stick only belong to one path, so, we can them to construct upper bound of NG..
Formally speaking, we can construct an upper bound NG., with corollary C.8 Specifically, when Xm+1 is blocked, number of g g-equivalent classes in the induced MAG will be less than 2 when X is open, the induced MAG G' will be eventually a lollipop with n covariates, as shown in Fig. 12 (c). So, have the following inequation: NG., < 2mh No:
Recursively performing the analysis on G²,-11 we have NG. S 2m 2 + No S S 2m m 1). As is constant number, have No. O(n). n
2m
2m 2m
30

---

Under review as a conference paper at ICLR 2023
Am
Am
An
Am
WM
M
8m
Hm
B
Bm
Bn
2+
(a) The original DAG
(b) The induced MAG when Am-1
open
(c)
The induced MAG when Amo 8mss are both open
Figure 13: Recovering g g-equivalence in 2-1ollipop.
Example 7 (k-lollipop). A lollipop with &s sticks called a k-lollipop. For any k-lollipop, we have No.. = (nn²).
Proof. Let's firstly look at the number of g-equivalent classes in 2-1ollipop, as shown in Fig. 13 (a). Following Example we can construct an upper bound of NG., by performing Corollary C.8 on [Am+1, 1m+1). There are 2 blocked, number of g-equivalent classes in the induced MAG will be less than 2tt (ii,iii) when one of covariate in {Am+1. mm+1 } is blocked, the other is open, the induced MAG will be eventually a /-lollipop, as shown in Fig. 13 (b). so the number of equivalent classes in this situation will be bounded by ( ): (iv) when both Am++ and Bm++1 are open, the induced MAG G' is a 2-1ollipop
situations, specifically, (i) when Am+1 and Bm+1 are both
with n 2 covariates, as shown in Fig. 13 (c). To conclude, we have the following inequation: Ng.. I 2 + ( NG
Recursively performing the analysis on G we have Ng. 2 +0( No 2 + 2
2 +
0("22) +NN s... < 2" + + 2" +(5) + 0 mm= + ((1) 0(2 +0
As m a constant number, we have No., (n²))
Inspired by this observation, we can analyze Na., in k-lollipop by induction. Formally speaking, Base. For 2-lollipop, NG., On²) holds.
Induction Hypotheses. For k-lollipop, No.. 0nn)) holds
construct the upper bound of NG., by performing Corollary C.8 on (X) (m=+1.-.. x m where X X is the left-most covariate the i-th stick. There are 2* situations: (i) when at least one of X +I... Xh+1 is blocked, the induced MAG will be a lollipop with less than equal to k sticks. According the induction hypotheses, number of g-equivalent classes X
Step. For k+ /-lllolppp
x
in these sub-grat aub-graphs will be at most (n*)). (ii) When all covariates of X1 m++1--- 1--1- m+11 open,
the induced G' will be eventually k /loollppp with k covariates. So, we have the following inequation: NG." S Onr*) + No..."
Recursively performing the analysis on G/ have No., I O(n²) + +No 50(nn)+0(( k)²) +... 0(1) = +1)
Xn-1 X Xa
Y
Xn-z
Xix
X1s+1
Y
Xis-z
Xn
Xn
(a) A skip chain with m
(b) Induced MAG when Xa-1 is blocked
(c) Induced MAG when X is open
Figure 14: Recovering g-equivalence in skip-chain with m
Example (Skip Chain). A skip chain graph is constructed by adding skip connections among Y and m covariates a chain. For any m-skip chain with n covariates, NG = O(n 2m).
Proof. Let's firstly look at the skip chain graph with m = 1, an example of which is shown in Fig. 14 (a). As we can this example is constructed by adding skip connections between Y and X
Following Example we can construct an upper bound of NG.,, by performing Corollary C.8 on
X When X blocked, the induced MAG is branches graph, as show in Fig. 14
31

---

Under review as a conference paper at ICLR 2023
X
X
Xt1-2
Xp-1
Xn
X1141
Xiss1
y
X12-1
X12-1
Y
X
X12
X22+1
X
Xi2-1
X
(a) A skip chuin with
(b) A constructed upper bound for (a)
(c) The induced MAG when X blocked
Figure 15: Recovering g-equivalence in skip chain with m = 2.
(b). So, the number of g-equivalent classes in this sub-graph is O(n²). When X is open, the
induced MAG G' will be eventually a skip chain with m 1 and covariates, as shown in Fig. 14
(c). So, we have the following inequation: No.,, I O ) + No. Recursively performing this analysis, we have No... S <O(n²) + No.,,.. S I O( nn O(( + 0
Now that we have analyzed the skip chain graph when m 1, let's look the situation when m an example of which is shown in Fig. 15 (a). Firstly, we can construct upper bound for any skip chain graph with m = 2 by adding extra connections among Y, X1,, X12, until these three vertices form complete connection,
shown in Fig. 15 (b). Then, again perform Corollary C.8 on X11-1. When X, is blocked, the induced MAG contains two disjoint branches, one which is a chain, the other is eventually skip chain graph with classes this sub-graph (nn)). When X is open, the induced MAG will be a skip chain with 2 and covariates. So, similarly by the inequation when m 1, we have No... O(
1. So, we have number of g-equivalent
Following the spirit of the above analysis, we show NGm.. (m²) by induction.
Base. When No... (n²)) holds.
Induction Hypotheses. Suppose for any skip chain graph with m, we have No mm.)).
Step. For skip chain graph with skip connections among Y and XaI,.., X Xim.y firstly construct an upper bound by adding extra connections among Y and Xı1, ..., X Xim+1 until these vertices form complete connection. Then, in the resulted graph, perform Corollary C.8 on X When X blocked, the induced MAG contains two disjoint branches, one of which eventually skip chain graph with less than or equal to m. So, we have number of g-equivalent classes in this sub-graph Onn2++1) When Xa11 open, the induced MAG will be skip chain with m and n covariates. So, similarly by the inequation when m = 1, we have NGm,, = (n 2m)
chain, the other is
Lemma D.2 (Property of Tree). A tree is an undirected graph in which any two vertices are connected by exactly one path. If there are dy leaves and dz3 vertices of degree at least three the tree, then dy 2 d23 + 2.
Proof. Denote number of all vertices in the tree as dr, then by the handshaking lemma,
dr
d, + 2(dT -d, dza) + 3d23 <
deg(V,) 2(ddr - 1),
which indicates dr 2 d23 + +
Proposition D.3. Complexity of Alg. i O(Nc).
Proof. Treat each call of the MAG(.) function in Alg. 7 as unit operation.
1. In the recursion tree of Alg. 6 number all vertices dy the complexity of Alg. while number of leaves d, is Ng-
32

---

Under review as a conference paper at ICLR 2023
2. Each interval vertices in the recursion tree has degree at least three (one parent vertex in the tree and at least two children vertices in the tree). By Lemma D.2 dr+1
S dr. S dr, which indicates
dr = (dd)) and thus complexity of Alg. is ((Na).
E APPENDIX FOR SEC. 3.4 SPARSE MIN-MAX OPTIMIZATION
In this section, introduce theoretical analysis of the following empirical min-max optimization problem and besides, a more efficient algorithm called Linearized Bregman Iteration (LBI).
Σ 1=1
(1= s, sr.,rr) + ||
min max 0,8 0
(11)
2n
where T1, En midd pl y 2m Jo(pa((xs)) For simplicity, use denote p(æ, y/am Js(pa(z)))) in the rest of this paper.
E.I STATISTICAL CONSISTENCY
Denote a, a, 0: x. :::y)==> (y Sa (ar58, ax)) Suppose we can obtain the worst-case risk with 0" == arg maxe Eple(S, a, 0: 0; y)]. Then we denote
1 (,3,a,0:00) vi), (a", arg min EpC(
a(3,a)
2n
0,
Based on this, we denote A == supp( In this regard, the optimization with respect to (a, 3) is Lasso with general loss. Since goal is to recover the optimal subset A and the predictor with (a", BB*), we are interested the model selection consistency and fz-consistency properties:
Model Selection Consistency: lim P(A.) = P(A), where A, =- supp(3n).
fg-Consistency: lim, Ik5. 11 where Ç = {a TT)
Here, we denote Č. == arg min, C(C) A |||| The model selection consistency can ensure us to find the optimal subset and the fg-consistency further guarantees the optimality of learned predictor. In the following, we discuss two settings: i) fixed when |S| = d is fixed: ii) high-dimensional increases with n. We first introduce some assumptions, which are commonly made in Lasso Zhao Yu (2066;; Negahban et al. (2012); Rejchel
Assumption E. (Restricted Strongly Convexity (RSC)), We assume that convex; L and Q := [, are twice differentiable and satisfies H := V2L(a", 38*) 2 * and H = (aQ
y * for some 7 > 0.
Assumption E.2 Sgguare-ntegrrability of the gradient). We assume E [il(a, oo for each (ar, () in some neighborhood of (a*, (a*, Assumption E.3 (Irrepresentable condition). We assume that
(sign((r)
0 sign((8)) )l <1
HAs H1 (aaA))
Remark E.4. The restricted strongly convexity condition has been widely assumed in variable selection Negahban et al. (2012); Zhao & Yu (2006), especially in high-dimensional statistics to ensure the identifiability of the oracle parameter. The irrepresentable condition was almost necessary to recover the true signal set. For regularity, it was needed in general convex loss Niemiro (992); Rejchel (2016) to ensure asymptotic normality.
Now ready to introduce results. We first introduce the model selection consistency with fixed setting. Before that, we first introduce two lemmas in Rejchel (2016).
Lemma E.5 (Corollary 2.3 in Rejchel 2016)). Under assumptions
E. E. and set Xn such that
lim, Au 0 and lim,
Aa. we have X, (Cn, (*) +p 5D arg mine V(C) with
Vn
V(C) 1 HC TTC EGsigg(())) Σ Gl-
jEA
JEA
33

---

Under review as a conference paper at ICLR 2023
Lemma E.6 (Theorem 2.3 in Rejchel (0166).
Under the same conditions in Lemma. E.3
we have
aC(0)
ac(s*)
Σ
(1(555)) +p O.
o
aç
ac
With this lemma, have the following model selection consistency results:
Theorem E.7. Under the same conditions in Lemma. E.5 and additionally assumption E.3 we have that lim, P(An)) P(A).
The proof is very similar to Corollary 2.4 in Rejchel 2016). We include it here for completeness.
Proof. Denote Lx,, (C) == ((() + /n/n. Note that if there exists E AA then we have P(j E A) P(zétan(j) = 0) + O according Lemma, E.5 Thus we have P(A C A.,) + 1. Next, we show that P(An C AA + Otherwise, Vn > 0, there exists j E An, but not belong to A.. Recall that So. minimizes Lx, (S), have that
aC(Sn) 08,
An 013,1 = 0.
Since 3, # 0. we have that
ac(Sn)
n
= 1.
An
083
Besides, we have that
n aC(Gn) n
(OC(GN) a3,
aC((") 035
(Š. -8')
n OC(S)) + An a3j
n H(S. 5*). An
+
An 03,
An
According to Lemma. E.6 we have that the first term converges to O in probability; besides, due to square-integrability and central limit theorem, we have that the second term also converges to in probability. From Lemma. B.5 the third term converges to HÇ in probability. Note that ( satisfies that:
(a T. (TTTT = =00 (0T, (sign(33 (2)))
H (ass(an)) TT TTT 2 |l
0
Therefore, we have [1 since the irrepresentable condition holds. In this regard, we have that
n aC(G.) An 05,
<1,
which contradicts to the fact that j E An.
Next we show that in both fixed and sigh-dimensionee settings, we have the following eg-consistency, which is natural conclusion applying the results in M-estimator Negahban et al. (2012):
Theorem E.8. Under assumptions and suppose An 2 2n||V [l(*))) Then we have
(
X2
Ir -
(IA) dim(a:)
n2n2
Proof. According to theorem Negahban et al. 2012)) we have that
X 2
ll - 3 0
n²-,²
if A, 2 2nR (VL(())) for some regularization funtion with R denoting the conjugate function of R and V/ = Enyo | The proof is completed by noting that R² 20.
R(y) In our setting, R(3) = || Therefore, we have v < VTA + dim(a).
34

---

Under review as a conference paper at ICLR 2023
Remark E.9. According to square-integratiliy and the large law number theorem, have that VC(() ta.s. O and thus [VC(()) -ta.a. O. Therefore, as long as A satisfies conditions in Thm. E.7 An satisfy Xn VnI[V )l In this regard, both model selection consistency and fz-consistency in fixed setting can hold; while in high-dimensional setting, we have f2-consistency, which is our ultimate goal, i..,, identifying the minimax optimal predictor.
E.2 LINEARIZED BREGMAN ITERATION
In this section, we introduce an alternative algorithm, namely Linearized Bregman Iteration (L.B) to replace the minimization step via Lasso. I.II was firstly proposed in Osher et al. (2005) in image denoising. In Osher et al. (2 Huang & Yao (018)). the authors established LBI's statistical model selection consistency from the perspective of differential inclusion. Such consistency holds under nearly the same condition in linear model; however additionally requires restricted strongly convexity to hold for each solution in the path, under general convex loss.
More importantly, LBI enjoys more efficiency in implementationn, compared to Lasso. Specifically, the condition on A, for model selection consistency and & fg-consistency is in an asymptotic form. sequence of hyperparameters and run
In practice, to select the optimal Lasso has optimization algorithm for each hyperparameter. In contrast, I.II can generate a whole regularization solution path, with each iteration corresponding solution in the pat. Motivated by this property, we proposed to replace the minimization step via LBI, which is composed of a gradient descent followed by settthhresholing step.
Combined with the gradient ascent step, the algorithm is showed as follows:
Maximization step:
B 0% + dVVl(Bk, ag. 0x), gradient ascent w.r.t. 0
Linearized Bregman Iteration:
a+++ O.k aovve((8 az 0x), gradient descent w.r.t.
2k+1 = 21 -V (Bk, ak Oz), gradient descent w.r.t. |
B&t1 (isign(2++1) max (0, |2k41| 1). soft- soft-thresholding to obtain B
Here, the step size, = == ||3|| + 2r | and k > 0 denotes the damping factor which is trade-off between efficiency and statistical properties. Specifically, Inverse Scale Space (ISS) which can return unbiased solutions, is the limit of LBI as и - oc; however, large x leads to computation inefficiency by noticing that 5 and к should satisfy 5K < 1/Amax (V2)) where Amax (A) denotes the maximal eigenvalue of A. Instead of running an optimization algorithm in the minimization step via Lasso, it only spends gradient descent and soft-thresholding steps, which is much more efficient for implementation. However, a vinadvantage of using LBI to replace Lasso lies in the lack of statistical consistency guarantees, the LBI alternates with the gradient ascent w.r.t. O.
35

---

Under review as a conference paper at ICLR 2023
F APPENDIX FOR SEC. 4 EXPERIMENT
Implementation of Baselines. Vanilla uses [x] to predict Y and is implemented by the same neu- ral network as fs_ (which will be introduced later). Other baselines are implemented by the authors' official codes. Specifically, ICP
IC (https)
Anchor regression
IRM (https:/7
HRM (https 7/github. com wwwmeosttu/ IB-IRM (https: wwttps::/giituubbb..... com/ jack///aaaa////aoaaaaaaaaaaaa/aa IB- IRM;; As
the Surgery Estimator did not provide official codes, we implement it following settings of our method.
F.1 SIMULATION
Implementation Details. In all three settings, SGD is used for optimization. In Setting-1,2 the structural equation æ1 f- 91 (41, y) + -1 estimated by a one-layer fully-connected neural network (FC), with training iterations set to 1000, the learning rate set to 0.01. In Setring-3, the equation is estimated by a two-layers FC with a sigmoid activation function in the hidden layer, with training iteration set to 1000, the learning rate to 0.01. In Setting-I: fs_ is parameterized by one-layer FC, with training iterations set to 2000, the learning rate set to 0.001. Jo is parameterized by the same structure, with training iterations parameterized by one-layer FC, with training iterations set to 1000, the learning rate set to 0.001. Jg parameterized by the same structure, with training iterations set to 5000, the learning rate to 0.05. In Se:ting-3: 3: fs_ is parameterized by a two-layers FC with sigmoid activation function in the hidden layer, with training iterations set to 5000, the learning rate set to 0.01. Jo is parameterized by the same structure, with training iterations set to 2000, the learning rate set to 0.01. The codes are implemented with Py Torch 1.10 and run on server with an Intel Xeon E5-2699A @ 2.40GZZ CPU.
2000, the learning rate set 0.05. In Setting-2: f's is
Additional Results on Causal Discovery. We randomly generate DAGs according to the Erdos- Renyi model Erdös al. (1960). We consider three low dimensional settings of nodes number {6, 8, 10} and a high dimensional setting with 100 nodes. For the low dimensional settings, we generate 10 domains, where the number of mutable variables is set to {2, 3). and the sample size n. is set to 200 for each domain. For the high dimensional setting, the generated graphs are sparse. We generate 20 domains, where the number of mutable variable is set 20, and the sample size n. is set to 500 for each domain. For the low dimensional settings, we implement the PC et al. (2000) algorithm learn the undirected skeletons. For the high dimensional setting, PC-stable Colombo et al. (0144) is used. Our algorithm is then used to determine local components. To remove the effect of randomness, we repeat for 40 times. We report the F1 score, precision, and recall in Tab, 2 As we can see, when the causal graphs are more complicated, our discovery algorithm can still give accurate results, which further validates its effectiveness and stability.
Table 2: Performance of Causal Discovery.
Nodes
8
10
100
Metries
F1
0.98
0.96
0.94
0.88
precision
0.97
0.96
0.94
0.86
recall
0.99
0.96
0.95
0.90
Comparison with Baselines. We report the maximal mean square error (max MSE) over the test sets for our method and baselines in Tab. 3 Besides, we in Tab. 4 report the standard deviation of mean square error (std. of MSE) over the test sets as a measure of transferring stability. As we can see, the maximum and standard deviation of MSE of our method are both low. For example, max MSE is 0.0075, and std. of MSE is 0.0006 in 2. This verifies that our method is both robust and stably transferable to distributional shifts. Besides, method has large improvement over baselines in
36

---

Under review as a conference paper at ICLR 2023
the highly non- non-linear setting-3, As for the slight improvements over the baseline, it may be due the simulation settings being simple enough for the vanilla method to only exploit X2 for prediction.
F.2 ALZHEIMER'S ' DISEASE DIAGNOSIS
Implementation Details. The imaging data are acquired from structural Magnetic Resonance Imaging (sMRI) scan, After prte-prepooeessing via Dartel VBM (Ashburner 2007) and Statistical Parametric Mapping (SPM) for segmentation, we partition the whole brain into 9 brain regions according to Tab. and Tab. 7 Data normalization (w.r.t. mean and standard deviation) is used. All structural equations are estimated by a two-layers FC with a sigmoid activation function in the hidden layer. For structural equations generating X3, the training takes 5000 iterations, with the learning rate set to 0.1. For those generating X4, Xs, X6. Xy, the training takes 2000 iterations, with the learning rate to 0.1. fs_ is parameterized by two-layers FC with sigmoid activation function in the hidden layer. with training iterations set to 5000, the learning rate set to 0.25 (decrease to 0.1 at iteration 4000). Jg is parameterized by the same structure, with training iterations set to 2000, the learning rate set to 0.25. SGD is used for optimization. For the sparsity-based optimization, we set the training iterations to 350, with learning rate set 0.05, and penalty weight to 2. Adam is used for optimization.
We pick four domains with more than 40 patients as the training domains and test on the rest three domains. To remove the effect of randomness, we replicate over all the 15 possible train-test splits.
Additional Results. We firstly report std. of MSE over the test sets for our method and baselines in Tab. 5 As we can see, our method outperforms other baselines by significant margin. This result demonstrates the utility of our method in learning stably transferable predictors. Then, we compare the performance of all g-equivalent classes with more than one member in Fig. [8] As we can see, most equivalent classes have similar performance (small deviations). As for the several classes with large deviations, it may be due to the approximation error incurred during inferring the causal graph. Next, we show the optimization curve of JJ) and max MSE for 100 randomly picked subsets S_ C S, in Fig. 16 and Fig. [7] As we can see, the optimization over Jo is well converged, and the performance of different subsets is consistent with our expectations. This observation again suggests the utility of Thm. 3.5 in finding the optimal predictor. Finally, we show the loss curve of sparsity-based optimization in Fig. 20 As we can see, the optimization over h is well converged.
F.3 GENE FUNCTION PREDICTION
Implementation Details. Data normalization (w..t. mean and standard deviation) used. Structural equations generating X2, X4 are estimated by two-layers FC with sigmoid activation function in the hidden layer, with training iterations set 5000, the learning rate set to 0.01. fs is parameterized by two-layers FC with sigmoid activation function in the hidden layer, with training iterations set to 20000, the Icarning rate set to 0.01. 19 is paramcterized by the same structure, with training iterations set to 120000, the learning rate set to 0.05. SGD is used for optimization. For the sparsity- based optimization, we set the training iterations to 350, with learning rate set to 0.05, and penalty weight set 0.1. Adam is used for optimization.
We use the wide-type mice and three kinds of gene knockouts in the training domains. To remove the effect of randomness, we generate 45 replications, with each trial appending out of the remaining 10 gene knockouts to the training domains and testing on the rest 8 gene knockouts.
Firstly, we report std. of MSE over the test sets for our method and baselines in Tab. 6 Similarly, our method outperforms other baselines by a significant margin, which together with the Alzheimer's 's disease experiment, shows the utility of our method in learning stably transferable predictors. Then, we show the optimization curve of h(S Jo) and the loss curve of sparsity- based optimization Fig. 19 and Fig. 2I respectively. As we can see, the optimizations both well converge.
Additional Results.
37

---

Under review as a conference paper at ICLR 2023
Table 3: Maximal MSE comparison on simulation data.
Vanilla
ICP
IC
IRM
AncReg
HRM
IB-IRM
Sung
Ours
setting-1 setting-2
1.90
2.17
1.68
1.38
1.34
2.69
1.58
1.184.06
1.18 00
007.00
.17..51
064.02
064.00
0071
33
.2
2005
.073
setting-3
1.72
1.61
1.544.62
2.98
2.34
1.75
1.71
1.10
1.10
Table 4: Mean (over randomization) of std. (over test domains) of MSE on simulation data.
Vanilla
ICP
IC
IRM
AncReg
HRM
IB-IRM
Surg
Ours
setting-1
.36
.22
.27
.01
.16
.56
.14
.10
.10
setting-2
.0057
.0111
.0051
.0034
.0005
.0820
.0142
.0006
.0006
setting-3 .35
.14
.27
.59
.40
.24
.20
.09
.09
Table 5: Mean (over randomization) of std. (over test domains) of MSE on ADNI dataset.
Vanilla
ICP
IC
IRM
AncReg
HRM
IB-IRM
Ours
0.267
0.270
0.252
0.166
0.161
0.294
0.244
0.018
Table 6: Mean (over randomization) std (over test domains) of MSE on IMPC gene dataset.
Vanilla
ICP
IRM
AneReg
HRM
IB-IRM
Ours
0.257
0.274
0.302
0.319
0.275
0.278
0.259
0.017
Table 7: Brain regions partition.
Brain Region
AAL Index
Frontal lobe (X1)
34..67789110.112211114.5555
Medial temporal lobe (X2)
85,86,87,88
Parietal lobe (Xx)
59,60,61,62
Occipital lobe (Xa)
1,,50.5,,,2,53.54
Cingulum (Xs)
31.32.33.34.35366
Insula (Xe)
29,30
Amygdala (X7)
41,42
Hippocampus (Xg)
37,38
Pallidum (Xg)
75,76
38

---

Under review as a conference paper at ICLR 2023
Table 8: Automatic Anatomical Labeling (AAL) indices for brain regions.
Brain Region
AAL Index
Brain Region
AAL Index
Precentral L
Precentral
2
Frontal_Sup_L.
Frontal_Supp SanntllSuupRRee
4
Frontal_ SitiallSup_OHae 1uup_rrbLL
Soonaa__Sup__a Orb_R
6
Frontal_Min
Frontal __d__Mid____ R
8
Frontal __d_Mid__00 Orb.
Frontal Mid Orb
Frontal Inf Oper L
11
Frontal Inf Oper R
12
Frontal_ Inf Tri Tri__.
13
Frontal Loontal_Inf__e Tii_Tri__
14
Frontal Inf St_rrb_l
15
Frontal_ Inf_Orb________
16
Rolandic Oper L
17
Rolandic Oper
18
Supp_Motor_Arca
Supp_ Hopp_MotorAAree
Olfactory_L.
21
Olfactory_R
Frontal Mud_Me
Medial L
23
Frontal Sup_ hup_Medial_ R
24
Frontal Montal_Mid___ Mid_Orb. I.
25
Frontal
26
Rectus_
27
Rectus R
28
Insula_ I.
29
Insula R
30
Cingulum Ant
31
Cingulum Ant
Cingulum Mid L
33
Cingulum Mid R
34
Cingulum Post L
35
Cingulum Post R
36
Hippocampus I.
37
Hippocampus R
38
Purallippoaampal L
39
Parallippocampal
40
Amygdala
41
Amygdala R
Calcarine_L
43
Calcarine_R
Cuncus_ L
45
Cuncus_R __
46
Lingual L
47
Lingual R
48
Occipital Sup
49
Occipital Sup_R
50
Occipital Mid
51
Occipital Mid
Occipital Inf
53
Occipital Inf_R R
Fusiform I.
55
Fusiform
56
Postcentral L
57
Postcentral R
58
Parictal_Stte Leteta__Sup____
59
Parietal_Suu S_tr__Su__R__
60
Parictal_Inf
61
Sorictal_Inf_RR
62
SupraMarginal
63
SupraMarginal R
Angular
65
Angular R
66
Precuncus
67
Precuneus_ R
68
Paracentral Lobule
69
Paracentral
Lobuul Lobule R
70
Caudate L
71
Caudate R
72
Putamen_1 L
73
Putamen_F R
Pallidum L
75
Pallidum R
76
Thalamus_1.
77
Thalamus
78
Heschl L
79
Heschl R
80
Temporal Sup
81
Temporal Sup R
82
Temporal_Pes Pole_Poole____e Sup____p____
83
Temporal_Pos Laclll_Po__SSSeeee SSp___Su__R__
84
Temporal Mid
85
Temporal Mid_ R
Temporal Pole Mid L
87
Temporal Pole Mid R
88
Temporal_ Inf
89
Temporal_Ing Inf_R
90
Cerebelum Crus
91
Cerebelum 1arrollum_Cruss1111 1_uus_RR
92
Cerehelum Crus2
93
Cerebelum 2tcbellum_Crus2_11
94
Cerehelum 3
95
Cerebelum 3
Cerebelum 4 L
97
Cerebelum
98
Ccrehclum_61 ___
99
Cerebelum 6_R
100
Cerehelum 7b
_b_
101
Cerebelum _a______RRR_____
102
Cerebelum 8
103
Cerebelum 8
104
Cerebelum
105
Cerebelum
Cerebelum 1
107
Cerebelum 10_R
108
Vermis 1
109
Vermis 3
110
Vermis 4
111
Vermis 6
112
Vermis 7
113
Vermis 8
114
Vermis
115
Vermis
39

---

Under review as a conference paper at ICLR 2023
12.00
X3 K:6
20
M.25 XE,10
3
15
22.5
15
22. 83,66A13
22.0 3425.8
24,85 X31
10
X3 20.
M4XT
5-
25
356
20
280
355
17.5
25
10 D.36X16 X13
32 RLA3AA99211 X21
15.0
X2 X2
12.5
10.0
10
7.5
29.96
72.45
XAX1L
XIXUX
Interations
(c)
3030X 212.011
20
X2 320A31
X5
10
K2.54
X0 21,66433
57.813
X7
26
120
100
15c
0-
200
16c
199
25
20 20
17.5
20
X0LXX
6 12.5
0 15
MAXT
7.5
332
10:
22 30.300403 206.11
X3 36.15,11.1188
180
10c
(h)
Figure 16: Optimization curves of h(S_, on ADNI dataset.
40

---

Under review as a conference paper at ICLR 2023
25
22.3 32.40
X5,07.100
12.0 2125
22,6
30351A12
22.00 SCOALL
15
10.6 XLXTAI
25,013X11
5-
22,00 22.00.07110..2
Interations
Interations
(a)
(b)
10.35 31.07
17.5
A5,01.000
5
81
12.5
26.0 X31
10.0
22,04 XX.XXX1
7.5
X3.84 X7.XX
5
2.5
Interations
Interations
(c)
(d)
XL
17.5
17.5
22.22.3.85
0
12.5
X5.33XX11
35,36,37300 200
10.0
12.5
7.5
7.5
5.0
X3.31
Bil
04
(c)
X01
XLL
25
15
15)
1:10 10
RL.X55XX X31
X8.KD
K1,85X11
the
0
120
63
176
Figure 17: Optimization curves of h(S_,, on ADNI dataset.
Max. Test MSE
5
10
15
20
25
30
35
Equivalent Classes
Figure 18: Performance of g-equivalent classes on ADNI.
41

---

Under review as a conference paper at ICLR 2023
4.0
empty
X.1
3.5
... 3
X_4
3.0
X,1.2
2
X_1XX4
2.5
X_2
$
2.0
X 2X
1.5
... 2X
X_1,X,2.044
1.0
X_2.X_3,__4
2000
4000
6000
8000
10000
Interations
Figure 19: Optimization curves of h(S-, Ja) on IMPC gene dataset.
20.0
h
17.5
15.0
12.5
Losses
10.0
Myd MM 1
7.5
5.0
2.5
50
100
150
200
250
300
350
Interations
Figure 20: Loss curve sparsity- -based optimization on ADNI dataset
4.0
h *
3.5
3.0
Losses
2.5
2.0
1.5
1.0
50
100
150 Interations
200
250
300
350
Figure 21: Loss curve of sparsity- based optimization on IMPC gene dataset
42


Reviews:
[{"review_id": "Xas6FrXJCxu", "reviewer": "Reviewer_ra9S", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "For a supervised learning task where data comes from multiple environments, invariant prediction relies on identifying a stable set of features that don't influence the shift in probability distributions over different training environments. This paper proposes a method to find the optimal sufficient set of stable features. The method is split into two cases depending on whether the optimal set of stable features is the entire set of stable features. The paper proposes a sufficient graphical condition for the latter that can be tested by causal discovery. If the condition fails, the paper parametrizes the dataset shift and does a brute-force search to get the min-max optimal subset. To reduce the search space, an equivalence class of the set of subsets is identified. While acknowledging that this brute force search might be inefficient, the authors frame an alternative sparse min-max optimization to solve the subset selection problem. The brute force search method is validated on synthetic data and both methods are validated on real world data.  ", "summary_of_the_review": "Overall, the authors propose two methods for an important problem - one of which potentially has an exponential search complexity and it's not clear to me how they parametrize the mutable shifts. The other method while having decent empirical performance is not addressed much in the theory portion of the paper. If the authors can address either of the two concerns, it would make the paper much stronger. ", "strength_and_weaknesses": "The paper’s strongest contribution is the empirical validation on the real world datasets considered albeit still being relatively low-dimensional. Although I didn’t check the proofs in detail, both Both the methods proposed seem theoretically sound and believable. The concept of equivalences is also novel and interesting in its own right. \n\nApart from the possibly exponential search complexity that the authors acknowledge, I have concerns about the complexity of searching over the parametrization of the shifts on mutable variables, I.e. estimation of h*(S_). How is this parametrized and why is the cost considered constant in the complexity analysis? If I am understanding the parametrization correctly, this should be exponential in the size of the mutable variables? \n\nGiven that the sparse optimization method is competitive relative to the min-max identification algorithm, it deserves more explanation than just a paragraph. \n\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "U8M_1-UwMsp", "reviewer": "Reviewer_s36p", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper concerns how to select/learn a predictive model from causal relations between variables in the distributional robustness setting. Specifically, the paper argues that given a DAG, a predictive model should be learned from a subset of the variables which are stable across a group of learning environments, while not using the information of those mechanisms which can vary between different environments. The idea is similar to that of Subbaswamy et al. (2019) except that the latter proposed to directly search for the optimal subset from the set of stable variables by comparing their validation errors. This paper provides a more complete analysis and argues that instead of using the validation error one can directly compare and minimize the worst-case risk over a class of distributions with varying mechanisms of how those variables in the mutable set are generated. In addition, in order to save the computational cost, the paper advocates to search equivalent classes of subsets instead of searching all subsets. In the case where the computation is still heavy, the paper suggests to adopt the lasso regression for subset selection. ", "summary_of_the_review": "This paper provides some analysis of the methods in Subbaswamy et al. (2019). It then proposed to compare the max (worst-case) risk directly, defining the equivalent class to save computation, and introduce lasso for consistent variable selection. Each of these could potentially be important but due to the amount of work that is presented in the paper, there are still questions on each of these. Some notations can be clarified at various places.", "strength_and_weaknesses": "+1: the paper is technically sound\n+2: the paper has both numerical implementations and theoretical justifications\n+3: the paper studies an important yet under-developed area, and is very timely.\n\n-1: the paper presents too many concepts but without digging up each one with sufficient depth. Theorem 3.1 is equivalent to the degeneration condition in Subbaswamy et al. (2019). Theorem 3.3 seems to be new. The equivalent class is interesting, but I wonder how much saving it will really lead to (considering that there are costs associated with identifying these equivalent classes, and there may also be errors in this identification.) The lasso part is somewhat disentangled from the rest of the topics and lasso is also not new. If the paper were to present fewer concepts but with enough discussion on each one to provide a better justification, it may work out better.\n-2: The main point of Theorem 3.3 is that the main source of variability of the distributions in $\\mathcal{P}$ is from the variety in $J$, which defines a mutable variable in terms of its parents using a definite function. I wonder how in practice the richness of $\\{P_J\\}$ aligns with the richness of $\\mathcal{P}$. While $\\mathcal{P}$ should contain the data distribution from all environments, how does one consider the maximization of all the possible $J$ functions? Note that all measurable functions $J$ are uncountable. If in practice one only uses a small parametric class of $J$ functions, then they may not be able to cover the base for all $\\mathcal{P}$. I do not really find how $\\{P_J\\}$ is defined at the implementation level from the paper. \n-3: As mentioned above, while the idea of the equivalent class could potentially save some training time, it also takes time to identify these equivalent classes. If the time for training each model is smaller than the time spent on the identification of these equivalent classes, then it may not be necessary to bother with these equivalent classes. \n-4: The section on lasso looks a bit disconnected from the rest of the paper. Here is why. The paper went through all the efforts with Theorem 3.3, dismissing the use of the validation error in Subbaswamy et al. (2019), and defining the equivalent classes. Then all of sudden, the paper says basically if this is too much then we can just use lasso, which is a regularized version of the empirical loss function. If lasso works then what's the point of all the previous discussions?\n-5: There is ambiguity as to what $f_\\alpha(x_S \\beta, x_M)$ is in the lasso formulation (3). What are the parameters $\\alpha$ and $\\beta$? Should $x_S \\beta$ be thought of as the inner product? Is $\\alpha$ then a parameter that governs how the single index $x_S \\beta$ and the mutable variable $x_M$ are used to form the prediction function $f$?\n-6: The variable selection consistency of lasso was only proved for the linear model by Zhao and Yu 2016. I am not sure if some results exist for an arbitrary non-linear model. I may be wrong about this.\n-7: I think the idea is that $x_M$ should not be used as part of the predictive model. That is what I saw in Algorithm 1. However, why x_M is included in the lasso model?\n-8: Some notations are confusing to me. I think D ⊥G Y | Z is defined as D and Y are d-separated by Z in the paper. But what is Y ⊥ G_{\\overline{X^0_M} K| X^0_M (this appears in Example 1). There is also Y ⊥ G_{\\overline{X} \\underline{Z}} Z | X (see top of page 16). What do these notations mean? Still d-separation?\n-9: please clarify what the two columns h*(S-) and max MSE mean respectively. As far as I can tell, h*(S-) is also a maximal MSE (over ...?) as it is defined. So it is unclear to me what is the difference. Is max MSE the test data maximal MSE over different observed environments?", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "GBfnqh0KnB", "reviewer": "Reviewer_AmVi", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper considers the problem of learning stable information to transfer to unseen environments. In particular, following the work of Subbaswamy et al. (2019), it investigates which subset of the whole stable information should the model transfer, in order to achieve optimal generalization ability. Authors propose to maximize over mutable mechanisms and to search over only equivalent classes in terms of worst-case risk, with certain theoretical characterization based on causal graph. The performance of the proposed method is also empirically validated.", "summary_of_the_review": "I think that the paper requires a major revision to make it readable and rigorous (e.g., clearly define the setting and notions of stable, unstable things, the do-operator; make the proof precise and give explanations where it is not very very obvious; give summary or high-level idea of your proof; etc.). Thus, I cannot recommend acceptance for the present version. I would like to re-evaluate the paper, and raise my score if authors revise their paper and could address my concerns/questions during the discussion period. I look forward to the upcoming discussions with authors.", "strength_and_weaknesses": "Pros:\n\n- the topic of investigating which subset of the whole stable information should the model transfer is very interesting \n- interesting theoretical charactizations, which lead to practically efficient algorithms\n- good empirical performance in the considered experiments\n\nCons and limitations:\n\n- the paper relies on causal discovery from observational data, which in practice cannot avoid estimating errors, particularly for large scale settings.\n- the proposed method seems not be able to work for image-like data where one has to learn representations.\n- writing is poor (see below),\n- minor one: causal sufficiency is required ", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "C2yo3tg3c6", "reviewer": "Reviewer_K5Kv", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper proposes a model selection procedure for choosing stable models that have minimax optimal risk across a family of distributions  induced by intervening on mutable variables. This work extends the shift-stable prediction principles established in Subbaswamy et al 2019, who established that one can construct a shift-stable predictor by estimating the intervened conditional expectation of Y given an intervention on mutable variables, and some subset of stable variables. This work addresses some practical gaps left in that paper: in particular, it is often ambiguous which set of stable variables one should condition on (in addition to intervening on mutable variables) from the perspective of minimax MSE, and cases where there exists no ambiguity are hard to identify. To address these gaps, the authors propose (1) a simple-to-confirm graphical criterion that is equivalent to the case where one should condition on all the stable variables to obtain a minimax predictor; (2) a selection criterion based on estimated minimax risk for selecting a minimax optimal subset of variables; and (3) algorithms for reducing the complexity of this subset selection. In addition, the authors propose using causal structure discovery from multiple environments (Huang et al 2020) to learn causal structure and mutable/immutable variables. The authors demonstrate that their algorithm reduces worst-case MSE in several cross-environment generalization experiments.", "summary_of_the_review": "The overall method here seems promising theoretically and empirically, and many of the claims are interesting. However, there are too many additional extraneous claims that are unsubstantiated or distracting.", "strength_and_weaknesses": "**Strengths:**\n + The authors identify several practical gaps in the proposed implementation of the Subbaswamy et al 2019 surgery estimator. The graphical criterion for conditioning on all stable variables is particularly nice, as it saves end-uses from having to run the ID algorithm for each DAG to understand whether the intervened distribution can be written as a conditional distribution.\n + The observation that the best subset of variables based on validation risk in the training environments does not correspond to the best subset for minimax risk is important. I wish the disconnect between these two were highlighted more prominently in a very simple toy example.\n + I like the claim that the minimax risk can be estimated by exploring the set of potential downstream distributions over $X_M$, although this was also not treated in enough detail.\n + The proposed application of LASSO style selection as an alternative to subset search is interesting, although this was also not treated in enough detail.\n + Experimental results seem promising and Figures 1, 6, and 7 are particularly compelling, although absolute numbers for MSE and comparisons of performance across different target distributions would be appreciated.\n\n**Weaknesses:**\n\n**Missing Technical Details**\n - The authors try to do too much in the paper, so many ideas are not treated thoroughly enough. While some results are shown with some rigor (e.g., the graphical criterion in Theorem 3.1, crucial details of other parts of the algorithm are missing. Several examples of missing details follow. I would suggest cutting down the claimed contributions of the paper, or expanding to a longer format venue.\n- How does one solve the maximization problem over $P_J$ that defines $h^*$ in practice? $P_J$ is a non-parametric class of functions that allows arbitrary dependence for the structural equation $Pa(X_M) \\rightarrow X_M$. This seems like a critical detail for the proposed algorithm to be usable. Does one need to make assumptions about $J$?\n - There appear to be missing assumptions for the structure identifiability claim in Proposition 3.4, and the proof for this proposition is not actually a proof, but the statement of an algorithm. At the very least, it seems important to have some assumption akin to positivity about the set of environments that is observed in training; as of now, there are no conditions stated about how the environments should be heterogeneous, or whether all environments need to change all mutable variables or something similar. I imagine there are many stronger conditions that are given in the Huang et al 2020 paper that is quoted extensively.\n - The g-equivalence algorithm seems quite critical given the claims in the paper, but is omitted from the main text. In particular, it seems important to understand when composing the g-equivalence algorithm with the reduced variable selection problem is actually more efficient than the original variable selection problem. I suspect it is, but it seems necessary to substantiate this claim. Meanwhile, the complexity analysis section goes into examples that don't feel particularly relevant.\n - It is not clear to me whether the Estimate $f_{S_-}$ section is correct, or why the algorithm needs to be so complicated; if $f_{S_-}$ is identifiable, shouldn't we be able to to estimate it by obtaining a functional from the ID algorithm? This seems more straightforward than estimating and generating from structural equations, which also requires additional assumptions for correctness. At any rate, it needs to be proved that the proposed algorithm actually estimates the appropriate interventional distribution.\n - The sparse min-max optimization section has far too little detail, and it is not clear what claims the authors want to make here. None of the cited work has results about doing variable selection with the inner maximization (rather they assume risk is minimized across the same distribution regardless of the $\\alpha, \\beta$ parameters), nor does it extend beyond linear models. If this is meant to be a main contribution, more details of the exact approach should be given in the main text. It is not necessary to prove anything about this practical approach, but the citations to selection consistency work suggest that the authors want to make stronger claims than they can substantiate.\n - For the actual implementations, no details are given about how the predictors are actually estimated. Is everything a linear model? How are hyperparameters set (for example, the very important robustness parameter in Anchor Regression)?\n - How are bi-directed edges handled by the algorithm? There is no mention of these, and some of the proofs assume a topological ordering exists, but the causal structure discovery algorithm returns bidirected edges in some cases.\n\n**Other issues**\n - In the set-up, there is no formal definition of a \"stable predictor\". Currently there is a sentence stating that such a predictor can be \"transferred to a broader family of environments without any adjustment\", but it is not clear what \"transferred\" means in this sentence. What is the actual invariance criterion, i.e., the property of f(X) that is the same across environments?\n - I think there is an error in Example 1. When the dashed arrow is absent, by the do calculus rules 1 then 3, $p(y \\mid k, do(x^0_m)) = p(y \\mid do(x^0_m)) = p(y)$. This is consistent, I believe, with Theorem 3.1, which reduces to the conditional $p(y \\mid k_2)$, but $k_2$ is the empty set in this example.\n - Notation in the Proof of Theorem 3.3 is different from the notation in the statement.\n\n**Suggestions:**\n - Can the proof of Theorem 1 be written using the rules of do calculus? It seems like this would be more straightforward than going back to conditional probability statements.\n - It would be nice to have a simple toy example showing how changes to the $X_M$ conditional distributions can induce changes in the risk of $f_{S_-}$, specifically showing that $f_{S_-}$ with the lowest in-distribution validation risk may not have the lowest worst-case validation risk.\n - The causal structure discovery portion of the algorithm seems like a distraction. Most of the result is quoting identification results from Huang et al and not a core contribution of this work. In addition, the assumptions for the structure discovery algorithms to work are stronger than the assumptions provided here. I would suggest dropping the structure discovery portion from the main development, and discussing structure discovery as a practical tool for implementation in the experiments section, as the contribution of this paper is mostly about doing the min-max optimization after causal structure has been ascertained, and there is no analysis of, say, the propagation of error from causal structure discovery to the final predictor. This would provide room for other important details, such as the algorithm for identifying g-equivalence classes or the sparse selection algorithm.\n - Comparisons of MSE across environments (not just worst-case MSE) would be useful here, since there is often a tradeoff between the two. It would also be useful to see how the variables selected by Surgery and the proposed algorithm differ. Finally, it would be useful to understand whether the MSE of these estimators is actually in a range that would be useful for prediction. Comparing to, e.g., an oracle estimator that is trained on data from the target environment, or a baseline that just takes the mean outcome would be useful.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}]

Author rebuttal:

Decision: Reject