Integer Programming for Learning Directed Acyclic Graphs from Continuous Data

TL;DR

This paper introduces a new mixed-integer quadratic optimization model for learning optimal directed acyclic graphs from continuous data, improving computational efficiency and scalability over existing methods.

Contribution

The paper proposes a layered network formulation that efficiently incorporates super-structures and regularizations, advancing the optimization approach for DAG learning.

Findings

The LN formulation outperforms existing models in computational time.

It scales better than algorithms using only  regularization.

The model effectively incorporates super-structures to reduce candidate DAGs.

Abstract

Learning directed acyclic graphs (DAGs) from data is a challenging task both in theory and in practice, because the number of possible DAGs scales superexponentially with the number of nodes. In this paper, we study the problem of learning an optimal DAG from continuous observational data. We cast this problem in the form of a mathematical programming model which can naturally incorporate a super-structure in order to reduce the set of possible candidate DAGs. We use the penalized negative log-likelihood score function with both $\ell_0$ and $\ell_1$ regularizations and propose a new mixed-integer quadratic optimization (MIQO) model, referred to as a layered network (LN) formulation. The LN formulation is a compact model, which enjoys as tight an optimal continuous relaxation value as the stronger but larger formulations under a mild condition. Computational results indicate that the proposed formulation outperforms existing mathematical formulations and scales better than available algorithms that can solve the same problem with only $\ell_1$ regularization. In particular, the LN formulation clearly outperforms existing methods in terms of computational time needed to find an optimal DAG in the presence of a sparse super-structure.

Figures

Tables

Table 1. Table 1 : The number of binary variables and the number of constraints

	Incomplete (moral) $ℳ$					Complete (moral) $ℳ$
	CP	LO	TO	LN		CP	LO	TO	LN
$\#$ Binary Vars (${\ell }_0$)	$\left|\overrightarrow{E}\right|$	$\left|\overrightarrow{E}\right|+\left(\frac{m}{2}\right)$	$m^2+\left|E\right|+\left|\overrightarrow{E}\right|$	$\left|E\right|+\left|\overrightarrow{E}\right|$		$2\left(\frac{m}{2}\right)$	$\left(\frac{m}{2}\right)$	$m^2+3\left(\frac{m}{2}\right)$	$3\left(\frac{m}{2}\right)$
$\#$ Binary Vars (${\ell }_1$)	$\left|E\right|$	$\left(\frac{m}{2}\right)$	$m^2+\left|E\right|$	$\left|E\right|$		$\left(\frac{m}{2}\right)$	$\left(\frac{m}{2}\right)$	$m^2+\left(\frac{m}{2}\right)$	$\left(\frac{m}{2}\right)$
$\#$ Constraints
(both ${\ell }_0$ and ${\ell }_1$)	Exp	2$\left(\frac{m}{3}\right)$	$\left|\overrightarrow{E}\right|+2m$	$\left|\overrightarrow{E}\right|$		2$\left(\frac{m}{3}\right)$	$2\left(\frac{m}{3}\right)$	$\left(\frac{m}{2}\right)+2m$	$\left(\frac{m}{2}\right)$

Table 2. Table 2 : Computational performance of LN versus A ⋆ -algorithm with ℓ 1 subscript ℓ 1 \ell_{1} regularization for λ ∈ { 0 , 0.1 } 𝜆 0 0.1 \lambda\in\{0,0.1\}

				Moral $\lambda =0$			Complete $\lambda =0$				Moral $\lambda =0.1$			Complete $\lambda =0.1$
Graphs (Data sets)	m	$\left|ℳ\right|$		A⋆-lasso	LN		A⋆-lasso	LN			A⋆-lasso	LN		A⋆-lasso	LN
dsep	6	16		0.017	0.378		0.0261	0.388			0.455	0.025		0.429	0.108
Asia	8	40		0.016	0.401		0.152	0.887			0.195	0.071		0.191	0.319
Bowling	9	36		0.018	0.554		0.467	1.31			0.417	0.225		0.489	0.291
Insurancesmall	15	76		391	2.719		547.171	613.543			2.694	1.135		3.048	0.531
Rain	14	70		119.15	1.752		101.33	246.25			51.737	0.632		69.404	3.502
Cloud	16	58		4433.29	1.421		18839	810.471			1066.08	0.426		2230.035	7.249
Funnel	18	62		6 hrs	1.291		6 hrs	6 hrs (.002)			6 hrs	0.395		6 hrs	3.478
Galaxy	20	76		6 hrs	1.739		6 hrs	6 hrs (.005)			6 hrs	0.740		6 hrs	9.615
Insurance	27	168		6 hrs	131.741		6 hrs	6 hrs (.162)			6 hours	12.120		6 hrs	6 hrs (.031)
Factors	27	310		6 hrs	6 hrs (.001)		6 hrs	6 hrs (.081)			6 hours	55.961		6 hrs	6 hrs (.01)

Equations

$$X_{k}=f_{k}\big{(}pa^{\mathcal{G}_{0}}_{k},\delta_{k}\big{)},\quad k=1,\dots,m,$$

$$X_k=\sum _{j\in p{a}_k^{{\mathcal{G}}_0}}{\beta }_{jk}{X}_j+{\delta }_k,k=1,\dots ,m,$$

$$\mathcal{X}=\mathcal{X}B+\Delta ,$$

$$l_n(\beta )=\frac{1}{2}\text{tr}\{(I-B)(I-B{)}^{T}S\},$$

$$\mathcal{M}$$

$$\mathcal{G}(B)\text{induces a DAG from}\overrightarrow{\mathcal{M}}.$$

$$pr({\hat{\mathcal{G}}}_n={\mathcal{G}}^0)\to 1(n\to \infty ),$$

$$\text{ℓ0-CP}\min$$

$$-M{g}_{jk}\le {\beta }_{jk}\le M{g}_{jk},$$

$$\sum _{(j,k)\in {\mathcal{C}}_A}{g}_{jk}\le |{\mathcal{C}}_A|-1,$$

$$g_{jk}\in \{0,1\},$$

$$g_{ij}+g_{jk}+g_{ki}\le 2,$$

$$g_{jk}+g_{kj}\le 1,$$

$$\text{ℓ0-LO}\min$$

$$-M{g}_{jk}\le {\beta }_{jk}\le M{g}_{jk},$$

$$g_{jk}\le w_{jk},$$

$$w_{jk}+w_{kj}=1,$$

$$w_{ij}+w_{jk}+w_{ki}\le 2,$$

$$w_{jk}\in \{0,1\},$$

$$g_{jk}\in \{0,1\},$$

$$\text{ℓ0-TO}\min$$

$$-M{g}_{jk}\le {\beta }_{jk}\le M{g}_{jk},$$

$$g_{jk}\le z_{jk},$$

$$z_{jk}+z_{kj}\le 1,$$

$$z_{jk}-m{z}_{kj}\le \sum _{s\in V}s(o_{ks}-o_{js}),$$

$$\sum _{s\in V}{o}_{rs}=1,$$

$$\sum _{r\in V}{o}_{rs}=1,$$

$$z_{jk}\in \{0,1\},$$

$$o_{rs}\in \{0,1\},$$

$$g_{jk}\in \{0,1\},$$

$$\min$$

$$-M{g}_{jk}\le {\beta }_{jk}\le M{g}_{jk},$$

$$g_{jk}\le z_{jk},$$

$$z_{jk}+z_{kj}=1,$$

$$z_{jk}-(m-1)z_{kj}\le {\psi }_k-{\psi }_j,$$

$$z_{jk}\in \{0,1\},$$

$$1\le {\psi }_k\le m,$$

$$g_{jk}\in \{0,1\},$$

$$\min$$

$$-M{z}_{jk}\le {\beta }_{jk}\le M{z}_{jk}$$

$$\text{eqref}LN-con4-\text{eqref}LN-con6.$$

$$z_{jk}-(m-1)z_{kj}\le \sum _{j\in V\setminus k}{z}_{jk}-\sum _{k\in V\setminus j}{z}_{kj}\forall j,k\in Vj\ne k.$$

$$1=z_{12}-m{z}_{21}\le {\psi }_2-{\psi }_1,$$

$$1=z_{23}-m{z}_{32}\le {\psi }_3-{\psi }_2,$$

$$⋮$$

$$1=z_{p-1,p}-m{z}_{p,p-1}\le {\psi }_p-{\psi }_1,$$

Full text

1 Introduction

The study of Probabilistic Graphical Models (PGMs) is an essential topic in modern artificial intelligence [23]. A PGM is a rich framework that represents the joint probability distribution and dependency structure among a set of random variables in the form of a graph. Once learned from data or constructed from expert knowledge, PGMs can be utilized for probabilistic reasoning tasks, such as prediction; see [23, 24] for comprehensive reviews of PGMs.

Two most common classes of PGMs are Markov networks (undirected graphical models) and Bayesian networks (directed graphical models). A Bayesian Network (BN) is a PGM in which the conditional probability relationships among random variables are represented in the form of a Directed Acyclic Graph (DAG). BNs use the richer language of directed graphs to model probabilistic influence among random variables that have clear directionality; they are particularly popular in practice, with applications in genetics [45], biology [27], machine learning [23], and causal inference [39].

Learning BNs is a central problem in machine learning. An essential part of this problem entails learning the DAG structure that accurately represents the (hypothetical) joint probability distribution of the BN. Although one can form the DAG based on expert knowledge, acquisition of knowledge from experts is often costly and nontrivial. Hence, there has been considerable interest in learning the DAG directly from observational data [3, 6, 9, 10, 17, 29, 34, 39].

Learning the DAG which best explains observed data is an NP-hard problem [5]. Despite this negative theoretical result, there has been interest in developing methods for learning DAGs in practice. There are two main approaches for learning DAGs from observational data: constraint-based and score-based. In constraint-based methods, such as the well-known PC-Algorithm [39], the goal is to learn a completed partially DAG (CPDAG) consistent with conditional independence relations inferred from the data. Score-based methods, including the approach in this paper, aim to find a DAG that maximizes a score that measures how well the DAG fits the data.

