CoLiDE: Concomitant Linear DAG Estimation

Authors’ response to Reviewer VEfR (Part 1)

Thanks for your time and effort in reviewing our manuscript, as well as for finding the problem highly relevant due to the importance of considering heteroscedastic settings in linear DAG learning. Moreover, we appreciate your valuable suggestions to improve our numerical experiments in several ways, and your request for further elaboration on model identifiability issues. Point-by-point responses to your constructive comments and associated requests for changes follow, which we believe have led to an improved revised paper. We strive to improve our paper and will be happy to continue the scholarly discussion if any lingering issues remain. **Identifiability.** In the second paragraph of Section 1 we discuss different choices of score functions (likelihood and regression-based) and their relative merits. But with regards to guarantees, you bring up an excellent point that certainly benefits from additional elaboration. In (Loh and Buhlmann, 2014), the overarching assumption is that *the exogenous noise covariance matrix is known* up to a common scale across variables (e.g., in the homoscedastic setting), while CoLiDE also endeavors to estimate the exogenous noise variances. In the homoscedastic setting, the weighted least squares (LS) criterion therein and the CoLiDE-EV score function are equivalent, since the latter is obtained from weighted LS via a constant scaling and shift that will not affect the optimal $\mathbf{W}$ solution. Hence, the identifiability results in (Loh and Buhlmann, 2014) will carry over for both linear Gaussian and non-Gaussian SEMs, again, provided the diagonal noise covariance matrix is known up to a constant factor. Now, the story is quite different in the more challenging heteroscedastic setting. As we now spell out in the introduction of the revised paper, for general linear Gaussian models the DAG structure is non-identifiable from observational data alone. Interestingly, just like GOLEM (Ng et al, 2020) and for general (non-identifiable) linear Gaussian SEMs, we can show that as $n\to\infty$ CoLiDE-NV probably yields a DAG that is *quasi-equivalent* to the ground truth graph. These guarantees for (heteroscedastic) linear Gaussian SEMs build on the interesting theoretical framework put forth in (Ghassami et al, 2020); see also Appendix C. Comments along these lines have been included in the revised manuscript; please check the paragraph immediately preceding Section 4.1. **Data standardization in heteroscedastic settings.** Following your suggestion, we have conducted an additional experiment involving standardized data in a heteroscedastic setting; please check Appendix E.4 in the revised manuscript. Despite an overall degradation in performance of all tested baselines resulting from such standardization (see Figure 9), CoLiDE-NV maintains its status as the superior method when compared to other state-of-the-art approaches. **Improved metrics over Markov equivalence classes.** This point is well taken. Following your suggestion, we have introduced a new metric termed SHD-C, defined in Appendix D.3 and computed as follows. We initially map both the estimated graph and the ground truth to their respective Completed Partially Directed Acyclic Graphs and subsequently calculate the SHD between them. Notably, SHD-C has been employed in (Ng et al, 2020); which also dealt with (non-identifiable) linear Gaussian SEMs. *SHD-C results are now reported in all relevant tables.* Furthermore, to complement the experiments presented in the main body of the paper, we have included Figures 5 and 8 in Appendix E. **Gaussian log-likelihood in the heteroscedastic setting.** Thanks again for this constructive feedback. We conducted new experiments specifically aimed at assessing the performance of the CoLiDE score functions in comparison with the Gaussian log-likelihood (as utilized in GOLEM), *by decoupling their effect from the need of a (soft or hard) DAG constraint*. Please check Appendix E.9 in the revised manuscript for a detailed description of the experimental setting. As an ablation study, we believe this is even more informative than augmenting both score functions with, say, the DAGMA acyclicity function. Our results indicate that for smaller DAGs (say with 10 nodes), the log-likelihood score functions outperform concomitant estimators. However, as the number of nodes increases, the score functions utilized in CoLiDE-EV and CoLiDE-NV exhibit superior performance (measured in terms of normalized SHD and SHD-C) when compared to the log-likelihood score functions of GOLEM-EV and GOLEM-NV.

