Venue: NeurIPS 2022 Conference

Reviews:

final_score: rating: 5: Borderline accept: Technically solid paper where reasons to accept outweigh reasons to reject, e.g., limited evaluation. Please use sparingly.
review_content: summary: The authors proposed a new DAG constraints based on the log-determinant of the adjacency matrices. Compare to the original polynomial based DAG constraints, the new DAG constraints do have some good properties and it achieves better performance on ER4 and SF4 graphs with node from 200 to 1000.
strengths_and_weaknesses: Pros:

1. The proof and derivation of the new DAG constraint is very clear.
2. The illustration of the behaviour of h_det is very good.


Cons:

1. The properties (ii) in Theorem 1 is actually as bad as the Exponential or Polynomial based DAG constraint. See [1,2] for more details. With  this property you will not be able to obtain an DAG unless you let mu in Algorithm 1 converge to zero. 

2. Corollary 2 may not be correct and Corollary 1 maybe meaningless. For Corollary 2, it only considers h_det, without consider the score function. By [1,2] without attain the limit point it would not be possible to obtain a DAG, and this will require infinite many steps of optimisation.

3. In the experiment part, the hyper parameter for GOLEM may not be quite right. I have tried GOLEM on SF{1,2,3,4}, and ER{1,2,3}, for numbers of nodes in {10, 20, 50, 100}. The performance of GOLEM performance is consistently better than NOTEARS in a large margin, which is quite different from the result reported in the paper. Also for all cases it would be good if the authors can provide results on graphs with nodes in  {10, 20, 50, 100}. This is because with the same number of expected edges, the more the nodes, the sparser the graph is. Also it would be good if the authors can provide results on ER2, ER3 and SF2, SF3, SF4 with all number of nodes from 10 to 1000. This can provide a fair comparison for different methods in different situations. I have also attached the GOLEM code I have used, you only need to replace the HTorch part with the exponential DAG constraints from NOTEARS.



[1] Wei, Dennis, Tian Gao, and Yue Yu. "DAGs with No Fears: A closer look at continuous optimization for learning Bayesian networks." Advances in Neural Information Processing Systems 33 (2020): 3895-3906.
[2] Ng, Ignavier, et al. "On the convergence of continuous constrained optimization for structure learning." International Conference on Artificial Intelligence and Statistics. PMLR, 2022.


The GOLEM code I have used.


"""

    import numpy as np
    import torch
    import torch.nn as nn
    from .dag import HTorch
    h_torch = HTorch.apply
    class GolemEVModel(nn.Module):
        """
        Set up the objective function of GOLEM.
        Hyperparameters
        (1) GOLEM-NV: equal_variances=False, lambda_1=2e-3, lambda_2=5.0.
        (2) GOLEM-EV: equal_variances=True, lambda_1=2e-2, lambda_2=5.0.
        """
        def __init__(self, d, lambda_1, lambda_2, eps=1e-6, h_type='exponential'):
            """
            Initialize self.
            Parameters
            d: int
                Number of nodes.
            lambda_1: float
                Coefficient of L1 penalty.
            lambda_2: float
                Coefficient of DAG penalty.
            equal_variances: bool
                Whether to assume equal noise variances
                for likelibood objective. Default: True.
            """
            super().__init__()

            self.d = d
            self.lambda_1 = lambda_1
            self.lambda_2 = lambda_2
            self.eps = eps
            self.h_type = h_type
            self._B = nn.Parameter(torch.zeros(self.d, self.d))

        def forward(self, cov_emp):
            # Placeholders and variables
            self.cov_emp = cov_emp
            self.B = self._preprocess(self._B)

            # Likelihood, penalty terms and score
            self.likelihood = self._compute_likelihood()
            self.L1_penalty = self._compute_L1_penalty()
            self.h = self._compute_h()
            self.score = (self.likelihood + self.lambda_1 * self.L1_penalty +
                      self.lambda_2 * self.h)

        def _preprocess(self, B):
            """
            Set the diagonals of B to zero.
            Parameters
            B: tf.Tensor
                [d, d] weighted matrix.

            Return
            torch.Tensor: [d, d] weighted matrix.
            """
            return (1. - torch.eye(self.d)) * B
    
        def _compute_likelihood(self):
            """
            Compute (negative log) likelihood in the linear Gaussian case.
            Return
            torch.Tensor: Likelihood term (scalar-valued).
            """
            I = torch.eye(self.d)
            return 0.5 * self.d * torch.log(
                torch.trace((I - self.B).T @ self.cov_emp
                        @ (I - self.B))) - torch.linalg.slogdet(I - self.B)[1]

        def _compute_L1_penalty(self):
            """
            Compute L1 penalty.
            Return
            tf.Tensor: L1 penalty term (scalar-valued).
            """
            return torch.norm(self.B, p=1)

        def _compute_h(self):
            """
            Compute DAG penalty.

            Return
            torch.Tensor: DAG penalty term (scalar-valued).
            """
            return h_torch(self.B * self.B, self.h_type, self.eps)
    def golem_ev(X,
             lambda_1=2e-2,
             lambda_2=5.0,
             learning_rate=1e-3,
             num_iter=2e+4,
             graph_thres=0.3,
             eps=1e-6,
             h_type='exponential'):
        n, d = X.shape
        cov_emp = np.cov(X.T, bias=True)
        cov_emp = torch.Tensor(cov_emp)
        model = GolemEVModel(d, lambda_1, lambda_2, eps, h_type)
        train_op = torch.optim.Adam(model.parameters(), lr=learning_rate)

        for i in range(int(num_iter)):
            model(cov_emp)
            score, likelihood, h, B_est = model.score, model.likelihood, model.h, model.B
            loss = score
            train_op.zero_grad()
            loss.backward()
            train_op.step()
   
        return B_est.detach().numpy()
"""
questions: Is the performance from different implementation? In NOTEARS, the authors does not decrease mu but increase the coefficients of h(W dot W). In GOLEM, the authors use fixed coefficient for L1. If put the DAG constraint in the same optimisation framework as NOTEARS, or if we put the NOTEARS in the same framework as the work, what is the performance?
limitations: It is discussed.
review_title: null

final_score: rating: 7: Accept: Technically solid paper, with high impact on at least one sub-area, or moderate-to-high impact on more than one areas, with good-to-excellent evaluation, resources, reproducibility, and no unaddressed ethical considerations.
review_content: summary: The paper introduces a new penalty for enforcing acyclicity of a weighted adjacency matrix for learning DAGs with differentiable optimization. The penalty is zero if and only if the adjacency matrix is acyclic. In contrast to previous proposals based on matrix polynomials, the introduced penalty is based on the log-determinant and is connected to the theory of M-matrices. The paper studies the properties of the new penalty and compares the penalty with previous proposals. From theoretical analysis and empirical observations, it is argued that the new penalty is more suited for learning DAGs / SEMs. Through numerical experiments, it is shown that the new penalty leads to speedups and higher precisions compared to previous proposals. 
strengths_and_weaknesses: ### Strengths
1. The paper makes progress on recent developments for differentiable learning of DAGs. 
2. The proposed penalty seems novel and has interesting connections to M-matrices. 
3. The paper compares the new penalty with previous ones in the literature from various aspects. 
4. The paper is very well-written. 

### Weaknesses
1. Several aspects in the comparison are based on heuristic arguments or empirical observations, which weakens the claim that the new penalty should be preferred over the earlier ones. A theoretical analysis is desirable to show that the advantage is real or "significant" in typical scenarios or when the problem size grows large. 


questions: 1. On discounting long cycles.

Argument (i) claims that the new penalty does not discount a long cycle in the same way some previous proposals do. From Figure 2, it looks persuasive that this could be a big advantage. I would like to see some asymptotic analysis that makes this comparison precise. 

2. On vanishing gradients. 

Lemma 5 seems insufficient to show that the new penalty does not suffer from vanishing gradients (when the other penalties do). Is there a more compelling argument?

3. On computational speedups. 

Despite the same computational complexity, argument (iii) curiously claims that the new penalty "can be computed in about an order of magnitude faster" than the other two in practice. Is there an explanation?

4. Minor issues

(1) line 42: missing power in the expression?

(2) Lemma 1: I would rather state (i) as "DAGs are in the interior of $\mathbb{W}^s$".
limitations: I do not foresee any potential negative societal impact of their work
review_title: null

final_score: rating: 7: Accept: Technically solid paper, with high impact on at least one sub-area, or moderate-to-high impact on more than one areas, with good-to-excellent evaluation, resources, reproducibility, and no unaddressed ethical considerations.
review_content: summary: * This paper considers the problem of learning a structure (DAG) over a set of observed variables with gradient-based methods. The approach designs a new continuous optimization problem capturing the acyclic property of DAGs, which defines the domain of the log-det function to be M-matrices. 
* The authors provide a detailed technical description of the approach as well as extensive experiments. 
* The authors describe technically that the proposed approach has more well behaved gradients (perhaps indicative of the improved performance empirically) as well as several other theoretical properties of the approach. 
* Empirically, the authors provide extensive experiments which demonstrate the efficiency and effectiveness of the proposed approach. 
strengths_and_weaknesses: Overall I believe that this is a strong paper with an interesting approach for what is becoming a more well studied area (continuous optimization approaches for structure learning). 

In particular strengths include: 

* **Technical Contribution** - I believe the authors have clearly explained the problem they are solving, why their methodical approach solves the problem, and what are the challenges and insights that are required in the contribution. I find the contribution to be interesting and meaningful. 
* **Empirical benefits** - The authors have demonstrated the performance characteristics of their approach with extensive experiments, which would allow a practitioner to consider such results and then determine which methods are most suitable for the application (based on metrics such as time, SHD, TPR, etc). And most often, DAGMA would be a top selection, especially in terms of efficiency compared to methods like NOTEARS. 

I have a few questions which indicate places for improvement regarding setup of experiments and measurements, please refer to questions section for detailed description of these weaknesses. 


questions: * Rather than sampling a DAG structure by ordering the nodes in an undirected graph, why not sample DAGs from something like Price's model? I realize that your experiment setup is following earlier work however. I am curious what if any properties we should be aware of for your generative model of dags. How do in-degree / out-degree distributions look like? How do they differ from degree distribution from undirected graph?
* When / why would quality $Q(\cdot)$ make sense / not-make sense for an evaluation measure of the DAGs? Could it be evaluated on held-out data of the same underlying distribution? 
* When / would it make sense to evaluate not just wall clock time, but number of samples $n$? 
* It seems like hyperparameter tuning could play a large role in the performance of these approaches. Do I understand correctly that previous work has tuned hyperparameters for these same exact datasets in the cases of re-use? Is there clarity on how all hyperparemters were selected? I wonder what you think about evaluating methods using bayesian hyperparameter optimization to find, for each dataset (e.g., setting of generative model hyperparameters) best setting of hyperparameters. 
* There is a note about the proposed approach being better at detecting large cycles. For give me if I have missed something, but was there analysis in the experiments section, which noted how much of the improved performance (e.g. SHD) was due to better recovery of such large cycles? 
* If these methods were used as part of some downstream approach, do you think the empirical analysis here is predictive of performance? Why/why not / what other measurements would be interesting to consider?
limitations: Yes, the authors have done a very nice job with this.
review_title: null