Existing DAG learning algorithms can also be divided into methods for discrete and continuous data. Score-based methods for learning DAGs from discrete data typically involve a two-stage learning process. In stage 1, the score for each candidate parent set (CPS) of each node is computed. In stage 2, a search algorithm is used to maximize the global score, so that the resulting graph is acyclic. Both of these stages require exponential computation time [42]. For stage 2, there exist elegant exact algorithms based on dynamic programming [13, 21, 22, 30, 32, 36, 44], A⋆ algorithm [43, 44], and integer-programming [2, 3, 7, 8, 9, 10]. The A⋆ algorithm identifies the optimal DAG by solving a shortest path problem in an implicit state-space search graph, whereas integer programming (IP) directly casts the problem as a constrained optimization problem. Specifically, the variables in the IP model indicate whether or not a given parent set is assigned to a node in the network. Hence, the problem involves $m2^{m-1}$ binary variables for $m$ nodes. To reduce the number of binary variables, a common practice is to limit the cardinality of each parent set [3, 9], which can lead to suboptimal solutions.

A comprehensive empirical evaluation of A*⋆* algorithm and IP methods for discrete data is conducted in [26]. The results show that the relative efficiency of these methods varies due to the intrinsic differences between them. In particular, state-of-the-art IP methods can solve instances up to 1,000 CPS per variable regardless of the number of nodes, whereas A*⋆* algorithm works for problems with up to 30 nodes, even with tens of thousands of CPS per node.

While statistical properties of DAG learning from continuous data have been extensively studied [15, 25, 34, 37, 41], the development of efficient computational tools for learning the optimal DAG from continuous data remains an open challenge. In addition, despite the existence of elegant exact search algorithms for discrete data, the literature on DAG learning from continuous data has primarily focused on approximate algorithms based on coordinate descent [14, 17] and non-convex continuous optimization [46]. To our knowledge, [29] and [42] provide the only exact algorithms for learning medium to large DAGs from continuous data. An IP-based model using the topological ordering of variables is proposed in [29]. An A*⋆-lasso algorithm for learning an optimal DAG from continuous data with an $\ell_{1}$ regularization is developed in [42]. A⋆-lasso incorporates the lasso-based scoring method within dynamic programming to avoid an exponential number of parent sets and uses the A⋆* algorithm to prune the search space of the dynamic programming method.

Given the state of existing algorithms for DAG learning from continuous data, there is currently a gap between theory and computation: While statistical properties of exact algorithms can be rigorously analyzed [25, 41], it is much harder to assess the statistical properties of approximate algorithms [1, 14, 17, 46] that offer no optimality guarantees [21]. This gap becomes particularly noticeable in cases where the statistical model is identifiable from observational data. In this case, the optimal score from exact search algorithms is guaranteed to reveal the true underlying DAG when the sample size is large. Therefore, causal structure learning from observational data becomes feasible [25, 33].

In this paper, we focus on DAG learning for an important class of BNs for continuous data, where causal relations among random variables are linear. More specifically, we consider DAGs corresponding to linear structural equations models (SEMs). In this case, network edges are associated with the coefficients of regression models corresponding to linear SEMs. Consequently, the score function can be explicitly encoded as a penalized negative log-likelihood function with an appropriate choice of regularization [29, 35, 41]. Hence, the process of computing scores (i.e., stage 1) is completely bypassed, and a single-stage model can be formulated [42]. Moreover, in this case, IP formulations only require a polynomial, rather than exponential, number of variables, because each variable can be defined in the space of arcs instead of parent sets. Therefore, cardinality constraints on the size of parent sets, which are used in earlier methods to reduce the search space and may lead to suboptimal solutions, are no longer necessary.

Contributions

In this paper, we develop tailored exact DAG learning methods for continuous data from linear SEMs, and make the following contributions:

• –

We develop a mathematical framework that can naturally incorporate prior structural knowledge, when available. Prior structural knowledge can be supplied in the form of an undirected and possibly cyclic graph (super-structure). An example is the skeleton of the DAG, obtained by removing the direction of all edges in a graph. Another example is the moral graph of the DAG, obtained by adding an edge between pairs of nodes with common children and removing the direction of all edges [39]. The skeleton and moral graphs are particularly important cases, because they can be consistently estimated from observational data under proper assumptions [20, 25]. Such prior information limits the number of possible DAGs and improves the computational performance.

• –

We discuss three mathematical formulations, namely, cutting plane (CP), linear ordering (LO), and topological ordering (TO) formulations, for learning an optimal DAG, using both $\ell_{0}$ and $\ell_{1}$-penalized likelihood scores. We also propose a new mathematical formulation to learn an optimal DAG from continuous data, the layered network (LN) formulation, and establish that other DAG learning formulations entail a smaller continuous relaxation feasible region compared to that of the continuous relaxation of the LN formulation (Propositions 3 and 4). Nonetheless, all formulations attain the same optimal continuous relaxation objective function value under a mild condition (Propositions 5 and 6). Notably, the number of binary variables and constraints in the LN formulation solely depend on the number of edges in the super-structure (e.g., moral graph). Thus, the performance of the LN formulation substantially improves in the presence of a sparse super-structure. The LN formulation has a number of other advantages; it is a compact formulation in contrast to the CP formulation; its relaxation can be solved much more efficiently compared with the LO formulation; and it requires fewer binary variables and explores fewer branch-and-bound nodes than the TO formulation. Our empirical results affirm the computational advantages of the LN formulation. They also demonstrate that the LN formulation can find a graph closer to the true underlying DAG. These improvements become more noticeable in the presence of a prior super-structure (e.g., moral graph).

• –

We compare the IP-based method and the A*⋆-lasso algorithm for the case of $\ell_{1}$ regularization. As noted earlier, there is no clear winner among A⋆-style algorithms and IP-based models for DAG learning from discrete data [26]. Thus, a wide range of approaches based on dynamic programming, A⋆* algorithm, and IP-based models have been proposed for discrete data. In contrast, for DAG learning from continuous data with complete super-structure, the LN formulation remains competitive with the state-of-art A*⋆-lasso algorithm for small graphs, whereas it performs better for larger problems. Moreover, LN performs substantially better when a sparse super-structure is available. This is mainly because the LN formulation directly defines the variables based on the super-structure, whereas the A⋆*-lasso algorithm cannot take advantage of the prior structural knowledge as effectively as the LN formulation.

In Section 2, we outline the necessary preliminaries, define the DAG structure learning problem, and present a general framework for the problem. Mathematical formulations for DAG learning are presented in Section 3. The strengths of different optimization problems are discussed in Section 4. Empirical results are presented in Sections 5 and 6. We end the paper with a summary and discussion of future research in Section 7.

2 Penalized DAG Estimation with Linear SEMs

The causal effect of (continuous) random variables in a DAG $\mathcal{G}_{0}$ can be described by SEMs that represent each variable as a (nonlinear) function of its parents. The general form of these models is given by [31]

[TABLE]

where $X_{k}$ is the random variable associated with node $k$; $pa^{\mathcal{G}_{0}}_{k}$ denotes the parents of node $k$ in $\mathcal{G}_{0}$, i.e., the set of nodes with arcs pointing to node $k$; $m$ is the number of nodes; and latent random variables, $\delta_{k}$ represent the unexplained variation in each node.

An important class of SEMs is defined by linear functions, $f_{k}(\cdot)$, which can be described by $m$ linear regressions of the form

[TABLE]

where $\beta_{jk}$ represents the effect of node $j$ on $k$ for $j\in pa^{\mathcal{G}_{0}}_{k}$. In the special case where the random variables are Gaussian, Equations (1) and (2) are equivalent, in the sense that $\beta_{jk}$ are coefficients of the linear regression model of $X_{k}$ on $X_{j}$, $j\in pa^{\mathcal{G}_{0}}_{k}$, and $\beta_{jk}=0$ for $j\notin pa^{\mathcal{G}_{0}}_{k}$ [35]. However, estimation procedures proposed in this paper are not limited to Gaussian random variables and apply more generally to linear SEMs [25, 35].

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$. From $\mathcal{M}$, generate a bi-directional graph $\overrightarrow{\mathcal{M}}=(V,\overrightarrow{E})$ where $\overrightarrow{E}=\{(j,k),(k,j)|(j,k)\in E\}$. Throughout the paper, we refer to directed edges as arcs and to undirected edges as edges.

Consider $n$ i.i.d. observations from the linear SEM (2). 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, and $m$ columns representing random variables. The linear SEM (2) can be compactly written as

[TABLE]

where ${B}=[\beta]\in{\mathbb{R}}^{m\times m}$ is a matrix with $\beta_{kk}=0$ for $k=1,\dots,m$ and $\beta_{jk}=0$ for all $(j,k)\notin\overrightarrow{E}$; $\Delta$ is the $n\times m$ noise matrix. More generally, $B$ defines a directed graph $\mathcal{G}(B)$ on $m$ nodes such that arc $(j,k)$ appears in $\mathcal{G}(B)$ if and only if $\beta_{jk}\neq 0$.

Let $\sigma^{2}_{k}$ denote the variance of $\delta_{k};k={1,2,,\dots,m}$. We assume that all noise variables have equal variances, i.e., $\sigma_{k}=\sigma$. This condition implies the identifiability of DAGs from Gaussian [33] and non-Gaussian [25] data. Under this condition, the negative log likelihood for linear SEMs with Gaussian or non-Gaussian noise is proportional to

[TABLE]

where $S=\mathcal{X}^{T}\mathcal{X}$ and $I$ is the identity matrix [25].

