Consistent Second-Order Conic Integer Programming for Learning Bayesian Networks

A Bayesian network (BN) is a probabilistic graphical model consisting of a labeled directed acyclic graph (DAG) $\mathcal{G}=(V,E)$, in which the vertex set $V=\{V_{1},\dots,V_{m}\}$ corresponds to $m$ random variables, and the edge set $E$ prescribes a decomposition of the joint probability distribution of the random variables based on their parents in $\mathcal{G}$. The edge set $E$ encodes Markov relations on the nodes in the sense that each node is conditionally independent of its non-descendents given its parents. BNs have been used in knowledge discovery [49, 9], classification [1], feature selection [22], latent variable discovery [29] and genetics [37]. They also play a vital part in causal inference [40].

In this paper, we propose reformulations of the mixed-integer quadratic programs (MIQP) for learning the optimal DAG structure of BNs given $n$ continuous observations from a system of linear structural equation models (SEMs). While there exist exact integer-programming (IP) formulations for learning DAG structure with discrete data [11, 12, 27, 50, 4, 34, 35, 5, 13, 14], the development of tailored computational tools for learning the optimal DAG structure from continuous data has received less attention. In principle, exact methods developed for discrete data can be applied to continuous data. However, such methods result in exponentially sized formulations in terms of the number of binary variables. A common practice to circumvent the exponential number of binary variables is to limit the in-degree of each node [12, 14, 5]. But, this may result in sub-optimal solutions. On the contrary, MIQP formulations for learning DAGs corresponding to linear SEMs require a polynomial number of binary variables. This is because for BNs with linear SEMs, the score function — i.e., the penalized negative log-likelihood (PNL) — can be explicitly written as a function of the coefficients of linear SEMs [45, 52, 38, 32]. In contrast to the existing MIQPs [38, 32], our reformulations exploit the convex quadratic objective and the relationship between the continuous and binary variables to improve the strength of the continuous relaxations.

Continuous BNs with linear SEMs have witnessed a growing interest in the statistics and computer science communities [52, 43, 31, 23, 47]. In particular, it has been shown that the solution obtained from solving the PNL augmented by $\ell_{0}$ regularization, which introduces a penalty on the number of non-zero arc weights in the estimated DAG, achieves desirable statistical properties [41, 52, 31]. Moreover, if the model is identifiable [41, 31], that is when the true causal graph can be identified from the joint distribution, then such a solution is guaranteed to uncover the true causal DAG when the sample size $n$ is large enough. However, given the difficulty of obtaining exact solutions, existing approaches for learning DAGs from linear SEMs have primarily relied on heuristics, using techniques such as coordinate descent [21, 2, 26] and non-convex continuous optimization [61]. Unfortunately, these heuristics are not guaranteed to achieve the desirable properties of the global optimal solution. Moreover, it is difficult to evaluate the statistical properties of a sub-optimal solution with no optimality guarantees [28]. To bridge this gap, in this paper we develop mathematical formulations for learning optimal BNs from linear SEMs using a PNL objective with $\ell_{0}$ regularization. By connecting the optimality gap of the mixed-integer program to the statistical properties of the solution, we also establish an early stopping criterion under which we can terminate the branch-and-bound procedure and attain a solution which asymptotically recovers the true parameters with high probability.

Our work is related to recent efforts to develop exact tailored methods for DAG learning from continuous data. [57] show that $A^{\ast}$-lasso algorithm tailored for DAG structure learning from continuous data with $\ell_{1}$-regularization, which introduces a penalty on the sum of absolute values of the arc weights, is more effective than the previous approaches based on dynamic programming [e.g., 46] that are suitable for both discrete and continuous data. [38] develop a mathematical program for DAG structure learning with $\ell_{1}$ regularization. [32] improve and extend the formulation by [38] for DAG learning from continuous data with both $\ell_{0}$ and $\ell_{1}$ regularizations. The numerical experiments by [32] demonstrate that as the number of nodes grows, their MIQP formulation outperforms $A^{\ast}$-lasso and the existing IP methods; this improvement is both in terms of reducing the IP optimality gap, when the algorithm is stopped due to a time limit, and in terms of computational time, when the instances can be solved to optimality. In light of these recent efforts, the current paper makes important contributions to this problem at the intersection of statistics and optimization.

