Markov Equivalence and Consistency in Differentiable Structure Learning

We would like to express our gratitude to the reviewer for their time, effort, and valuable suggestions. We would like to take this opportunity to address all the concerns raised in the reviews.

The method can be regarded as a direct combination of several works. For example, the log-likelihood loss was also proposed in GOLEM [43]. The MCD penalty was first proposed in [72].

Thanks for reviewer’s questions on our contribution. We would like to emphasize that our paper is not merely a combination of these previous works.

(1) New Identifiability Result: We utilize our log-likelihood loss with a sparsity regularizer (MCP) to prove an identifiability result under common assumptions from the optimization literature. Unlike GOLEM [1], which has specific assumptions (such as the Triangle Assumption [1]) and relies on a linear model, our method does not require these conditions. Moreover, our proof is purely based on the properties of the score function and does not rely on any assumptions about the underlying graph structure, unlike GOLEM. Additionally, as far as we know, this is the first time MCP has been introduced in the continuous DAG learning literature. Previously, $\ell_1$ was the commonly used sparsity regularizer, which led to biased estimation of the underlying parameters and could negatively affect the recovered structure. By utilizing quasi-MCP, we overcome these drawbacks of $\ell_1$ and achieve unbiased estimation. Moreover, quasi-MCP also introduces significant technical challenges that our analysis handles (see also Common Concern (3)).

(2) Scale-invariant property: We address a common concern for many DAG learning algorithms regarding scale invariance, wherein the algorithm may output different DAGs if the data is rescaled. To address this, we demonstrate that the log-likelihood score is scale-invariant under Gaussian noise. By showing this result, we aim to correct the common misunderstanding that the algorithm itself is affected by the scale of the data. Instead, it is the use of an incorrect score function that is sensitive to the scale.

We believe these points distinguish our contributions from previous works, and are significant.

When applied to nonlinear models, the author used $\frac{1}{2n}\sum_{i=1}^d \log (\\|x_i -\hat{f}_i(x)\\|^2)$ as the log-likelihood loss. This seems only valid for homogeneous Gaussian noise where the variance in the denominator can be canceled out. However, in the experimental setting, the noise variance is heteroscedastic, making the likelihood improper.

The log-likelihood for nonlinear models seems improper due to the heteroscedastic noise. Either the setting or the likelihood should be corrected before rebooting the experiments.

Thank you for this important question! We are happy to clarify the misunderstanding and will include the full derivations in the paper for clarity.

When applied to a nonlinear model with heteroscedastic noise, the negative log-likelihood is:

$$\frac{1}{2n}\sum\_{i=1}^d \log (\\|x_i -\hat{f}_i(x)\\|^2)$$

However, when the noise is homogeneous, the negative log-likelihood becomes:

$$\frac{d}{2}\log \left(\sum\_{i=1}^d \\|x_i - \hat{f}_i(x)\\|^2\right)$$

Although these formulations appear similar, they are fundamentally different. With heteroscedastic noise, the loss is the log of the sum of least squares. In contrast, with homogeneous noise, the loss is the sum of the logs of least squares. The full derivation can be found in [1] (which covers both heteroscedastic and homogeneous noise in the linear case, with $\hat{f}_i(x) = B_i^\top x^{(i)}$) and [2] (which uses heteroscedastic noise in their setting).

Thanks for this question! We address this point later in the Question “The log-likelihood….before rebooting the experiments.”

Real-world datasets can be added to validate the proposed method.

It would be interesting to compare the proposed method with the existing ones on real datasets.

We include the real-data application in attached pdf (See Common Concern (2)). These results will be updated into our paper.

Code is not available.

Thank you for bringing this up, the code will be released in the camera ready.

The conclusion section is missing.

Thanks for pointing it out. Please refer to Common Concern (4) for conclusion and limitation.

Section D.3.2, time spent on the simulation. The plots are about SHD rather than about the running time.

