Interventional Causal Structure Discovery over Graphical Models with Convergence and Optimality Guarantees

Causal structure learning is defined as learning a DAG, represented as $\mathbb{G}(V,E)$, over the data $\boldsymbol{X}\in\mathbb{R}^{N\times D}$ sampled from a joint distribution $\mathbb{P}(X)$. In graph $\mathbb{G}$, each node $i\in V$ corresponds to a random variable $X_{i}\in\left\{X_{1},\ldots,X_{D}\right\}$, and edge $(i,j)\in E$ denotes a direct causal relationship from the variable $X_{i}$ to $X_{j}$, i.e., $X_{i}\rightarrow X_{j}$. The distribution of variable $X_{i}$ is $\mathbb{P}_{i}\left(X_{i}\mid\mathrm{Pa}\left(X_{i}\right)\right)$, where $\mathrm{Pa}\left(X_{i}\right)$ denotes the parent set of node $X_{i}$. Intervention on variable $X_{i}$ is defined as replacing the conditional distribution $\mathbb{P}\left(X_{i}\mid\mathrm{Pa}\left(X_{i}\right)\right)$ with a new distribution $\mathbb{P}_{\mathrm{new}}$, including perfect and imperfect interventions. Perfect intervention refers to removing the effects of all parent variables, i.e., $\mathbb{P}_{\mathrm{new}}=\widetilde{\mathbb{P}}(X_{i})$; while imperfect intervention replaces the original conditional distribution with a new conditional distribution, i.e., $\mathbb{P}_{\mathrm{new}}=\widetilde{\mathbb{P}}(X_{i}|\mathrm{Pa}(X_{i}))$.

There are many causal structure learning methods, which can be broadly categorized into FCMs, constraint-based and score-based methods. FCMs aim to learn and represent the causal relationship between cause and effect variables using a predefined function containing independent noise terms. Typical FCMs include LiNGAM [29], PNL [30] and ANM [31]. Constraint-based methods leverage conditional independence (CI) tests between the variables to identify causal relationships, such as PC and FCI[8]. COmbINE [32] and HEJ [33] support the introduction of interventional data, and they typically rely on Boolean satisfiability solvers. Mooij et al. [34] proposed the joint causal inference framework that can handle unknown interventions with by coupling with many constraint-based algorithms. Score-based methods aim to obtain a DAG by optimizing a score function. GES [11] searches for the highest scoring graphs from a discrete space $\mathbb{G}$ by iteratively adding, removing and flipping edges. GIES [35] and GNIES [36] are the variants of GES that can be used for interventional data. IGSP is a hybrid method [24]. However, these methods typically search in discrete spaces.

A bilevel optimization problem is an optimization problem with two levels, each of which has its own objective function and constraints. The general form can be denoted as:

where $\mathcal{F}$ and $\mathcal{G}$ are the upper- and lower-level objective function, $\boldsymbol{x}$ and $\boldsymbol{y}$ are decision variables, $\mathcal{f}_{i}$ is the $i_{th}$ upper-level constraint, and $\mathcal{g}_{j}$ denotes the $j_{th}$ lower-level constraint.

There are many approaches for bilevel optimization [45], here we focus on the method based on optimal conditions. This method aims to replace the lower-level optimization problem with the optimal conditions (e.g., Karush-Kuhn-Tucker (KKT) conditions) [46] and reformulates it as a Mathematical Problem with Complimentary Constraints, which can be expressed as follows.

Moreover, the problem is shown to be equivalent to the original bilevel optimization problem when Slater’s condition (Definition 1) is satisfied [47].

The general form of the polynomial optimization problem (POP) denotes as follows:

The objective function $\mathcal{F}$ and constraints $\mathcal{f}_{i}$ are assumed to be polynomials in the variables $\boldsymbol{x}=(x_{1},\dots,x_{N})$, expressed as $\mathcal{f}(\boldsymbol{x})=\sum_{\boldsymbol{\alpha}\in A}f_{\alpha}\boldsymbol{x}^{\alpha}$, where $f_{\alpha}\in\mathbb{R}$ represents the coefficients, and $A\subseteq\mathbb{N}^{N}$. The monomial $\boldsymbol{x}^{\boldsymbol{\alpha}}$ denotes $x_{1}^{\alpha_{1}}\dots x_{N}^{\alpha_{N}}$, where $\boldsymbol{\alpha}=(\alpha_{1},\dots,\alpha_{N})$ is the order of each variable $x_{n}$. The degree of monomial $\boldsymbol{x}^{\boldsymbol{\alpha}}$ is $\operatorname{deg}(\boldsymbol{x}^{\boldsymbol{\alpha}})=\sum_{i=1}^{N}\alpha_{i}$. The support set of a polynomial function $\mathcal{f}$ is defined as $\operatorname{supp}(\mathcal{f}):=\left\{\boldsymbol{\alpha}\in A\mid f_{\alpha}\neq 0\right\}$, and the degree of a polynomial function $\mathcal{f}$ is $\operatorname{deg}(\mathcal{f})=\max\{\operatorname{deg}(\boldsymbol{x}^{\boldsymbol{\alpha}}):\boldsymbol{\alpha}\in\operatorname{supp}(\mathcal{f})\}$. Let $\mathbb{R}\left[x\right]$ be the set of real $n$-variate polynomials, and the set of polynomials of degree no more than $2d$ is represented as $\mathbb{R}_{2d}[x]$. For $r\in\mathbb{Z}_{+}$, let $\mathbb{S}^{r\times r}$ (resp. $\mathbb{S}_{+}^{r\times r}$) denote the set of symmetric matrices (resp. positive semidefinite matrices).

Traditional approaches [19] achieve solutions by constructing a sum-of-squares (SOS) hierarchy for Eq. (4) and solving a sequence of corresponding SDP problems. The polynomial is expressed as $f(\boldsymbol{x})=(\boldsymbol{x}^{\mathbb{N}_{d}^{N}})^{T}\boldsymbol{Q}(\boldsymbol{x}^{\mathbb{N}_{d}^{N}})$, where $\boldsymbol{Q}\in\mathbb{S}_{+}^{r\times r}$ is the Gram Matrix [48]. Here, $\boldsymbol{x}^{\mathbb{N}_{d}^{N}}$ represents the standard monomial basis of $\boldsymbol{x}$ with a degree not exceeding $d$, denoted as $\mathbb{N}_{d}^{N}$ in the sequel. Let $\boldsymbol{y}$ be the univariate vector corresponding to $\mathbb{N}{d}^{N}$, then the monomial variable $\boldsymbol{x}^{\boldsymbol{\alpha}}$ can be replaced by a real variable $y_{\boldsymbol{\alpha}}$. Define $\mathcal{L}_{\boldsymbol{y}}:\mathbb{R}\left[\boldsymbol{x}\right]\rightarrow\mathbb{R}$ as the linear function:

Let $\boldsymbol{M}_{d}(\boldsymbol{y})$ be the $d$-order moment matrix associated with $\boldsymbol{y}$. The items in the matrix can be obtained as follows:

where $1\leq a,b\leq|\boldsymbol{y}|$ denote indexes, and $|\boldsymbol{y}|$ denotes the dimension of $\boldsymbol{y}$.

For the constraint function $\mathcal{f}_{i}(\boldsymbol{x})$, the elements in the $d$-order localizing matrix $\boldsymbol{M}_{d}(\mathcal{f}_{i}\boldsymbol{y})$ are presented as:

where $y_{\alpha}=\mathcal{L}_{y}(\boldsymbol{x}^{\alpha})$. Let $d_{i}=\lceil\deg(\mathcal{f}_{i})/2\rceil(i=1,...,I)$ and $d_{\mathrm{min}}=\max\{\lceil\mathrm{deg}(\mathcal{F})\rceil/2,d_{1},...,d_{l}\}$.

By introducing Eq. (6) and (7), the problem Eq. (4) can be relaxed into a SDP problem (with relaxation order $d\geq d_{\mathrm{min}}$):

However, it is more complicated to solve the above problem directly due to the scale of the issue [49].

We aim to learn a DAG of an causal graphical model given observational and intervention data. Here, we first make the causal faithful and causal sufficient assumptions for the causal model, i.e., the observed data is consistent with the real causal relationship, and all variables for inferring the causal relationship are included. However, in the following, we relax this assumption appropriately and experimentally discuss the causal structure discovery performance of the proposed method in scenarios with potentially confounding causal variables.