• •

  The statistical properties of optimal PNL with $\ell_{0}$ regularization have been studied extensively [31, 52]. However, it is often difficult to obtain an optimal solution and no results have been established on the statistical properties of approximate solutions. In this paper, we give an early stopping criterion for the branch-and-bound process under which the approximate solution gives consistent estimates of the true coefficients of the linear SEM. Our result leverages the statistical consistency of the PNL estimate with $\ell_{0}$ regularization [52, 41] along with the properties of the branch-and-bound method wherein both lower and upper bound values on the objective function are available at each iteration. By connecting these two properties, we obtain a concrete early stopping criterion, as well as a proof of consistency of the approximate solution. To the best of our knowledge, this result is the first of its kind for DAG learning.

• •

  In spite of recent progress, a key challenge in learning DAGs from linear SEMs is enforcing bounds on arc weights. This is commonly modeled using the standard “big-$M$ constraint” approach [38, 32]. As shown by [32], this strategy leads to poor continuous relaxations for the problem, which in turn results in slow lower bound improvement in the branch-and-bound tree. In particular, [32] establish that all existing big-$M$ formulations achieve the same continuous relaxation objective function under a mild condition (see Proposition 4). To circumvent this issue, we present a mixed-integer second-order cone program (MISOCP), which gives a tighter continuous relaxation than existing big-$M$ formulations under certain conditions discussed in detail in Section 5.3. This formulation can be solved by powerful state-of-the-art optimization packages. Our numerical results show the superior performance of MISOCP compared to the existing big-$M$ formulations in terms of improving the lower bound and reducing the optimality gap. We also compare our method against the state-of-the-art benchmarks [9, 24] both for identifiable and non-identifiable instances, and show that our method provides the best estimation among all methods in most of the networks, especially for the non-identifiable cases.

The rest of the paper is organized as follows. In Section 2, we define the DAG structure learning problem corresponding to linear SEMs, and give a general framework for the problem. In Section 3, we present our early stopping criterion and establish the asymptotic properties of the solution obtained under this stopping rule. We review existing mathematical formulations in Section 4, and present our proposed mathematical formulations in Section 5. Results of comprehensive numerical studies are presented in Section 6. We end the paper with a summary in Section 7.

Let $\mathcal{M}=(V,E)$ be an undirected and possibly cyclic super-structure graph with node set $V=\{1,2,\dots,m\}$ and edge set $E\subseteq V\times V$; let $\overrightarrow{\mathcal{M}}=(V,\overrightarrow{E})$ be the corresponding bi-directional graph with $\overrightarrow{E}=\{(j,k),(k,j)|(j,k)\in E\}$. We refer to undirected edges as edges and directed edges as arcs.

Let $\mathcal{X}=(\mathcal{X}_{1},\dots,\mathcal{X}_{m})$ be the $n\times m$ data matrix with $n$ rows representing i.i.d. samples from each random variable, and $m$ columns representing random variables $X_{1},\ldots,X_{m}$. The linear SEM (1) can be compactly written in matrix form as $\mathcal{X}=\mathcal{X}{B}+\mathcal{E}$, where ${B}=[\beta]\in\mathbb{R}^{m\times m}$ is a matrix with $\beta_{kk}=0$ for $k=1,\dots,m$, $\beta_{jk}=0$ for all $(j,k)\notin E$, and $\mathcal{E}$ is the $n\times m$ ‘noise’ matrix. Then, $\mathcal{G}(B)$ denotes the directed graph on $m$ nodes such that arc $(j,k)$ appears in $\mathcal{G}(B)$ if and only if $\beta_{jk}\neq 0$. Throughout the paper, we will use $B$ and $\beta$ to denote the matrix of coefficients and its vectorized version.

A key challenge when estimating DAGs by minimizing the loss function (2) is that the true DAG is generally not identifiable from observational data. However, for certain SEM distributions, the true DAG is identifiable from observational data; that is when the true causal graph can be identified from the joint distribution. Two important examples are linear SEMs with possibly non-Gaussian homoscedastic noise variables [41], as well as linear SEMs with unequal noise variances that are known up to a constant [31]. In these special cases, the true DAG can be identified from observational data, without requiring the (strong) ‘faithfulness’ assumption, which is known to be restrictive in high dimensions [51, 48]. Given these important implications, in this paper we focus on learning Bayesian networks corresponding to the above identifiable linear SEMs, i.e., settings where the error variances are either equal, or known up to a constant.