Thanks for bring up this typo. We put the wrong name in the title and y-axis, it should be title: “$d$ vs Time for methods”, y-axis: "Time (second)”. The value in the plot is correct. This will be update in the paper correspondingly.

[72] in the reference should be cited correctly.

Thanks; this reference will be updated in the final version.

Can you provide further examples/illustrations on whether BIC or $\ell_1$ penalty can recover minimal models? I think this will strengthen the proposed method when faithfulness fails.

Thanks for the suggestion. Please refer to Common Concern (2) for more examples.

[1] Ng, Ignavier, AmirEmad Ghassami, and Kun Zhang. "On the role of sparsity and dag constraints for learning linear dags." Advances in Neural Information Processing Systems 33 (2020): 17943-17954.

[2] Bühlmann, Peter, Jonas Peters, and Jan Ernest. "CAM: Causal additive models, high-dimensional order search and penalized regression." (2014): 2526-2556.

We sincerely thank the reviewer for the insightful critiques and comprehensive understanding of our work, and for providing such useful feedback. We will try our best to address the reviewer’s concern.

Comparison to previous work

We thank the reviewer for highlighting this related work; we will certainly update our paper to cite it and provide a discussion. Although our work and [1] prove similar results on identifiability through optimizing the log-likelihood with sparsity regularization, there are important differences. Most importantly, the use of $\ell_0$ significantly simplifies the analysis (similar identifiability results as Theorem 1, 2, and 4 are easily obtained, see Common Concern (3) for details) and ultimately does not lead to a differentiable program. Our focus is on a fully differentiable formulation. In fact, this is precisely what leads to Assumptions A and B, which are not needed when $\ell_0$-regularization is used. This is not surprising since $\ell_0$ leads to a combinatorial optimization problem which we are trying to avoid. Again, see Common Concern (3) for details.

Moreover, our results rely on slightly weaker assumptions to get identifiability results, i.e., Sparsest Markov Representation (SMR) assumption, which has been shown to be weaker than faithfulness (with only observational data, $\mathcal{I}^*$-faithfulness assumption is regular faithfulness assumption [1]). Our proof is entirely different and more straightforward, relying solely on the score function, without the consideration of different skeletons and immoralities [1]. Our paper also includes new results on scale-invariance to address current concerns with many DAG learning algorithms.

Experiments investigate only the linear case

Originally, due to space limits, we put the nonlinear experiments results in Appendix D.3.3. (Page 31), which may have been easy to miss. Nonetheless, we have added more nonlinear experiments; please see Common Concern (2).

Limitations of the theory

Thanks for these useful comments! Assumptions are needed because the our penalty term is quasi-MCP, introducing fundamental difficulties compared to previous work [1] using $\ell_0$. If penalty term is $\ell_0$, all these assumptions can be removed, and similar results hold. Please refer to Common Concern (3) for more details.

Also, can you give a concrete example of function B(\psi) that is Lipschitz.

First, in a (general) linear model, $\psi = B$. In this sense, $B(\psi) = \psi = B$, and $B(\psi)$ is 1-Lipschitz. Another example can be found in Common Concern (1). In this case, $\psi = \\{\beta^S\\}_{S\subseteq [p]}$, $[B(\psi)]\_{ji} = \sum\_{\{S:j\in S, i\in S, S\subseteq A\}}(\beta^{S})^2$, which is Lipschitz.

Moreover, [2] uses a neural network to model $f_j$ in Equation (2). Let $\psi_j$ be the weights that connect the input layer and the first hidden layer in the $j$-th neural network. In this case, $\psi = (\psi_1, \ldots, \psi_p)$. It turns out that $[B(\psi)]\_{ij} = \\|\text{i-th-column}(\psi_j)\\|_2$. As a consequence, $\\|B(\psi)\\|\_2 = \sum\_{i,j} \\|\text{i-th-column}(\psi_j)\\|_2 \leq p \sum_j \\|\psi_j\\|_2 \leq p^2 \\|\psi\\|_2$. Therefore, it is $p^2$-Lipschitz. In general, $B(\psi)$ is a function of $\psi$, and usually, the nonzero part of $\psi$ indicates that a certain entry of $B(\psi)$ is nonzero. $\ell_1$ or $\ell_2$ norms are used to ensure nonzeros are mapped to nonzeros and zeros are mapped to zeros, which makes $B(\psi)$ inherit the Lipschitz property of $\ell_1$ or $\ell_2$.