In practice, we are often interested in learning sparse DAGs. Thus, a regularization term is used to obtain a sparse estimate. For the linear SEM (2), the optimization problem corresponding to the penalized negative log-likelihood with super-structure $\mathcal{M}$ (PNL$\mathcal{M}$) for learning sparse DAGs is given by

[TABLE]

The objective function (5a) consists of two parts: the quadratic loss function, $l_{n}(\beta)$, from (4), and the regularization term, $\phi({B})$. Popular choices for $\phi(B)$ include $\ell_{1}$-regularization or lasso [40], $\phi({B})=\sum_{(j,k)\in E^{\rightarrow}}|\beta_{jk}|$, and $\ell_{0}$-regularization, $\phi({B})=\sum_{(j,k)\in\overrightarrow{E}}\mathbf{1}(\beta_{jk})$, where $\mathbf{1}(\beta_{jk})=1$ if $\beta_{jk}\neq 0$, and 0 otherwise. The tuning parameter $\lambda$ controls the degree of regularization. The constraint (5b) stipulates that the resulting directed subgraph (digraph) has to be an induced DAG from $\overrightarrow{\mathcal{M}}$.

When the super-structure $\mathcal{M}$ is a complete graph, PNL$\mathcal{M}$ reduces to the classical PNL. In this case, the consistency of sparse PNL for DAG learning from Gaussian data with an $\ell_{0}$ penalty follows from an analysis similar to [41]. In particular, we have

[TABLE]

where $\hat{\mathcal{G}}_{n}$ is the estimate of the true structure ${\mathcal{G}^{0}}$. An important advantage of the PNL estimation problem is that its consistency does not require the (strong) faithfulness assumption [33, 41].

The mathematical model (5) incorporates a super-structure $\mathcal{M}$ (e.g., moral graph) into the PNL model. When $\mathcal{M}$ is the moral graph, the consistency of sparse PNL$\mathcal{M}$ follows from the analysis in [25], which studies the consistency of the following two-stage framework for estimating sparse DAGs: (1) infer the moral graph from the support of the inverse covariance matrix; and (2) choose the best-scoring induced DAG from the moral graph. The authors investigate conditions for identifiability of the underlying DAG from observational data and establish the consistency of the two-stage framework.

While PNL and PNL$\mathcal{M}$ enjoy desirable statistical properties for linear SEMs with Gaussian [41] and non-Gaussian [25] noise, the computational challenges associated with these problems have not been fully addressed. This paper aims to bridge this gap.

3 Mathematical Formulations

Prior to presenting mathematical formulations for solving PNL$\mathcal{M}$, we discuss a property of DAG learning from continuous data that distinguishes it from the corresponding problem for discrete data. To present this property in Proposition 1, we need a new definition.

Definition 1**.**

A tournament is a directed graph obtained by specifying a direction for each edge in the super-structure graph $\mathcal{M}$ (see Figure 1).

Proposition 1**.**

There exists an optimal solution $\mathcal{G}(B)$ to PNL$\mathcal{M}$ (5) with an $\ell_{0}$ (or an $\ell_{1}$) regularization that is a cycle-free tournament.

All proofs are given in Appendix I. Proposition 1 implies that for DAG learning from continuous variables, the search space reduces to acyclic tournament structures. This is a far smaller search space when compared with the super-exponential O\Big{(}m!2^{m\choose 2}\Big{)} search space of DAGs for discrete variables. However, one has to also identify the optimal $\beta$ parameters. This search for optimal $\beta$ parameters is critical, as it further reduces the super-structure to the edges of the DAG by removing the edges with zero $\beta$ coefficients.

A solution method based on brute force enumeration of all tournaments requires $m!\gamma$ computational time when $\mathcal{M}$ is complete, where $\gamma$ denotes the computational time associated with solving PNL$\mathcal{M}$ given a known tournament structure. This is because when $\mathcal{M}$ is complete, the total number of tournaments (equivalently the total number of permutations) is $m!$. However, when $\mathcal{M}$ is incomplete, the number of DAGs is fewer than $m!$. The topological search space is $m!$ regardless of the structure of $\mathcal{M}$ and several topological orderings can correspond to the same DAG. The TO formulation [29] is based on this search space. In Section 3.2, we discuss a search space based on the layering of a DAG, which uniquely identifies a DAG, and propose the corresponding Layered Network (LN) formulation, which effectively utilizes the structure of $\mathcal{M}$. We first discuss existing mathematical formulations for the PNL$\mathcal{M}$ optimization problem (5) in the next section. Given the desirable statistical properties of $\ell_{0}$ regularization [41] and the fact that existing mathematical formulations are given for $\ell_{1}$ regularization, we present the formulations for $\ell_{0}$ regularization. We outline the necessary notation below.

**Indices

$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}$

$\lambda:$ tuning parameter (penalty coefficient)

Continuous optimization variables

$\beta_{jk}$: weight of arc $(j,k)$ representing the regression coefficients $\forall(j,k)\in\overrightarrow{E}$

Binary optimization variables

$z_{jk}=1\,\,\text{if arc}\,\,(j,k)\,\text{exists in a DAG};\text{otherwise}\,0,\,\forall(j,k)\in\overrightarrow{E}$

$g_{jk}=1\,\,\text{if}\,\,\beta_{jk}\neq 0;\,\text{otherwise}\,0,\,\forall(j,k)\in\overrightarrow{E}$

3.1 Existing Mathematical Models

The main source of difficulty in solving PNL$\mathcal{M}$ is due to the acyclic nature of DAG imposed by the constraint in (5b). A popular technique for ensuring acyclicity is to use cycle elimination constraints, which were first introduced in the context of the Traveling Salesman Problem (TSP) in [11].

Let $\mathcal{C}$ be the set of all possible cycles and $\mathcal{C}_{A}\in\mathcal{C}$ be the set of arcs defining a cycle and define F(\beta,g):=n^{-1}\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\sum_{(j,k)\in\overrightarrow{E}}g_{jk}. Then, the $\ell_{0}$-PNL$\mathcal{M}$ model can be formulated as

[TABLE]

Following [29], the objective function (6a) is an expanded version of $l_{n}(\beta)$ in PNL$\mathcal{M}$ (multiplied by $2n^{-1}$) with an $\ell_{0}$ regularization. The constraints in (6b) stipulate that $\beta_{jk}\neq 0$ only if $z_{jk}=1$, where $M$ is a sufficiently large constant. The constraints in (6c) rule out all cycles. Note that for $|\mathcal{C}_{A}|=2$, constraints in (6c) ensure that at most one arc exists among two nodes. The last set of constraints specifies the binary nature of the decision vector $g$. Note that $\beta$ variables are continuous and unrestricted; however, in typical applications, they can be bounded by a finite number $M$. This formulation requires $|\overrightarrow{E}|$ binary variables and an exponential number of constraints. A cutting plane method [28] that adds the cycle elimination inequalities as needed is often used to solve this problem. We refer to this formulation as the cutting plane (CP) formulation.

Remark 1*.*

For a complete super-structure $\mathcal{M}$, it suffices to impose the set of constraints in (6c) only for cycles of size 2 and 3 given by

[TABLE]

In other words, the CP formulation (both with $\ell_{0}$ and $\ell_{1}$ regularizations) needs a polynomial number of constraints for complete super-structure $\mathcal{M}$.

The second formulation is based on a well-known combinatorial optimization problem, known as linear ordering (LO) [16]. Given a finite set $S$ with $q$ elements, a linear ordering of $S$ is a permutation $\mathcal{P}\in S_{q}$ where $S_{q}$ denotes the set of all permutations with $q$ elements. In the LO problem, the goal is to identify the best permutation among $m$ nodes. The “cost” for a permutation $\mathcal{P}$ depends on the order of the elements in a pairwise fashion. Let $p_{j}$ denote the order of node $j\in V$ in permutation $\mathcal{P}$. Then, for two nodes $j,k\in\{1,\dots,m\}$, the cost is $c_{jk}$ if the order of node $j$ precedes the order of node $k$ ($p_{j}\prec p_{k}$) and is $c_{kj}$ otherwise ($p_{j}\succ p_{k}$). A binary variable $w_{jk}$ indicates whether $p_{j}\prec p_{k}$. Because $w_{jk}+w_{kj}=1$ and $w_{jj}=0$, one only needs ${m\choose 2}$ variables to cast the LO problem as an IP formulation [16].

The LO formulation for DAG learning from continuous data has two noticeable differences compared with the classical LO problem: (i) the objective function is quadratic, and (ii) an additional set of continuous variables, i.e., $\beta$s, is added. Cycles are ruled out by directly imposing the linear ordering constraints. The PNL$\mathcal{M}$ can be formulated as (8).

[TABLE]

The interpretation of constraints (8b)-(8d) is straightforward. The constraints in (8c) imply that if node $j$ appears after node $k$ in a linear ordering ($w_{jk}=0$), then there should not exist an arc from $j$ to $k$ ($g_{jk}=0$). The set of inequalities (8e) implies that if $p_{i}\prec p_{j}$ and $p_{j}\prec p_{k}$, then $p_{i}\prec p_{k}$. This ensures the linear ordering of nodes and removes cycles.

The third approach for ruling out cycles is to impose a set of constraints such that the nodes follow a topological ordering. A topological ordering is a linear ordering of the nodes of a graph such that the graph contains an arc $(j,k)$ if node $j$ appears before node $k$ in the linear order. Define decision variables $o_{rs}\in\{0,1\}$ for all $r,s\in\{1,\dots,m\}$. This variable takes value 1 if topological order of node $r$ (i.e., $p_{r}$) equals $s$, and 0, otherwise. If a topological ordering is known, the DAG structure can be efficiently learned in polynomial time [35], but the problem remains challenging when the ordering is not known. The topological ordering prevents cycles in the graph. This property is used in [29] to model the problem of learning a DAG with $\ell_{1}$ regularization. We extend their formulation to $\ell_{0}$ regularization. The topological ordering (TO) formulation is given by