The negative log likelihood for an identifiable linear SEM (1) with equal noise variances is proportional to

here $\widehat{\Sigma}=n^{-1}\mathcal{X}^{\top}\mathcal{X}$ is the empirical covariance matrix, and $I$ is the identity matrix [45, 52].

To learn sparse DAGs, van de Geer and Bühlmann [52] propose to augment the negative log likelihood with an $\ell_{0}$ regularization term. Given a super-structure $\mathcal{M}$, the optimization problem corresponding to this penalized negative log-likelihood (PNL$\mathcal{M}$) is given by

where the tuning parameter $\lambda_{n}$ controls the degree of the $\ell_{0}$ regularization

where $\mathbbm{1}(\beta_{jk})$ is an indicator function with value one if $\beta_{jk}\neq 0$, and 0 otherwise. The constraint (3b) stipulates that the resulting directed subgraph is a DAG induced from $\overrightarrow{\mathcal{M}}$. When $\mathcal{M}$ corresponds to a complete graph, PNL$\mathcal{M}$ reduces to the original PNL of van de Geer and Bühlmann [52].

The choice of $\ell_{0}$ regularization in (3) is deliberate. Although $\ell_{1}$ regularization has attractive computational and statistical properties in high-dimensional regression [8], many of these advantages disappear in the context of DAG structure learning [21, 2]. By considering $\ell_{0}$ regularization, [52] establish the consistency of PNL under appropriate assumptions. More specifically, for a Gaussian SEM, they show that the estimated DAG has (asymptotically) the same number of edges as the DAG with minimal number of edges (minimal-edge I-MAP), and establish the consistency of PNL for learning sparse DAGs. These results are formally stated in Proposition 1 in the next section.

In our earlier work [32], we observe that existing mathematical formulations are slow to converge to a provably optimal solution, $\beta^{\star}$, of (3) using the state-of-the-art optimization solvers. Therefore, the solution process needs to be terminated early to yield a feasible solution, $\hat{\beta}$ with a positive optimality gap, i.e., a positive difference between the upper bound on $\operatorname{\mathcal{L}}(\beta^{\star})$ provided by $\operatorname{\mathcal{L}}(\hat{\beta})$ and a lower bound on $\operatorname{\mathcal{L}}(\beta^{\star})$ provided by the best continuous relaxation obtained by the branch-and-bound algorithm upon termination. However, statistical properties of such a sub-optimal solution are not well-understood. Therefore, there exists a gap between theory and computation: while the optimal solution has nice statistical properties, the properties of the solutions obtained from approximate computational algorithms are not known. Moreover, due to the non-convex and complex nature of the problem, characterizing the properties of the solutions provided by heuristics is especially challenging. In the next section, we bridge this gap by developing a concrete early stopping criterion and establishing the consistency of the solution obtained using this criterion.

In this section, we establish a sufficient condition for the approximate solution of PNL$\mathcal{M}$, $\hat{\beta}$ to be consistent for the true coefficients, $\beta^{0}$; that is $\|\beta^{0}-\hat{\beta}\|_{2}^{2}=\mathcal{O}\left(s^{0}\log(m)/n\right)$, where $s^{0}$ is the number of arcs in the true DAG, and $x\asymp y$ means that $x$ converges to $y$ asymptotically. This result is obtained by leveraging an important property of the branch-and-bound process for integer programming that provides both lower and upper bounds on the objective function $\operatorname{\mathcal{L}}(\beta^{\star})$ upon early stopping, as well as the consistency results of the PNL estimate with $\ell_{0}$ regularization. Using the insight from this new result, we then propose a concrete stopping criterion for terminating the branch-and-bound process that results in consistent parameter estimates.

Let $\mathrm{LB}$ and $\mathrm{UB}$, respectively, denote the lower and upper bounds on the optimal objective function value (3a) obtained from solving (3) under an early stopping criterion (i.e., when the obtained solution is not necessarily optimal). We define the difference between the upper and lower bounds as the absolute optimality gap: $\mathrm{GAP}=\mathrm{UB}-\mathrm{LB}$. Let $\hat{\mathcal{G}}$ and $\hat{\beta}$ denote the structure of the DAG and coefficients of the arcs from optimization model (3) under the early stopping condition with sample size $n$ and regularization parameter $\lambda_{n}$. Let ${\mathcal{G}^{\star}}$ and $\beta^{\star}$ denote the DAG structure and coefficients of arcs obtained from the optimal solution of (3), and $\mathcal{G}^{0}$ and $\beta^{0}$ denote the true DAG structure and the coefficient of arcs, respectively. We denote the number of arcs in $\hat{\mathcal{G}}$, $\mathcal{G}^{0}$, and ${\mathcal{G}^{\star}}$ by $\hat{s}$, $s^{0}$, and $s^{\star}$, respectively. The score value in (3a) of each solution is denoted by $\operatorname{\mathcal{L}}(\phi)$ where $\phi\in\{\beta^{\star},\hat{\beta},\beta^{0}\}$.