Unclear limit statements in Theorem 1 to 4

Thank you for the suggestion, and we apologize for any confusion. Indeed, our results apply to the population case, which we indicate by using $n \rightarrow \infty$. For Theorems 1 to 4, we are stating that when considering the population case (infinite number of samples), certain equality relationships hold. We will revise these statement to make this more clear.

Minor

(1) Originally, we include the plot of quasi-MCP, but it is removed due to space limit. But it seems it will be better idea to get it back again.

(2) Yes, this is a typo, it should be $\mathcal{M}(G^0) = \mathcal{G}(\mathcal{E}_{\min}(\psi^0,\xi^0))$

(3) Good suggestion! it will be fixed.

(4) Sparsest Markov representation(SMR) is not unique itself in general, but it is unique up to Markovian class.

Limitation

Thanks for bring it up. We provide conclusions and limitations in Common Concern (4).

[1] Brouillard, Philippe, et al. "Differentiable causal discovery from interventional data." Advances in Neural Information Processing Systems 33 (2020): 21865-21877.

[2] Zheng, Xun, et al. "Learning sparse nonparametric dags." International Conference on Artificial Intelligence and Statistics. Pmlr, 2020.

We would like to express our gratitude to the reviewer for acknowledging the value of our contributions and clarity of presentation. All the points addressed below will be updated in our paper.

The section on scale invariance …. or remark.

We will move these details to the appendix, and shorten the discussion in the main paper.

The authors state that the finite parameter…..mean on l297-298?

First, consider simple case of Gaussian models, and then below we discuss how to generalize. Let $X \sim \mathcal{N}(0, \Sigma)$ and $\Sigma\succ0$. It belongs to the exponential family and is full-rank. The covariance $\Sigma$ is identifiable, however, we are interested in the matrix $B$ in the SEM $X = B^\top X + N$. As in Section 4.1, we know that $\Sigma = (I - B)^{-\top} \Omega (I - B)^{-1}$. As we discuss in Section 4.2, there are at most $p!$ different $B$ that satisfy this condition, making our problem unidentifiable with respect to $B$, but with a finite number of equivalent parameters.

As for unidentifiable exponential families in [4], we can always reduce it to an identifiable model by projecting canonical statistics to a subspace. For the same example, if $\Sigma$ is singular, we can decompose $X = (X_a, X_b)$ such that $\text{Cov}(X_a)\succ 0$ and has the same rank as $\Sigma$, and $X_b$ is a linear combination of $X_a$. As discussed before, there is a finite equivalence class for $X_a$, and $X_b$ has a deterministic relationship to $X_a$.

Another explicit example (GLM with binary data) can be found in the Common Concern (1). This also extends to more general exponential families. When exponential families are full-rank, although their parameters are identifiable, the parameters $(\psi, \xi)$ as well as the SEM $B(\psi)$ may not be identifiable, which are our interest. As before, parameter is $\Sigma$, $\psi = B, \xi = \Omega$, and $\Sigma$ is function of $(B,\Omega)$. Moreover, in many cases, the size of $\mathcal{E}(\psi^0, \xi^0)$ is closely related to $p!$, which is finite since there are $p!$ topological sorts $\pi$. For each $\pi$, it is usually associated with a pair $(\psi(\pi), \xi(\pi))$ in a certain way such that $P(x; \psi(\pi), \xi(\pi)) = P(x; \psi^0, \xi^0)$. This is why we claim Assumption A(1) is reasonable.

Regarding Assumption A and B, more discussion can be found in the Common Concern (3).

The authors do not provide…scores that are useful for causal discovery?.