[TABLE]

In this formulation, $z_{jk}$ is an auxiliary binary variable which takes value 1 if an arc exists from node $j$ to node $k$. Recall that, $g_{jk}=1$ if $|\beta_{jk}|>0$. The constraints in (9c) enforce the correct link between $g_{jk}$ and $z_{jk}$, i.e., $g_{jk}$ has to take value zero if $z_{jk}=0$. The constraints in (9d) imply that there should not exist a bi-directional arc among two nodes. This inequality can be replaced with equality (Corollary 1). The constraints in (9e) remove cycles by imposing an ordering among nodes. The set of constraints in (9f)-(9g) assigns a unique topological order to each node. The last two sets of constraints indicate the binary nature of decision variables $o$ and $z$.

Corollary 1, which is a direct consequence of Proposition 1, implies that we can use $z_{jk}=1-z_{kj}\,\,$ for all $j<k$ and reduce the number of binary variables.

Corollary 1**.**

The constraints in (9d) can be replaced by $z_{jk}+z_{kj}=1,\,\,\forall(j,k)\in\overrightarrow{E}$.

In constraints (9b)-(9i), both variables $z$ and $g$ are needed to correctly model the $\ell_{0}$ regularization term in the objective function (see Figure 2(a)). This is because the constraints (9e) satisfy the transitivity property: if $z_{ij}=1$ for $(i,j)\in\overrightarrow{E}$ and $z_{jk}=1$ for $(j,k)\in\overrightarrow{E}$, then $z_{ik}=1$ for $(i,k)\in\overrightarrow{E}$, since $z_{ij}=1$ implies that $\sum_{s\in V}s(o_{js}-o_{is})\geq 1$; similarly, $z_{jk}=1$ implies $\sum_{s\in V}s\,(o_{ks}-o_{js})\geq 1$. If we sum both inequalities, we have $\sum_{s\in V}s\,(o_{ks}-o_{is})\geq 2$, which enforces $z_{ik}=1$. Such a transitivity relation, however, need not hold for the decision vector $g$. In other words, the decision variable $g_{jk}$ is used to keep track of the number of non-zero weights $\beta_{jk}$ associated with the arc $(j,k)$, and the decision vector $z$ is used to remove cycles via the set of constraints in (9e) by creating an acyclic tournament on the super-structure $\mathcal{M}$. A tournament on super-structure $\mathcal{M}$ assigns a direction for each edge in an undirected super-structure $\mathcal{M}$. In other words, if we were to let $g_{ij}=z_{ij}$ for $(i,j)\in\overrightarrow{E}$ and hence use the decision variable $z$ in the objective, then we would be counting the number of edges, equal to $|E|$, instead of number of non-zero $\beta$ values.

3.2 A New Mathematical Model: The Layered Network (LN) Formulation

As an alternative to the existing mathematical formulations, we propose a new formulation for imposing acyclicity constraints that is motivated by the layering of nodes in DAGs [18]. More specifically, our formulation ensures that the resulting graph is a layered network, in the sense that there exists no arc from a layer $v$ to layer $u$, where $u<v$. Let $\psi_{k}$ be the layer value for node $k$. One may interpret $\psi_{k}$ as $\sum_{s=1}^{m}{s\,o_{ks}}$ for all $k\in V$, where the variables $o_{ks}$ are as defined in the TO formulation. However, note that the notion of $\psi_{k}$ is more general because $\psi_{k}$ need not be integer. Figure 3 depicts the layered network encoding of a DAG. With this notation, our layered network (LN) formulation can be written as

[TABLE]

The interpretation of the constraints (10b)-(10c) is straightforward. The constraints in (10e) ensure that the graph is a layered network. The last set of constraints indicates the continuous nature of the decision variable $\psi$ and gives the tightest valid bound for $\psi$. It suffices to consider any real number for layer values $\psi$ as long as layer values of any two nodes differ by at least one if there exists an arc between them. Additionally, LN uses a tighter inequality compared to TO, by replacing $m$ with parameter $m-1$ in (10e). This is because the difference between the layer values of two nodes can be at most $m-1$ for a DAG with $m$ nodes. The next proposition establishes the validity of the LN formulation.

Proposition 2**.**

An optimal solution to (10) is an optimal solution to (5).

The LN formulation highlights a desirable property of the layered network representation of a DAG in comparison to the topological ordering representation. Let us define nodes in layer 1 as the set of nodes that have no incoming arcs in the DAG, nodes in layer 2 as the set of nodes that have incoming arcs only from nodes in the layer 1, layer 3 as the set of nodes that have incoming arcs from layer 2 (and possibly layer 1), etc. (see Figure 3). The minimal layer number of a node is the length of the longest directed path from any node in layer 1 to that node. For a given DAG, there is a unique minimal layer number, but not a unique topological order. As an example, Figure 2(a) has three valid topological orders: (i) 1,2,4,3, (ii) 1,4,2,3, and (iii) 4,1,2,3. In contrast, it has a unique layer representation, 1,2,3,1.

There is a one-to-one correspondence between minimal layer numbering and a DAG. However, the solutions of the LN formulation, i.e., $\psi$ variables (layer values), do not necessarily correspond to the minimal layer numbering. This is because the LN formulation does not impose additional constraints to enforce a minimal numbering and can output solutions that are not minimally numbered. However, because branch-and-bound does not branch on continuous variables, alternative (non-minimal) feasible solutions for the $\psi$ variables do not impact the branch-and-bound process. On the contrary, we have multiple possible representations of the same DAG with topological ordering. Because topological ordering variables are binary, the branch-and-bound method applied to the TO formulation explores multiple identical DAGs as it branches on the topological ordering variables. This enlarges the size of the branch-and-bound tree and increases the computational burden.

Layered network representation also has an important practical implication: Using this representation, the search space can be reduced to the total number of ways we can layer a network (or equivalently the total number of possible minimal layer numberings) instead of the total number of topological orderings. When the super-structure $\mathcal{M}$ is complete, both quantities are the same, and equal to $m!$. Otherwise, a brute-force search for finding the optimal DAG has computational time $\mathcal{L}{\gamma}$, where $\mathcal{L}$ denotes the total number of minimal layered numberings, and $\gamma$ is the computational complexity of solving PNL$\mathcal{M}$ given a known tournament structure.

We close this section by noting that a set of constraints similar to (10e) was introduced in [7] for learning pedigree. However, the formulation in [7] requires an exponential number of variables. In more recent work on DAGs for discrete data, Cussens and colleagues have focused on a tighter formulation for removing cycles, known as cluster constraints [8, 9, 10]. To represent the set of cluster constraints, variables have to be defined according to the parent set choice leading to an exponential number of variables. Thus, such a representation is not available in the space of arcs.

3.2.1 Layered Network with $\ell_{1}$ regularization

Because of its convexity, the structure learning literature has utilized the $\ell_{1}$-regularization for learning DAGs from continuous variables [17, 29, 35, 42, 46]. The LN formulation with $\ell_{1}$-regularization can be written as

[TABLE]

Remark 2*.*

For a complete super-structure $\mathcal{M}$, $\psi_{k}=\sum_{j\in V\setminus k}z_{jk}\,\forall k\in V$. Thus, the LN formulations (both $\ell_{0}$ and $\ell_{1}$) can be encoded without $\psi$ variables by writing (10e) as

[TABLE]

Remark 3*.*

CP and LO formulations reduce to the same formulation for $\ell_{1}$ regularization when the super-structure $\mathcal{M}$ is complete by letting $w_{ij}=g_{ij}$ in formulation (8) for all $(j,k)\in\overrightarrow{E}$.

An advantage of the $\ell_{1}$-regularization for DAG learning is that all models (CP, LO, TO and LN) can be formulated without decision variables $g_{jk}$, since counting the number of non-zero $\beta_{jk}$ is no longer necessary.

Table 1 shows the number of binary variables and the number of constraints associated with cycle prevention constraints in each model. Evidently, $\ell_{0}$ models require additional binary variables compared to the corresponding $\ell_{1}$ models. Note that the number of binary variables and constraints for the LN formulation solely depend on the number of edges in the super-structure $\mathcal{M}$. This property is particularly desirable when the super-structure $\mathcal{M}$ is sparse. The LN formulation requires the fewest number of constraints among all models. The LN formulation also requires fewer binary variables than the TO formulation. More importantly, different topological orders for the same DAG are symmetric solutions to the associated TO formulation. Consequently, branch-and-bound requires exploring multiple symmetric formulations as it branches on fractional TO variables. As for the LO formulation, the number of constraints is $O(m^{3})$ which makes its continuous relaxation cumbersome to solve in the branch-and-bound process. The LN formulation is compact, whereas the CP formulation requires an exponential number of constraints for incomplete super-structure $\mathcal{M}$. The CP formulation requires fewer binary variables for $\ell_{0}$ formulation than LN; both formulations need the least number of binary variables for $\ell_{1}$ regularization.

In the next section, we discuss the theoretical strength of these mathematical formulations and provide a key insight on why the LN formulation performs well for learning DAGs from continuous data.

4 Continuous Relaxation

One of the fundamental concepts in IP is relaxations, wherein some or all constraints of a problem are loosened. Relaxations are often used to obtain a sequence of easier problems which can be solved efficiently yielding bounds and approximate, not necessarily feasible, solutions for the original problem. Continuous relaxation is a common relaxation obtained by relaxing the binary variables of the original mixed-integer quadratic program (MIQP) and allowing them to take real values. Continuous relaxation is at the heart of branch-and-bound methods for solving MIQPs. An important concept when comparing different MIQP formulations is the strength of their continuous relaxations.

Definition 2**.**

A formulation $A$ is said to be stronger than formulation $B$ if $\mathcal{R}(A)\subset\mathcal{R}(B)$ where $\mathcal{R}(A)$ and $\mathcal{R}(B)$ correspond to the feasible regions of continuous relaxations of $A$ and $B$, respectively.

Proposition 3**.**

*The LO formulation is stronger than the LN formulation, that is,

$\mathcal{R}(LO)\subset\mathcal{R}(LN)$.*

Proposition 4**.**

When the parameter $m$ in (9e) is replaced with $m-1$, the TO formulation is stronger than the LN formulation, that is, $\mathcal{R}(TO)\subset\mathcal{R}(LN)$.

These propositions are somewhat expected because the LN formulation uses the fewest number of constraints. Hence, the continuous relaxation feasible region of the LN formulation is loosened compared to the other formulations. The next two results justify the advantages of the LN formulation.

Proposition 5**.**

Let $\beta^{\star}_{jk}$ denote the optimal coefficient associated with an arc $(j,k)\in\overrightarrow{E}$ from (5). For both $\ell_{0}$ and $\ell_{1}$ regularizations, the initial continuous relaxations of the LN formulation attain as tight an optimal objective function value as the LO, CP, TO formulations if $M\geq 2\underset{(j,k)\in\overrightarrow{E}}{\max}\,|\beta^{\star}_{jk}|$.

Proposition 5 states that although the LO and TO formulations are tighter than the LN formulation with respect to the feasible region of their continuous relaxations, the continuous relaxation of all models attain the same objective function value (root relaxation).

Proposition 6**.**

For the same variable branching in the branch-and-bound process, the continuous relaxations of the LN formulation for both $\ell_{0}$ and $\ell_{1}$ regularizations attain as tight an optimal objective function value as LO, CP and TO, if $M\geq 2\underset{(j,k)\in\overrightarrow{E}}{\max}\,|\beta^{\star}_{jk}|$.

Proposition 6 is at the crux of this section. It shows that not only does the tightness of the optimal objective function value of the continuous relaxation hold for the root relaxation, but it also holds throughout the branch-and-bound process under the specified condition on $M$, if the same branching choices are made. Thus, the advantages of the LN formulation are due to the fact that it is a compact formulation that entails the fewest number of constraints, while attaining the same optimal objective value of continuous relaxation as tighter models.

In practice, finding a tight value for $M$ is difficult. Our computational results show that the approach suggested in [29] to obtain a value of $M$, which is explained in Section 5 and used in our computational experiments, always satisfies the condition in Proposition 6 across all generated instances.

5 Comparison of MIQP Formulations

We present numerical results comparing the proposed LN formulation with existing approaches. Experiments are performed on a cluster operating on UNIX with Intel Xeon E5-2640v4 2.4GHz. All MIQP formulations are implemented in the Python programming language. Gurobi 8.0 is used as the MIQP solver. A time limit of $50m$ (in seconds), where $m$ denotes the number of nodes, is imposed across all experiments after which runs are aborted. Unless otherwise stated, an MIQP optimality gap of $0.001$ is imposed across all experiments; the gap is calculated by $\frac{UB-LB}{UB}$ where UB denotes the objective value associated with the best feasible integer solution (incumbent) and LB represents the best obtained lower bound during the branch-and-bound process.

For CP, instead of incorporating all constraints given by (6c), we begin with no constraint of type (6c). Given an integer solution with cycles, we detect a cycle and impose a new cycle prevention constraint to remove the detected cycle. Depth First Search (DFS) can detect a cycle in a directed graph with complexity $O(|V|+|E|)$. Gurobi Lazy Callback is used, which allows adding cycle prevention constraints in the branch-and-bound algorithm, whenever an integer solution with cycles is found. The same approach is used by [29]. Note that Gurobi solver follows a branch-and-cut implementation and adds many general-purpose and special-purpose cutting planes.

To select the $M$ parameter in all formulations we use the proposal of [29]. Specifically, given $\lambda$, we solve each problem without cycle prevention constraints. We then use the upper bound $M=2\underset{(j,k)\in\overrightarrow{E}}{\max}\,|\beta_{jk}|$. The results provided in [29] computationally confirm that this approach gives a large enough value of $M$. We also confirmed the validity of this choice across all our test instances.

5.1 Synthetic datasets

We use the R package pcalg to generate random Erdős-Rényi graphs. Firstly, we create a DAG using randomDAG function and assign random arc weights (i.e., $\beta$) from a uniform distribution, $\mathcal{U}[0.1,1]$. This ground truth DAG is used to assess the quality of estimates. Next, the resulting DAG and random coefficients are input to the rmvDAG function, which uses linear regression as the underlying model, to generate multivariate data (columns of matrix $\mathcal{X}$) with the standard normal error distribution.

We consider $m\in\{10,20,30,40\}$ nodes and $n\in\{100,1000\}$ samples. The average outgoing degree of each node, denoted by $d$, is set to 2. We generate 10 random graphs for each setting ($m$, $n$, $d$). The raw observational data, $\mathcal{X}$, for the datasets with $n=100$ is the same as first 100 rows of the datasets with $n=1000$.

We consider two types of problem instances: (i) a set of instances for which the moral graph corresponding to the true DAG is available; (ii) a set of instances with a complete undirected graph, i.e., assuming no prior knowledge. The first class of problems is referred to as moral instances, whereas the second class is called complete instances. The raw observational data, $\mathcal{X}$, for moral and complete instances are the same. The function moralize(graph) in the pcalg R-package is used to generated the moral graph from the true DAG. The moral graph can also be (consistently) estimated from data using penalized estimation procedures with polynomial complexity [20, 25]. However, since the quality of the moral graph equally affects all optimization models, the true moral graph is used in our experiments.

We use the following IP-based metrics to measure the quality of a solution: Optimality gap (MIQP GAP), computation time in seconds (Time), Upper Bound (UB), Lower Bound (LB), computational time of root continuous relaxation (Time LP), and the number of explored nodes in the branch-and-bound tree.

We also evaluate the quality of the estimated DAGs by comparing them with the ground truth. To this end, we use the average structural Hamming distance $(\mathrm{SHD})$, as well as true positive (TPR) and false positive rates (FPR). These criteria evaluate different aspects of the quality of the estimated DAGs: $\mathrm{SHD}$ counts the number of differences (addition, deletion, or arc reversal) required to transform predicted DAG to the true DAG; TPR is the number of correctly identified arcs divided by the total number of true arcs, $P$; FPR is the number of incorrectly identified arcs divided by the total number of negatives (non-existing arcs), $N$. For brevity, TPR and FPR plots are presented in Appendix II.

5.2 Comparison of $\ell_{0}$ formulations

Figure 5 reports the average metrics across 10 random graphs for $\ell_{0}$ formulations with $n=1000$. The LO formulation fails to attain a reasonable solution for one graph (out of 10) with $m=40$ and $\lambda\in\{0.1,1\}$. This is due to the large computation time for solving its continuous relaxation. We excluded these two instances from LO results.

Figure 5(a) shows that the LN formulation outperforms other formulations in terms of the average optimality gap across all number of nodes $m\in\{10,20,30,40\}$ and regularization parameters, $\lambda\in\{0.1,1\}$. The difference becomes more pronounced for moral instances. For moral instances, the number of binary variables and constraints for LN solely depends on the size of moral graph. Figure 5(b) also indicates that the LN formulation requires the least computational time for small instances, whereas all models hit the time limit for larger instances.

Figures 5(c)-(d) show the performance of all methods in terms of their upper and lower bounds. For easier instances (e.g., complete instances with $m\in\{10,20\}$ and moral instances), all methods attain almost the same upper bound. Nonetheless, LN performs better in terms of improving the lower bound. For more difficult instances, LN outperforms other methods in terms of attaining a smaller upper bound (feasible solution) and a larger lower bound.

Figures 5(e)-(f) show the continuous relaxation time of all models, and the number of explored nodes in the branch-and-bound tree, respectively. The fastest computational time for the continuous relaxation is for the TO formulation followed by the LN formulation. However, the number of explored nodes provides more information about the performance of mathematical formulations. In small instances, i.e., $m=10$, where an optimal solution is attained, the size of the branch-and-bound tree for the LN formulation is smaller than the TO formulation. This is because the TO formulation has a larger number of binary variables, leading to a larger branch-and-bound tree. On the other hand, for large instances, the number of explored nodes in the LN formulation is larger than the TO formulation. This implies that the LN formulation explores more nodes in the branch-and-bound tree given a time limit. This may be because continuous relaxations of the LN formulation are easier to solve in comparison to the continuous relaxations of the TO formulation in the branch-and-bound process. As stated earlier, the branch-and-bound algorithm needs to explore multiple symmetric formulations in the TO formulation as it branches on fractional topological ordering variables. This degrades the performance of the TO formulation. The LO formulation is very slow because its continuous relaxation becomes cumbersome as the number of nodes, $m$, increases. Thus, we can see a substantial decrease in the number of explored nodes in branch-and-bound trees associated with the LO formulation. The CP formulation is implemented in a cutting-plane fashion. Hence, its number of explored nodes is not directly comparable with other formulations.

Figures 5(a)-(f) show the importance of incorporating available structural knowledge (e.g., moral graph). The average optimality gap and computational time are substantially lower for moral instances compared to complete instances. Moreover, the substantial difference in the optimality gap elucidates the importance of incorporating structural knowledge. Similar results are obtained for $n=100$ samples; see Appendix II.

We next discuss the performance of different methods in terms of estimating the true DAG. The choice of tuning parameter $\lambda$, the number of samples $n$, and the quality of the best feasible solution (i.e., upper bound) influence the resulting DAG. Because our focus in this paper is on computational aspects, we fixed the values of $\lambda$ for a fair comparison between the formulations, and used $\lambda=0.1$ based on results in preliminary experiments. Thus, we focus on the impact of sample size as well as the quality of the feasible solution in the explanation of our results.