Authors’ response to Reviewer VEfR (Part 2)

**Comparison with GOLEM.** Certainly, we always followed the guideline of initializing GOLEM-NV with the output of GOLEM-EV, as suggested in (Ng et al, 2020). To avoid any ambiguity, we have spelled this out in the revised manuscript; please check Appendix D.5. Our GOLEM experiments are based on the profiled likelihood approach in (Ng et al, 2020), and we used Ng et al’s own implementation. We carefully checked their paper and could not find a version of their method without profiling. **Hyperparameter tuning.** Our optimization methodology closely follows DAGMA. Given that the experimental settings are similar, for fair comparison we selected $\lambda=0.05$ and all other hyperparameters as in (Bello et al, 2022). To assess the robustness of said choice, we tested several other $\lambda$ values and found $0.05$ to be preferable. Importantly, CoLiDE decouples $\lambda$ from $\sigma$, eliminating the need for recalibration when there are variations in the noise levels. Clarifications along these lines have been included in the revised manuscript; please check Section 4.1. For other baselines, we chose the recommended hyperparameters as outlined in Appendix C.5. In closing, we would like to clarify that while some experimental settings considered here are different from those in the GOLEM and DAGMA papers, a few (see e.g., Table 1) are purposely identical to those therein; and CoLiDE’s advantage over its competitors is apparent. **Decoupling of sparsity regularization parameter and noise level.** We apologize for the lack of clarity in this key concept, which was probably due to our suboptimal exposition. We clarify this here. Optimal rates for lasso hinge on selecting $\lambda \asymp \sigma \sqrt{\log d / n}$; see e.g, (Li et al, 2020). However, the exogenous noise variance $\sigma^2$ is rarely known in practice. Interestingly, by virtue of the scaled residuals in concomitant lasso estimators then minimax optimality requires $\lambda\asymp \sqrt{\log d/n}$ (Li et al, 2020). This strategy renders $\lambda$ independent of $\sigma$, which is exactly what we mean when we say `decoupling the sparsity parameter from the exogenous noise levels’. Clarifying comments along these lines were included in the revised manuscript; see Sections 2, 4, and Appendix A.1. Thanks again for your valuable feedback that has contributed to an improved revised manuscript, and we look forward to addressing any further concerns that may arise.

Authors’ response to Reviewer wJd3

Thanks for your time and effort in reviewing our manuscript, as well as for finding the problem of broad interest to the ICLR community, the empirical analysis extensive, and the paper well written. We are glad to hear you believe this work’s impact can be potentially high, with benefits that could permeate to other continuous optimization approaches to DAG learning. Even more, we believe our ideas could also contribute to order-based methods, and in the revised paper we conduct a preliminary test based on valuable feedback provided by Reviewer xkAR. Point-by-point responses to your comments and associated requests for changes follow. We strive to improve our paper and will be happy to continue the scholarly discussion if any outstanding issues remain. **Decomposability of scale estimator.** This is an excellent point. It is true that (unlike other approaches) the adjacency matrix $\mathbf{W}$ can be efficiently updated with stochastic gradient steps using individual samples or mini-batches, but *the same holds true for the noise scale estimators* in (4) and (6). This follows because these estimators depend on $\operatorname{cov}(\mathbf{X}) := \frac{1}{n} \mathbf{X} \mathbf{X}^{\top}$, and the sample covariance matrix can be updated *recursively* from individual samples or mini-batches. In this sense, the noise scale estimators are decomposable as well and can be implemented from partial or streaming data using, e.g., online algorithms that we sketch in the revised manuscript. Please check the `Online updates’ sections under Appendices B.1 and B.2. A full exploration of this exciting direction is left as future work. Moreover, in the revised paper we explicitly state that CoLiDE’s score function is fully decomposable; see the discussion following (3). **Noise-dependent term in the score function.** The origin of concomitant scale estimators can be traced back to (Huber, 1981), who proposed *convex* criteria to jointly estimate linear regression coefficients along with the noise standard deviation $\sigma$. A first key ingredient is to note that for convex loss $\rho(u)$ (e.g., $\rho(u)=u^2$ in ordinary least squares), then it follows that $\rho(u/\sigma)\times \sigma$ is jointly convex in $(u,\sigma)$. Second, Huber pointed out that the inclusion of the (convexity preserving) linear term $d\sigma/2$ as in (2), yields an estimator $\hat{\sigma}^2$ that is consistent under Gaussianity [it is actually the maximum likelihood estimator; see (18)]. A clarifying comment on the reason behind the addition of this linear term in the objective function has been included in the revised paper; please check the discussion following (2). It should thus be clear that the sparsity regularization (which only acts on $\mathbf{W}$) should not affect the *form* of the noise scale estimator $\hat{\sigma}$. It may affect noise estimation *performance* though, since $\hat{\sigma}$ depends on $\hat{\mathbf{W}}$ [cf. (4)]. One could certainly use the edge cardinality function $||\mathbf{W}|||_0$ as a penalty, but we opted against it so as to retain convexity of the score function $\mathcal{S}(\mathbf{W},\sigma)$. **Hyperparameter tuning.** Our optimization methodology closely follows DAGMA. Given that the experimental settings are similar, for fair comparison we selected $\lambda=0.05$ as in (Bello et al, 2022). To assess the robustness of said choice, we tested several other $\lambda$ values and found $0.05$ to be preferable. Importantly, CoLiDE decouples $\lambda$ from $\sigma$, eliminating the need for recalibration when there are variations in the noise levels. Clarifications along these lines have been included in the revised manuscript; please check Section 4.1. **SortNRegress in Figure 1.** Following your suggestion, we have incorporated results for SortNRegress into Figure 1. We agree that plotting these results helps better convey the non-trivial nature of the task at hand. We have also fixed the typo in the gradient expression (16). Good catch! Thanks again for your valuable feedback that has contributed to an improved revised manuscript, and we look forward to addressing any further concerns that may arise. We appreciate your willingness to increase your rating.

Authors’ response to Reviewer xkAR (Part 1)

Thanks for your time and effort in reviewing our manuscript, as well as for recognizing the contributions and CoLiDE’s improved performance over state-of-the-art approaches. We are glad to hear you found the paper to be clearly written. Moreover, we appreciate your valuable suggestions to further improve DAG structure estimation, by bringing to bear exciting optimization advances along with the CoLiDE score function. Point-by-point responses to your comments and associated requests for changes follow. We strive to improve our paper and will be happy to continue the scholarly discussion if any lingering issues remain. **Significance.** This point is well taken. Indeed, as even the paper title indicates the scope of this study is admittedly limited to DAG structure learning from observational data adhering to *linear* SEMs. We are the first ones to acknowledge this as one of the limitations of this work (see Section 6). While it is clear that CoLiDE *can* be extended to non-linear SEMs, we feel it is not prudent to include additional experiments in what is already a fully-packed paper. This issue notwithstanding, we contend that the linear case is still a significant and challenging problem, particularly in the heteroscedastic setting dealt with where. CoLiDE, our novel score-based approach exhibits significant improvements relative to existing state-of-the-art methods, across a broad swath of metrics and settings spanning different graph ensembles, exogenous noise distributions (with and without equal variances across variables), and even a real dataset that has been widely used for benchmarking. In terms of potential broader impacts, our novel convex score function incorporates concomitant estimation of scale and has wide applicability when combined with (or used in place of) conventional loss functions in other state-of-the-art DAG learning methods. Finally, while this work is framed within the continuous constrained relaxation paradigm to linear DAG learning, preliminary experiments with TOPO (thanks for this valuable suggestion!) show that CoLiDE also benefits order-based methods. **Novelty.** This work contributes several innovations to the timely field of learning DAG structure from observational data. First, we propose a new convex score function for sparsity-aware learning of linear DAGs, which incorporates concomitant estimation of scale parameters to enhance DAG topology inference using continuous first-order optimization. To the best of our knowledge, this is the first time that ideas from concomitant scale estimation permeate benefits to DAG learning. Second, we develop scalable (inexact) block-coordinate descent iterations to jointly estimate the DAG topology along with the exogenous noise levels – both in homoscedastic and heteroscedastic settings and with a marginal increase in complexity over structure estimation alone. Unlike some other existing approaches, we also show our score function is fully decomposable across samples, thus facilitating mini-batch based optimization if needed for scalability. Third, in existing methods the sparsity regularization parameter depends on the unknown exogenous noise levels, making the calibration challenging. CoLiDE effectively removes this coupling, leading to minimum (or no) recalibration effort across diverse problem instances – hence a valuable form of robustness that leads to tangible computational savings. Finally, through a comprehensive experimental evaluation protocol we demonstrate CoLiDE outperforms state-of-the-art methods without incurring added complexity, especially when the DAGs have several hundreds of nodes and the noise level profile is heterogeneous.

Authors’ response to Reviewer xkAR (Part 2)

**Identifiability.** This is an excellent point that certainly benefits from additional elaboration. In (Loh and Buhlmann, 2014), the overarching assumption is that *the exogenous noise covariance matrix is known* up to a common scale across variables (e.g., in the homoscedastic setting), while CoLiDE also endeavors to estimate the exogenous noise variances. In the homoscedastic setting, the weighted least squares (LS) criterion therein and the CoLiDE-EV score function are equivalent, since the latter is obtained from weighted LS via a constant scaling and shift that will not affect the optimal $\mathbf{W}$ solution. Hence, the identifiability results in (Loh and Buhlmann, 2014) will carry over for both linear Gaussian and non-Gaussian SEMs, again, provided the diagonal noise covariance matrix is known up to a constant factor. For the aforementioned reasons, we respectfully believe there will be no added benefit in combining CoLiDE-EV with the weighted LS loss in (Loh and Buhlmann, 2014). Now, the story is quite different in the more challenging heteroscedastic setting. As we now spell out in the introduction of the revised paper, for general linear Gaussian models the DAG structure is non-identifiable from observational data alone. Interestingly, just like GOLEM (Ng et al, 2020) and for general (non-identifiable) linear Gaussian SEMs, we can show that as $n\to\infty$ CoLiDE-NV probably yields a DAG that is *quasi-equivalent* to the ground truth graph. These guarantees for (heteroscedastic) linear Gaussian SEMs build on the interesting theoretical framework put forth in (Ghassami et al, 2020); see also Appendix C. Comments along these lines have been included in the revised manuscript; please check the paragraph immediately preceding Section 4.1. **Combining the CoLiDE score function with TOPO.** Thanks for bringing both relevant works by Deng et al to our attention, which we now cite and discuss in the revised manuscript. Following your excellent suggestion, we have implemented a new method that combines the convex CoLiDE score functions along with TOPO. Our encouraging preliminary results in Table 4 show that utilizing TOPO as an optimizer enables CoLiDE-NV to attain better performance, with up to 20% improvements in terms of SHD and for the particular setting tested (10 different ER graphs with 50 nodes from a linear Gaussian SEM, where the noise variance of each node was randomly selected from the range $[0.5, 10]$). Additional metrics are presented in Table 4 for a comprehensive evaluation. Please check Appendix E.8 for additional details of this new experiment. We have appended our CoLiDE-TOPO implementation to the anonymized code that we made publicly available as supplementary material accompanying this submission. In closing, we note that the selection of the solver depends on the problem at hand. While discrete optimizers excel over smaller graphs, continuous optimization tends to be preferred for larger DAGs (with hundreds of nodes). Our experiments show that CoLiDE offers advantages in both of these settings, which reinforces the significance of the proposed score function for linear DAG learning. Thanks again for your valuable feedback that has contributed to an improved revised manuscript, and we look forward to addressing any further concerns that may arise.

Thanks

We are glad to hear our responses clarified your questions and we agree the paper is stronger now. We appreciate your feedback and willingness to adjust your score.