Next, we present our main result. Our result extends van de Geer and Bühlmann’s result on consistency of PNL$\mathcal{M}$ for the optimal, but computationally unattainable, estimator, $\beta^{\star}$ to an approximate estimator, $\hat{\beta}$, obtained from early stopping. In the following (including the statement of our main result, Proposition 2), we assume that the super-structure $\mathcal{M}$ is known a priori. The setting where $\mathcal{M}$ is estimated from data is discussed at the end of the section. We begin by stating the key result from [52] and the required assumptions. Throughout, we consider a Gaussian linear SEM of the form (1). We denote the variance of error terms, $\epsilon_{j}$, by $\sigma_{jj}^{2}$ and the true covariance matrix of the set of random variables, $(X_{1},\ldots,X_{m})$ by the $m\times m$ matrix $\Sigma$.

Here, $\lambda=\lambda_{n}/n$, because the loss function (2) is that of [52] scaled by the sample size $n$. The next result establishes the consistency of the approximate estimator, $\hat{\beta}$, obtained using our proposed early stopping strategy.

Proposition 2 suggests that the branch-and-bound algorithm can be stopped by setting a threshold $c^{\star}n\lambda s^{0}$ on the value of $\mathrm{GAP}=|\mathrm{UB}-\mathrm{LB}|$ for a constant $c^{\star}>0$, say $c^{\star}=1$. Such a solution will then achieve the same desirable statistical properties (in terms of parameter consistency) as the optimal solution $\beta^{\star}$. While $\lambda$ can be chosen data-adaptively (as discussed in Section 6), both of these choices depend on the value of $s^{0}$, which is not known in practice. However, one can find an upper bound for $s^{0}$ based on the number of edges in the super-structure $\mathcal{M}$. In particular, if $\mathcal{M}$ is the moral graph [40] with $s_{m}$ edges, then $s^{0}\leq s_{m}$. However, while in many applications $s_{m}\asymp s^{0}$, this is not always guaranteed. Thus, to ensure consistent estimation when replacing $s^{0}$ with $s_{m}$ and setting $c^{\star}=1$ in practice, we use the more conservative threshold of $\lambda s^{0}\asymp s^{0}\log(m)/n$. With this choice the first and third terms in (3) would be of the same (vanishing) order, and the consistency rate would be driven by the convergence rate of $\|\beta^{\star}-\beta^{0}\|_{2}^{2}$. We investigate the performance of this choice in Section 6.4.

The above results, including the specific choice of early stopping criterion, are also valid if the super-structure $\mathcal{M}$ corresponding to the moral graph is not known a priori. That is because the moral graph can be consistently estimated from data using recent developments in graphical modeling; see Drton and Maathuis [16] for a review of the literature. While some of the existing algorithms based on $\ell_{1}$-penalty require an additional irrepresentability condition [33, 44], this assumption can be relaxed by using instead an adaptive lasso penalty or by thresholding the initial lasso estimates [8].

In light of Proposition 2, it is of great interest to develop algorithms that converge to a solution with a small optimality gap expeditiously. To achieve this, one approach is to obtain better lower bounds using the branch-and-bound process from strong mathematical formulations for (3). To this end, we next review existing formulations of (3).

In this section, we review known mathematical formulations for DAG learning with linear SEMs. We first outline the necessary notation below.