Figures 6(a)-(b) show the SHDs for all formulations for $n=1000$ and $n=100$, respectively. Comparing Figure 6(a) with Figure 6(b), we observe that the SHD tends to increase as the number of samples decreases. As discussed earlier, when $n\rightarrow\infty$, penalized likelihood likelihood estimate with an $\ell_{0}$ regularization ensures identifiability in our setting [33, 41]. However, for a finite sample size, identifiability may not be guaranteed. Moreover, the appropriate choice of $\lambda$ for $n=100$ may be different than the corresponding $\lambda$ for $n=1000$.

Figure 6(a) shows that all methods learn the true DAG with $\lambda=0.1$, and given a moral graph for $m\in\{10,20,30\}$. In addition, SHD is negligible for LN and CP formulations for $m=40$. However, we observe a substantial increase in SHD (e.g., from 0.2 to near 10 for LN) for complete graphs. These figures indicate the importance of incorporating available structural knowledge (e.g., a moral graph) for better estimation of the true DAG.

While, in general, LN performs well compared with other formulations, we do not expect to see a clear dominance in terms of accuracy of DAG estimation either due to finite samples or the fact that none of the methods could attain a global optimal solution for larger instances. On the contrary, we retrieve the true DAG for smaller graphs for which optimal solutions are obtained. As pointed out in [29], a slight change in the objective function value could significantly alter the estimated DAG. Our results corroborate this observation.

5.3 Comparison of $\ell_{1}$ formulations

Figure 7 shows various average metrics across 10 random graphs for $\ell_{1}$ regularization with $n=1000$ samples. Figure 7(a) shows that the LN formulation clearly outperforms other formulations in terms of average optimality gap across all number of nodes, $m\in\{10,20,30,40\}$, and regularization parameters, $\lambda\in\{0.1,1\}$. Moreover, Figure 7(b) shows that the LN formulation requires significantly less computational time in moral instances, and in complete instances with $m\in\{10,20\}$ compared to other methods. In complete instances, all methods hit the time limit for $m\in\{30,40\}$. Figures 7(c)-(f) can be interpreted similar to the Figures 5(c)-(f) for $\ell_{0}$ regularization. Similar to $\ell_{0}$ regularization, Figures 7(a)-(b) demonstrate the importance of incorporating structural knowledge (e.g., a moral graph) for $\ell_{1}$ regularization. Similar results are observed for $n=100$ samples; see Appendix II.

As expected, the DAG estimation accuracy with $\ell_{1}$ regularization is inferior to the $\ell_{0}$ regularization. This is in part due to the bias associated with the $\ell_{1}$ regularization, which could be further controlled with, for example, adaptive $\ell_{1}$-norm regularization [47]. Nonetheless, formulations for $\ell_{1}$ regularization require less computational time and are easier to solve than the corresponding formulations for $\ell_{0}$ regularization.

6 Comparison with the A*⋆*-lasso algorithm

In this section, we compare the LN formulation with A*⋆*-lasso [42], using the MATLAB code made available by the authors. For this comparison, the same true DAG structures are taken from [42] and the strength of arcs ($\beta$) are chosen from $\mathcal{U}[-1,-0.1]\cup\mathcal{U}[0.1,1]$. The number of nodes in the 10 true DAGs varies from $m=6$ to $m=27$ (see Table 2). The true DAG and resulting random $\beta$ coefficients are used to generate $n=500$ samples for each column of data matrix $\mathcal{X}$.

A time limit of six hours is imposed across all experiments after which runs are aborted. In addition, for a fairer comparison with A*⋆*-lasso, we do not impose an MIQP gap termination criterion of 0.001 for LN and use the Gurobi default optimality gap criterion of 0.0001.

Similar to synthetic data described in Section 5.1, we consider two cases: (i) moral instances and (ii) complete instances. For the former case, the moral graph is constructed from the true DAG as done in Section 5.1. The raw observational data (i.e., $\mathcal{X}$) for moral and complete instances are the same.

We compare the LN formulation with A*⋆-lasso using $\ell_{1}$ regularization. Note that A⋆-lasso cannot solve the model with $\ell_{0}$ regularization. Furthermore, the original A⋆-lasso algorithm assumes no super-structure. Therefore, to enhance its performance, we modified the MATLAB code for A⋆*-lasso in order to incorporate the moral graph structure, when available. We consider two values for $\lambda\in\{0,0.1\}$ for our comparison. As $\lambda$ decreases, identifying an optimal DAG becomes more difficult. Thus, it is of interest to evaluate the computational performance of both methods for $\lambda=0$ (i.e., no regularization) to assess the performance of these approaches on difficult cases (see, e.g., the computational results in [46] and the statistical analysis in [25] for $\lambda=0$). We note that model selection methods (such as Bayesian Information Criterion [38]) can be used to identify the best value of $\lambda$. However, in this section, our focus is to evaluate the computational performance of these approaches for a given $\lambda$ value.

Table 2 shows the solution times (in seconds) of A*⋆-lasso versus the LN formulation for complete and moral instances. For the LN formulation, if the algorithm cannot prove optimality within the 6-hour time limit, we stop the algorithm and report, in parentheses, the optimality gap at termination. For complete instances with $\lambda=0$, the results highlight that for small instances (up to 14 nodes) A⋆-algorithm performs better, whereas the LN formulation outperforms A⋆-lasso for larger instances. In particular, we see that the LN formulation attains the optimal solution for the Cloud data set in 810.47 seconds and it obtains a feasible solution that is provably within 99.5% of the optimal objective value for Funnel and Galaxy data sets. For moral instances, we observe significant improvement in the computational performance of the LN formulation, whereas the improvement in A⋆-lasso is marginal in comparison. This observation highlights the fact that dynamic programming-based approaches cannot effectively utilize the super-structure knowledge, whereas an IP-based approach, particularly the LN formulation, can significantly reduce the computational times. For instance, LN’s computational time for the Cloud data reduces from $\sim 810$ seconds to less than two seconds when the moral graph is provided. In contrast, the reduction in computational time for A⋆*-algorithm with the moral graph is negligible. For $\lambda=0.1$, the problem is easier to solve. Nevertheless, LN performs well in comparison to A⋆-lasso and the performance of the LN formulation improves for moral instances.

For DAG learning from discrete data, an IP-based model (see e.g., [19]) outperforms A⋆ algorithms when a cardinality constraint on the number of the parent set for each node is imposed; A⋆ tends to perform better if such constraints are not enforced. This is mainly because an IP-based model for discrete data requires an exponential number of variables which becomes cumbersome if such cardinality constraints are not permitted. In contrast, for DAG learning from continuous data with linear SEMs, our results show that an IP-based approach does not have such a limitation because variables are encoded in the space of arcs (instead of parent sets). That is why LN performs well even for complete instances (i.e., no restriction on the cardinality of parent set).

There are several fundamental advantages of IP-based modeling, particularly the LN formulation, compared to A*⋆-lasso: (i) The variables in IP-based models (i.e., TO, LN, and CP) depend on the super-structure. Therefore, these IP-based models can effectively utilize the prior knowledge to reduce the search space, whereas A⋆-lasso cannot utilize the super-structure information as effectively. This is particularly important for the LN formulation, as the number of variables and constraints depend only on the super-structure; (ii) all IP-based models can incorporate both $\ell_{0}$ and $\ell_{1}$ regularizations, whereas A⋆-lasso can solve the problem only with $\ell_{1}$ regularization; (iii) all IP-based methods in general enjoy the versatility to incorporate a wide variety of structural constraints, whereas A⋆-lasso and dynamic programming approaches cannot accommodate many structural assumptions. For instance, a modeler may prefer restricting the number of arcs in the DAG, this is achievable by imposing a constraint on an IP-based model, whereas one cannot impose such structural knowledge on A⋆-lasso algorithm; (iv) A⋆-lasso is based on dynamic programming; therefore, one cannot abrupt the search with the aim of achieving a feasible solution. On the other hand, one can impose a time limit or an optimality gap tolerance to stop the search process in a branch-and-bound tree. The output is then a feasible solution to the problem, which provides an upper bound as well as a lower bound which guarantees the quality of the feasible solution; (v) algorithmic advances in integer optimization alone (such as faster continuous relaxation solution, heuristics for better upper bounds, and cutting planes for better lower bounds) have resulted in 29,000 factor speedup in solving IPs [4] using a branch-and-bound process. Many of these advances have been implemented in powerful state-of-the-art optimization solvers (e.g., Gurobi), but they cannot be used in dynamic programming methods, such as A⋆*-lasso.

Figures 9(a)-(b) illustrate the progress of upper bound versus lower bound in the branch-and-bound process for the Factor dataset and highlight an important practical implication of an IP-based model: such models often attain high quality upper bounds (i.e., feasible solutions) in a short amount of time whereas the rest of the time is spent to close the optimality gap by increasing the lower bound.

The results for the moral graph in Figure 9(b) highlight another important observation. That is, providing the moral graph can significantly accelerate the progress of the lower bound in the branch-and-bound process. In other words, information from the moral graph helps close the optimality gap more quickly.

7 Conclusion

In this paper, we study the problem of learning an optimal DAG from continuous observational data using a score function, where the causal effect among the random variables is linear. We cast the problem as a mathematical program and use a penalized negative log-likelihood score function with both $\ell_{0}$ and $\ell_{1}$ regularizations. The mathematical programming framework can naturally incorporate a wide range of structural assumptions. For instance, it can incorporate a super-structure (e.g., skeleton or moral graph) in the form of an undirected and possibly cyclic graph. Such super-structures can be estimated from observational data. We review three mathematical formulations: cutting plane (CP), topological ordering (TO), and Linear Ordering (LO), and propose a new mixed-integer quadratic optimization (MIQO) formulation, referred to as the layered network (LN) formulation. We establish that the continuous relaxations of all models attain the same optimal objective function value under a mild condition. Nonetheless, the LN formulation is a compact formulation in contrast to CP, its relaxation can be solved much more efficiently compared to LO, and enjoys a fewer number of binary variables and traces a fewer number of branch-and-bound nodes than TO.

Our numerical experiments indicate that these advantages result in considerable improvement in the performance of LN compared to other MIQP formulations (CP, LO, and TO). These improvements are particularly pronounced when a sparse super-structure is available, because LN is the only formulation in which the number of constraints and binary variables solely depend on the super-structure. Our numerical experiments also demonstrate that the LN formulation has a number of advantages over the A*⋆*-lasso algorithm, especially when a sparse super-structure is available.

At least two future research avenues are worth exploring. First, one of the difficulties of estimating DAGs using mathematical programming techniques is the constraints in (6b). The big-$M$ constraint is often very loose, which makes the convergence of branch-and-bound process slow. It is of interest to study these constraints in order to improve the lower bounds obtained from continuous relaxations. Second, in many real-world applications, the underlying DAG has special structures. For instance, the true DAG may be a polytree [12]. Another example is a hierarchical structure. In that case, it is natural to learn a DAG such that it satisfies the hierarchy among different groups of random variables. This problem has important applications in discovering the genetic basis of complex diseases (e.g., asthma, diabetes, atherosclerosis).

Appendix I

PROOF OF PROPOSITION 1

Let $({\hat{\mathbf{\beta}}},\hat{z})$ be an optimal solution for (5) with an optimal objective value $F(\hat{\beta})$. Let us refer to the DAG structure corresponding to this optimal solution by $DAG(V,\hat{E}^{\rightarrow})$. Suppose that for some $(j,k)$, we have $\hat{z}_{jk}+\hat{z}_{kj}=0$. To prove the proposition, we construct an optimal solution which satisfies $z_{jk}+z_{kj}=1$ for all pairs of $(j,k)$ and meets the following conditions: (i) this corresponding DAG (tournament) is cycle free (ii) this tournament has the same objective value, i.e., $F(\hat{\beta})$, as an optimal DAG.

Select a pair of nodes, say $p,q\in\mathcal{M},p\neq q$ from $DAG(V,\hat{E})$ for which $\hat{z}_{pq}+\hat{z}_{qp}=0$. If there is a directed path from $p$ to $q$ (respectively $q$ to $p$), then we can add the following arc $(p,q)$ (respectively, $(q,p)$). This arc does not create a cycle in the graph. If there is no directed path between $p$ and $q$, we can add an arc in either direction. In all cases, set $\beta$ value corresponding to the added arc to zero. We repeat this process for all pairs of nodes with $\hat{z}_{pq}+\hat{z}_{qp}=0$.

We can add such arcs without creating any cycle. This is because if we cannot add an arc in either direction, it implies that we should have a directed path from $p$ to $q$ and a directed path from $q$ to $p$ in graph $DAG(V,\hat{E})$ which is a contradiction because it implies a directed cycle in an optimal DAG. Note that in each step, we maintain a DAG. Hence, by induction we conclude that condition (i) is satisfied. The pair of nodes can be chosen arbitrarily.

Since in the constructed solution we set $\beta$ for the added arcs as zero, the objective value does not change. This satisfies condition (ii), and completes the proof. $\square$

PROOF OF PROPOSITION 2

First we prove that (10e) removes all cycles. Suppose, for contradiction, that a cycle of size $p\geq 2$ is available and represented by $(1,2,\dots,p,1)$. This implies $z_{j+1,j}=0\,\text{and}\,z_{j,j+1}=1,\,\forall j=\{1,\dots,p-1\}$, and $z_{p,1}=1,z_{1,p}=0$. Then,

[TABLE]

We sum the above inequalities and conclude $p\leq 0$, a contradiction.

To complete the proof, we also need to prove that any DAG is feasible for the LN formulation. To this end, we know that each DAG has a topological ordering. For all the existing arcs in the DAG, substitute $z_{jk}=1$ and assign a topological ordering number to the variables $\psi_{k},k\in V$ in LN. Then, the set of constraints in (10e) is always satisfied. $\square$

PROOF OF PROPOSITION 3

The LO formulation is in $w$-space whereas LN is in $(z,\psi)$-space. Hence, we first construct a mapping between these decision variables.

Given a feasible solution $w_{jk}$ for all $j,k\in V,\,j\neq k$ in the LO formulation, we define $z_{jk}=w_{jk}$, $(j,k)\in\overrightarrow{E}$ and $\psi_{j}=\sum_{\ell\in V\backslash{\{j\}}}w_{\ell j}$, $j\in V$. Let $\textbf{1}({x\geq 0})$ be a function which takes value 1 if $x\geq 0$ and 0 otherwise. Given $z_{jk}$ for all $(j,k)\in\overrightarrow{E}$ and $\psi_{j}$ for $j\in V$ in the LN formulation, we map $w_{jk}=z_{jk}$ for all $(j,k)\in\overrightarrow{E}$ and $w_{jk}=\textbf{1}(\psi_{k}-\psi_{j}+1\geq 0)$ for all $(j,k)\notin\overrightarrow{E}$. Note that $w$-space is defined for all pair of nodes whereas $z$-space is defined for the set of arcs in $\overrightarrow{E}$. In every correct mapping, we have $w_{jk}=z_{jk}\text{,}\,\,\forall(j,k)\in\overrightarrow{E}$.

Fixing $j$ and $k$ for each $(j,k)\in\overrightarrow{E}$ and summing the left hand-side of inequalities (8e) over $i\in V\setminus\{j,k\}$ we obtain

[TABLE]

where the equivalences follow from constraints (8d). Given our mapping, $z_{jk}=w_{jk}$ for all $(j,k)\in\overrightarrow{E}$ and $\psi_{j}=\sum_{\ell\in V\backslash{\{j\}}}w_{\ell j}$ for $j\in V$, the above set of constraints can be written as

[TABLE]

which satisfies (10e) in the LN formulation. This implies $\mathcal{R}(LO)\subseteq\mathcal{R}(LN)$ for $\ell_{1}$-regularization. For $\ell_{0}$-regularization, we need to add that $g_{jk}\leq w_{jk}=z_{jk}\text{,}\,\forall(j,k)\in E^{\rightarrow}$. This implies $\mathcal{R}(LO)\subseteq\mathcal{R}(LN)$ for $\ell_{0}$ regularization.

To show strict containment, we give a point that is feasible to the LN formulation that cannot be mapped to any feasible point in the LO formulation. Consider $m=3$, $z_{13}=1,z_{31}=0,z_{32}=0.5+\epsilon,z_{23}=0.5-\epsilon,z_{12}=z_{21}=0.5$, $\psi=(1,2,2)$, for $0<\epsilon<\frac{1}{6}$, with an appropriate choice of $\beta$. It is easy to check that this is a feasible solution to the LN formulation. Because we must have $w_{ij}=z_{ij},\forall(i,j)\in\overrightarrow{E}$, we have $w_{13}+w_{32}+w_{21}>2$. Therefore, the corresponding point is infeasible to the LO formulation and this completes the proof.

$\square$

PROOF OF PROPOSITION 4

This proof is for TO formulation when the parameter $m$ on (9e) is replaced with $m-1$.

In the TO formulation, define the term $\sum_{s\in V}so_{ks}$ as $\psi_{k}$ and the term $\sum_{s\in V}so_{js}$ as $\psi_{j}$. Further, remove the set of constraints in (9f), (9g), and (9i). This implies that $\mathcal{R}(TO)\subseteq\mathcal{R}(LN)$. To see strict containment, consider the point described in the proof of Proposition 3, which is feasible to the LN formulation. For this point, there can be no feasible assignment of the decision matrix $o$ such that $\psi_{j}=\sum_{s\in V}so_{js}$, hence $\mathcal{R}(TO)\subset\mathcal{R}(LN)$.

$\square$

PROOF OF PROPOSITION 5

Let $\bar{F}(\beta_{X}^{\star})$ denote the optimal objective value associated with the continuous relaxation of model $X\in\{CO,LO,LN\}$.

Part A. $\bar{F}(\beta_{LO}^{\star})=\bar{F}(\beta_{LN}^{\star})$.

Case 1. $\ell_{1}$ regularization

Suppose (${\beta_{LN}^{\star}},z^{\star}$) is an optimal solution associated with the continuous relaxation of the LN formulation (10a)-(10g) and (${\beta_{LO}^{\star}},w^{\star}$) is an optimal solution associated with continuous relaxation of the LO formulation with $\ell_{1}$-regularization.

Given Proposition 3, we conclude that $\bar{F}(\beta_{LN}^{\star})\leq\bar{F}(\beta_{LO}^{\star})$. We prove that $\bar{F}(\beta_{LN}^{\star})\geq\bar{F}(\beta_{LO}^{\star})$ also holds in an optimal solution. To this end, we map an optimal solution in continuous relaxation of the $\ell_{1}$-LN formulation to a feasible solution in continuous relaxation of the $\ell_{1}$-LO formulation with the same objective function. This implies $\bar{F}(\beta_{LN}^{\star})\geq\bar{F}(\beta_{LO}^{\star})$.

Given an optimal solution (${\beta_{LN}^{\star}},z^{\star}$) to the continuous relaxation of the $\ell_{1}$-LN formulation, we construct a feasible solution $(\beta_{LO},w)$ to the continuous relaxation of the LO formulation as

$w_{jk}=w_{kj}=\begin{cases}\frac{1}{2},\quad z_{jk}^{\star}>0,\,(i,j)\in\overrightarrow{E},\\ 0,\quad\text{otherwise},\end{cases}$

and let $\beta_{LO}=\beta_{LN}^{*}$.

