Revisiting Differentiable Structure Learning: Inconsistency of $\ell_{1}$ Penalty and Beyond

In the context of DAG learning, the DAG constraint, denoted as $h(B)$, ensures that the learned structure is a DAG when $h(B)=0$. NOTEARS (Zheng et al., 2018) employs a hard DAG constraint, whereas GOLEM (Ng et al., 2020) introduces a soft DAG constraint.

In structure learning for linear Gaussian cases, score-based methods, specifically with the BIC score (Schwarz, 1978; Chickering, 2002), often aim to recover the sparsest underlying DAG that best explains the observed data (Singh & Moore, 2005; Cussens, 2011; Yuan & Malone, 2013; Chickering, 2002). In the asymptotic case where the population covariance matrix, denoted by $\Sigma^{*}$, is available, this can, loosely speaking, be formulated as

$$\min_{B,\Omega}\;\|B\|_{0}\quad\text{subject to}\quad(I-B)^{-\top}\Omega(I-B)^{-1}=\Sigma^{*}\quad\text{and}\quad B\text{ is a DAG}.$$

Recall that $B$ is the weighted adjacency matrix representing the structure of the DAG and $\Omega$ is the diagonal matrix of noise variances. That is, the goal is to minimize the $\ell_{0}$ norm of $B$, i.e., the number of edges, while maximizing the likelihood by satisfying the covariance constraint and ensuring that $B$ is a valid DAG. In other words, the objective is to recover the sparsest DAG, $\hat{B}$, along with its corresponding $\hat{\Omega}$, that can generate the observed covariance matrix $\Sigma^{*}$ (i.e., $(I-\hat{B})^{-\top}\hat{\Omega}(I-\hat{B})^{-1}=\Sigma^{*}$). Under the sparsest Markov faithfulness assumption (Raskutti & Uhler, 2018), the estimated $\hat{B}$ will be Markov equivalent to the true graph $B^{*}$.

Many previous work, including GOLEM, replace the $\ell_{0}$ penalty with the more tractable $\ell_{1}$ penalty. The corresponding optimization problem becomes

$$\min_{B,\Omega}\;\|B\|_{1}\quad\text{subject to}\quad(I-B)^{-\top}\Omega(I-B)^{-1}=\Sigma^{*}\quad\text{and}\quad B\text{ is a DAG}.$$

(1)

The $\ell_{1}$ penalty encourages smaller edge weights but introduces a key inconsistency: it doesn’t guarantee true sparsity in the resulting structure. Minimizing the $\ell_{1}$ norm may lead to a denser structure with more edges than the solution from minimizing the $\ell_{0}$ norm. This occurs because $\ell_{1}$ favors edges with small absolute values, even if they represent spurious edges. In some cases, the sum of the absolute values of more edges can be smaller than that of fewer, larger-magnitude edges, which leads $\ell_{1}$-based methods to include unnecessary edges. As a result, the learned structure may deviate from the true DAG or its Markov equivalence class. Taking GOLEM as an example, even if the covariance constraint $(I-B)^{-\top}\Omega(I-B)^{-1}=\Sigma^{*}$ is satisfied in both $\ell_{0}$- and $\ell_{1}$-based formulations, the structural properties of the solution can differ significantly.

In Section 3.2, we demonstrate this with a large number of experiments, showing that a large proportion of DAGs $\tilde{B}$ satisfying the covariance constraint $(I-\tilde{B})^{-\top}\Omega(I-\tilde{B})^{-1}=\Sigma^{*}$ have smaller $\ell_{1}$ norms than the true graph $B^{*}$. Moreover, in these DAGs satisfying the covariance constraint and having smaller $\ell_{1}$ norms than the true graph $B^{*}$, we can always find ones with much larger $\ell_{0}$ values than $B^{*}$, meaning they do not correspond to the true graph $B^{*}$ or its Markov equivalence class, with the structural Hamming distance (SHD) computed over the true and estimated CPDAG (Completed Partially Directed Acyclic Graph), which we refer to as SHD of CPDAG in this paper, far from zero.

In this section, we demonstrate and evaluate the inconsistency of the $\ell_{1}$ penalty in likelihood-based GOLEM through experiments. We generated 1,000 true DAGs $B^{*}$, compute their corresponding covariance matrices $\Sigma^{*}$ under infinite sample conditions. For each $B^{*}$, we use Cholesky decomposition to generate large number of DAGs $\tilde{B}$ that can generate the same $\Sigma^{*}$. Our goal is to identify the $\tilde{B}$ with the minimum $\ell_{1}$ norm among these $\tilde{B}$s, denoted as ${B}_{\ell_{1}}$, and compare it with $B^{*}$ in terms of $\ell_{1}$ norm, $\ell_{0}$ norm (i.e., edge count), and record its SHD of CPDAG. Additionally, we also record the proportion of $\tilde{B}$s that satisfy the covariance constraint (i.e., generate the same $\Sigma^{*}$) but have a smaller $\ell_{1}$ norm than $B^{*}$, to give an intuitive sense of the extent of $\ell_{1}$ inconsistency.

In this section, we present a 3-node counterexample to demonstrate the inconsistency of the $\ell_{1}$ penalty in differentiable structural learning. Specifically, we compare a true weighted adjacency matrix $B^{*}$ with an estimated adjacency matrix $\tilde{B}$, and show that although both matrices can generate the same covariance matrix, their $\ell_{0}$ norm (edge count), $\ell_{1}$ norm, and structural differences, measured by SHD of CPDAG, reveal the inconsistency of the $\ell_{1}$ penalty.

The true adjacency matrix $B^{*}$ and its corresponding noise covariance matrices $\Omega^{*}$ are given as:

$$B^{*}=\begin{bmatrix}0&\frac{1}{2}&0\\ 0&0&0\\ 0&-1&0\end{bmatrix},\quad\Omega^{*}=\begin{bmatrix}16&0&0\\ 0&4&0\\ 0&0&1\end{bmatrix}.$$

The estimated adjacency matrix $\tilde{B}$ and its corresponding noise covariance matrices $\tilde{\Omega}$ are:

$$\tilde{B}=\begin{bmatrix}0&\frac{1}{2}&\frac{1}{10}\\ 0&0&-\frac{1}{5}\\ 0&0&0\end{bmatrix},\quad\tilde{\Omega}=\begin{bmatrix}16&0&0\\ 0&5&0\\ 0&0&\frac{4}{5}\end{bmatrix}.$$

Both matrices $B^{*}$ and $\tilde{B}$, along with their respective noise covariance matrices, generate the same covariance matrix:

$$\Sigma^{*}=\begin{bmatrix}16&8&0\\ 8&9&-1\\ 0&-1&1\end{bmatrix}.$$

We have $\|B^{*}\|_{0}=2$ and and $\|\tilde{B}\|_{0}=3$, indicating that $B^{*}$ is sparser than $\tilde{B}$.

However, when considering the $\ell_{1}$ norm, we observe that: $\|B^{*}\|_{1}=\frac{3}{2}\quad\text{and}\quad\|\tilde{B}\|_{1}=\frac{4}{5}.$ That is, although $\tilde{B}$ has a higher $\ell_{0}$ norm, it achieves a lower $\ell_{1}$ norm, highlighting the inconsistency between the two norms. Therefore, the optimization problem in Eq. equation 1 may return $\tilde{B}$, which is clearly not Markov equivalent to $B^{*}$.

This counterexample demonstrates the inconsistency of the $\ell_{1}$ penalty: it may lead to solutions with smaller total edge weights (resulting in a lower $\ell_{1}$ norm), but these solutions may still have more edges (a higher $\ell_{0}$ norm) and deviate from the true DAG structure and its Markov equivalence class, even if these solutions can generate the same covariance matrix as the ground truth DAG.

Across all experiments in section 4, we simulate Erdös–Rényi graphs (ERDdS & R&wi, 1959) with $kd$ edges, denoted as ERk graphs, with edge weights uniformly sampled from $[-2,-0.5]\cup[0.5,2]$. For all experiments, the data is generated with a fixed noise ratio of 16. Specifically, the variances of two randomly selected noise variables are set to 1 and 16, respectively, while the variances of the remaining noise variables are sampled uniformly from the range $[1,16]$. This setting ensures a realistic variation in noise across the variables, aligning with the assumptions of non-equal noise variances (NV). Further implementation details of our experiments in section 4 are in Appendix A. From this point forward, unless otherwise specified, we use ”CALM” to refers to the specific version of the method where the $\ell_{0}$-penalty is approximated using Gumbel Softmax.

The following four sub-sections evaluate CALM’s performance: first, we compare soft vs. hard DAG constraints and moral graph vs. no moral graph; second, we assess the effect of data standardization; third, we test different $\ell_{0}$ approximation methods. Finally, we compare CALM with baseline approaches. For each scenario, we conducted 10 experiments and calculated the mean SHD of CPDAG, precision of skeleton, recall of skeleton, and their standard errors.

Recall that NOTEARS (Zheng et al., 2018) adopts a hard DAG constraint while GOLEM (Ng et al., 2020) uses a soft DAG constraint. Here, we evaluate the impact of incorporating the moral graph and using either soft or hard DAG constraints on the results of the $\ell_{0}$-penalized likelihood optimization. We consider linear Gaussian models with 50 variables and ER1 graphs. Here, we focus on the nonconvexity aspect of the optimization problem, and thus set the sample size to infinity to eliminate finite sample errors (the way we achieve infinite samples is in Appendix A.4). The experiment results for finite samples are included in Section 4.6. Furthermore, following Reisach et al. (2021); Kaiser & Sipos (2022), we standardize the data. All experiments were conducted with the Gumbel-Softmax approach to approximate $\ell_{0}$ penalty. The implementation details of Gumbel-Softmax-based $\ell_{0}$ penalty and the hard DAG constraints is in Appendix A.2. The implementation details of soft DAG constraints is in Appendix A.3.

Ng et al. (2024) previously pointed out that the original GOLEM-NV-$\ell_{1}$ formulation performed poorly both before and after data standardization. To further evaluate the robustness of CALM, we conduct experiments to compare its performance with (CALM-Standardized) and without data standardization (CALM-Non-Standardized). Here, we use infinite samples to eliminate finite sample error and consider a 50-node linear Gaussian model with ER1 graphs. The implementation details of Gumbel-Softmax-based $\ell_{0}$ penalty and the hard DAG constraints for CALM is in Appendix A.2. In Table 3, CALM shows consistently low SHD of CPDAG before and after data standardization, demonstrating its stability and robustness across both standardized and non-standardized data. Interestingly, this is in constrast with the observation by Reisach et al. (2021); Kaiser & Sipos (2022) that differentiable structure learning methods do not perform well after data standardization, which further validate the robustness of our method.

We compare our method with three $\ell_{0}$ approximations (CALM, CALM-STG, CALM-Tanh) and the original GOLEM-NV-$\ell_{1}$. For all methods here, We used infinite samples to eliminate finite sample error and considered 50- node linear Gaussian model with ER1 graphs. Additional results for ER4 with 50 nodes and ER1 with 100 nodes are presented in Appendix B. Here, CALM maintains our definition, specifically referring to our method that employs the Gumbel-Softmax approximation for the $\ell_{0}$ penalty. CALM-STG (Stochastic Gates) refers to our method that uses stochastic gates to approximate $\ell_{0}$ penalty (Yamada et al., 2020), while CALM-Tanh (Hyperbolic Tangent) refers to our method that employs the smooth hyperbolic tangent function to approximate $\ell_{0}$ penalty (Bhattacharya et al., 2021). The key parameter selection and implementation details for STG and Tanh in approximating the $\ell_{0}$ penalty can be found in Appendix A.5.

Table 4 shows that CALM-STG achieves competitive results after data standardization but underperforms without standardization. CALM-Tanh shows the weakest performance among the three $\ell_{0}$ approximations methods. Only CALM performs well both with and without data standardization, highlighting its robustness, which is why we ultimately selected Gumbel-Softmax as the $\ell_{0}$ approximation method as our representative implementation. Additionally, all three methods outperform GOLEM-NV-$\ell_{1}$, underscoring the importance of $\ell_{0}$ penalty approximation, as $\ell_{1}$ penalty suffers from inconsistency (as discussed in Section 3). Moreover, this also shows that the effectiveness of combining $\ell_{0}$ approximation with moral graphs and hard constraints.

We finally compare the performance of CALM against several baseline methods, including the original GOLEM-NV-$\ell_{1}$, NOTEARS, PC (Spirtes & Glymour, 1991), and FGES (Ramsey et al., 2017) (see Appendix A.6 for baseline methods’ implemention details). We evaluated the methods at two different sample sizes: $n=1000$ and $n=10^{6}$. We considered a 50-node linear Gaussian model with ER1 and ER4 graphs, as well as a 100-node linear Gaussian model with ER1 graphs. Here, the moral graph in CALM is estimated from finite samples (see Appendix A.1 for how we estimate the moral graph). Specifically, for the 1000-sample experiments, the moral graph is estimated from 1000 samples, and for the $10^{6}$-sample experiments, the moral graph is estimated from $10^{6}$ samples. All results presented here are from experiments with standardized data. The additional experimental results for data without standardization are presented in Appendix C.

Table 5 summarizes the performance comparison between CALM and the baseline methods. The results clearly demonstrate that CALM consistently outperforms both NOTEARS and the original GOLEM-NV-$\ell_{1}$ across all graph structures and sample sizes. This highlights the effectiveness and robustness of incorporating the Gumbel-Softmax approximation to $\ell_{0}$, moral graph, and hard DAG constraints. It is observed that the results of CALM are not as competitive as those obtained by the PC and FGES methods for sparse graphs such as ER1 graphs. This outcome is expected, given that the continuous optimization in linear likelihood-based formulation struggles with such high levels of nonconvexity.

However, it is worth noting that for ER1 graphs, in the case of 1000 samples, the results of CALM are nearly identical to those of PC. This indicates that in practical scenarios with relatively small sample sizes, even in sparse graphs, CALM is able to compete with the performance of discrete methods like PC. This represents a significant breakthrough.

Furthermore, in more dense graphs, the 50-node ER4 graphs, CALM demonstrates superior performance compared to the PC and FGES methods. This result suggests that in higher-density graphs, CALM enables continuous optimization methods to outperform discrete methods.