Furthermore, we assume that interventions are created by affecting only a single variable, and each intervention is independent of the other. The scope of this intervention closely follows Ke et al. [53]. Besides, we assume that all interventions are perfect, i.e., the distribution of the intervened variable is independent of the parent variables. And we assume interventions are conducted on each variable, and all samples, including the intervention target, are provided. In next section, we also try to relax this assumption and experimentally discuss the performance of the proposed method in the case of imperfect intervention.

Pearl’s Causal Hierarchy theory [1] posits that causal inference consists of three ascending levels: association, intervention, and counterfactual. Each level cannot provide higher-level information. In the first level, association problem is typically addressed by using passively observed data. On the second level, interventional data can offer more information of directed causality. However, learning causal structures relying solely on observational data is typically challenging. This is because, under the faithfulness assumption, the model may only be able to learn the Markov equivalence class of the true graph, whereas interventional data can effectively enhance identifiability.

Therefore, we can define two parameters to jointly determine the graph structure, respectively: $\boldsymbol{P}\in[-1,1]^{D\times D}$ denotes the existence of undirected edges, $\boldsymbol{Q}\in[0,1]^{D\times D}$ the direction of the edges. In particular, the diagonal entries of both $\boldsymbol{P}$ and $\boldsymbol{Q}$ are $0$, i.e., $P_{ii}$ (or $Q_{ii}$ )=$0$, $i\in[D]{:=}\{1,\ldots,D\}$, besides $P_{ij}=P_{ji}$, $Q_{ij}+Q_{ij}=1$. Then the weighted adjacency matrix of the DAG can be determined by $\boldsymbol{W}=\boldsymbol{P}\odot\boldsymbol{Q}$. When $W_{ij}=P_{ij}\cdot Q_{ij}=0$, it means there is no directed causal relationship between variable $i$ and $j$. By decomposing the adjacency matrix into these two parameters, we can model the causal structure learning with observational and interventional data as a bilevel optimization problem thereby improving the identifiability of the structure. However, this modelling may entail extensive matrix operations, particularly with acyclic constraints, resulting in heightened complexity and slow convergence. Let $\boldsymbol{p}\in[-1,1]^{D_{2}}$ and $\boldsymbol{q}\in[0,1]^{D_{1}}$ be the concatenation of the non-diagonal elements from the upper triangular matrices of the matrices $\boldsymbol{P}$ and $\boldsymbol{Q}$, respectively, where $\overline{D}=D_{1}=D_{2}=\frac{D^{2}-D}{2}$.Thus, the elements in $\boldsymbol{P}$ and $\boldsymbol{Q}$ can be easily expressed by $\boldsymbol{p}$ and $\boldsymbol{q}$. And we can model the above problem with the following bilevel polynomial optimization formulation.

where $\boldsymbol{X}_{obs}\in\mathbb{R}^{N_{obs}\times D}$ is the observational data with $D$ representing the number of random variables. $\boldsymbol{X}_{int}^{I_{t}}\in\mathbb{R}^{N_{int}\times D}$ represents the interventional data, where interventions are performed on variable $t\subset[D]$. And the interventional data under different variables is denoted as $\boldsymbol{X}_{int}=\{\boldsymbol{X}_{int}^{I_{t}}\}_{t=1}^{D}$. $\mathcal{F}$ and $\mathcal{f}_{i}$ denote the polynomial objective function and polynomial constraints of the upper-level problem, respectively. The function $\mathcal{F}$ is defined as the least squares loss, $\mathcal{F}=\mathcal{L}(\boldsymbol{q},\boldsymbol{p}|\boldsymbol{X}_{int},\boldsymbol{X}_{obs})=\sum_{t=1}^{D}\mathcal{L}(\boldsymbol{q},\boldsymbol{p}|\boldsymbol{X}_{int}^{I_{t}})+\alpha\mathcal{L}(\boldsymbol{q},\boldsymbol{p}|\boldsymbol{X}_{obs})$ ($\alpha$ is an empirical parameter, relevant details are given in the Appendix B), and the function $\mathcal{f}_{i}=\mathcal{h}_{i}(q,p)$ represents an acyclic constraint with different step lengths. Here, we use both observational and interventional data to learn variable $\boldsymbol{q}$ in the upper-level function, which not only fully leverages the information from observational data, but also improves the accuracy of results by incorporating interventional data. When calculating the upper-level objective function, we masked the intervened variable because the parent set of it changes during perfect intervention. $\mathcal{G}$ and $\mathcal{g}_{j}$ represent the polynomial objective function and polynomial constraints of the lower-level problem respectively, with $\mathcal{G}={\cal L}(\boldsymbol{q},\boldsymbol{p}|\boldsymbol{X}_{obs})+\lambda_{sp}{\cal L}_{sp}(\boldsymbol{p})$ and $\mathcal{g}_{j}(\boldsymbol{p})=-p_{j}^{2}(1-p_{j}^{2})$. ${\cal L}_{sp}$ is the regularization function. Moreover, due to the sparsity of the causal structure, in most scenarios, we can simplify $\boldsymbol{q}$ and $\boldsymbol{p}$ by preprocessing or even introducing expert knowledge. For example, when there is obviously no causal relationship between two variables, the corresponding edge $p_{i}=0$.

The proposed method can effectively utilize observational and interventional data, and improve the accuracy of causal structure learning. For some scenarios where there is no or limited intervention data, we can utilize data augmentation techniques to simulate the interventions to improve the usability of our method. M. Ilse [54] introduced Intervention-augmentation equivalence (IAE), demonstrating the feasibility of simulating interventional data through data augmentation using only observational data under IAE conditions. In an IAE causal process $f_{X}{:}\,{\mathcal{D}}\times\mathcal{Y}\to{\mathcal{X}}$, every stochastic data augmentation transformation $\operatorname{aug}(\cdot)$ on $x\in\mathcal{X}$ is equivalent to a corresponding noise intervention $do(\cdot)$ on $d\in\mathcal{D}$ such that: $\operatorname{aug}(f_{X}(d,y))=f_{X}(do(d),y)$. Hence, by verifying if the causal relationship between variables meets the IAE condition, an appropriate data augmentation model can be trained for each variable and simulate interventions through the addition of noise. More details can be found in the original publication.

Reformulation based on optimal conditions. Although Eq. (11) is a bilevel polynomial optimization problem with both upper and lower constraints, fortunately, the lower optimization problem can be proved convex because of the operations that preserve convexity [55]. According to Proposition 1, since the lower-level problem satisfies the Slater’s condition (Definition 1), Eq. (11) can be transformed into a single-level polynomial optimization problem.

We denote the inequality constraints and equality constraints as follows:

$\hat{\mathcal{f}}_{k_{1}}(\boldsymbol{q},\boldsymbol{p},\boldsymbol{\lambda})=\\ \quad\quad\begin{cases}-\mathcal{f}_{k_{1}}(\boldsymbol{q},\boldsymbol{p}),&k_{1}=1,...,I\\ \lambda_{k_{1}-1},&k_{1}=I+1,...,I+J+1\\ -\mathcal{g}_{j}(\boldsymbol{p}),&k_{1}=I+J+2,...,I+2J+1\end{cases}$ and

$\widehat{\mathcal{g}}_{k_{2}}(\boldsymbol{q},\boldsymbol{p},\boldsymbol{\lambda})=\\ \quad\quad\begin{cases}\left(\lambda_{0}\nabla_{\boldsymbol{p}}\mathcal{G}(\boldsymbol{q},\boldsymbol{p}|\boldsymbol{X}_{obs})+\sum_{j=1}^{J}\lambda_{j}\nabla_{\boldsymbol{p}}\mathcal{g}_{j}(\boldsymbol{p})\right)_{k_{2}},\\ \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad k_{2}=1,...,\overline{D}\\ \lambda_{j}\mathcal{g}_{j}(\boldsymbol{p}),\quad\quad\quad\quad\quad\quad\quad\quad k_{2}=\overline{D}+1,...,\overline{D}+J{,}\end{cases}$

then the problem Eq. (12) can be expressed as

Proof. In order to prove the above proposition, we need to show that for any $\boldsymbol{q}$, there exists $\boldsymbol{\lambda}\in\mathbb{R}^{j+1}$ such that the feasible set $\boldsymbol{p}\in\mathbb{P}\left(\boldsymbol{q}\right)$ of the lower optimization problem is equivalent to the following condition that

According to the Theorem 2.1 [56], under the slater condition, a point $\boldsymbol{p}$ is the global optimal point of the lower-level problem if and only if $\boldsymbol{p}$ is the KKT point, i.e., there is $a_{j}\geq 0,j=1,...,J$ such that

Then we have $\lambda_{0}=\frac{1}{\sqrt{1+\sum_{j=1}^{J}a_{j}^{2}}}$ and $\lambda_{j}=\frac{a_{j}}{\sqrt{1+\sum_{j=1}^{J}a_{j}^{2}}}$ ($j=1,...,J$), where $a_{j}>0$. $\square$

In this subsection, we investigate how to solve the above single-level polynomial optimization problem Eq. (13), which can be outlined into three steps. First, the variables are partitioned according to the correlative sparsity; then, the term sparsity graph is constructed according to the correlation between the monomials; finally, the SDP relaxation problem of the POP problem is provided and solved.

However, due to privacy concerns, local datasets are not allowed to be uploaded to a central server in some scenarios [65, 44, 43]. Therefore, we attempt to further extend the Bloom algorithm to distributed settings, enabling it to learn the graph structure from distributed data without sharing locally stored data.

Distributed Data. Let $C=\{c_{1},c_{2},...,c_{M}\}$ be the client set which includes $M$ different clients, and $S$ represent the central server. The dataset $X^{c_{m}}=\{\boldsymbol{X}^{c_{m}}_{obs},\boldsymbol{X}^{c_{m}}_{int}\}$ represent the local data owned by the client $c_{m}$. The dataset $X=\{X^{c_{1}},X^{c_{2}},...,X^{c_{M}}\}$ is called distributed dataset. And we define $X$ as a homogeneous distributed dateset, which means they are sampled from an identical distribution.

To learn causal structures from distributed data, the distributed Bloom solves each subproblem Eq. (19) by distributing it across all local clients. Since data is not shared between clients and the server, data privacy is significantly ensured. During training, the server and clients will exchange updated variables in each communication round to facilitate coordinated joint learning of the causal structure.

Client Update. Essentially, solving the optimization problem on each client can be seen as an independent process. Therefore, each client will calculate sub-problems with local dataset, following the steps in Algorithm 1. After each communication round, the client will receive the updated variables, which will be used as the initial values for the next round of solving problem.

Server Update. After clients completed local updates, the server randomly selects $r$ clients and collects the learned variables to the set $\widehat{x}_{r}=\{\widehat{\boldsymbol{x}}^{c_{1}},\widehat{\boldsymbol{x}}^{c_{2}},...,\widehat{\boldsymbol{x}}^{c_{r}}\}$. By averaging the collected variables, we will get $\widehat{\boldsymbol{x}}_{new}$, which is then distributed to all clients.

The specific steps are shown in the Algorithm 2.

Privacy Protection. To avoid leakage of raw data in the client, distributed Bloom only exchanges the learned variables during the training process. Therefore, we consider the information leakage of local data to be relatively limited. For the graph structure information that may be contained in the transmission, we can address it by selecting a client as a proxy server [44]. It is worth mentioning that our work only provides a possible distributed extension of Bloom. Regarding further privacy protection efforts, we could introduce more advanced privacy protection techniques [66], which will be the focus of future research.

Communication Cost. Distributed Bloom only requires exchanging parameters between the server and clients during communication rounds. Despite introducing some communication costs, we consider this to be acceptable with relatively low communication overhead. In each communication round, servers only need to collect learned parameters from selected clients and distribute updated variables to each client. Furthermore, the trade-off between performance and communication costs can also be controlled by choosing the number of selected clients $r$.

Additionally, for some large-scale structure learning problems, solving each subproblem still imposes significant demands on the computational resources of individual clients [42]. Fortunately, due to the sparsity of the constructed SDP problems Eq.(18), we can apply many existing distributed solving methods [67, 68, 69]. Although this is not the main focus of this article, we also present a feasible approach in Appendix B.

We conduct experiments on synthetic and real data respectively to demonstrate its effectiveness by comparing the proposed method with some typical methods. All experiments are implemented on a server with Intel(R) Xeon(R) CPU E5-2690 v3 CPUs.

We tested the overall performance of all methods on two synthetic datasets and real datasets respectively.