A Skewness-Based Criterion for Addressing Heteroscedastic Noise in Causal Discovery

The author mainly used ER1 to evaluate the methods on HSNMs without latent confounders, which is relatively sparse.

Thank you for your feedback. We have conducted experiments on ER2 graphs without latent confounders. The results are presented in the table below, where SkewScore consistently outperforms the other methods.

Table 1: Topological order divergence with different Dimensionalities: ER2

I think the justification of Equation (6) was a bit missing. Proposition 5 justified that Assumption 1 is mild, and therefore Equation (5) in the bivariate setting holds. However, there is no such explanation for Equation (6). Can you provide some concrete examples where Equation (6) is reduced to some known results, or simplified by specific function and noise classes?

Thank you for your question. The multivariate extension of Proposition 5 is natural, and we provided the discussion in Lines 778-787 of the manuscript. To enhance clarity, we revise the previous discussion and provide the following explanation:

Similar to Proposition 5, the set of functions $(f_i,\sigma_i)$ that violate Eq. (6) has a property analogous to having Lebesgue zero measure. Specifically, we parameterize the causal model by parameterizing functions $f_i^\theta$ and $\sigma_i^\theta$, where $\theta$ belongs to the parameter space $\Theta \subset \mathbb{R}^K$. The left-hand side of Eq. (6) defines a function of the parameter $\theta$, and we write it as $F(\theta)$. With this parameterization, the model that violates Eq. (6) corresponds to the solution of the equation $F(\theta)=0$. Under regular assumptions on $f_i$, $\sigma_i$, $p_X$, and $p_{N_i}$, $F$ is a real-analytic function. Then by applying (Mityagin, 2015, Proposition 0), the violated set $S_{\text{nc}}(\Theta) = \\{\theta \in \Theta : F(\theta) = 0\\}$ has zero measure as long as the function $F$ is not a zero constant function. It is worth noting, however, that in the homoscedastic Gaussian linear model, $F$ is identically zero, and thus (Mityagin, 2015, Proposition 0) no longer applies. We have included this discussion in Lines 857-866 in the revision. Thanks again for your suggestion.

We provide an explicit example involving 3 nodes. The model is given by: $X = f_ 2(Z) + N_ 1$, $Y = f_ 1(Z) + f_ 3(X) + N_ 2$, where $Z,N_ 1, N_ 2$ are independent Gaussian variables with distribution $\mathcal{N}(0,1)$, and $f_ i(x) = a_ ix^2+b_ ix+c_ i$ are quadratic functions. In this case, Equation (6) can be simplified into: $a_ 2 b_ 2 + a_ 3 b_ 3\neq 0$, $a_ 3\neq 0$, which are non-restrictive and can be easily satisfied.

Which ER graph model is used in the experiments, $G(d, p)$ or $G(d, m)$?

$G(d, m)$ is used. Thanks for the question.

Can you explain in more detail why $c \neq 0$ in general in Remark 6?

Thank you for your question. As discussed in the first response to Reviewer 8vJA, a non-constant $\sigma(x)$ will make the summation term in Equation (*) more complicated and "non-degenerated". Intuitively, the complicated interactions between the link functions in the summation term are less likely to cancel out unless the model is designed artificially. Specifically, when $p_ X$ and $\sigma(x)$ are both symmetric function and their symmetric axis-es coincide, then $c=0$.

More rigorously, if we parameterize $\sigma$ as $\sigma(\theta)$, where $\theta$ belongs to the parameter space $\Theta$. Then similar to the proof of proposition 5, under smooth assumption, we could prove that $\\{\theta:c(\theta) = 0\\}$ has 0 Lebesgue measure. This illustrates that $c(\theta)\neq 0$ generally.

What is the runtime for the order learning and independent tests separately?

Thank you for your question. The runtime for order learning and independent tests separately is provided below. As shown in Table 2, the runtime for both order learning and DAG recovery increases with the number of variables $d$. However, the rate of increase for DAG recovery is higher compared to order learning. This suggests that while both phases are affected by dimensionality, DAG recovery tends to become more computationally expensive in higher-dimensional settings.

Table 2: Wall-clock runtime (in seconds) of SkewScore with 1000 samples. The experiment is run on 12 CPU cores with 24 GB RAM.

Dear Reviewer iLTQ,

Thank you for your time and effort in reviewing our paper. We have carefully addressed your concerns in the rebuttal and would greatly appreciate it if you could review our response. If you have further questions or comments after reading the rebuttal, we hope to have the opportunity to address them.

Thanks again for your time and consideration.

Best regards, Authors of Paper 9319

The generality of Assumption 1, which defines the applicable HSNMs, could be explored further, although the paper presents sufficient results for a conference submission.

Thank you for your suggestion. We further discuss the generality of Assumption 1 as follows. The score function of HSNMs given by Eq. (1) in the manuscript is given by (provided Eq. (14), Line 1063)

c(a,b)= -a^3 - \frac{3a}{4} - \frac{3a}{4} e^{-b} e^{\frac{a^2}{2}} (1 + 3a^2) - \frac{9a^3}{4} e^{-2b} e^{2a^2} + \frac{15a^3}{8} e^{-3b} e^{\frac{9a^2}{2}}.

How is Theorem 1 applied to establish the causal direction between variables $X$ and $Y$ in the experiments? Equations (4) and (5) represent the population mean, while in experiments, we can only obtain the sample mean.

Thank you for your question. Theorem 1 is to establish identifiability results in the infinite sample limit, which serves as a foundational step towards establishing guarantees on finite samples. We acknowledge that in practice, we can only compute the sample mean as an estimator of the population mean. Under standard regularity conditions, the sample mean is a consistent estimator of the population mean; as the sample size increases, the sample mean converges to the true population mean. Therefore, we use the sample estimates of skewness as proxies for their population counterparts in applying Theorem 1.

To assess the impact of using sample estimates, we conducted experiments examining the effect of sample size on our method’s performance (see Table 5 in Appendix J of our manuscript). The results show that our empirical skewness estimates become more accurate with larger sample sizes. Notably, even with smaller sample sizes, our method maintains strong performance, suggesting robustness to finite-sample variations.

Our identifiability results in the infinite sample limit demonstrate that the causal direction can be identified under appropriate assumptions. Extending our theoretical findings to finite samples is an important avenue for future work, including extending the results of [Zhu et al., 2024] to non-Gaussian scenarios. Additionally, to further enhance robustness in finite samples, future work will explore hypothesis testing approaches for detecting non-zero skewness. Since we do not assume a specific noise distribution, non-parametric methods such as bootstrapping and permutation tests can be employed to assess the significance of our skewness estimates.

In summary, although Theorem 1 is based on population means, we apply it in practice using sample estimates, and our experiments indicate that our method is effective even with finite samples.

The performance of SkewScore in terms of accuracy, as shown in Figure 3, is not consistently superior to other methods. Specifically, SkewScore outperforms others in the case of GP-sig Student's t in Figure 3(a), and in the cases of GP-sig Gaussian and GP-sig Student's t in Figure 3(b). What are the computation times associated with the different methods?

Thank you for your feedback. We would like to emphasize that SkewScore consistently achieves at least 95% accuracy across all eight settings in Figure 3, whereas other baselines drop below 85% in at least one setting. This highlights the robustness of SkewScore, which we have clarified in the revised manuscript by including detailed numbers in Appendix J for completeness. Additionally, while SkewScore shows comparable or superior accuracy in the bivariate cases presented in Figure 3, it provides a significant advantage in more complex multivariate tasks, as demonstrated in Figure 4.

Regarding computation times, we have summarized them in the following table for clarity. While SkewScore’s runtime is higher than HECI, it is still faster than LOCI, providing a balanced trade-off between efficiency and performance. SkewScore’s superior accuracy and robustness make it especially valuable for multivariate settings. Additionally, it is worth noting that methods like LOCI and HECI are specifically designed for two-variable scenarios, and extending them to multivariate cases presents considerable challenges.

Table: Computation times (in seconds) for different methods in the pairwise settings presented in the paper. Results are shown for GP-Sig (Gaussian) and GP-Sig (Student's t), both with and without latent confounders. Experiments are run on 12 CPU cores with 24 GB RAM.

References:

Zhu, Z., Locatello, F., & Cevher, V. (2024). Sample complexity bounds for score-matching: causal discovery and generative modeling. Advances in Neural Information Processing Systems, 36.

After obtaining the sample version of $*SkewScore*_ {x_ j}(p)$, we identify the variable $x_ j$ with the minimum estimated $*SkewScore*_ {x_ j}(p)$ to be the leaf node (sink node), as detailed in Line 6 of Algorithm 1 in the manuscript. This leaf node represents the variable that is most likely to be causally downstream, having the least skewness in the score's projection. We then remove this leaf node from the set of nodes (Lines 8-9 of Algorithm 1) and repeat the process for the remaining variables. By iteratively removing the node with the smallest $*SkewScore*$, we establish a topological order that reflects the causal direction among the variables.

To further clarify, we do not determine the causal direction based on whether the sample version of $*SkewScore*_ {x_ j}(p)$ is larger or smaller than a specific threshold. Instead, we employ the above process, which provides a systematic determination of the causal order. Thank you for your question.

Dear Reviewer 8P9r,

Thank you for your efforts in reviewing this paper. We have provided responses to your further comments. Could you please check whether the responses properly addressed your concerns, or if you have further comments? Your feedback would be appreciated. Thanks a lot for your time.

Best regards,

Authors of Paper 9319

Fact 2 is incorrect. Let $X$ follow a univariate exponential distribution with parameter $\lambda$, then $\mathbb{E}\left[\frac{d \log p(X)}{dx}\right] = -\lambda$. Its proof is also incorrect.

Fact 2 is correct under the assumption $p_ X(x)\in C^1(\mathbb{R}^d,\mathbb{R}_ +)$ stated in Line 150 of the original manuscript, where $C^1(\mathbb{R}^d, \mathbb{R}_ +)$ is the space of continuously differentiable functions from $\mathbb{R}^d$ to the positive real numbers $\mathbb{R}_ +$ (Line 117). The exponential distribution does not satisfy this assumption due to the singularity at $x=0$, where $p(0-)=0$ and $p(0+)=1$.

The proof of Fact 2 is based on the Newton-Leibniz formula. In the last step of the proof, since $p$ is continuously differentiable, by the Newton-Leibniz formula: $\int_ {\mathbb{R}}\frac{\partial}{\partial x_ i}p(x_ 1,x_ 2\cdots,x_ d) d x_ i = p(x_ 1,x_ 2\cdots,x_ d) |^{x_ i=+\infty}_ {x_ i = -\infty} = 0-0 = 0,$ where $\lim_ {x_ i\rightarrow+\infty}p(x_ 1,x_ 2\cdots,x_ d) =\lim_ {x_ i\rightarrow -\infty}p(x_ 1,x_ 2\cdots,x_ d) =0$ because $p$ is a probability density.

We also note that the assumption $p_ X(x)\in C^1(\mathbb{R}^d,\mathbb{R}_ +)$ can be relaxed. For example, Fact 2 holds for uniform distribution on $[-1,1]$ despite having singularities at $-1$ and $1$, as the densities at these two endpoints cancel out after the Newton-Leibniz step.

We have included these derivations and discussion in Appendix D of the revised manuscript. We have also added your example to help readers understand Fact 2 more comprehensively. Thanks for the question.

The proof of Theorem 1 heavily relies on the assumption of symmetric noise. Without this symmetry, skewness no longer provides useful implications for the causal direction, as the proof of Theorem 1 breaks down.

Thank you for highlighting the reliance of Theorem 1 on the assumption of symmetric noise. We would like to emphasize the following two points.

(1) We believe that the symmetric noise assumption is a reasonable generalization of the Gaussian noise assumption utilized in previous work on multivariate heteroscedastic causal models (Duong & Nguyen, 2023). Symmetric noise distributions cover a wide range of commonly used noise types, including Gaussian and Student’s t-distributions, as mentioned in Lines 139-142 of our manuscript.

(2) To evaluate the robustness of SkewScore under slightly skewed noise distributions, we conducted experiments using Gumbel noise in the bivariate heteroscedastic noise model (see Appendix J). The results indicate that while there is a decrease in accuracy, SkewScore’s performance degradation is less significant compared to HOST (Duong & Nguyen, 2023). This suggests that SkewScore may retain effectiveness even when the noise deviates from perfect symmetry.

We agree that extending our theoretical framework to accommodate non-symmetric noise distributions is both important and achievable. We will discuss this further in our response to your next question.

Is it possible to provide identification results such as "the skew score is larger under the wrong causal direction" without the symmetric noise assumption?

Thank you for raising this insightful question. Providing identification results without the symmetric noise assumption is indeed more challenging, but it is possible under certain conditions. Specifically, if the asymmetry introduced by the link functions $f$ (conditional mean) and $\sigma$ (conditional standard deviation) is sufficiently strong compared to the asymmetry of the exogenous noise, similar identification results can be achieved. One potential approach is to quantify the asymmetry of the exogenous noise using measures like skewness or the SkewScore. By doing so, we can establish sufficient conditions on the link functions that account for the noise’s skewness. This would likely involve developing an inequality version of Assumption 1 that incorporates both the properties of the link functions and the skewness of the exogenous noise. We are optimistic that with further research, we can relax the symmetric noise assumption and provide more general identification results.

Definition 1 needs detailed discussions. As stated, "this measure captures the asymmetry...," it would be helpful to demonstrate that it is zero when the distribution is symmetric.

Thank you for your suggestion. We have included additional discussions in Appendix A in the revision. We demonstrate that the SkewScore measure is zero when the conditional density is symmetric. Additionally, we provide two examples calculating the SkewScore for skewed distributions. Along with the original manuscript's discussion on the motivation for considering the SkewScore measure, we believe these enhancements will improve the reader's understanding of this definition.

Reference:

Duong, B., & Nguyen, T. (2023). Heteroscedastic Causal Structure Learning. In ECAI 2023 (pp. 598-605). IOS Press.

1. Thank you for this suggestion. At the same time, we would like to point out that $f\in C^1(\mathbb{R}^d, \mathbb{R}_ +)$ inherently implies that the $f$ is strictly positive on $\mathbb{R}^d$, because it means that $f$ maps every point in the domain $\mathbb{R}^d$ to the codomain $\mathbb{R}_ +=(0,\infty)$. The definitions of $\mathbb{R}_ +$ and $C^1(\mathbb{R}^d, \mathbb{R}_ +)$ are outlined in Lines 116-118 of the manuscript. For the exponential distribution example, we can state that $f_ {exp} \in C^1([0,\infty), \mathbb{R}_ +)$ or $f_ {exp} \in \textrm{piecewise} \ C^1(\mathbb{R}, \mathbb{R}_ + \cup \\{0\\})$, but $f_ {exp} \notin C^1(\mathbb{R}, \mathbb{R}_+)$. However, in light of your concern and to avoid any confusion, we have emphasized that "the density function is strictly positive on $\mathbb{R}^d$" in Line 151 of the revision. Thanks again for this comment, which helps improve the clarity.

2. Thank you for your suggestion. It is a good idea to check whether the model satisfies our assumption by checking if the smallest SkewScore is close to zero. Potential approaches for this verification are non-parametric hypothesis testing methods such as bootstrapping and permutation tests, since we do not assume a specific type of noise (e.g., Gaussian). We have included the discussion on checking the model assumption in Lines 277-280 of the revision. While we acknowledge the value of incorporating such a verification step, we believe this challenge is inherent to most causal discovery results and does not represent a major weakness or one specific to our paper.

Correctness guarantees for the multi-variable case

Thank you for you question. The guarantee for multivariate settings is provided in Corollary 2 (Lines 256-265 of the manuscript), which is an extension of Theorem 1.

Latent Confounders in multivariate settings: incorrect causal graphs. For $X_ 1 \to X_ 2 \to X_ 3$ with a latent confounder $L$ influencing both $X_ 1$ and $X_ 3$ ($X_ 1 \leftarrow L \to X_ 3$), the method would infer the structure $X_ 1 \to X_ 3$ along with a spurious direct link $X_ 1 \to X_ 3$.

Thank you for your insightful question on the latent confounders in multivariate settings. Our current procedure first infers the topological ordering with the SkewScore criterion, and then uses conditional independence tests to further remove edges and identify the DAG. We agree that the way we perform conditional independence tests will lead to spurious direct links when there are latent confounders in the multivariate setting. A potential way to mitigate this issue is to leverage the idea from FCI, which will involve more sophisticated ways of performing conditional independence tests (after obtaining the topological ordering using our SkewScore criterion) and output the Partial Ancestral Graph (PAG) instead of just a DAG. Note that this is an interesting direction and is beyond the scope of our work, which is left for future investigation. We have included this discussion in the revision (Lines 356-360) and hope it addresses your concern.

Scalability; complexity scales polynomially or exponentially with the number of variables.

Thank you for your question. The computational complexity is $O(d^2)$ with respect to the number of variables $d$. SkewScore can be divided into two phases: the causal order recovery with the SkewScore criterion and the DAG recovery with independence tests. The first phase includes $d$ iterations that perform SSM, where SSM scales to high dimensional data with Adam optimizer (Song et al., 2020). The second phase includes $O(d^2)$ conditional independence tests.

We acknowledge that further exploration is required to fully extend SkewScore to high-dimensional cases. Future work will focus on improving its scalability to enhance its applicability.

Relation between (1) constraints-based algorithms such as the PC algorithm, and (2) continuous optimization-based algorithms such as the one in https://arxiv.org/abs/1803.01422, which provides scalability?

(1) The PC algorithm runs in exponential time in $d$ (the number of variables) in the worst case, since it performs the CI tests for all subsets of variables in the skeleton recovery step. SkewScore reduces this computation complexity by first estimating the topological order, which constrains the DAG to be a sub-graph of a certain fully connected DAG, and then performs only $O(d^2)$ CI tests. The PC algorithm imposes fewer assumptions than SkewScore, as it does not rely on a functional causal model, and thus it can only identify the underlying structure up to Markov equivalence class but not the complete DAG.

(2) Continuous optimization based methods typically also make certain assumptions on the functional causal model, such as NOTEARS that assumes a linear structural equation model with equal noise variances, NOTEARS-MLP that assumes additive noise models with equal noise variances, and GraN-DAG that assumes additive noise models with Gaussian noises. Although these methods are relatively computationally efficient, recent studies have shown that they may exploit certain type of artifacts of simulated data (Reisach et al., 2021) and may be susceptible to the nonconvexity issue (Ng et al., 2024), thus limiting their applicability.

Summary: The PC algorithm imposes fewer assumptions but incurs higher computational costs and only identifies the Markov equivalence class. Continuous optimization based methods are computationally efficient but have several practice issues. SkewScore provides a balanced tradeoff, combining moderate computational complexity with broader applicability.

Small sample sizes

Thanks for your question. We conducted the experiment on the effect of sample size (see Table 5 in Appendix J of the manuscript). Our method becomes more accurate with a larger sample size. Even with smaller sample sizes, SkewScore shows strong performance, maintaining a reliable level of accuracy.

On the theoretical side, the asymptotic normality of the SSM estimator is established under appropriate assumptions in Theorem 3 of [Song et al., 2020]. However, we acknowledge that deriving nonasymptotic guarantees under general assumptions that cover heteroscedastic causal models remains an open challenge. A limitation is the Lipschitzness of the score function assumed in previous work [Zhu et al., 2024], which restricts the family of causal models SkewScore can address. Future work will aim to close this gap by extending existing results, such as those in [Zhu et al., 2024], to accommodate more general, non-Gaussian scenarios.

Confidence intervals or uncertainty quantification in measuring the skewness or determining causal directionality

Yes, we can perform hypothesis testing for skewness $\neq 0$. Since we do not assume a specific type of noise (e.g., Gaussian), non-parametric methods such as bootstrapping and permutation tests could be employed. Thank you for your question.

References:

Song, Y., Garg, S., Shi, J., & Ermon, S. (2020, August). Sliced score matching: A scalable approach to density and score estimation. In Uncertainty in Artificial Intelligence (pp. 574-584). PMLR.

Spirtes, P., & Glymour, C. (1991). An algorithm for fast recovery of sparse causal graphs. Social science computer review, 9(1), 62-72.

Zheng, X., Aragam, B., Ravikumar, P. K., & Xing, E. P. (2018). Dags with no tears: Continuous optimization for structure learning. Advances in neural information processing systems, 31.

Spirtes, P., Glymour, C., & Scheines, R. (2001). Causation, prediction, and search. MIT press.

Reisach, A., Seiler, C., & Weichwald, S. (2021). Beware of the simulated dag! causal discovery benchmarks may be easy to game. Advances in Neural Information Processing Systems, 34, 27772-27784.

Ng, I., Huang, B., & Zhang, K. (2024, March). Structure learning with continuous optimization: A sober look and beyond. In Causal Learning and Reasoning (pp. 71-105). PMLR.

Dear Reviewer h4sn,

Thank you for your time and effort in reviewing our paper. We have carefully addressed your concerns in the rebuttal and would greatly appreciate it if you could review our response. If you have further questions or comments after reading the rebuttal, we hope to have the opportunity to address them.

Thanks again for your time and consideration.

Best regards, Authors of Paper 9319