final_score: rating: 8: Strong Accept: Technically strong paper, with novel ideas, excellent impact on at least one area, or high-to-excellent impact on multiple areas, with excellent evaluation, resources, and reproducibility, and no unaddressed ethical considerations.
review_content: summary: This paper proposes a new formulation promoting acyclicity. Going beyond the widely used approaches based on the trace of polynomials of the hadamard square of the adjacency matrix, this paper rather takes a different angle by considering a reduced set of matrices, M-matrices. This M-matrix based acyclicity regularizer is thoroughly investigated so that its empirically superior behavior is also understood theoretically.
strengths_and_weaknesses: Strengths
* A novel approach based on M-matrix is proposed for an acyclicity regularizer free from optimization hassle in exiting approaches.
* The behavior of the log-det acyclicity regularizer -- the gradient and other regularization properties -- is analyzed to a quite detailed extent with supporting simulation results.
* DAG discovery with the proposed method outperforms others even with the advantage in the run-time.

Weaknesses
* For some theoretical results, it would be better if some technical detail is provided at a level to help to understand. 
  * In line 177, $(4) \Rightarrow (3)$ required a bit of thought to understand so some brief comments would help instantaneous understanding.
  * In Lemma 3, the statement that $\otimes$ is the Kronecker product may be added so that reader can look up its definition.
* It seems that sometimes theorems assumes slightly different definition for $\mathbb{W}^2$ with/without nonnegative entry condition (details below)
questions: * In Corollary 1, what does it mean by 'point towards the interior'. Is this term 'interior' is the topological term or a general term for somewhere inside? It seems that $\mathbb{W}^s$ is an open subset of $R^{d \times d}$. Does this mean that its direction is always directly toward a stationary point?

* It seems that in eq(5) the condition that $W \ge 0$, that is, all entries are nonnegative is missing. It seems that some Theorem may hold without this nonnegative entry condition but others do not. Theorems/Lemma critically relying on Proposition 1 should require the nonnegative entry condition. However, for example, Lemma 2 seems to consider the cases without nonnegative entry conditions, otherwise, $sign(W_{i,j})$ is always nonnegative. It would be better if the authors clarify this.


limitations: * Except for general limitations of differentiable approaches for DAG discovery, there is no other specific limitation. To discuss such such general limitations some other references can be added. In Differentiable Causal Discovery Under Unmeasured Confounding; Rohit Bhattacharya, Tushar Nagarajan, Daniel Malinsky, Ilya Shpitser; AISTATS 2021. there are many interesting discussions and some solutions on such general limitations.
  * Differentiable approach for ADMG for semi-Markovian case
  * How to find MEC not a single DAG
review_title: null

Meta-review: metareview: Overall, reviews for this paper are quite positive. The paper presents an interesting and effective new approach to incorporating a DAG constraint into an optimization problem by using a characterization of DAGs in terms of the logdet function.

During discussion, the reviewers raised several important questions/points for clarification, which the authors largely addressed in their responses. I encourage the authors to use these responses to guide editing of the paper for the final version.

There was some disagreement during author-reviewer discussion in regards to a comparison to GOLEM made by one of the reviewers. After the discussion period, the reviewer has provided some useful information on possible reasons for inconsistencies. I hope that the authors will investigate these points carefully and update empirical results/discussions as needed in the final version. 

From the reviewer:

I carefully compared my version of GOLEM with the version of GOLEM used by the author (released with the GOLEM paper), and there are several differences.
	1	My version is implemented in PyTorch, theirs is implemented in TensorFlow,
	2	In my version, a learning rate scheduler is used to apply smaller and smaller learning rates to solve the problem as it approaches a local optimum, and in theirs, a fixed learning rate is used.

Considering this, I think it is reasonable that the authors did not observe the same performance as I did. Here the learning rate scheduler might play an important role since I also observed that with some large learning rate GOLEM may not converge so that a fixed learning rate may finally converge to a bad solution, or fail to converge.
…
I think the paper can be further enhanced if the authors can replace the DAG constraint in GOLEM with theirs to obtain a new algorithm. From my experience, it is highly possible that with a proper optimization algorithm it can achieve far better performance than the current version.
confidence: Certain
Decision: Accept