Index Sets
$V=\{1,2,\dots,m\}$: index set of random variables;
$\mathcal{D}=\{1,2,\dots,n\}$: index set of samples.
Input
$\mathcal{M}=(V,E)$: an undirected super-structure graph (e.g., the moral graph);
$\overrightarrow{\mathcal{M}}=(V,\overrightarrow{E})$: the bi-directional graph corresponding to the undirected graph $\mathcal{M}$;
$\mathcal{X}=(\mathcal{X}_{1},\dots,\mathcal{X}_{m})$, where $\mathcal{X}_{v}=(x_{1v},x_{2v},\dots,x_{nv})^{\top}$ and $x_{dv}$ denotes $d$th sample ($d\in\mathcal{D}$) of random variable $X_{v}$; note $\mathcal{X}\in\mathbb{R}^{n\times m}$;
$\lambda_{n}:$ tuning parameter (penalty coefficient for $\ell_{0}$ regularization).

Let $F(\beta,g)=\sum_{k\in V}\sum_{d\in\mathcal{D}}\Big{(}x_{dk}-\sum_{(j,k)\in\overrightarrow{E}}\beta_{jk}x_{dj}\Big{)}^{2}+\lambda_{n}\sum_{(j,k)\in\overrightarrow{E}}g_{jk}$. The PNL$\mathcal{M}$ problem can be cast as the following optimization model:

There are two sources of difficulty in solving (13a)-(13d): (i) the acyclic nature of DAG imposed by the combinatorial constraints in (13b); (ii) the set of nonlinear constraints in (13c), which stipulates that $\beta_{jk}\neq 0$ only if there exists an arc $(j,k)$ in $\mathcal{G}(B)$. In Section 4.1, we discuss related studies to address the former, whereas in Section 4.2 we present relevant literature for the latter.

In this section, we discuss improved mathematical formulations for learning DAG structure of a BN based on convex (instead of linear) encodings of the constraints in (13c).

Problem (13) is an MIQP with non-convex complementarity constraints (13c), a class of problems which has received a fair amount of attention from the operations research community over the last decade [17, 18, 19, 20, 25, 30, 55, 56]. There has also been recent interest in leveraging these developments to solve sparse regression problems with $\ell_{0}$ regularization [42, 15, 58, 3, 54].

Next, we review applications of MIQPs with complementarity constraints of the form (13c) for solving sparse regression with $\ell_{0}$ regularization. [20] develop a so-called projected perspective relaxation method, to solve the perspective relaxation of mixed-integer nonlinear programming problems with a convex objective function and complementarity constraints. This reformulation requires that the corresponding binary variables are not involved in other constraints. Therefore, it is suitable for $\ell_{0}$ sparse regression, but cannot be applied for DAG structure learning. [42] show how a broad class of $\ell_{0}$-regularized problems, including sparse regression as a special case, can be formulated exactly as optimization problems. The authors use the Tikhonov regularization term $\mu\|\beta\|_{2}^{2}$ and convex analysis to construct an improved convex relaxation using the reverse Huber penalty. In a similar vein, [6] exploit the Tikhonov regularization and develop an efficient algorithm by reformulating the sparse regression mathematical formulation as a saddle-point optimization problem with an outer linear integer optimization problem and an inner dual quadratic optimization problem which is capable of solving high-dimensional sparse regressions. [58] apply the perspective formulation of sparse regression optimization problem with both $\ell_{0}$ and the Tikhonov regularizations. The authors establish that the continuous relaxation of the perspective formulation is equivalent to the continuous relaxation of the formulation given by [6]. [15] propose perspective relaxation for $\ell_{0}$ sparse regression optimization formulation and establish that the popular sparsity-inducing concave penalty function known as the minimax concave penalty [59] and the reverse Huber penalty [42] can be obtained as special cases of the perspective relaxation – thus the relaxations of formulations by [59, 42, 6, 58] are equivalent. The authors obtain an optimal perspective relaxation that is no weaker than any perspective relaxation. Among the related approaches, the optimal perspective relaxation by [15] is the only one that does not explicitly require the use of Tikhonov regularization.

The perspective formulation, which in essence is a fractional non-linear program, can be cast either as a mixed-integer second-order cone program (MISOCP) or a semi-infinite mixed-integer linear program (SIMILP). Both formulations can be solved directly by state-of-the-art optimization packages. [15] and [3] solve the continuous relaxations and then use a heuristic approach (e.g., rounding techniques) to obtain a feasible solution (an upper bound). In this paper, we directly solve the MISOCP and SIMILP formulations for learning sparse DAG structures.

Next, we present how perspective formulation can be suitably applied for DAG structure learning with $\ell_{0}$ regularization. We further cast the problem as MISOCP and SIMILP. To this end, we express the objective function (16a) in the following way: