Constraint-Free Structure Learning with Smooth Acyclic Orientations

We thank the reviewer for appreciating our main goal and the soundness of our proposal. We took particular care in considering their insightful comments, questions, and concerns, which we will now gladly address.

Contribution Claim. In line with yours and other reviewers' comments, we have adjusted our claim. Further, as stated in the Related Works, we only considered for our empirical evaluation approaches that ensure the acyclicity of the solution, while GOLEM requires to iteratively threshold the solution to ensure acyclicity. Given the reviewer's suggestion, we included in Appendix E.17 an empirical comparison with GOLEM to further motivate our claim. Concerning the comparison with DP-DAG, even if possibly less apt for practical applications, simpler and scalable optimization schemes like ours are more suitable to be integrated as components in continuous causal discovery approaches. Given the reviewer's concern, we further clarified the issue in the Related Works and reported additional experiments in Appendix E.19. The additional experiments report results for DAGuerrotype (Zantedeschi et al., 2023), a similar but faster method comparable to DP-DAG, which however already highlights significant scalability issues when compared to optimization-schemes such as COSMO.

AUC, NHD, and Acyclicity By following NOTEARS, most methods perform a thresholding step to prune off small weights and spurious arcs (Ng et al., 2023). To compare with them we adopted the same practice to report NHD, which can be computed only on binary matrices. Overall, we consider AUC to be a more suitable metric since it does not require thresholding the matrix and better captures the tradeoff between true positive and false positive ratio. Notably, the same reasoning underlies the metric chosen for validating DP-DAG in Charpentier et al. (2022). Concerning the reviewer's question, our solution is guaranteed to be acyclic before the thresholding step. To further back this claim, we updated the manuscript with a new section in Appendix G, which reports the trend of temperature and acyclicity and shows how the model reaches acyclicity zero even for non-zero temperatures.

Non-Linear. We agree with the reviewer that there is indeed a performance drop in the non-linear scenario, which is however shared by the baselines and mainly due to the increased difficulty of the task itself. We did not include NOTEARS-MLP for larger graphs because of the significantly higher computational time required for already 20 nodes in the non-linear scenario.

Regularization. L1 regularization favors feature selection and sparsity that, while fundamental for the overall weighted adjacency matrix $\mathbf{W}$, should not affect the priority vector. Therefore, we apply L1 regularization only on the directed matrix $\mathbf{H}$. Further, we employ L2 regularization to avoid breaking the ordering induced by the priority vector through the smooth orientation matrix. Given the concern expressed by the reviewer, we further motivated our choice by plotting in Appendix F the acyclicity during training for different values of L2 regularization. Finally, the priority vector regularization aims to avoid far apart priorities which could lead to zero-gradient issues. As we discuss at the end of Section 5.3 by reporting an ablation study, priority regularization is beneficial in particular for smaller graphs.

Space of DAGs. As stated in Theorem 1 and proved in Appendix A.2, the formulation of Equation 7 covers all and only weighted DAGs.

NOBEARS References. There was indeed an error in referencing Fang et al. which we fixed in the updated manuscript. Thanks for the catch!

We thank the reviewer for acknowledging the importance of the task, appreciating the scalability of our proposal, and most importantly for their insightful comments. We hope to have addressed all the raised concerns in our following comments and in the updated manuscript in which we highlighted in red all additions and modifications.

Related Works. Given your suggestions and the comments by other reviewers we have carefully enriched the related works section, which is now available in the updated version of our manuscript. Overall, we propose a parameterization and an optimization strategy to continuously model the space of DAGs, strictly following the line of work of NOTEARS. On the other hand, both DP-DAG (Charpentier et al. 2021) and BCD-Nets (Cundy et al. 2021) are Bayesian causal discovery algorithms that learn a distribution over DAGs given observational data. Similarly, DAGuerreotype (Zantedeschi et al. 2023) retrieves a DAG by fitting a distribution over permutations. We added and discussed the necessary references to clarify the scope of our work to the readers. For this reason, while our main empirical evaluation still focuses on other optimization schemes, we also added in Appendix E.19 empirical results of DAGuerreotype on the same experimental testbed. Our results highlight that DAGuerreotype does not scale on the size of the graph as simpler optimization schemes, which could be then more easily integrated as components in other continuous causal discovery approaches.

Evaluation. As we remarked in the main body, we are aware that continuous approaches are particularly susceptible to data normalization and might exploit variance ordering between variables, as shown by Reisach et al. 2021 and originally by Loh and Bühlmann (2014). Therefore, empirical results on simulated datasets that do not explicitly control this condition might not generalize to real-world scenarios. However, we propose a loss-agnostic optimization scheme that aims at being integrated as a building block for causal applications as it is the case for NOTEARS, DAGMA, or NOCURL. By comparing on the exactly the same test bed of these methods, we show that at parity of conditions with the existing literature on acyclic optimization we get comparable results for significantly less time. While surely interesting as a future direction, defining function classes, appropriate losses, and identifiability results are characteristics of a causal discovery method, which is out of the scope of the current work. However, the suggestion by the reviewer is indeed extremely valuable to match our remarks with experimental results. Therefore, we include in Appendix E.18 a comparison for our and the remaining models on variance normalized data. Our results highlight how performance degradation is mostly due to the choice of MSE and does not depend in general on the particular optimization strategy, as COSMO is affected comparably to all competing methods.

Temperature. We agree that temperature scheduling is fundamental in our solution. In practice, as we describe in the methodological section, we simply apply cosine annealing from $t=0.45$ to a relatively small temperature during training. In practice, because of the introduced regularization, the acyclicity of the solution consistently drops before reaching zero. For a visualization of this behavior please refer to the newly introduced Appendix G, where we report the trends of acyclicity and temperature. In particular, in the last phase of training, for non-zero temperature the acyclicity already reaches zero and both the directed and the smooth orientation are still optimized in this "acyclic regime".

Thresholding. While L1 regularization already encourages sparsity, all mentioned methods prune weights with absolute values smaller than 0.3, which is a de facto standard threshold used by NOTEARS, DAGMA, NOCURL, GOLEM, and others (Ng et al. 2023). While we do not delve into this matter, we report in our experimental comparison how existing work criticizes this practice, e.g., Xu et al. (2022). Further, we emphasized in the main body the AUC metric, since it does not depend on a fixed threshold and allows us to evaluate discovery for a varying threshold. Notably, DP-DAG shares the same choice of metric.

We thank the reviewer for the additional and valuable feedback.

Evaluation. COSMO does not employ early stopping or any other algorithmic conditional based on the convergence. Therefore, the running time of our parameterization is not affected by the "triviality" of the task. We produce a fair comparison by evaluating all methods in the exact same setting in which the corresponding papers defined them, where they all should be able to perform adequately. As suggested, we already provided the required results in Appendix E.18, which show how variance-normalization equally affects ours and the competing methods and does not add further information on their comparison. We will integrate results with "sortnregress" to better clarify the limitation, which we recall being an issue of the chosen loss function and not of ours or any particular structure learning scheme.

Related Works.

1. We thank the reviewer for the comment on our related works summary, which we propose to rephrase as "In the context of differentiable approaches, of particular interest are causal discovery methods that avoid acyclicity constraints by fitting a distribution over permutations parameterized either explicitly (Cundy et al., 2021; Charpentier et al.,2022) or by a topological ordering (Zantedeschi et al., 2022)."
2. As we report in the paper, "we propose a parameterization and an optimization scheme to fastly optimize acyclic adjacencies that could be then easily integrated in causal discovery methods" and we consequently focus on alternative optimization schemes instead of many of the available and more complex DAG discovery strategies. Given the reviewer's suggestion, we included results on DAGuerreotype in Appendix E.19 to back our point and highlight how scalability is a pressing issue compared to simpler optimization strategies.

We thank the reviewer for appreciating our contribution and our empirical evaluation. We have carefully considered their insightful questions and the identified weaknesses, which we will individually address below.

Novelty. Considering the weakness identified by the reviewer, we acknowledge in the main body the similarities between ours and NOCURL parameterization. However, given the theoretical guarantees ensured by our smooth orientation matrix (Theorems 1 and 2), we are able to define an optimization scheme significantly different from NOCURL. As we comment in Appendix D.2, NOCURL proposes an algorithm that extracts the ordering from a solution obtained through NOTEARS, which in turn requires evaluating acyclicity during training. On the other hand, COSMO jointly learns the ordering, is entirely differentiable, and does not evaluate acyclicity at any optimization step.

Epsilon. The reviewer is correct that with a sufficient small $\varepsilon$, the smooth orientation matrix $S_{t,\varepsilon}$ would be a fully connected matrix. However the sparsity of the overall solution is controlled also by the directed matrix $\mathbf{H}$, since the overall DAG is given by the product $\mathbf{W}=\mathbf{H}*S_{t,\varepsilon}(\mathbf{p})$. Therefore, even if an orientation would be complete, we can still learn a sparse graph given the sparsity of $\mathbf{H}$. Notably, $\varepsilon$ cannot be zero to avoid cycles between variables with the same priority.

L2 Regularization. We thank the reviewer for the suggestion on thresholding the smooth orientation matrix during training. We experimented with a similar approach during preliminary experiments, where we used a hard threshold to test the difference between priorities in the forward pass, while keeping the sigmoid for the backward pass. However, we found it beneficial to avoid thresholding in both passes. As a further reference, we included in Appendix F an experiment on how acyclicity evolves during training for different values of L2 regularization.

Temperature and Acyclicity. In our methodological section we report that we perform cosine annealing of the temperature during training, but we have added a pointer to the Appendix where the choice of temperature and other hyperparameters range is discussed. Further, the matrix is guaranteed to be acyclic without thresholding, which we only perform to remove spurious edges. We included in Appendix G a visualization of the optimization process that shows the trend of acyclicity and performance metrics against the temperature and highlights how COSMO reaches an acyclic solution even for non-zero temperatures.

Further Additions. We added in Appendix C.4 a table reporting the hyperparameter ranges for each baseline, including DAGMA. Also, we have updated the references according to yours and other reviewers comments, please refer to the extended discussion on related works in the updated draft.

We thank the reviewer for recognizing the novelty of our contribution but also for their constructive criticism, which led to a significant refinement of the paper and to the addition of a novel appendix section to the manuscript. We will now punctually address each concern.

Figure 3. The goal of the experiment reported in Figure 3 is to test the average training epoch time of each method. To this end, we limited the training of each method to a small and finite and small number of epochs and computed the average duration. The results are therefore not meaningful in terms of accuracy or performance.

Proof of Theorem 3. We thank the reviewer for the remark. There was indeed a typo which we have corrected in the rebuttal of the manuscript.

Contribution Claim. In Appendix D.2 we discuss the differences between our parameterization and the one used by NOCURL. The section was already referenced in the related works, but thanks to the reviewer suggestion we also added a pointer from the methodological section. There, we also changed the claim in “similar model [...] with a significantly different optimization scheme” to better highlight the differences. Also, we have adjusted our claim on the primacy of unconstrained optimization models according to the reviewer’s opinion.

Topology Change. Overall, we observed that the topology mostly changes at the beginning of the optimization process and at the end, when the temperature approaches zero. We also consistently notice that the model reaches acyclic solutions for non-zero temperatures. We back this claim by introducing Appendix G, where we report a visualization of the smooth orientation matrix, the directed matrix and the overall graph against the trends of temperature, dagness and NHD.

We thank the reviewer for having recognized the novelty of our method, for the insightful comments and the proposed suggestions to improve our manuscript. We hope to have addressed most of the identified weaknesses.

Figure 1. COSMO learns a directed acyclic graph as a function of a directed graph and a priority vector. The purpose of Figure 1 is to depict the parameterization that we also describe in intuitive terms in the Introduction section before summarizing our contributions. To clarify this point, in the rebuttal update of the manuscript, we remarked in the caption that the smooth acyclic orientation is a function of the priority vector and clarified the notation.

Computational Time. In terms of computational complexity, the construction of the directed acyclic graph given the priority vector $\mathbf{p}$ and the directed matrix $\mathbf{H}$ requires a quadratic number of steps in the number of nodes. Other than complexity results, we report in each table the computation time of each model on each dataset. In terms of sample complexity, we train all models on datasets containing 1000 realizations for each graph type, in close adherence with the NOTEARS experimental testbed.

Identifiability and Classic Methods. As we report in the background, the identifiability of the true underlying graph requires several assumptions on the function class and the loss function. What we propose is a general approach to optimize acyclic graphs, which our empirical results show to be faster while comparable in optimization with other continuous structure learning methods. The minimization of the MSE without further assumptions does not ensure per se particular results on the identifiability, which require a continuous structure learner, like ours, to be embedded in a more general causal discovery framework. As suggested, we remarked on this aspect in the Related Works, to better guide the reader in the context compared to classic structure learning or causal discovery methods.

Temperature and Initialization. We agree with the reviewer that the temperature and the initial conditions are indeed sensitive hyperparameters that might affect the method. In particular, we discuss the initialization strategy of the priority vector in Section 4. We introduced in Appendix G a visualization of the optimization process during temperature annealing, which highlights how COSMO converges to an acyclic solution before reaching zero temperature while continuing to optimize the acyclic matrix. Finally, concerning Figure 4, we gave priority to yours and other reviewers' comments and will consider moving it to the appendix according to the requested additions from the discussion.