On the Sparse DAG Structure Learning
Based on Adaptive Lasso

A DAG can be encoded as $G=(V,E)$, where $V$ refers to the set of nodes and $E$ is the set of directed edges. The nodes in $G$ represent variables under consideration. And directed edges correspond to the direct conditional relationships between parental nodes and child nodes. For example, $V_{1}\rightarrow V_{2}$, represents $V_{1}$ is the parent of $V_{2}$. Figure 2 illustrates three possible three-node DAGs with two edges.

The task of DAG structure learning is to seek a DAG from the DAGs space of $d$ nodes (discrete) that best fits the observational data $\mathbb{X}$. Let $\mathbb{X}\in\mathbb{R}^{n\times d}$ represent the observations consisting of $n$ $d$-dimensional vectors. We assume that these vectors are realizations of $n$ identically independently distributed (iid) random vectors with joint distribution $P(X)$, where $X=(X_{1},X_{2},...,X_{d})$. Importantly, with no latent variable and Markovian assumptions [40], the joint distribution $P(X)$ can be decomposed as follows:

where $Pa(X_{i})$ represents the parental set of $X_{i}$. Then, we model the relationships of variables in $X$ by a structure equation model (SEM), denoted as $X_{i}=f_{i}(Pa(X_{i}),\varepsilon_{i}),i=1,...,d$, where $\varepsilon_{i}$ refers to the uncovered variables. In this paper, we narrow our focus on the DAG with additive noise models (ANMs). By introducing the concept of weighted adjacency matrix $W\in\mathbb{R}^{d\times d}$, we can establish the linear ANMs as $X=XW+N$, where $N$ is a $d$-dimensional random vector, denoted as $N=(N_{1},N_{2},...,N_{d})$. Each $N_{i}$ represents the noise term added to $X_{i}$. In this paper, we assume $N_{1},N_{2},...,N_{d}$ are independent and identically distributed. The information about the structure of a DAG is completely involved in $W$, which is defined in the following way:

We note $G_{W}$ as the graph induced by $W$. To derive $W$, we can maximize the likelihood of observations $\mathbb{X}$. Benefiting from the well-known probability density transformation method [35], this target can be equally achieved by minimizing the least-square loss function $\frac{1}{2n}||X-XW||_{2}^{2}$ over the search space. A regularization term $\lambda||W||_{1}$ is added to deduce a sparse solution. Then, the overall score function is

Previous score-based methods suffer from searching the discrete DAGs space and lead to low efficiency. To overcome this problem, Zheng et al.[1] proposed the NOTEARS method, which seeks to use a function $h$ satisfying these four desiderata: i) $h(W)=0$ is a necessary and sufficient condition for $G_{W}$ to be a DAG; ii) its value reflects the level of acyclicity of the graph; iii) smooth; iv) its first derivative exists and is easy to compute. Based on the method in [41], which starts to use the matrix algebra to investigate the number of circles in a graph, NOTEARS verifies that $h(W)=\text{Tr}(e^{W\circ W})-d$ is a proper function to characterize the acyclicity with the derivative $\nabla h(W)=(e^{W\circ W})^{T}\circ 2W$. After determining the smooth equality constraint, the problem can be written as the following optimization program:

In order to solve the program (4), this constrained optimization program is converted to a sequence of unconstrained problems by the augmented Lagrangian method [42]. The procedure is shown as follows.

For each iteration, the first step can be solved by the L-BFGS proximal quasi-Newton method [43, 44] with the derivative as $-\frac{1}{n}X^{T}(X-XW)+\lambda\mathbbm{1}_{d\times d}+\left(\rho|h(W)|+\alpha_{k}\right)\nabla h(W).$ To overcome the problem caused by the non-smooth point zero in the penalty term, i.e., $W_{ij}=0,\quad i,j=1,..,d$, the original weighted adjacency matrix is cast into the difference between two matrices $(W^{1},W^{2})$ with non-negative elements [45].

Remark After obtaining the initial estimate, NOTEARS zeros out small coefficients if their absolute values are less than the predefined threshold to reduce the number of false positive discovery edges. However, this strategy leads to more missing edges with small weights and is hard to adapt to different scenarios. For example, if the weights of a DAG are generated from a uniform distribution over $(0,2)$, by setting $w_{t}=0.3$, the estimates of small weights, which are located in the interval $(0,0.3)$, have a high risk to be forced to $0$ in the post-processing step. Since the simulation experiments in previous work commonly set a gap to zero, for instance, in the simulation part of [1], the weights of edges commonly sampled uniformly from $[-2\alpha,-0.5\alpha]\cup[0.5\alpha,2\alpha]$, where $\alpha$ is a weight scale, this side effect is not obvious. The specific algorithm of NOTEARS with a fixed predefined threshold is presented in Algorithm 1.

In recent decades, variable selection has been studied extensively. In particular, Lasso [46] was introduced to improve the prediction accuracy and, simultaneously, reduce overfitting. In the linear model, the lasso estimates are obtained by optimizing the following program:

where $\lambda_{n}$ is related to the sample size $n$. Although Lasso variable selection has been shown to be consistent under some conditions [47], Zou[48] illustrated certain scenarios where the lasso is an inconsistent procedure in terms of variable selection. Zou[48] argued that the setup in lasso is somewhat unfair since it treats all the coefficients with indistinctive penalty levels in the sparsity term. To fix this problem, Zou[48] proposed a new methodology called adaptive lasso. In adaptive lasso, the weights of penalty levels are totally data-driven and cleverly chosen. The adaptive lasso estimates are obtained by optimizing the following program:

where $c_{i}$ is a predefined penalty weight. For example, we can use the ordinary least square estimates $\widehat{\beta_{ols}}$ to define $c_{i}$. Choose a $\gamma>0$, and define the penalty weight $c_{i}=\frac{1}{|\widehat{\beta_{ols,i}}|^{\gamma}}$, where $\gamma$ commonly be $0.5$, $1$, $2$ or determined by the cross-validation method.

From Figure 1, we can find that when the penalty level is relatively low, the learned false edges cannot be zeroed out. On the other hand, assigning a larger weight to the regularization term may lead to even worse results. This implies that Lasso put much more penalty on large coefficients rather than false positive ones, which is not consistent with our requirement. Motivated by this observation, to improve the criterion (4) on variable selection, one intuition is to apply the adaptive penalty levels to different coefficients, which has been developed in [39]. Following this motivation, we modify Eq. (4) to

where $c_{ij}$ is the specified penalty for $W_{ij}$. Expecting the minor false positive edges can be shrunk to zeros and reserve the true edges at the same time, we would like the corresponding $c_{ij}$ to be large for the former and small for the latter. Hence, the minor false edges are heavily penalized, and the true positive edges are lightly penalized. More details about how to choose $c_{ij}$ are given in Section III-C. With this modification of the score function, no extra predefined threshold is needed to post-processing the initial estimate anymore.

In this section, we show that with a proper choice of $c_{ij}$ and under some mild conditions, our method possesses the oracle properties when the sample number $n$ goes to infinity. Similar to procedure (5), via the augmented Lagrangian method, a dual ascent process can solve the ECP (8). The first step of each iteration can be regarded as finding the local minimum of the Lagrangian with a fixed Lagrange multiplier $\alpha$, that is,

Let $W^{*}$ denote the underlying true adjacency matrix. And then, we define

which means that $\mathcal{A}$ collects the indices for terms whose true parameters are nonzero, and $\mathcal{A}_{c}$ contains the indices for terms that do not exist in the underlying true model. Here, we first make two conditions:

1. 1.

   $X_{i}=X_{i}W^{*}+N_{i},i=1,2,...,n$, where $N_{1},N_{2},...,N_{n}$ are identically independent distributed with mean $\mathbb{0}$ and variance $\mathbb{\sigma}^{2}<\infty$.

2. 2.

   $\frac{1}{n}X^{T}X$ is finite and converges to a positive definite matrix $D\in\mathbb{R}^{d\times d}$.

Under these two conditions, the oracle properties of NOTEARS-AL are described in the following theorems. Rigorous proofs are given in Appendix VI.

In condition 1), we already assume $\frac{1}{n}X^{T}X\to D$ as $n\to\infty$. Furthermore, without loss of generality, we assume $\mathcal{A}_{i}=\{(i,1),(i,2),...,(i,d_{i})\}$. Let

where $D_{i0}$ is a $d_{i}\times d_{i}$ matrix corresponding to the coefficients with index in $\mathcal{A}_{i}$.

Firstly, for the zero coefficients part, Theorem 2 shows that given $a_{n}=o_{p}(n^{q})$ for some $q\leq-\frac{1}{2}$ and $\sqrt{n}b_{n}\rightarrow\infty$, NOTEARS-AL can consistently remove all irrelevant variables with probability tending to 1. Secondly, by magnifying the difference by $\sqrt{n}$ for the nonzero estimators, we can see that the pattern turns out to be a normal distribution. Then, there is only one question left. That is how to determine $c_{ij}$’s such that the nonzero related part $a_{n}=o_{p}(n^{q})$ for some $q\leq-\frac{1}{2}$ and zero related part $\sqrt{n}b_{n}\rightarrow\infty$. Here, using the fact that ordinary least square (OLS) estimates $\widehat{W}_{ols}$ are $\sqrt{n}$-consistent estimates of $W$ (the process of proof is consistent with Theorem1), suppose that $\lambda_{n}\sqrt{n}\rightarrow 0$ and $\lambda_{n}n^{\frac{\gamma+1}{2}}\rightarrow\infty$, then, we can define $c_{ij}$’s by $\widehat{W}_{ols}$ with some specific power $\gamma$ so that the above two requirements are met. Besides, specific details to determine $\lambda_{n}$ are given in III-E.

The algorithm can be divided into two parts. The first part is to find the solution $\widehat{W}_{ols}$ for Eq. (10), and the second part is plugging in the adaptive penalty parameters to Eq. (9) and solving it. These equality-constrained programs (ECPs) are then solved by the augmented Lagrangian method, i.e., using the Lagrangian method to handle the hard DAG constraint $h(W)=0$. The first part follows the same procedure as the original NOTEARS without the thresholding step. And then, we obtain a matrix $C\in\mathbb{R}^{d\times d}$ with element $c_{ij}=\frac{1}{|\widehat{W}_{ols_{ij}}|}$. Let $W_{C}\vcentcolon=C\circ W$, where $\circ$ is the Hadamard product. Thus, we have $W_{ij}=W_{C,ij}/c_{ij}$. In matrix notation, it is $W=W_{C}\varoslash C$, where $\varoslash$ refers to Hadamard division. Then, ECP (4) can be transformed as

where $\widehat{W}_{als}$ is the final estimate for $W$. Then, similar to (5), via the augmented Lagrangian method, the above ECP (11) can be solved by a dual ascent process.

To complete the above optimization, we need to derive a new form of the gradient for the first step in Eq. (12). Here we also reformulate the original $W$ as the difference between two matrices with positive elements. For the Loss function:

For the acyclicity term $h$:

To sum up, the algorithm proceeds as Algorithm 2.

Choosing the tuning parameter $\lambda_{n}$ is an important issue. We found that the results of NOTEARS-based methods are sensitive to the value of hyper-parameter, thus it is essential to choose a good one. Our goal is to choose the optimal model $s_{n}$ which contains the true variables as many as possible. Here, following Zou [48], we utilize cross-validation to find an optimal $\lambda_{n}$. From NOTEARS-AL, we can obtain a set of candidate models $S=\{s_{n},n=1,2,...,N\}$, where $s_{n}=\{(i,j)\in(1:d)\times(1:d):\widehat{W}_{n,ij}\neq 0\}$. The cross-validation procedure we use shows as Algorithm 3.

Another key modification is that we need a more sufficient validation set to find the best model to achieve the restricted model-selection consistency, especially when the candidate model set is large[49]. Motivated by this, instead of using the traditional $K$-fold cross-validation method, whose validation set is only $1/K$ of data, we suggest swapping the proportions of the validation set and training set.

Zheng et al.[27] also extend the idea of NOTEARS to generalized linear models (GLMs). Traditionally, a GLM is formulated as follows:

where $g_{i}$ is a link function used in GLMs, and $d_{i}\in\mathbb{R}^{d}$. To extend the DAG constraint beyond the linear setting, the key point is to find a sensible way to redefine the weighted adjacency matrix $W$ since there is no explicit form of $W$ in the GLMs setting. Based on the idea in [50], which encodes the importance of each variable considered via utilizing corresponding partial derivatives, it is then easy to show that the dependence among variables can be precisely measured by the $L^{2}$-norm of corresponding partial derivatives [27]. This implies the weighted adjacency matrix $W$ in GLMs can be defined as follows:

Thus, for linear mean functions, i.e. $f_{i}(X)=d_{i}^{T}X$, $W_{ki}=0$ is equivalent to $d_{ik}=0$. We can simplify the DAG constraint $h(W)=0$ by plugging in $W=({d_{1}}^{T},{d_{2}}^{T},...,{d_{d}}^{T})$. To construct the score function for GLMs, we need a more suitable loss function here. In this paper, we focus on one specific model in GLMs, the logistic regression model, where $X_{i}\in\{0,1\},i=0,...,d$. For logistic regression, the least-square loss is not appropriate anymore since it will result in a non-convex optimal function. Here we introduce the concept named Log Loss, which is commonly used as the loss function for logistic regression, defined as follows:

where $t_{i}=x_{i}\times({d_{1}}^{T},{d_{2}}^{T},...,{d_{d}}^{T})$.

Compared to the original weight-generating interval of NOTEARS, which is Uniform $([-2.0,-0.5]\cup[0.5,2.0])$, there are two points we consider to improve in our method. The first one is about the upper and lower bound, which are set to be $\pm 2\alpha$. We expect our method can extend NOTEARS to a wider range with a better estimate. The second one is the coefficients near $0$, which cannot be learned by the fixed-thresholding post-processing step. For the above two improvement points, we first perform the structure learning study in Gaussian-weight cases, where the weights of existing edges are sampled from a normal distribution with mean $0$ and variance $2^{2}$. In this way, we can investigate the performance of NOTEARS-AL under a wider range of coefficients.

In our experiments, we generate {ER$1$, ER$2$, ER$4$} graphs with different numbers of nodes $d=\{10,20,50\}$, $15$ graphs for each type. And then, we simulate $n=1000$ samples with iid Gaussian noises. The results are shown in Figure 3. For better visualization, only the average values of metrics are reported. Consistent with previous work GES, it shows a significant advantage in the ER$1$ case, no matter which metric we evaluate. But the performances of GES deteriorate rapidly when the level of sparsity decreases. Also, as we can see from Figure 3, for ER$4$ with $50$ nodes, the algorithm of GES takes too long to obtain the results because of the large searching space and low searching efficiency. Not surprisingly, the performance of NOTEARS and NOTEARS-AL show a positive relationship. This is natural since NOTEARS-AL extracts the weights of penalty level from conventional NOTEARS, which means that the performance of NOTEARS-AL is also sensitive to the results of the conventional one. Nonetheless, our approach performs uniformly better than NOTEARS across different graphs in terms of all metrics.

Figure 4 and Figure 5 show structure recovery results for SEM with exponential noise and Gumbel noise, respectively, where the weights of existing edges are sampled from a normal distribution with mean $0$ and variance $2^{2}$.

As we mentioned, one side effect of the fixed-threshold strategy is that small coefficients in the true model may be deleted indiscriminately. Since the simulation experiments in previous work commonly set a gap to zero, which is $\alpha[-2,-0.5]\cup\alpha[0.5,2]$, this side effect is not obvious. Hence, we investigate the value of missing edges when weights are generated from a range that covers zero, e.g., Uniform $(-5,5)$, Uniform $(-10,10)$, Gaussian $(0,2^{2})$ and Gaussian $(0,3^{2})$. For each kind of distribution, we generated $200$ DAG structures from the ER$1$ model with $10$ nodes and then simulated $1000$ samples. Based on the different methodologies of these two methods, one with a hard threshold of $0.1$ and the other with adaptive penalty levels, our method is more flexible and suitable to handle the coefficients around zero. Indeed, Figure 6 confirms this intuition. We plot the histogram of missing edges. In all three cases, our method (blue bar) loses fewer edges than NOTEARS with hard thresholding (yellow bar), especially the area around zero. Moreover, as the range becomes wider or the variance increases, the gap between these two methods becomes more and more significant.

To investigate the performance of our method for estimating DAG from high-dimensional graphical models, we increase the number of nodes to $100$ and compare the results to those of fixed-thresholding NOTEARS. Here, to illustrate the validity of the evaluation criteria, we draw error bars with sample standard deviations. The results are shown in Figure 7. We can find that under nearly the same value of SHD, the true positive rate of NOTEARS-AL is significantly higher than that of NOTEARS with a fixed threshold. Besides, in terms of FDR (lower is better), NOTEARS-AL has lower values for all graphs. In general, NOTEARS-AL consistently outperforms the original fixed-threshold one, even in high-dimensional settings.

For GLMs, we concentrate on logistic regression for dichotomous data $X_{i}\in\{0,1\}$, $i=1,...,d$ with link function $g(t)=e^{t}/(e^{t}+1)$. As we introduced in Section III-F, for linear mean functions, the elements of $W$ in DAG constraint $h(W)=0$ can be substituted by $d_{ik}$’s, the coefficients in linear mean functions. We compare the performance of fixed-thresholding NOTEARS to NOTEARS-AL with edge weights following Gaussian $(0,2^{2})$.

Figure 8 shows the structure recovery results for $n=1000$. By comparison, We can find the TPR of NOTEARS-AL is significantly higher than that of NOTEARS with a fixed threshold. At the same time, in terms of the other two metrics, i.e., SHD and FDR, both methods keep at a similar level. In general, NOTEARS-AL consistently outperforms the original NOTEARS.

Finally, consistent with Zheng et al. [1], we compare our method to NOTEARS on the same real dataset, which records the expression levels of $11$ phosphoproteins and phospholipids in human cells. This dataset is commonly used in the literature on Bayesian network structure inference as it provides a network that is widely accepted by the biological community. Based on $d=11$ cell types and $n=911$ observational samples, the ground truth structure given in [5] contains $17$ edges. In our experiments, NOTEARS with $0.1$ threshold estimates $22$ edges with an SHD of $20$, compared to $20$ estimates for NOTEARS-AL with an SHD of $19$.

Let $a_{n,ij}=\lambda_{n}c_{ij}$, $(i,j)\in\mathcal{A}$, we can derive

where we define

For $\@slowromancap i@$,

For $\@slowromancap ii@$,

For $\@slowromancap iii@$,

where $A=\max\{a_{n,ij}/n^{q}:(i,j)\in\mathcal{A}\}$, and $A\in O(1)$.