The point is that any likelihood can be used in conjunction with the quasi-MCP. Two examples of non-Gaussian likelihoods are logistic (binary) models and additive noise models. A third example is any model for the joint distribution with GLMs for CPDs. All told, these comprise a very rich class of models with explicitly derivable score functions. Unfortunately, we have not made this point explicit and will certainly do so in the camera ready.

For additive noise models, the log-likelihood function is $\frac{1}{2n}\sum_{i=1}^d \log (\\|x_i -\hat{f}\_i(x)\\|^2)$[1], which is used in our experiments (see Line 855 in Appendix). It is different from MSE loss $\frac{1}{2n}\sum_{i=1}^d\\|x_i - \hat{f}_i(x)\\|_2^2$ used in NOTEARS [2]. For the improvement of the experiments, it is based on the correct use of loss function and penalty, thus it leads to better performance. For other examples, please refer to Common Concern (1).

You mention that $\ell_1$ loss ... when replaced with $\ell_1$?

Using $\ell_1$ as a penalty can lead to biased parameter estimates because it generally shrinks coefficients to zero. In contrast, MCP provides unbiased estimates. Specifically, with MCP, $\mathcal{O}_{n,\lambda,\delta} = \mathcal{E}\_{\min}(\Theta)$ when population case is considered. However, this is not true for $\ell_1$.

Consider $X_1\sim \mathcal{N}(0,1)$, $X_2\sim X_1+\mathcal{N}(0,1)$, with the known order ($X_1\rightarrow X_2$), so that the log-likelihood with $\ell_1$ penalty is

$$\log((1-a)^2+1)+\lambda|a|$$

It is easy to see $a = 1$ is not minimizer. This indicates that $\ell_1$ leads to biased estimates.

Although the likelihood is "structurally"….. in terms of thresholding?

Good point! Although scaling does not affect which entries in the adjacency matrix are zero, it changes the magnitude of the entries. With a finite number of samples, this alters the thresholding we should use, which significantly affect the recovered structure in practice.

With respect to MSE loss, it depends: It is scale-invariant for a fixed topological sort, but this is not true in general when the sort is unknown. Let us clarify:

1. For fixed topo sort, the MSE in NOTEARS is scale-invariant. Since the matrix $B$ can be recovered by linear regression, and changing the scale of the data $X$ does not change the positions of the nonzero entries in $B$, leading to invariance. This is similar to the intuition of Theorem 3.
2. However, the MSE loss in NOTEARS is not scale-invariant when the sort is unknown. For instance, consider $X = (X_1,X_2)$ where $X_1 = N_1, X_2 =aX_1+ N_2, N_1, N_2\overset{i.i.d}{\sim}\mathcal{N}(0,1)$. Then the optimal for MSE loss is

$$B^* = \begin{bmatrix} 0 & a \\\\ 0 & 0 \end{bmatrix}$$

However, if you standardize $X$ to get $Z =(X_1,X_2/\sqrt{1+a^2})$ and optimal are

$$B^*_1 = \begin{bmatrix}0&a/\sqrt{a^2+1} \\\\ 0&0\end{bmatrix}\qquad B^*_2 = \begin{bmatrix}0&0 \\\\ a/\sqrt{a^2+1}&0\end{bmatrix}$$

So MSE is not scale-invariant in the sense that global optimal solution can not have the same structure for different scale of data. However, this is not the case for likelihood (Theorem 3).

Limitation and conclusion

Please refer to Common Concern (4). More discussions on Assumption A1 can be found in previous response. We will add this to the paper with the extra page in the camera ready.

[1] Bühlmann, et al "CAM: Causal additive models, high-dimensional order search and penalized regression" (2014)

[2] Zheng, et al. "Learning sparse nonparametric dags" (2020)

Thank you for the reviewer’s further comments, and we apologize for the unclear part. We are happy to address any concerns you have and ensure a clear understanding of our work. Please rest assured, we will incorporate all the suggestions into the paper, including moving Section C.2 to the main paper to provide readers with more background knowledge.

Is it the case that for a given sparsity pattern (i.e., DAG), only one $B$ setting in the Gaussian example will yield the same distribution?

You are correct! From my understanding, if a sparsity pattern here refers to a topological sort $\pi$. For each $\pi$, there will be only one $B$ that is consistent with $\pi$ and can generate the distribution.

Consider a simple three-node example. Let $X \sim \mathcal{N}(0, \Sigma)$, where $\Sigma \succ 0$:

$$\Sigma = \begin{bmatrix} 1/2&1&1\\\\1&3&3\\\\1&3&4 \end{bmatrix}$$

We would like to recover $B$ (DAG) such that $X = B^\top X + N$ and $X \sim \mathcal{N}(0, \Sigma)$. In this case, there are $p! = 3! = 6$ different topological sorts, and at most six different $B$s can generate $X$. Based on Section C2, it is worth noting that aside from these different $B$s, no other $B$ can generate the distribution $X$ [1]. These $B$s and $\pi$s are:

$$B(\pi_1) = \begin{bmatrix} 0&0&0\\\\1/3&0&0\\\\0&3/4&0 \end{bmatrix}\qquad \pi_1 = [2,1,0]$$ $$B(\pi_2) = \begin{bmatrix} 0&0&0\\\\1/3&0&1\\\\0&0&0 \end{bmatrix}\qquad \pi_2 = [1,2,0]$$ $$B(\pi_3) = \begin{bmatrix} 0&1&0\\\\0&0&0\\\\1/4&1/2&0 \end{bmatrix}\qquad \pi_3 = [2,0,1]$$ $$B(\pi_4) = \begin{bmatrix} 0&1&2\\\\0&0&0\\\\0&1/2&0 \end{bmatrix}\qquad \pi_4 = [0,2,1]$$ $$B(\pi_5) = \begin{bmatrix} 0&0&0\\\\1/3&0&1\\\\0&0&0 \end{bmatrix}\qquad \pi_5 = [1,0,2]$$ $$B(\pi_6) = \begin{bmatrix} 0&2&0\\\\0&0&1\\\\0&0&0 \end{bmatrix}\qquad \pi_6 = [0,1,2]$$

Here $B(\pi_5) = B(\pi_2)$, although $\pi_2\ne \pi_5$. This is why we say it at most $p!$ different $B$.

For unidentifiable exponential families, consider the case where $X \sim \mathcal{N}(0,\Sigma)$ and $\Sigma$ is singular. For example,

$$\Sigma = \begin{bmatrix} 1&1&1\\\\1&1&1\\\\1&1&2 \end{bmatrix}$$

We can reduce this to a lower dimension. Note that here $X_1 = X_2$, so we only need to consider $(X_2, X_3)$:

$$\text{Cov}(X_2,X_3) = \begin{bmatrix} 1&1\\\\1&2 \end{bmatrix}\succ0$$

Then it reduces to previous case.

As for more general model, for each $\pi$, it is usually associated with only one pair $(\psi(\pi), \xi(\pi))$ in a certain way such that $P(x; \psi(\pi), \xi(\pi)) = P(x; \psi^0, \xi^0)$, similar to Linear Guassian case. Another example in Common concern (1) follows such property.

I'm wondering if you mean C.2? I can't actually find this statement anywhere.

Apologies for the misunderstanding. You are correct; it should be Section C.2. We will include more details in Section 4.2 for clarity!

As stated in Section C.2, for any topological sort $\pi$, it corresponds to a pair $(\tilde{B}(\pi), \tilde{\Omega}(\pi)) \in \mathcal{E}(\Theta)$. For all DAGs with $p$ nodes, there are at most $p!$ different topological sorts. This explains why there are at most $p!$ different $B$ matrices that satisfy this condition.

[1] Aragam, Bryon, and Qing Zhou. "Concave penalized estimation of sparse Gaussian Bayesian networks." The Journal of Machine Learning Research 16.1 (2015): 2273-2328.