We now show that this mapping is always valid for the $\ell_{1}$-LO formulation. Recall that for the $\ell_{1}$ regularization, we do not have the decision vector $g$ in the formulation. Three set of constraints have to be satisfied in the $\ell_{1}$-LO formulation.

$|\beta_{jk}|\leq Mw_{jk},\quad\forall(j,k)\in\overrightarrow{E}$  (8b)

$w_{jk}+w_{kj}=1,\quad\forall(j,k)\in\overrightarrow{E}$ (8d)

$w_{ij}+w_{jk}+w_{ki}\leq 2,\quad\forall(i,j),(j,k),(k,i)\in\overrightarrow{E}\,i\neq j\neq k.$  (8e)

The set of constraints (8b) is trivially satisfied because we set $M\geq 2\underset{(j,k)\in\overrightarrow{E}}{\max}\,\quad|\beta_{jk}^{\star}|$. The set of constraints (8d) is trivially satisfied. The set of constraints (8e) is satisfied because the left hand side of inequality can take at most $\frac{3}{2}$ given this mapping. Therefore, $\bar{F}(\beta_{LO}^{\star})\leq\bar{F}(\beta_{LO})=\bar{F}(\beta_{LN}^{\star})$. This completes this part of the proof.

Case 2. $\ell_{0}$ regularization

Suppose (${\beta_{LN}^{\star}},g_{LN}^{\star},z^{\star}$) is an optimal solution associated with a continuous relaxation of the $\ell_{0}$-LN formulation, and (${\beta_{LO}^{\star}},g_{LO}^{\star},w^{\star}$) is an optimal solution associated with a continuous relaxation of the $\ell_{0}$-LO formulation.

Given Proposition 3, $\bar{F}(\beta_{LN}^{\star},g_{LN}^{\star})\leq\bar{F}(\beta_{LO}^{\star},g_{LO}^{\star})$. We now prove that, in an optimal solution $\bar{F}(\beta_{LN}^{\star},g_{LN}^{\star})\geq\bar{F}(\beta_{LO}^{\star},g_{LO}^{\star})$ also holds.

Given an optimal solution (${\beta_{LN}^{\star}},g_{LN}^{\star},z^{\star}$) for the continuous relaxation of the $\ell_{0}$-LN formulation, we construct a feasible solution $(\beta_{LO},g_{LO},w)$ for the continuous relaxation of the $\ell_{0}$-LO formulation as

$w_{jk}=w_{kj}=\begin{cases}\frac{1}{2},\quad z_{jk}^{\star}>0,\,(j,k)\in\overrightarrow{E},\\ 0,\quad otherwise,\end{cases}$

and let $\beta_{LO}=\beta_{LN}^{*}$ and $g_{LO}=g_{LN}^{\star}$.

We now show that this mapping is always valid for the LO formulation. Four sets of constraints have to be satisfied in the LO formulation.

$|\beta_{jk}|\leq Mg_{jk},\quad\forall(j,k)\in\overrightarrow{E}$, (8b)

$g_{jk}\leq w_{kj},\quad\forall(j,k)\in\overrightarrow{E}$ (8c)

$w_{jk}+w_{kj}=1,\quad\forall(j,k)\in\overrightarrow{E}$ (8d)

$w_{ij}+w_{jk}+w_{ki}\leq 2,\quad\forall i,j,k\in V,i\neq j\neq k,$  (8e)

The set of constraints in (8d) is satisfied similar to $\ell_{1}$ case. The proof that constraints (8b), (8c), and (8e) are also met is more involved. If the $\ell_{0}$-LN formulation attains a solution for which $g_{ij}^{\star}=g_{jk}^{\star}=g_{ki}^{\star}=1$ (we dropped subscript LN), then our mapping leads to an infeasible solution to the $\ell_{0}$-LO formulation, because it forces $w_{ij}+w_{jk}+w_{ki}\geq 2$ for $\ell_{0}$-LO. Next we show that this will not be the case and that our mapping is valid. To this end, we show that in an optimal solution for the continuous relaxation of $\ell_{0}$-LN, we always have $g^{\star}_{jk}\leq\frac{1}{2}\text{,}\,\forall(j,k)\in\overrightarrow{E}$. Note that in the LN formulation, we have $|\beta_{jk}|\leq Mg_{jk}\text{,}\,\forall(j,k)\in E^{\rightarrow}$ and $g_{jk}\leq z_{jk}$ (we dropped the subscript LN). Suppose that $g^{\star}_{jk}\geq\frac{1}{2}$. In this case, the objective function forces $g^{\star}_{jk}$ to be at most $\frac{1}{2}$. Note that the objective function can reduce $g^{\star}_{jk}$ up to $\frac{1}{2}$ and decreases the regularization term without any increase on the loss function. This is because $\beta_{jk}\leq Mg_{jk}\text{,}\,\forall(j,k)\in E^{\rightarrow}$ can be replaced by $\beta_{jk}\leq M\frac{1}{2}\text{,}\,\forall(j,k)\in E^{\rightarrow}$ without any restriction on $\beta$ because $M\geq 2\underset{(j,k)\in\overrightarrow{E}}{\max}\,|\beta^{\star}_{jk}|$. Therefore, our mapping is valid and implies that $\bar{F}(\beta_{LO}^{\star},g^{\star}_{LO})\leq\bar{F}(\beta_{LO},g_{LO})=\bar{F}(\beta_{LN}^{\star},g^{\star}_{LN})$.

Part B. $\bar{F}(\beta_{TO}^{\star})=\bar{F}(\beta_{LN}^{\star})$.

Case 1. $\ell_{1}$ regularization

Given Proposition 4, we conclude that $F(\beta_{LN}^{\star})\leq F(\beta_{TO}^{\star})$. We now prove that $\bar{F}(\beta_{LN}^{\star})\geq\bar{F}(\beta_{TO}^{\star})$ also holds in an optimal solution. We map an optimal solution in continuous relaxation of the $\ell_{1}$-LN formulation to a feasible solution in continuous relaxation of the $\ell_{1}$-TO formulation with the same objective value. This implies $\bar{F}(\beta_{LN}^{\star})\geq\bar{F}(\beta_{TO}^{\star})$.

Next we construct a feasible solution $(\beta_{TO},z_{TO},o)$ to the TO formulation. Given an optimal solution (${\beta_{LN}^{\star}},z_{LN}^{\star}$, $\psi^{\star}$) for a continuous relaxation of the LN formulation, rank the $\psi^{\star}_{j}$ in non-descending order. Ties between $\psi$ values can be broken arbitrarily in this ranking. Then, for each variable $j\in\{1,\dots,m\}$, $o_{jr}=1$ where $r$ denotes the rank of $\psi_{j}$ in a non-descending order (the first element is ranked 0). This mapping satisfies the two assignment constraints (9f)-(9g). Let $\beta_{TO}=\beta_{LN}^{*}$, $z_{TO}=z_{LN}^{*}$. This gives a feasible solution for the TO formulation. Thus, $\bar{F}(\beta_{LN}^{\star})\geq\bar{F}(\beta_{TO}^{\star}).$

Case 2. $\ell_{0}$ regularization

Define $g^{\star}_{LN}=g_{{TO}}$. The rest of the proof is similar to the previous proofs.

Part C. $\bar{F}(\beta_{CP}^{\star})=\bar{F}(\beta_{LN}^{\star})$.

To prove this, we first prove the following Lemma.

Lemma. The LO formulation is at least as strong as the CP formulation, that is $\mathcal{R}(LO)\subseteq\mathcal{R}(CP)$.

Proof. Consider (8d)-(8e) in the linear ordering formulation given by

[TABLE]

Consider the set of constraints in CP given by (6c) as

[TABLE]

Consider the same mapping as in Proposition 3. Select an arbitrary constraint from (6c). Without loss of generality, consider $g_{1,2}+g_{2,3}+\dots+g_{p-1,p}+g_{p,1}\leq p-1$. We arrange the terms in (8d)-(8e) as

[TABLE]

Summing the above set of inequalities and substituting $w_{ij}=1-w_{ji}$, when appropriate, gives $w_{1,2}+w_{2,3}+w_{3,4}+\dots+w_{p,1}\leq p-1$. Thus, we conclude that $\mathcal{R}(LO)\subseteq\mathcal{R}(CP)$. $\square$

Given this lemma and Proposition 3, we conclude that $\bar{F}(\beta_{LN}^{\star})\geq\bar{F}(\beta_{TO}^{\star})$. $\square$

PROOF OF PROPOSITION 6

This is a generalization of Proposition 5. Suppose we have branched on variable $w_{jk}$ in the LO formulation or correspondingly on variable $z_{jk}$ in the LN formulation. If $w_{jk}=0$, then $\beta_{jk}=0$. In this case, it is as if we now need to solve the original model with $(j,k)\notin\overrightarrow{E}$. Thus, Proposition 5 (Part A) implies that both models attain the same continuous relaxation. On the other hand, if $w_{jk}=1$ in LO (or correspondingly $z_{jk}=1$ in LN), we define our mapping as

$w_{jk}=w_{kj}=\begin{cases}\frac{1}{2},\quad z_{jk}^{\star}>0,\,(j,k)\in\overrightarrow{E},\\ 1,\quad(j,k)\in\mathcal{B}\\ 0,\quad\text{otherwise},\end{cases}$

where $\mathcal{B}$ is the set of $(j,k)$ for which $z_{jk}=1$. The rest of the proof follows from Proposition 5 (Part A).

The proofs for the TO and CP formulations are almost identical with Proposition 5. Suppose we have branched on variable $z_{jk}$ on any of the models (CP, LO, TO). If $z_{jk}=0$, then $\beta_{jk}=0$. In this case, it is as if we now need to solve the original model with $(j,k)\notin\overrightarrow{E}$. If $z_{jk}=1$, one defines the mapping in Proposition 5. The proof follows from Proposition 5 parts B and C, respectively. $\square$