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

## Abstract

Abstract Real-world data often violates the equal-variance assumption (homoscedasticity), making it essential to account for heteroscedastic noise in causal discovery. In this work, we explore heteroscedastic symmetric noise models (HSNMs), where the effect Y 𝑌 Y italic_Y is modeled as Y = f ⁢ ( X ) + σ ⁢ ( X ) ⁢ N 𝑌 𝑓 𝑋 𝜎 𝑋 𝑁 Y=f(X)+\sigma(X)N italic_Y = italic_f ( italic_X ) + italic_σ ( italic_X ) italic_N , with X 𝑋 X italic_X as the cause and N 𝑁 N italic_N as independent noise following a symmetric distribution. We introduce a novel criterion for identifying HSNMs based on the skewness of the score (i.e., the gradient of the log density) of the data distribution. This criterion establishes a computationally tractable measurement that is zero in the causal direction but nonzero in the anticausal direction, enabling the causal direction discovery. We extend this skewness-based criterion to the multivariate setting and propose SkewScore , an algorithm that handles heteroscedastic noise without requiring the extraction of exogenous noise. We also conduct a case study on the robustness of SkewScore in a bivariate model with a latent confounder, providing theoretical insights into its performance. Empirical studies further validate the effectiveness of the proposed method.

## 1 Introduction

Discovering causal relationships among variables from data, known as causal discovery, is of great interest in many fields such as biology (Sachs et al.,, ) and Earth system science (Runge et al.,, ). The primary approaches to causal discovery include constraint-based methods (Spirtes and Glymour,, ; Spirtes et al.,, ), which rely on conditional independence tests, and score-based methods (Heckerman et al.,, ; Chickering,, ; Huang et al.,, ), which optimize a certain objective function. These methods often identify the causal structure only up to Markov equivalence (Spirtes et al.,, ; Glymour et al.,, ), since the causal directions in general cannot be determined without prior knowledge or additional assumptions (Pearl,, ).

Another line of approaches, based on functional causal models, impose restrictions on the causal mechanisms to make the causal directions identifiable. The key idea is to show that, under these restrictions, no forward and backward models that result in independent noise can co-exist, which leads to an asymmetry that helps identify the causal direction. Examples include linear non-Gaussian models (Shimizu et al.,, ), nonlinear additive noise models (), and post-nonlinear causal models (Zhang and Hyvärinen,, ). Most of these models assume that the noise terms are homoscedastic, i.e., they have constant variances across different samples. This assumption may limit their applicability because the noise terms are often heteroscedastic in real-world data, such as those in environmental science (Merz et al.,, ) and robotics (Kersting et al.,, ).

To address this limitation, there is growing interest in developing more general functional causal models capable of inferring causal directions despite the presence of heteroscedastic noise. Previous approaches for handling heteroscedastic noise typically involve extracting noise terms through nonlinear regressions to estimate conditional means and variances, followed by independence tests on the residuals (Strobl and Lasko,, ; Immer et al.,, ). This procedure needs to be repeated for each variable, which can be computationally intensive. Alternative strategies that avoid extracting noise terms include the likelihood-based approach, which often requires the noise to follow Gaussian distributions (Immer et al.,, ; Khemakhem et al.,, ; Duong and Nguyen,, ), and the kernel-based criterion proposed by Mitrovic et al., (), whose theoretical guarantees for handling heteroscedastic noise may not be clear.

In this work, we introduce a novel criterion to address symmetric heteroscedastic noise in causal discovery. This is motivated by the fact that symmetric noise encompasses a wide range of widely used assumptions in noise distributions, such as Gaussian, centered uniform, Laplace, and Student’s t distributions, with applications in, e.g., finance (Anderson and Arnold,, ), environmental science (Damsleth and El-Shaarawi,, ), and psychology (McGill,, ).
Specifically, we leverage the skewness of the score (the gradient of the log density) of the data distribution to establish a measurement that is zero in the causal direction but positive in the anticausal direction, thereby uncovering the causal direction.
We extend this skewness-based criterion to the multivariate setting and propose SkewScore, an algorithm that effectively manages heteroscedastic noise without needing to extract the exogenous noise. This approach adopts the two-phase ordering-based search framework for DAGs introduced by Teyssier and Koller, (). First, we estimate the topological order by iteratively identifying the sink node of the causal graph using the skewness-based criterion. Then, we prune the directed acyclic graph (DAG) associated with this topological order using classical methods, such as conditional independence tests (Zhang et al.,, ). This reduces the computational complexity, lowering the number of conditional independence tests from exponential to polynomial in the number of vertices.

Furthermore, we conduct a case study to evaluate the robustness of SkewScore in a bivariate model with a latent confounder. Existing approaches based on functional causal models typically do not allow latent confounders, as their presence will violate the independent noise condition required for causal direction identification. However, our theoretical insights show that in this bivariate model with a latent confounder, when the causal effect between the observed variables dominates that of the latent confounder, SkewScore can still correctly identify the causal direction. Additionally, our method enables the quantification of the latent confounder’s impact within this triangular model.

## 2 Related Work

Heteroscedastic Noise Models (HNMs). Causal discovery with HNMs has gained significant interest in recent years.
HNMs relax the assumption of homoscedastic noise, which was generally assumed in previous studies.
Xu et al., () adopted a modeling choice where the variance of the scaled noise variable is a piecewise constant function of its parents, which may restrict the possible mechanisms of the variance.
Tagasovska et al., (); Strobl and Lasko, (); Immer et al., () considered a more general class of HNMs, where the conditional variance is a deterministic function of its parents. The works by Tagasovska et al., () and Immer et al., () focus primarily on identifying pairwise causal relationships. Both Strobl and Lasko, () and Duong and Nguyen, () include a noise-extraction step to estimate the sink node, with Strobl and Lasko, () using mutual information and Duong and Nguyen, () using normality as criteria.
Yin et al., () propose a two-phase continuous optimization framework to learn the causal graph.

Latent Confounder. Causal discovery aims to identify causal relationships from observational data. However, traditional methods typically assume the absence of latent confounders in the causal graph, an assumption that often does not hold in real-world scenarios. To address this challenge, extensive research has focused on learning causal structures that allow latent variables. These approaches include methods based on conditional independence tests (Spirtes et al.,, ; Colombo et al.,, ; Akbari et al.,, ), over-complete ICA-based techniques (; Salehkaleybar et al.,, ), Tetrad condition (Silva et al.,, ; Kummerfeld and Ramsey,, ), high-order moments (Shimizu et al.,, ; Cai et al.,, ; Xie et al.,, ; Adams et al.,, ; Chen et al.,, ), matrix decomposition techniques (Anandkumar et al.,, ), mixture oracles (Kivva et al.,, ) and rank constraints (Huang et al.,, ; Dong et al.,, ). Nevertheless, these methods either recover the causal graph only up to a broad equivalence class or assume that each latent confounder has at least two pure children that are not adjacent, which forbids the existence of triangle structure involving latent confounders, that is $X\leftarrow Z\rightarrow Y$ and $X\rightarrow Y$, where $Z$ is the latent variable.

Causal Discovery with Score Matching. A recent class of causal discovery algorithms employs constraints on the score, i.e., the gradient of the log density function, to uniquely identify the causal graph in additive noise models (ANMs). This line of thought determines the topological order by iteratively identifying sink nodes (i.e., leaf nodes) using the score (Rolland et al.,, ; ; Sanchez et al.,, ; ). Interestingly, empirically indicates that the presence of latent confounders may not significantly disrupt the inference of the topological ordering; score-matching-based approaches still provide reliable ordering. Effective gradient estimation techniques are crucial in connecting the gradient of the log-likelihood to successful causal discovery. Various methods have been proposed to estimate the score from samples generated by an unknown distribution, including score matching (Hyvärinen and Dayan,, ; Vincent,, ; Song and Ermon,, ; Song et al.,, ), and kernel score estimators based on Stein’s methods (Li and Turner,, ; Shi et al.,, ; Zhou et al.,, ). These methods have been successfully applied to generative modeling, representation learning, and addressing intractability in approximate inference algorithms.

## 3 Problem Setup and Preliminaries

Symbols and notations. Let $\mathbb{R}_{+}=(0,\infty)$. We denote by $C^{1}(\mathbb{R}^{k})$ the set of all continuously differentiable functions on $\mathbb{R}^{d}$, and by $C^{1}(\mathbb{R}^{k},\mathbb{R}_{+})$ the subset of $C^{1}(\mathbb{R}^{k})$ consisting of functions that take only positive values. We denote by $L^{\infty}(\mathbb{R}^{k})$ the set of all essentially bounded functions on $\mathbb{R}^{k}$.
Furthermore, let $L^{\infty}:=\cup_{k=1}^{\infty}L^{\infty}(\mathbb{R}^{k})$, and $C^{1}=\cup_{k=1}^{\infty}C^{1}(\mathbb{R}^{k})$. For brevity, we write $\frac{\partial}{\partial x_{i}}$ as $\partial_{x_{i}}$. For a random variable $X=(X_{1},X_{2},\ldots,X_{d})$ in $\mathbb{R}^{d}$, we represent $(X_{1},\ldots,X_{i-1},X_{i+1},\ldots,X_{d})$ as ${X_{-i}}$, and similarly, $(x_{1},\ldots,x_{i-1},x_{i+1},\ldots,x_{d})$ as $x_{-i}$. We use $p_{X}(x)$, or simply $p(x)$, to denote the probability density of $X$, and $p_{X_{i}|X_{-i}}(x_{i}|x_{-i})$ as $p(x_{i}|x_{-i})$ for conditional densities. We denote by $\mathrm{pa}(X_{i})$ the set of parents (direct causes) of $X_{i}$ in a directed acyclic graph (DAG), and by $\mathrm{child}(X_{i})$ the set of children (direct effects) of $X_{i}$. A node with no children, i.e., no outgoing edges, is called a sink or leaf node.

### 3.1 Heteroscedastic Symmetric Noise Models (HSNM)

We consider a structural causal model over $d$ random variables $\{X_{i}\}_{i=1}^{d}$, where each variable $X_{i}$ corresponds to a node in the directed acyclic graph (DAG) $G$. The Heteroscedastic Symmetric Noise Model (HSNM) is defined as

$$ $X_{i}=f_{i}(\mathrm{pa}(X_{i}))+\sigma_{i}(\mathrm{pa}(X_{i}))N_{i},\quad i=1, \dots,d,$ (1) $$

where $N_{i}$ is an exogenous noise variable that follows a symmetric distribution with mean zero, i.e., $p_{N_{i}}(-n_{i})=p_{N_{i}}(n_{i})$. The noise terms $\{N_{i}\}_{i=1}^{d}$ are jointly independent, and each $N_{i}$ is independent of $\mathrm{pa}(X_{i})$ when $\mathrm{pa}(X_{i})\neq\emptyset$.
We assume that $p_{X}\in C^{1}(\mathbb{R}^{d},\mathbb{R}_{+})$, and for each $i\in\{1,2,\dots,d\}$, $p_{N_{i}}\in C^{1}(\mathbb{R},\mathbb{R}_{+})$, with $f_{i},\sigma_{i}\in C^{1}$ and $\sigma_{i}\in L^{\infty}$. Additionally, all partial derivatives $\partial_{x_{j}}\sigma_{i}\in L^{\infty}$. We also assume there exists a constant $r>0$ such that $\inf_{i\in\{1,\dots,d\}}\inf_{x}\sigma_{i}(x)\geq r$. Furthermore, we assume that $p_{X}$, $p_{N}$, $f_{i}$, and $\sigma_{i}$ satisfy weak regularity conditions to ensure that the derivatives and expectations involved in our analysis are well-defined.

It is worth noting that symmetric noise encompasses a wide range of widely used assumptions in noise distributions, such as Gaussian, centered uniform, Laplace, and Student’s t distributions, with applications in, e.g., finance (Anderson and Arnold,, ), environmental science (Damsleth and El-Shaarawi,, ), and psychology (McGill,, ).

### 3.2 Score Function, Score Matching, and Skewness

In this paper, the *score* function refers to the gradient of the log density, $\nabla p(x)$, which can be estimated using various score matching methods Hyvärinen and Dayan, (); Vincent, (); Song and Ermon, (); Song et al., (); Li and Turner, (); Shi et al., (); Zhou et al., (). We introduce two important properties of the score function and define the skewness of the score.

#### Facts

Let $p(x_{1},x_{2},\cdots,x_{d})\in C^{1}(\mathbb{R}^{d},\mathbb{R}_{+})$ be a density function. The following two facts hold for every $i=1,2\cdots,d$:

###### Fact 1 .

$\partial_{x_{i}}\log p(x_{i}|x_{1},\cdots,x_{i-1},x_{i+1},\cdots,x_{d})=
\partial_{x_{i}}\log p(x_{1},x_{2},\cdots,x_{d}).$

###### Fact 2 .

$\mathbb{E}_{X\sim p}\left[\partial_{x_{i}}\log p(X)\right]=0.$

The proofs of Fact 1 and Fact 2 are provided in Appendix D. These two facts are crucial for deriving Theorem 1, a skewness-based criterion for HSNMs. Fact 1 allows us to work with derivatives of the log joint density rather than directly estimating conditional densities, while Fact 2 shows that the expectation of the score is a zero vector.

Skewness measures the asymmetry of a probability distribution. In particular, Pearson’s moment coefficient of skewness is defined as the normalized third central moment. For a real-valued random variable $W$, the unnormalized skewness is given by $\textit{Skew}(W)=\mathbb{E}\left[\left(W-\mathbb{E}[W]\right)^{3}\right]$.
Inspired by this, we introduce the notion of skewness for the score function, noting that by Fact 2, the score function has zero mean.

###### Definition 1 .

Given a probability density $p$ for a random variable $X=(X_{1},X_{2}\cdots,X_{d})\in\mathbb{R}^{d}$, we define its skewness of the score in the $x_{i}$-coordinate as:

$$ $\textit{SkewScore}_{x_{i}}(p):=\left|\mathbb{E}_{X\sim p}\left[\left(\partial_ {x_{i}}\log p(X)\right)^{3}\right]\right|.$ (2) $$

This measure captures the asymmetry within each coordinate of the score, which plays a key role in identifying causal directions in HSNMs.

###### Fact 1 .

###### Fact 2 .

###### Definition 1 .

## 4 Causal Discovery with the Skewness-Based Criterion

In this section, we focus on identifying causal relationships within heteroscedastic symmetric noise models (HSNM), as formalized in Eq. (1). We introduce a novel identification criterion, the skewness-of-score criterion, which captures the asymmetry between the causal and anti-causal directions to determine the causal order in the infinite sample limit. We begin by providing intuitive explanations of this criterion, present theoretical guarantees for the bivariate causal model, and extend the results to multivariate cases. Finally, we present an algorithm based on this criterion in Section 4.2.

### 4.1 Identification with the Skewness-Based Criterion

We begin with the bivariate case to build intuition for the skewness-based criterion. The extension to multivariate models follows naturally, as shown in Corollary 2.

Consider the following bivariate heteroscedastic symmetric noise model (HSNM):

$$ $Y=f(X)+\sigma(X)N,$ (3) $$

where $N$ is symmetric exogenous noise independent of $X$ with mean zero, i.e., $p_{n}(-n)=p_{n}(n)$, and $X\perp\!\!\!\perp N$.
We assume the same regularity conditions stated in Section 3.1, specified as Assumption E and E in Appendix E.

The key insight for the identifiability of the heteroscedastic symmetric noise model (HSNM), as given in Eq. (3), is that the conditional distribution $P(Y\mid X)$ is symmetric, while $P(X\mid Y)$ is generally asymmetric under mild assumptions on $f$ and $\sigma$. This contrast is illustrated in Figure 1. To quantify the asymmetry of these conditional distributions, a typical measure is the skewness, defined as the standardized central third moment. However, directly estimating the skewness of the conditional distributions from data is challenging due to the complexity of the conditional expectation.

To circumvent this challenge, we instead turn to an auxiliary random variable — the *score function* $\nabla\log p(X,Y)$, where $(X,Y)$ follows the data distribution $p$. Specifically, Fact 1 allows us to convert the partial derivatives of the log conditional density into the gradient of the log joint density, which is the score function. This transformation reframes the problem of analyzing the conditional distribution as one involving the score function, thereby simplifying the estimation process. Consequently, this leads to the identifiability criterion for heteroscedastic symmetric noise models, formally stated in Theorem 1, with the model restriction given in Assumption 4.1.

Figure: Figure 1: Identifiability of the heteroscedastic symmetric noise model (HSNM). For the causal direction $X\rightarrow Y$ (left), the conditional distribution $p(y|x)$ is symmetric for any $x$ (top right), and the variance of $p(y|x)$ varies with different values of $x$ due to heteroscedasticity. In contrast, $p(x|y)$ is asymmetric for some $y$ (bottom right).
Refer to caption: x1.png

- Assumption 1 (Model Restriction).
Denote $A(x)=\frac{p^{\prime}_{x}(x)}{p_{x}(x)}-\frac{\sigma^{\prime}(x)}{\sigma(x)}$, $B(x)=-\frac{f^{\prime}(x)}{\sigma(x)}$, and $C(x)=\frac{\sigma^{\prime}(x)}{\sigma(x)^{2}}$, where $p^{\prime}_{x}$, $\sigma^{\prime}$, and $f^{\prime}$ represent the derivatives of $p_{x}$, $\sigma$, and $f$ respectively. Assume the following inequality holds:
$\displaystyle\int\left[\left(A(x)+C(x)\tfrac{p_{n}^{\prime}(u)u}{p_{n}(u)}
\right)^{3}+3\left(B(x)\tfrac{p_{n}^{\prime}(u)}{p_{n}(u)}\right)^{2}\left(A(x
)+C(x)\tfrac{p_{n}^{\prime}(u)u}{p_{n}(u)}\right)\right]p_{n}(u)p_{x}(x)
\mathrm{d}x\mathrm{d}u\neq 0.$

Assumption 4.1 excludes the scenario where the skewnesses of the score’s projections are both zero in the causal and anti-causal directions. In particular, it rules out linear Gaussian additive models. In Proposition 5 in Appendix B, we demonstrate that Assumption 4.1 holds generically by showing that the set of models violating this assumption has a property analogous to having a zero Lebesgue measure. More discussion is provided in Appendix B.

We now present a key result that shows how skewness of the score of the data distribution provides useful information that helps uncover the causal direction.

###### Theorem 1 .

Let $(X,Y)$ follows the heteroscedastic symmetric noise model given by Eq. (3).
Let $p(x,y)$ denote the joint distribution of $(X,Y)$. The following asymmetry holds.

- 1.
The score function’s component in the $Y$-coordinate (the effect direction) follows an unskewed distribution:
$\textit{SkewScore}_{y}(p):=\left|\mathbb{E}_{X,Y}\left[\left(\partial_{y}\log p
(X,Y)\right)^{3}\right]\right|=0.$
(4)
- 2.
Moreover, under Assumption 4.1, the score function’s component in the $X$-coordinate (the cause direction) follows a skewed distribution:
$\textit{SkewScore}_{x}(p):=\left|\mathbb{E}_{X,Y}\left[\left(\partial_{x}\log p
(X,Y)\right)^{3}\right]\right|\neq 0.$
(5)

The proof is provided in Appendix E. The asymmetry in the above theorem can be used to determine the causal direction between $X$ and $Y$. Specifically, one could estimate the score function $\nabla\log p(x,y)$ and determine the causal direction based on the magnitude of the skewness of its $x$- and $y$-coordinate.

Theorem 1 shows that in bivariate HSNMs, only the score’s component only along the sink-node coordinate is unskewed, under weak assumptions. This result naturally generalizes to the multivariate case for $d\geq 2$, as stated in the following corollary.

We first clarify the notations: Let $X=(X_{1},X_{2},\cdots X_{d})\in\mathbb{R}^{d}$ follow the multivariate heteroscedastic symmetric noise model given by Eq. (1), and $p(x)=p(x_{1},x_{2},\cdots,x_{d})$ denote the joint distribution of $X$. For simplicity, denote $\mathrm{d}x=\mathrm{d}x_{1}\mathrm{d}x_{2}\cdots\mathrm{d}x_{d}$, $f_{i}=f_{i}(\mathrm{pa}_{i}(x))$, $\sigma_{i}=\sigma_{i}(\mathrm{pa}_{i}(x))$, $p_{i}=p_{N_{i}}\left(\frac{x_{i}-f_{i}}{\sigma_{i}}\right)$ and $p^{\prime}_{i}=p^{\prime}_{N_{i}}\left(\frac{x_{i}-f_{i}}{\sigma_{i}}\right)$.

###### Corollary 2 .

Assume that the following condition holds for every node $k$ with direct causes (i.e., $\mathrm{child}(k)\neq\emptyset$):

$$ $\int_{\mathbb{R}^{d}}\bigg{[}\frac{p^{\prime}_{k}}{\sigma_{k}p_{k}}-\sum_{i\in \mathrm{child}(k)}\left(\frac{p^{\prime}_{i}}{p_{i}}\cdot\frac{\sigma_{i} \partial_{x_{k}}f_{i}+(x_{i}-f_{i})\partial_{x_{k}}\sigma_{i}}{\sigma_{i}^{2}} +\frac{\partial_{x_{k}}\sigma_{i}}{\sigma_{i}}\right)\bigg{]}^{3}\prod_{j=1}^{ d}\frac{p_{j}}{\sigma_{j}}\mathrm{d}x\neq 0.$ (6) $$

Then, for every node $j$, it is a sink if and only if the score’s component in the $x_{j}$-coordinate follows an unskewed distribution, i.e.,

$$ $\textit{SkewScore}_{x_{j}}(p)=0,$ $$

where $\textit{SkewScore}_{x_{j}}(p)$ is defined in Definition 1.

The proof is provided in Appendix F. This corollary enables us to find the sink nodes of a DAG with the knowledge of the score function. By iteratively applying this process, we can determine the topological order of the graph. In Section 4.2, we provide a practical algorithm that leverages this Corollary for identifying the causal order.

###### Theorem 1 .

###### Corollary 2 .

### 4.2 Algorithm

We apply Corollary 2 and propose an algorithm to determine the causal ordering, named Skewscore (Algorithm 1). Our approach follows the topological order search framework of Rolland et al., (), which determines the topological order by iteratively identifying sink nodes. Different from their method, we identify the node with the minimum skewness of score as the sink node in each iteration. Compared to existing methods, Skewscore circumvents the need to extract the exogenous noise (Immer et al.,, ; Duong and Nguyen,, ; ) or to estimate the diagonal of the Hessian of the log density (Rolland et al.,, ; Sanchez et al.,, ).

Once the topological order is obtained using Algorithm 1, the complete DAG can be constructed by adding edges with conditional independence tests (Zhang et al.,, ) associated with the topological order. In this order-based search, $\mathcal{O}(d^{2})$ CI tests are performed, compared to exponential complexity in classical methods (Spirtes and Glymour,, ) in worst-case scenarios.

Figure: Algorithm 1 SkewScore

## 5 Case Study: Robustness to the Latent Confounder in a Bivariate Setting

It has been shown that various types of assumption violations may affect the performance of causal discovery (Reisach et al.,, ; ; Ng et al.,, ).
In this section, we conduct a case study to analyze the robustness of Algorithm 1 to the latent confounder in a pairwise causal model given by Eq. (7). We provide theoretical insights into this robustness in Proposition 3, with experimental results presented in Section 6.

We focus on a three-variable model consisting of two observed variables, $X$ and $Y$, and one latent variable, $Z$, as depicted in Figure 2. The model is formally expressed as:

$$ $\displaystyle X$ $\displaystyle=\phi_{0}(Z)+N_{0},$ (7) $\displaystyle Y$ $\displaystyle=\lambda\tilde{f}(X)+\phi_{1}(Z)+\sigma(X)N_{1},$ $$

where $N_{i}\sim q_{i}$ with $q_{i}(n)=q_{i}(-n)$ and $q_{i}\in C^{1}(\mathbb{R},\mathbb{R}_{+})$, for $i=0,1$, and $Z\sim p_{z}(z)$. In this model, the effect of the latent confounder has an additive structure, and the conditional standard deviation $\sigma(X)$ only depends on the observed variable $X$, without interaction with the latent variable. The function $f$ is replaced with $\lambda\tilde{f}$, where $\lambda>0$ scales the effect of $X$ on $Y$. To analyze the confounding effects of $Z$, we vary $\lambda$ while keeping $\tilde{f},\sigma,\phi_{i},q_{i}$, and $p_{z}(z)$ fixed. As in Section 3.1, we assume $\tilde{f},\sigma\in C^{1}(\mathbb{R})$, $\sigma,\sigma^{\prime}\in L^{\infty}(\mathbb{R})$, and that there exists a constant $r>0$ such that $\sigma(x)\geq r$.

Figure: Figure 2: Three-variable HSNM model with latent confounder given by Eq. (7).

The criterion in Theorem 1 identifies the causal order without requiring exogenous-noise extraction. This simplicity allows us to quantify the impact of latent confounders. Notably, to recover the causal direction in the infinite sample limit, Line 6 in Algorithm 1 is equivalent to observing that the score’s component in the effect coordinate is less skewed than in the cause. The following proposition shows that in Model (7) with $\lambda$ large enough, the causal direction $X\rightarrow Y$ can be identified by Algorithm 1 in the infinite sample limit, with only mild regularity assumptions on the function $\tilde{f}$, $\phi_{i}$, $q_{i}$ for $i\in\{0,1\}$, and $p_{z}$.

###### Proposition 3 .

Suppose two variables $X$ and $Y$ follow the causal models in Eq.(7). Given $\tilde{f}$, $\sigma$, $\phi_{i}$, $q_{i}$ for $i\in\{0,1\}$, and $p_{z}$,
let $p_{\lambda}(x,y)$ denote the joint distribution of $(X,Y)$. Assume the regularity conditions are satisfied
(^*^**This assumption is specified by Assumption H in Appendix H.). Then there exists $\lambda^{*}\in\mathbb{R}_{+}$, such that for every $\lambda\geq\lambda^{*}$, we have

$$ $\textit{SkewScore}_{x}(p_{\lambda})>\textit{SkewScore}_{y}(p_{\lambda}),$ $$

where $\textit{SkewScore}_{x}(\cdot),\textit{SkewScore}_{y}(\cdot)$ is defined in Definition 1.
This indicates the effect $Y$ can be identified in Line 6 in Algorithm 1 in the infinite sample limit.

Mild regularity conditions on $\tilde{f}$, $\phi$, and $p_{z}$, that are sufficient for this analysis, are discussed in Appendix H, where we provide the proof.

Roughly speaking, when $\lambda$ is sufficiently large, similar to a high signal-to-noise ratio (SNR), the effect of the latent variable $Z$ can mostly be treated as noise. In this setting, $\lambda$ serves as an indicator of how much the observed variables dominate over the latent confounder. The threshold $\lambda^{*}$ is the critical value beyond which the causal influence of $X$ on $Y$ becomes strong enough that the skewness-based measure of the causal direction $X\rightarrow Y$ consistently surpasses that of the reverse direction $Y\rightarrow X$. For $\lambda<\lambda^{*}$, this skewness-based measure may not clearly distinguish the causal direction. However, once $\lambda$ exceeds $\lambda^{*}$, the log-density derivatives of $X$ show much higher skewness, allowing for the reliable identification of $Y$ as the effect variable, even when there is a latent confounder.

Finally, we present an explicit analysis of a simplified homoscedastic case of the triangular model defined by Eq. (7), where the exogenous noises and the cause follow Gaussian distributions, and the effects of the latent variable $Z$ on both $X$ and $Y$ are assumed to be linear.

###### Example 4 .

Consider the following model:

$$ $\displaystyle\begin{cases}X&=Z+N_{0},\\ Y&=Z+f(X)+N_{1},\end{cases}$ $$

where $Z,N_{i}$ for $i\in\{0,1\}$ are independent variables, all distributed as $\mathcal{N}(0,1)$. Then, we have:

$$ $\displaystyle\begin{cases}\textit{Skew}_{x}&=\frac{2\sqrt{2\pi}}{\pi}\left| \int_{\mathbb{R}}xf^{\prime}(x)(1+f^{\prime}(x))e^{-\frac{x^{2}}{2}}\mathrm{d} x\right|,\\ \textit{Skew}_{y}&=0.\end{cases}$ (8) $$

The explicit calculation of Eq. (8) can be found in Appendix I. If $f$ is a linear function, then $\textit{Skew}_{x}=0$ and the model is unidentifiable by Algorithm 1. However, in most cases, $\textit{Skew}_{x}>0$, making the model identifiable by Algorithm 1. For example, for a quadartic function $f(x)=ax^{2}+bx+c$ with $a\neq 0$, $\textit{Skew}_{x}>0$ if and only if $b\neq\frac{1}{2}$. Similarly, for almost all polynomials with a degree greater than $2$, we have $\textit{Skew}_{x}>0$, indicating that Algorithm 1 can effectively determine the causal direction in the infinite sample limit.

###### Proposition 3 .

###### Example 4 .

## 6 Experiments

Setup. We compare SkewScore to state-of-the-art causal discovery methods to demonstrate the effectiveness of our proposed methods. We conduct experiments using synthetic data generated from HSNMs in three settings: (1) bivariate HSNMs; (2) latent-confounded triangular HSNM setting (bivariate model with a latent confounder, as discussed in Section 5); and (3) multivariate HSNMs. For each scenario, noise is sampled from two types of symmetric distributions: Gaussian distribution and Student’s t distribution. The degree of freedom of the Student’s t distribution is uniformly sampled from $\{2,3,4,5\}$.

For the multivariate setting, the causal conditional mean function $f$ is modeled by a Gaussian process with a radial basis function (RBF) kernel (unit bandwidth), while the conditional standard deviation, $\sigma$, is parameterized using a sigmoid function. In both the bivariate and triangular settings, two data generation procedures are used: one with the Gaussian process-sigmoid setup (denoted as “GP-sig”), the same procedure as the multivariate setting, and another with $f$ modeled as an invertible sigmoid function and $\sigma$ set to the absolute function (denoted as “Sig-abs”). In latent-confounded triangular setting, the latent effects $\phi_{0}$ and $\phi_{1}$ are set as linear functions with slopes uniformly sampled from $(-2,-0.5)\cup(0.5,2)$. The data is then standardized. The causal graphs are constructed via the Erdös-Rényi (ER) model, where for a given number of nodes $d$, the average number of edges is also set to $d$, for $d\geq 2$. The sample size for all experiments is $5000$. The score estimator employs sliced score matching (SSM) with a $3$-layer MLP. We use a batch size of 128 and configure the hidden dimensions to 128 for $d<10$ and 512 for $d\geq 10$. Optimization is performed using the Adam optimizer with a learning rate of $10^{-3}$, and we subsample 1000 points for conditional independence testing after obtaining the topological order.

Comparisons. We empirically evaluate SkewScore against state-of-the-art bivariate causal discovery methods designed for heteroscedastic noise models: HECI (Xu et al.,, ); LOCI (Immer et al.,, ) and HOST (Duong and Nguyen,, ). Given that Since HOST can also handle multivariate settings, we include it in those comparisons. We also test DiffAN (Sanchez et al.,, ), which uses the second-order derivative of the log-likelihood, estimated through diffusion models. Additionally, we evaluate NoTEARS (Zheng et al.,, ), which uses continuous optimization for linear DAG learning, and the MLP version of DAGMA (Bello et al.,, ), which handles nonlinear DAG learning. For methods that output a DAG, we compute the corresponding topological order. For baseline methods designed only for cause-effect pairs, we perform comparisons using data generated from pairwise causal relations with $d=2$.

Metrics. For the pairwise causal discovery settings, we report the accuracy of the estimated causal directions. For multivariate cases, we compute two quantities: the topological order divergence (Rolland et al.,, ), and the structural Hamming distance (SHD) between the estimated and the ground truth graphs. The topological order divergence evaluates how well the estimated topological order aligns with the ground truth. Given an ordering $\pi$ and a binary ground-truth adjacency matrix $A$, the topological order divergence is defined as $D_{top}(\pi,A)=\sum_{i=1}^{d}\sum_{j:\pi_{i}>\pi_{j}}A_{ij}$.
Lower values of $D_{\text{top}}(\pi,A)$ indicate higher accuracy in the estimated topological order, reflecting more precise causal discovery. The structure Hamming distance (SHD) measures the difference between two directed acyclic graphs (DAGs) by counting the number of edge additions, deletions, or reversals required to transform one into the other.

Results. The results of the synthetic experiments for the bivariate, latent-confounded-triangular, and multivariate scenarios are shown in Figure 3(a), Figure 3(b) and Figure 4, respectively. For the pairwise causal discovery task, the results are averaged over 100 independent runs, while for the multivariate task, they are averaged over 10 runs. Our approach, SkewScore, consistently outperforms or matches other methods across the bivariate, latent-confounded triangular, and multivariate settings.

In pairwise causal discovery tasks (Figure 3), methods designed for heteroscedastic noise models (HNMs) generally perform well, but tend to make more errors when confounders are present or when the noise follows a Student’s *t* distribution. In contrast, methods like NoTEARS, DAGMA, and DiffAN—which are not specifically designed for HNMs—face performance limitations in both pairwise and multivariate scenarios.
SkewScore demonstrates consistently strong performance in the latent-confounded triangular setting, supported by the theoretical insights discussed in Section 5. All HNM-based methods perform well in the latent-confounded setting with the Sig-abs structural causal model formulation, where LOCI handles both Gaussian and Student’s t-distributed noise effectively, and HOST performs well for the Gaussian noise case. However, when the latent-confounded data is generated using the GP-sig formulation, these HNM baselines do not match the performance of SkewScore.
In more challenging multivariate scenarios, our method outperforms the others, especially as the dimensionality increases. HOST depends on noise extraction and normality tests, which can be difficult with complex data. DiffAN performs well with a small number of variables, but its performance drops as dimensionality grows. This is likely due to its non-retraining approach and the score derived from ANM, which may lead to error accumulation over iterations.
Additionally, all experiments are run on CPUs, though our method can easily be deployed on GPUs. This suggests the potential for further scalability by leveraging modern computational advancements. More experiments and discussions are provided in Appendix J.

Figure: (a) Bivariate heteroscedastic noise model.
Refer to caption: x4.png

Figure: Figure 4: Topological order divergence and structural Hamming distance (SHD) across different number of variables (the dimension $d$). Lower values indicate better performance for both metrics.
Refer to caption: x6.png

## 7 Conclusion

In this study, we introduced a novel criterion—the skewness-of-score criterion—for causal discovery with heteroscedastic symmetric noise models. Leveraging this criterion, we developed the algorithm SkewScore, which effectively identifies causal directions without the need for extracting exogenous noise. We also conduct a case study on the robustness of our algorithm in a triangle latent-confounded model, providing theoretical insights into its performance under this setup. Theoretical guarantees and empirical validations support our methodology.

Future work will focus on extending the method to extremely high-dimensional cases, where the number of variables and the complexity of their interactions significantly increase. Additionally, future research will explore causal discovery in causal models involving multiple latent confounders and more general latent structures.

## Acknowledgements

The research is partially supported by the NSF awards: SCALE MoDL-2134209, CCF-2112665 (TILOS), DMS-2428058.
It is also supported by the U.S. Department of Energy, the Office of Science, the Facebook Research Award, as well as CDC-RFA-FT-23-0069 from the CDC’s Center for Forecasting and Outbreak Analytics. YL thanks Sumanth Varambally and Shibo Hao for the discussion. YH thanks Xuyang Lin for the discussion.

## Appendix A Visualization of the Score Function

We provide the following illustration for score functions on data points.
Figure 1 illustrates that to capture the skewness of the residuals, we can consider the skewness of conditional distribution. The following Figure 5 illustrates how to quantify the skewness of conditional distributions $p(y|x)$ and $p(x|y)$ in Figure 1. Rather than working directly with the conditional random variables $X|Y$ and $Y|X$, which are often difficult to handle, we introduce an auxiliary random variable, $\nabla\log p(X,Y)$, where $(X,Y)$ follows the data distribution $p$. This two-dimensional variable, induced by the score function on the data distribution, is represented as the black directed line segment in Figure 5. This auxiliary random vector’s first element, $\partial_{x}\log p(X,Y)$, is visualized as the orange directed line segment, and the second element, $\partial_{y}\log p(X,Y)$, is shown in blue. Fact 1 connects the conditional distribution $p(y|x)$ and $p(x|y)$ and this auxiliary random vector, enabling us to transform the conditional-distribution-estimation problem to a joint-distribution one. Theorem 1 shows that under Assumption 4.1, the orange directed segment (representing $\partial_{x}\log p(X,Y)$) follows an asymmetric distribution, whereas the blue directed segment (representing $\partial_{y}\log p(X,Y)$) follows a symmetric distribution. We summarize the motivations in Figure 6.

Figure: Figure 5: Score function (a vector field) visualized on data points with projections on the $x$-axis and $y$-axis. For clarity, only the projections for one data point are shown (highlighted in orange and blue). Theorem 1 examines the skewness of these projections across the data distribution.
Refer to caption: x8.png

Figure: Figure 6: Conceptual flow from skewness of exogenous noise to skewness of the score.

## Appendix B Discussion on Model

Assumption 4.1

The following proposition demonstrates that Assumption 4.1 is not overly restrictive and can generally be satisfied; specifically, the set of models that violate this assumption is relatively small. To establish this, we fix $p_{x}$, $p_{n}$, and $\sigma(x)$, and consider $f$ within a finite-dimensional space with a sufficiently large dimension to capture the features of the data. We then demonstrate that the set of functions $f(x)$ that do not satisfy the assumption has zero measure.

###### Proposition 5 .

Fix $p_{x}$, $p_{n}$, and $\sigma(x)$. Consider $f$ within an $M$-dimensional linear space $V$. In this space $V$, define the non-compliant set as

$$ $S_{\mathrm{nc}}(V):=\left\{f\in V:f\text{ does not satisfy Assumption \ref{ assmp:E_partial_x_neq_0}}\right\}.$ $$

With the notation in Assumption 4.1, define

$$ $\displaystyle\begin{cases}h(x):=3\left(\int_{\mathbb{R}}\frac{p_{n}^{\prime}(u )^{2}}{p_{n}(u)}\mathrm{d}u\right)\frac{A(x)p_{x}(x)}{\sigma(x)^{2}}+3\left( \int_{\mathbb{R}}\frac{p_{n}^{\prime}(u)^{3}u}{p_{n}(u)^{2}}\mathrm{d}u\right) \frac{C(x)p_{x}(x)}{\sigma(x)^{2}},\\ c:=\int\left(A(x)+C(x)\tfrac{p_{n}^{\prime}(u)u}{p_{n}(u)}\right)^{3}p_{n}(u)p _{x}(x)\mathrm{d}x\mathrm{d}u,\end{cases}$ $$

which do not involve $f$.

Since $\mathrm{dim}V=M$, $V$ is isomorphic to $\mathbb{R}^{M}$ and thus naturally inherits the Lebesgue measure on $\mathbb{R}^{M}$. This induced measure does not depend on the choice of basis and is denoted by $m$. Assume there exists $x_{0}\in\mathbb{R}$ such that $h(x_{0})\neq 0$, then the following holds.

- 1.
If $c\neq 0$, then the non-compliant set has zero measure, i.e.,
$m(S_{\mathrm{nc}}(V))=0.$
- 2.
If $c=0$, assume $V$ is the space of polynomial functions of degree at most $M-1$, i.e., $V=\left\{\sum_{i=0}^{M-1}a_{i}x^{i}:a_{0},a_{2},\dots,a_{M-1}\in\mathbb{R}\right\}$. Assume $\int_{\mathbb{R}}\left(p_{x}(x)^{2}+p^{\prime}_{x}(x)^{2}\right)e^{\mu x^{2}}
\mathrm{d}x<\infty$, for some $\mu\in(0,+\infty)$. Then there exists $M_{0}\in\mathbb{N}$, such that for every $M\geq M_{0}$, the non-compliant set has zero measure, i.e.,
$m(S_{\mathrm{nc}}(V))=0.$

This proposition indicates that Assumption 4.1 is generally satisfied for most models within a sufficiently large finite-dimensional space, as the set of functions violating the assumption has zero measure. The proof is provided in Appendix G.

###### Remark 6 .

We note that when $\sigma(x)$ is not constant (i.e., in the heteroscedastic model), $c\neq 0$ in most cases. However, if $\sigma(x)$ is constant and $p_{x}$ is symmetric, then $c=0$. In this scenario, if $X$ is a Gaussian variable, the assumption in the second case is always satisfied; that is, there exists some $\mu\in(0,+\infty)$ such that $\int_{\mathbb{R}}\left(p_{x}(x)^{2}+p^{\prime}_{x}(x)^{2}\right)e^{\mu x^{2}}
\,\mathrm{d}x<\infty$.

For the multi-variate case, we provide the following intuition that the assumption in Corollary 2 is generally satisfied. Similar to proposition 5, the set of functions $f$ that violate Eq. (6) has a property analogous to a Lebesgue zero measure set. Specifically, we assume that $f_{i}$ is generated from the linear function space $V=\{\sum_{j=0}^{m}a_{j}g_{i}:a_{j}\in\mathbb{R}\}$, where $\{g_{i}(x_{1},x_{2},\cdots x_{d})\}_{i=1}^{m}$ are the bases. Denoting by $f_{i}=\sum_{j=0}^{m}a_{i,j}g_{i}$, then $\{f_{i}\}_{i=1}^{d}$ violating (6) is equivalent to the parameters $a=\{a_{i,j}\}_{i\in\{1,2,\cdots,d\},j\in\{1,2,\cdots,m\}}$ satisfying the equation $F(a)=0$, for a fixed multivariate function $F$. Under regular assumption of $g_{i},\sigma_{i}(\text{pa}_{i}(x))$, $p_{X}(x)$, $p_{N_{i}}(x)$, $F$ is a real-analytic function, and thus we can apply (, Proposition 0), which confirm that the solution set $\{a:F(a)=0\}$ has zero measure. It is also noteworthy that in the Gaussian linear model, such $F$ is identically zero, so the theorem from no longer applies.

We provide the following examples to give a more intuitive understanding of Assumption 4.1.

###### Example 7 .

In this example we consider a Homoscedastic model, where $\sigma(x)=\sigma$ is a constant. Let $f$ be a two-layer Neural Network with Softplus activation function:

$$ $f(x)=w_{2}\log(1+e^{(w_{1}x+b_{1})})+b_{2},\ w_{1},w_{2},b_{1},b_{2}\in\mathbb {R}.$ $$

Assuming $p_{x}$ is an unimodal distribution symmetric with respect to 0 0, then Assumption 4.1 is violated if and only if $w_{1}=0$ or $w_{2}=0$, i.e., when $f$ is a constant function.

- Proof
The LHS of Assumption 4.1 equals to
$\frac{3w_{1}^{2}w_{2}^{2}}{\sigma^{2}}\int_{\mathbb{R}}\frac{1}{(1+e^{-(w_{1}x
+b)})^{2}}p^{\prime}_{x}(x)dx\int_{\mathbb{R}}\frac{(p^{\prime}_{N}(u))^{2}}{p
_{N}(u)}du.$
Notice that $\frac{(p^{\prime}_{N}(u))^{2}}{p_{N}(u)}$ is always non-negative, so $\int_{\mathbb{R}}\frac{(p^{\prime}_{N}(u))^{2}}{p_{N}(u)}du>0$. We only need to prove $\int_{\mathbb{R}}\frac{1}{(1+e^{-(w_{1}x+b)})^{2}}p^{\prime}_{x}(x)dx\neq 0$ when $w_{1}\neq 0$.
Notice that $p_{x}(x)$ is an even function thus $p^{\prime}_{x}(x)$ is an odd function, thus
$\displaystyle\int_{\mathbb{R}}\frac{p^{\prime}_{x}(x)}{(1+e^{-(w_{1}x+b)})^{2}
}dx$
$\displaystyle=\int_{0}^{\infty}\frac{p^{\prime}_{x}(x)}{(1+e^{-(w_{1}x+b)})^{2
}}p^{\prime}_{x}(x)dx+\int_{-\infty}^{0}\frac{p^{\prime}_{x}(x)}{(1+e^{-(w_{1}
x+b)})^{2}}dx$
$\displaystyle=\int_{0}^{\infty}\left(\frac{1}{(1+e^{-(w_{1}x+b)})^{2}}-\frac{1
}{(1+e^{-(-w_{1}x+b)})^{2}}\right)p^{\prime}_{x}(x)dx$
$\displaystyle=\int_{0}^{\infty}\frac{(2+e^{-(w_{1}x+b)}+e^{-(-w_{1}x+b)})e^{-b
}(e^{w_{1}x}-e^{-w_{1}x})}{(1+e^{-(w_{1}x+b)})^{2}(1+e^{-(-w_{1}x+b)})^{2}}p^{
\prime}_{x}(x)dx$
Since $w_{1}\neq 0$, and $p_{x}$ is an unimodal distribution symmetric with respect to 0 0. The integrand is always of the same sign (always strictly greater than 0 or strictly less than 0) over the integration interval $(0,\infty)$. Therefore, the above integral is non-zero.
$\blacksquare$

###### Example 8 .

In this example we consider the Heteroscedastic case with $\sigma(x)=e^{\lambda(x-\mu)^{2}}$, where $\mu\neq 0,~{}\lambda>0$. Let $f(x)=ax+b$, $X\sim\mathcal{N}(0,1)$, $N\sim\mathcal{N}(0,1)$. Then
$f$ violates Assumption 4.1 if and only if

$$ $a=\pm\sqrt{-\frac{c}{\gamma}},$ $$

where $\gamma=3\mu\lambda\bigg{[}\frac{4e^{\frac{\mu^{2}}{2(4\lambda+1)}}}{(4\lambda+
1)\sqrt{4\lambda+1}}+\frac{6e^{\frac{\mu^{2}}{2(6\lambda+1)}}}{(6\lambda+1)
\sqrt{6\lambda+1}}\bigg{]}$ and $c$ is defined in Proposition 5. Especially, when $\frac{c}{\gamma}>0$, $f$ always satisfy Assumption 4.1.

It is evident that the non-compliant set $S_{\mathrm{nc}}(V)=\{\sqrt{-\frac{c}{\gamma}}\}\times\mathbb{R}\cup\{-\sqrt{-
\frac{c}{\gamma}}\}\times\mathbb{R}$ always has 0 measure (it is an empty set when $\frac{c}{\gamma}<0$), providing an example for proposition 5. Additionally, due to the simple structure, we already have $m(S_{\mathrm{nc}}(V))=0$ when $c=0$. Therefore, the discussion of the second case in Proposition 5 is unnecessary. The proof of this example is provided as follows.

- Proof
Recall in the proof of proposition 5 (Eq. (16)), we proved that $f$ violates Assumption 4.1 equals to
$\int_{\mathbb{R}}f^{\prime}(x)^{2}h(x)\mathrm{d}x=-c.$
As $f(x)=ax+b$ is linear, it equals to
$a^{2}\int_{\mathbb{R}}h(x)\mathrm{d}x=-c.$
On the other hand, we can simplify $h(x)$ as follows:
$h(x)=\frac{3p_{X}^{\prime}(x)}{\sigma(x)^{2}}-\frac{9p_{X}(x)\sigma(x)}{\sigma
(x)^{4}}=-\frac{3x}{\sqrt{2\pi}}e^{-\frac{4\lambda(x-\mu)^{2}+x^{2}}{2}}-\frac
{18\lambda(x-\mu)}{\sqrt{2\pi}}e^{-\frac{6\lambda(x-\mu)^{2}+x^{2}}{2}}.$
Thus we can compute that
$\int_{\mathbb{R}}h(x)\mathrm{d}x=3\mu\lambda\bigg{[}\frac{4e^{\frac{\mu^{2}}{2
(4\lambda+1)}}}{(4\lambda+1)\sqrt{4\lambda+1}}+\frac{6e^{\frac{\mu^{2}}{2(6
\lambda+1)}}}{(6\lambda+1)\sqrt{6\lambda+1}}\bigg{]}.$
Then we know that $\int_{\mathbb{R}}h(x)\mathrm{d}x\neq 0$. Therefore, we conclude that $f$ violates Assumption 4.1 if and only if
$a^{2}=-\frac{c}{\int_{\mathbb{R}}h(x)\mathrm{d}x},$
which equals to
$a=\pm\sqrt{-\frac{c}{\gamma}}.$
$\blacksquare$

In the following example, we discuss a case where $f$ does not satisfy Assumption 4.1.

###### Example 9 .

Let $\sigma(x)=\sigma$ be a constant and $f$ be linear, then $f$ does not satisfy Assumption 4.1 if and only if:

$$ $\int_{\mathbb{R}}\frac{p^{\prime}_{x}(x)^{3}}{p_{x}(x)^{2}}\mathrm{d}x=0.$ (9) $$

When $X$ follows a symmetric distribution, Eq. (9) holds. Specifically, when $X$ follows a Gaussian distribution, Eq. (9) holds, implying that $f$ fails to satisfy Assumption 4.1.
This indicates that the Gaussian Linear model is unidentifiable by our algorithm, consistent with previous research such as (). However, as illustrated in example 8, if $\sigma(x)$ is replaced by a non-constant function, the model becomes identifiable for most cases. The proof of this example is provided as follows.

- Proof
We prove it by directly simplifying the formula in Assumption 4.1. Notice that when $\sigma(x)=\sigma$ is a constant, $f(x)=ax+b$ is linear, $C(x)=0$, $B(x)=-\frac{a}{\sigma}$, $A(x)=\frac{p^{\prime}_{x}(x)}{p_{x}(x)}$. So $f$ does not satisfy Assumption 4.1 equals to
$\int_{\mathbb{R}}\frac{p^{\prime}_{x}(x)^{3}}{p_{x}(x)^{2}}\mathrm{d}x+\frac{3
a^{2}}{\sigma^{2}}\int_{\mathbb{R}}\frac{p_{n}^{\prime}(u)^{2}}{p_{n}(u)}\int_
{\mathbb{R}}p^{\prime}_{x}(x)\mathrm{d}x=0.$
Notice that $\int_{\mathbb{R}}p^{\prime}_{x}(x)\mathrm{d}x=0$, it then equals to
$\int_{\mathbb{R}}\frac{p^{\prime}_{x}(x)^{3}}{p_{x}(x)^{2}}\mathrm{d}x=0,$
thus completing the proof.
$\blacksquare$

###### Proposition 5 .

###### Remark 6 .

###### Example 7 .

###### Example 8 .

###### Example 9 .

## Appendix C Discussion on the Skewness Definition

To enhance the flexibility of our algorithm, we redefine the skewness of a random variable $W$ as $Skew_{\psi}(W)=\mathbb{E}[\psi(W-\mathbb{E}W)]$, employing a nonlinear odd test function $\psi$ that satisfies $\psi(-s)=-\psi(s)$. Specifically, in Theorem 1, the chosen function is $\psi(s)=s^{3}$. Note that for any symmetric random variable $W$, $Skew_{\psi}(W)=0$ for any odd function $\psi$, thus upholding the validity of Eq. (4) under this generalized definition. To ensure that a $\psi$-analog of Eq. (5) also holds in the anti-causal direction, it is essential to verify the $\psi$-analog of Assumption 4.1. A necessary condition for $\psi$ is nonlinearity due to Fact 2.

## Appendix D Proofs of supporting lemmas and facts

- Proof of Fact 1
Without loss of generality, it suffices to provide the proof when $p\in C^{1}(\mathbb{R}^{2})$. Notice that the derivative of the log conditional density translates to the derivative of the log joint distribution, then
$\displaystyle\frac{\partial}{\partial x_{i}}\log p(x_{i}|x_{-i})$
$\displaystyle=\frac{\partial}{\partial x_{i}}\log\left(\frac{p(x_{1},x_{2},
\cdots,x_{d})}{p(x_{-i})}\right)$
$\displaystyle=\frac{\partial}{\partial x_{i}}\left(\log p(x_{1},x_{2},\cdots,x
_{d})-\log p(x_{-i})\right)$
$\displaystyle=\frac{\partial}{\partial x_{i}}\log p(x_{1},x_{2},\cdots,x_{d}).$
$\blacksquare$

- Proof of Fact 2
$\displaystyle\mathbb{E}\left[\frac{\partial}{\partial x_{i}}\log p(x_{1},x_{2}
,\cdots,x_{d})\right]$
$\displaystyle=\int_{\mathbb{R}^{d}}\frac{\frac{\partial}{\partial x_{i}}p(x_{1
},x_{2},\cdots,x_{d})}{p(x_{1},x_{2},\cdots,x_{d})}p(x_{1},x_{2},\cdots,x_{d})
\mathrm{d}x_{-i}\mathrm{d}x_{i}$
$\displaystyle=\int_{\mathbb{R}^{d-1}}\int_{\mathbb{R}}\frac{\partial}{\partial
x
_{i}}p(x_{1},x_{2},\cdots,x_{d})\mathrm{d}x_{i}\mathrm{d}x_{-i}$
$\displaystyle=0.$
$\blacksquare$

## Appendix E Proof of Theorem 1

- Assumption 2 (Regularity assumption for $p_{x},p_{n}$). $p_{n},p_{x}\in C^{1}$,
1.
$\int_{\mathbb{R}}\frac{\left|(p_{n}^{\prime}(u))^{3}\right|}{p_{n}^{2}(u)}du<\infty$,
2.
$\int_{\mathbb{R}}\left|\frac{p_{x}^{\prime}(x)^{3}}{p_{x}(x)^{2}}\right|dx<\infty$, $\int_{\mathbb{R}}\frac{\left|(p_{n}^{\prime}(u))^{3}u^{3}\right|}{p_{n}^{2}(u)
}du<\infty$.

- Assumption 3 (Regularity assumption for $f$). $f\in C^{1}$,
$\int_{\mathbb{R}}\left|f^{\prime}(x)^{3}\right|p_{x}(x)\mathrm{d}x<\infty.$

###### Remark 10 .

Assumption E and E ensures all the integrals we discuss here are well-defined.

- Proof of Theorem 1
We first prove Eq. (4). As $Y=f(X)+\sigma(X)N$, the density of $(X,Y)$ should be
$p(x,y)=p_{x}(x)p(y|x)=\frac{p_{x}(x)}{\sigma(x)}p_{n}\left(\frac{y-f(x)}{
\sigma(x)}\right).$
$\displaystyle\mathbb{E}_{X,Y}\left(\partial_{y}\log p(X,Y)\right)^{3}$
$\displaystyle=\int_{\mathbb{R}^{2}}\frac{p_{x}(x)}{\sigma(x)^{4}}\frac{\left(p
_{n}^{\prime}\left(\frac{y-f(x)}{\sigma(x)}\right)\right)^{3}}{p_{n}^{2}\left(
\frac{y-f(x)}{\sigma(x)}\right)}\mathrm{d}x\mathrm{d}y$
$\displaystyle=\int_{\mathbb{R}^{2}}\frac{p_{x}(x)}{\sigma(x)^{3}}\frac{\left(p
_{n}^{\prime}(u)\right)^{3}}{p_{n}^{2}(u)}\mathrm{d}x\mathrm{d}u$
$\displaystyle=\int_{\mathbb{R}}\frac{p_{x}(x)}{\sigma(x)^{3}}\mathrm{d}x\int_{
\mathbb{R}}\frac{(p_{n}^{\prime}(u))^{3}}{p_{n}^{2}(u)}\mathrm{d}u$
$\displaystyle=0.$
For Eq. (5), we first expand it explicitly:
$\displaystyle\mathbb{E}_{X,Y}\left(\partial_{x}\log p(X,Y)\right)^{3}$
$\displaystyle=\int_{\mathbb{R}^{2}}\bigg{[}\frac{p^{\prime}_{x}(x)}{p_{x}(x)}-
\frac{\sigma^{\prime}(x)}{\sigma(x)}-\frac{p_{n}^{\prime}(u)}{p_{n}(u)}\left(
\frac{f^{\prime}(x)}{\sigma(x)}+\frac{\sigma^{\prime}(x)u}{\sigma(x)^{2}}
\right)\bigg{]}^{3}p_{x}(x)p_{n}(u)\mathrm{d}x\mathrm{d}u$
(10)
$\displaystyle=\int_{\mathbb{R}}p_{x}(x)\int_{\mathbb{R}}\bigg{[}A(x)+B(x)\frac
{p_{n}^{\prime}(u)}{p_{n}(u)}+C(x)\frac{p_{n}^{\prime}(u)u}{p_{n}(u)}\bigg{]}^
{3}p_{n}(u)\mathrm{d}u\mathrm{d}x.$
(11)
where $A(x)=\frac{p^{\prime}_{x}(x)}{p_{x}(x)}-\frac{\sigma^{\prime}(x)}{\sigma(x)}$, $B(x)=-\frac{f^{\prime}(x)}{\sigma(x)}$ and $C(x)=\frac{\sigma^{\prime}(x)}{\sigma(x)^{2}}$, as is denoted in Assumption 4.1.
As $p_{n}(u)=p_{n}(-u)$, we derive $\int_{\mathbb{R}}\frac{p_{n}^{\prime}(u)^{3}}{p_{n}(u)^{2}}du=\int_{\mathbb{R}
}p_{n}^{\prime}(u)du=0$, so we can can expand Eq. (11) and get
$\displaystyle\mathbb{E}_{X,Y}\left(\partial_{x}\log p(X,Y)\right)^{3}$
$\displaystyle=\int_{\mathbb{R}}p_{x}(x)\int_{\mathbb{R}}\Bigg{[}\left(A(x)+C(x
)\frac{p_{n}^{\prime}(u)u}{p_{n}(u)}\right)^{3}\mathrm{d}x$
(12)
$\displaystyle~{}~{}+3\left(B(x)\frac{p_{n}^{\prime}(u)}{p_{n}(u)}\right)^{2}
\left(A(x)+C(x)\frac{p_{n}^{\prime}(u)u}{p_{n}(u)}\right)\Bigg{]}p_{n}(u)
\mathrm{d}u.$
(13)
By assumption, $\mathbb{E}_{X,Y}\left(\partial_{x}\log p(X,Y)\right)^{3}\neq 0$, thus completing the proof.
$\blacksquare$

###### Remark 10 .

## Appendix F Proof of corollary 2

- Proof of corollary 2
As in prior work, we assume all integrals below are well-defined. We first compute that: for every $k\in\{1,2,\cdots d\}$,
$\partial_{x_{k}}\log p(x)=\frac{p^{\prime}_{k}}{\sigma_{k}p_{k}}-\sum_{i\in
child
(k)}\left(\frac{p^{\prime}_{i}}{p_{i}}\cdot\frac{\sigma_{i}\partial_{x_{k}}f_{
i}+(x_{i}-f_{i})\partial_{x_{k}}\sigma_{i}}{\sigma_{i}^{2}}+\frac{\partial_{x_
{k}}\sigma_{i}}{\sigma_{i}}\right).$
(14)
It suffices to notice that the joint distribution of $X$ can be written as:
$\displaystyle p(x)$
$\displaystyle=\prod_{i=1}^{d}p\left(x_{i}|\mathrm{pa}_{i}(x)\right)$
$\displaystyle\log p(x)$
$\displaystyle=\sum_{i=1}^{d}\log p\left(x_{i}\mid\mathrm{pa}_{i}(x)\right)$
$\displaystyle=\sum_{i=1}^{d}\log\left({p_{N_{i}}\left(\frac{x_{i}-f_{i}\left(
\mathrm{pa}_{i}(x)\right)}{\sigma_{i}(\mathrm{pa}_{i}(x))}\right)}\right)-\sum
_{i=1}^{d}\log{\sigma_{i}(\mathrm{pa}_{i}(x))}.$
We then differentiate with respect to $x_{k}$ and directly derive Eq. (14).
We are now prepared to prove the first assertion. Considering a sink $X_{k}$, we know that $\text{child}(k)=\emptyset$. Therefore, by the formula in the first line of the proof:
$\partial_{x_{k}}\log p(x)=\frac{p^{\prime}_{k}}{\sigma_{k}p_{k}}.$
Then we can compute:
$\displaystyle\mathbb{E}_{X}\left(\partial_{x_{k}}\log p(X)\right)^{3}$
$\displaystyle=\int_{\mathbb{R}^{d}}\left[\frac{p^{\prime}_{N_{k}}\left(\frac{x
_{k}-f_{k}\left(\mathrm{pa}_{k}(x)\right)}{\sigma_{k}(\mathrm{pa}_{k}(x))}
\right)}{\sigma_{k}(\mathrm{pa}_{k}(x))\cdot p_{N_{k}}\left(\frac{x_{k}-f_{k}
\left(\mathrm{pa}_{k}(x)\right)}{\sigma_{k}(\mathrm{pa}_{k}(x))}\right)}\right
]^{3}\prod_{i=1}^{d}p\left(x_{i}|\mathrm{pa}_{i}(x)\right)\mathrm{d}x$
$\displaystyle=\int_{\mathbb{R}^{d}}\frac{p(\mathrm{pa}_{k}(x))}{\sigma_{k}(
\mathrm{pa}_{k}(x))^{4}}\cdot\frac{\left(p^{\prime}_{N_{k}}\left(\frac{x_{k}-f
_{k}\left(\mathrm{pa}_{k}(x)\right)}{\sigma_{k}(\mathrm{pa}_{k}(x))}\right)
\right)^{3}}{p^{2}_{N_{k}}\left(\frac{x_{k}-f_{k}\left(\mathrm{pa}_{k}(x)
\right)}{\sigma_{k}(\mathrm{pa}_{k}(x))}\right)}\prod_{i\neq k}p\left(x_{i}|
\mathrm{pa}_{i}(x)\right)\mathrm{d}x$
$\displaystyle=\int_{\mathbb{R}^{d}}\frac{p(\mathrm{pa}_{k}(x))}{\sigma_{k}(
\mathrm{pa}_{k}(x))^{3}}\cdot\frac{\left(p_{N_{k}}^{\prime}(u)\right)^{3}}{p_{
N_{k}}^{2}(u)}~{}\mathrm{d}u~{}\prod_{i\neq k}(p\left(x_{i}|\mathrm{pa}_{i}(x)
\right)\mathrm{d}x_{i})$
$\displaystyle=\int_{\mathbb{R}^{d-1}}\frac{p(\mathrm{pa}_{k}(x))}{\sigma_{k}(
\mathrm{pa}_{k}(x))^{3}}\cdot\prod_{i\neq k}(p\left(x_{i}|\mathrm{pa}_{i}(x)
\right)\mathrm{d}x_{i})\int_{\mathbb{R}}\frac{\left(p_{N_{k}}^{\prime}(u)
\right)^{3}}{p_{N_{k}}^{2}(u)}~{}\mathrm{d}u$
$\displaystyle=0.$

We obtain the last equality because $\int_{\mathbb{R}}\frac{\left(p_{N_{k}}^{\prime}(u)\right)^{3}}{p_{N_{k}}^{2}(u
)}~{}du=0$. More explicitly, due to the assumption that $p_{N_{k}}^{\prime}(u)$ is symmetric, $\frac{\left(p_{N_{k}}^{\prime}(u)\right)^{3}}{p_{N_{k}}^{2}(u)}$ is a odd function and thus its integral is 0 0.

We then prove the second assertion. Recall that

$$ $\partial_{x_{k}}\log p(x)=\frac{p^{\prime}_{k}}{\sigma_{k}p_{k}}-\sum_{i\in child (k)}\left(\frac{p^{\prime}_{i}}{p_{i}}\cdot\frac{\sigma_{i}\partial_{x_{k}}f_{ i}+(x_{i}-f_{i})\partial_{x_{k}}\sigma_{i}}{\sigma_{i}^{2}}+\frac{\partial_{x_ {k}}\sigma_{i}}{\sigma_{i}}\right).$ $$

Therefore,

$$ $\mathbb{E}_{X}\left(\partial_{x_{i}}\log p(X)\right)^{3}=\int_{\mathbb{R}^{d}} \bigg{[}\frac{p^{\prime}_{k}}{\sigma_{k}p_{k}}-\sum_{i\in child(k)}\left(\frac {p^{\prime}_{i}}{p_{i}}\cdot\frac{\sigma_{i}\partial_{x_{k}}f_{i}+(x_{i}-f_{i} )\partial_{x_{k}}\sigma_{i}}{\sigma_{i}^{2}}+\frac{\partial_{x_{k}}\sigma_{i}} {\sigma_{i}}\right)\bigg{]}^{3}\prod_{j=1}^{d}\frac{p_{j}}{\sigma_{j}}\mathrm{ d}x.$ $$

Thanks to our assumption (6), we derive that:

$$ $\mathbb{E}_{X}\left(\partial_{x_{i}}\log p(X)\right)^{3}\neq 0,$ $$

thus completing the proof.
$\blacksquare$

## Appendix G Proof of Proposition 5

###### Lemma 11 .

If $\int_{\mathbb{R}}\left(p_{x}(x)^{2}+p^{\prime}_{x}(x)^{2}\right)e^{\mu x^{2}}
\mathrm{d}x<\infty$ for some $\mu\in(0,+\infty)$, recalling

$$ $h(x):=3\left(\int_{\mathbb{R}}\frac{p_{n}^{\prime}(u)^{2}}{p_{n}(u)}\mathrm{d} u\right)\frac{A(x)p_{x}(x)}{\sigma(x)^{2}}+3\left(\int_{\mathbb{R}}\frac{p_{n} ^{\prime}(u)^{3}u}{p_{n}(u)^{2}}\mathrm{d}u\right)\frac{C(x)p_{x}(x)}{\sigma(x )^{2}},$ $$

then we have :

$$ $\int_{\mathbb{R}}h(x)^{2}e^{\mu x}\mathrm{d}x<\infty.$ $$

- Proof
As $A(x)=\frac{p^{\prime}_{x}(x)}{p_{x}(x)}-\frac{\sigma^{\prime}(x)}{\sigma(x)}$ and $C(x)=\frac{\sigma^{\prime}(x)}{\sigma(x)^{2}}$ and $\sigma,\sigma^{\prime}$ are bounded and $\sigma(x)\geq r>0$, there exists a constant $b>0$ s.t. $\left|A(x)\right|\leq\left|\frac{p^{\prime}_{x}(x)}{p_{x}(x)}\right|+b$, $\left|C(x)\right|\leq b$. Notice that $\int_{\mathbb{R}}\frac{p_{n}^{\prime}(u)^{2}}{p_{n}(u)}\mathrm{d}u$ and $\int_{\mathbb{R}}\frac{p_{n}^{\prime}(u)^{3}u}{p_{n}(u)^{2}}\mathrm{d}u$ are all constant, then there exists constant $c_{1}>0$ s.t.
$\left|h(x)\right|\leq c_{1}\left(p_{x}(x)+\left|p_{x}^{\prime}(x)\right|\right).$
Therefore
$\displaystyle\left|h(x)\right|^{2}$
$\displaystyle\leq 2c_{1}^{2}\left(p_{x}(x)^{2}+p_{x}^{\prime}(x)^{2}\right),$
$\displaystyle\implies\int_{\mathbb{R}}h(x)^{2}e^{\mu x}\mathrm{d}x$
$\displaystyle\leq\int_{\mathbb{R}}\left(p_{x}(x)^{2}+p^{\prime}_{x}(x)^{2}
\right)e^{\mu x^{2}}\mathrm{d}x<\infty.$
$\blacksquare$

Defining

$$ $L^{2}(\mathbb{R},e^{\mu x^{2}}):=\left\{h\in L^{2}(\mathbb{R}):\int_{\mathbb{R }}h(x)^{2}e^{\mu x}\mathrm{d}x<\infty\right\},$ $$

lemma 11 equals to saying that $h\in L^{2}(\mathbb{R},e^{\mu x^{2}})$. The following lemma is interesting as the integral is taken over $\mathbb{R}$ instead of a finite interval, and thus the Stone Weierstrass theorem cannot be applied.

###### Lemma 12 .

For $h\in L^{2}(\mathbb{R},e^{\mu x^{2}})$, if $\int_{\mathbb{R}}x^{n}h(x)dx=0$ for every $n\in\mathbb{N}$, then $h=0$ almost everywhere.

- Proof
Without loss of generality, we assume $\mu=1$. It is well known that $L^{2}(\mathbb{R},e^{\mu x^{2}})$ is a Hilbert space. Then by knowledge of Hermite polynomials, there exist a series of polynomials $\left\{H_{n}\right\}_{n=1}^{\infty}$ which is a orthogonal basis of $L^{2}(\mathbb{R},e^{\mu x^{2}})$. $\int_{\mathbb{R}}x^{n}h(x)dx=0$ for every $n\in\mathbb{N}$ implies that $\int_{\mathbb{R}}H_{n}(x)h(x)dx=0$ for every $n\in\mathbb{N}$, thus proving $h=0$.
$\blacksquare$

- Proof of Proposition 5
For every finite dimensional subspace $V\subset V_{0}$, $f\in S_{\text{nc}}(V)$ equals to
$\displaystyle 0=\int\left[\left(A(x)+C(x)\tfrac{p_{n}^{\prime}(u)u}{p_{n}(u)}
\right)^{3}+3\left(B(x)\tfrac{p_{n}^{\prime}(u)}{p_{n}(u)}\right)^{2}\left(A(x
)+C(x)\tfrac{p_{n}^{\prime}(u)u}{p_{n}(u)}\right)\right]p_{n}(u)p_{x}(x)
\mathrm{d}x\mathrm{d}u.$
(15)
In the above equation, only $B(x)$ involves $f(x)$. So we can simplified Eq. (15) into
$\int_{\mathbb{R}}f^{\prime}(x)^{2}h(x)\mathrm{d}x+c=0,$
(16)
where
$\displaystyle\begin{cases}h(x)=3\left(\int_{\mathbb{R}}\frac{p_{n}^{\prime}(u)
^{2}}{p_{n}(u)}\mathrm{d}u\right)\frac{A(x)p_{x}(x)}{\sigma(x)^{2}}+3\left(
\int_{\mathbb{R}}\frac{p_{n}^{\prime}(u)^{3}u}{p_{n}(u)^{2}}\mathrm{d}u\right)
\frac{C(x)p_{x}(x)}{\sigma(x)^{2}},\\
c=\int\left(A(x)+C(x)\tfrac{p_{n}^{\prime}(u)u}{p_{n}(u)}\right)^{3}p_{n}(u)p_
{x}(x)\mathrm{d}x\mathrm{d}u.\end{cases}$
Let $f_{i}$ ($1\leq i\leq d$) be a basis of $V$, $f=\sum_{i}a_{i}f_{i}$, then Eq. (16) equals to
$\sum_{1\leq i,j\leq d}a_{i}a_{j}\int_{\mathbb{R}}f_{i}^{\prime}(x)f_{j}^{
\prime}(x)h(x)\mathrm{d}x+c=0.$
Let $J_{V}(a_{1},a_{2}\cdots,a_{n}):=\sum_{1\leq i,j\leq d}a_{i}a_{j}\int_{\mathbb{
R}}f_{i}^{\prime}(x)f_{j}^{\prime}(x)h(x)\mathrm{d}x+c$. $J$ is a polynomial and thus a real analytic function. Define
$KerJ_{V}:=\left\{(a_{1},a_{2},\cdots,a_{d}):J_{V}(a_{1},a_{2},\cdots,a_{n})=0\right\}$,
then
$\displaystyle f\in S_{\text{nc}}(V)\iff(a_{1},a_{2},\cdots,a_{n})\in KerJ_{V}.$
We first consider the second case where $\int_{\mathbb{R}}\left(p_{x}(x)^{2}+p^{\prime}_{x}(x)^{2}\right)e^{\mu x^{2}}
\mathrm{d}x<\infty$ for some $\mu\in(0,+\infty)$. By lemma 11, $\int_{\mathbb{R}}h(x)^{2}e^{\mu x}\mathrm{d}x<\infty$. If $\int_{\mathbb{R}}x^{n}h(x)dx=0$ for every $n\in\mathbb{N}$, then by lemma 12, $h=0$ almost everywhere. Notice that $h$ is also continuous, implying that $h(x)=0$ for every $x\in\mathbb{R}$. This contradicts with our assumption on $h$. So we know that there exist $M_{0}$, s.t. $\int_{\mathbb{R}}x^{M_{0}}h(x)\mathrm{d}x\neq 0$. Therefore, for $M\geq M_{0}+1$, if we take $V$ to be the space of polynomials with degree smaller than $M$ and $f_{i}(x)=x^{i-1}$, then $J_{V}$ is a non-zero function. By basic knowledge of real analysis [Proposition 0], $KerJ_{V}$ has zero measure. So $m(S_{\text{nc}}(V))=m(KerJ_{V})=0$
For the first case, if $c\neq 0$, then $J_{V}(0,0,\cdots 0)=c\neq 0$. So $J_{V}$ is a non-zero function for every finite dimensional space $V\subset C^{1}(\mathbb{R})$. Similar to the analysis above, we derive $m(S_{\text{nc}}(V))=0$ and thus completing the proof.
$\blacksquare$

###### Lemma 11 .

###### Lemma 12 .

## Appendix H Proof of Proposition 3

Let

$$ $h(x,v):=\int_{\mathbb{R}}p(z)q_{0}(x-\phi_{0}(z))\sigma(x)^{-1}q_{1}\left( \frac{v-\phi_{1}(z)}{\sigma(x)}\right)dz.$ $$

Under our assumption in 5, it is easy to verify that $h$ is well-defined and belong to $C^{1,1}(\mathbb{R}^{2})\cap L^{\infty}$.
Formally, we denote

$$ $\displaystyle\begin{cases}A_{3}&=-\int_{\mathbb{R}^{2}}\frac{f^{\prime}(x)^{3} \partial_{v}h(x,v)^{3}}{h(x,v)^{2}}dxdv,\\ A_{2}&=3\int_{\mathbb{R}^{2}}\frac{f^{\prime}(x)^{2}\partial_{v}h(x,v)^{2} \partial_{x}h(x,v)}{h(x,v)^{2}}dxdv,\\ A_{1}&=-3\int_{\mathbb{R}^{2}}\frac{f^{\prime}(x)\partial_{v}h(x,v)\partial_{x }h(x,v)^{2}}{h(x,v)^{2}}dxdv,\\ A_{0}&=\int_{\mathbb{R}^{2}}\frac{\partial_{x}h(x,v)^{3}}{h(x,v)^{2}}dxdv.\end {cases}.$ $$

The well-defineness of $A_{i}$ ($0\leq i\leq 3$) will be verified in the following lemma.

- Assumption 4 (Restriction on Model with a Latent Confounder).
1.
It exists an open subset $U\subset(0,+\infty)$ s.t. $\mathbb{E}_{X,Y}\left|\left(\partial_{i}\log p_{\lambda}(X,Y)\right)^{3}\right
|<\infty$ for $i\in\left\{x,y\right\}$ and every $\lambda\in U$,
2.
$A_{1}^{2}+A_{2}^{2}+A_{3}^{2}\neq 0$.

###### Lemma 13 .

Under assumption H, $A_{i}$ ($i=0,1,2,3$) are all well-defined integral and

$$ $\mathbb{E}_{X,Y}\left(\partial_{x}\log p_{\lambda}(X,Y)\right)^{3}=A_{3} \lambda^{3}+A_{2}\lambda^{2}+A_{1}\lambda+A_{0}.$ $$

- Proof
$\displaystyle\mathbb{E}_{X,Y}\left(\partial_{x}\log p_{\lambda}(X,Y)\right)^{3}$
$\displaystyle=\int_{\mathbb{R}^{2}}\left(\partial_{x}\log p_{\lambda}(x,y)
\right)^{3}p_{\lambda}(x,y)dxdy$
$\displaystyle=\int_{\mathbb{R}^{2}}\frac{\left(\partial_{x}p_{\lambda}(x,y)
\right)^{3}}{p_{\lambda}(x,y)^{2}}dxdy$
$\displaystyle=\int_{\mathbb{R}^{2}}\frac{\left(\partial_{x}h(x,v)-\lambda f^{
\prime}(x)\partial_{v}h(x,v)\right)^{3}}{h(x,v)^{2}}dxdv$
,
Consequently we can write the integrand as:
$\displaystyle\frac{\left(\partial_{x}h(x,v)-\lambda f^{\prime}(x)\partial_{v}h
(x,v)\right)^{3}}{h(x,v)^{2}}=\sum_{i=0}^{3}g_{i}(x,v)\lambda^{i},$
where
$\displaystyle\begin{cases}g_{3}(x,v)&=-\frac{f^{\prime}(x)^{3}\partial_{v}h(x,
v)^{3}}{h(x,v)^{2}},\\
g_{2}(x,v)&=3\frac{f^{\prime}(x)^{2}\partial_{v}h(x,v)^{2}\partial_{x}h(x,v)}{
h(x,v)^{2}},\\
g_{1}(x,v)&=-3\frac{f^{\prime}(x)\partial_{v}h(x,v)\partial_{x}h(x,v)^{2}}{h(x
,v)^{2}},\\
g_{0}(x,v)&=\frac{\partial_{x}h(x,v)^{3}}{h(x,v)^{2}}.\end{cases}$
As $h(x,v)>0$, $g_{i}$ are continuous. Taking $\lambda_{k}\in U$ ($0\leq k\leq 3$) s.t. $\lambda_{i}\neq\lambda_{j}$ for $i\neq j$. By assumption H, for $0\leq k\leq 3$, $\int_{\mathbb{R}^{2}}\left|\sum_{i=0}^{3}g_{i}(x,v)\lambda_{k}^{i}\right|dxdv<\infty$, i.e. $\sum_{i=0}^{3}\lambda_{k}^{i}g_{i}\in L^{1}(\mathbb{R}^{2})$. Notice that the matrix $A=(A_{ik})_{0\leq i,k\leq 3}$ such that $A_{ik}=\lambda_{i}^{k}$ is a Vandermonde matrix and thus invertible. Then since $L^{1}{(\mathbb{R}^{2})}$ is a linear space, we conclude that $g_{i}\in L^{1}{(\mathbb{R}^{2})}$ for $0\leq i\leq 3$. Consequently $A_{i}=\int_{\mathbb{R}^{2}}g(x,v)dxdv$ is well-defined and
$\mathbb{E}_{X,Y}\left(\partial_{x}\log p(X,Y)\right)^{3}=A_{3}\lambda^{3}+A_{2
}\lambda^{2}+A_{1}\lambda+A_{0}.$
$\blacksquare$

- Proof of Proposition 3

$$ $\displaystyle p(x,y,z)=p(z)p(x|z)p(y|x,z)=p(z)q_{0}(x-\phi_{0}(z))\sigma(x)^{- 1}q_{1}\left(\frac{y-\lambda f(x)-\phi_{1}(z)}{\sigma(x)}\right).$ $$

We write $p_{\lambda}(x,y)$ as $p(x,y)$ for short, therefore,

$$ $\displaystyle p(x,y)=\int_{\mathbb{R}}p(z)q_{0}(x-\phi_{0}(z))\sigma(x)^{-1}q_ {1}\left(\frac{y-\lambda f(x)-\phi_{1}(z)}{\sigma(x)}\right)dz.$ $$

Let $h(x,v):=\int_{\mathbb{R}}p(z)q_{0}(x-\phi_{0}(z))\sigma(x)^{-1}q_{1}\left(
\frac{v-\phi_{1}(z)}{\sigma(x)}\right)dz$. Notice that $h$ is independent of $\lambda$, $f$ and $h\in C^{1,1}(\mathbb{R}^{2})$ under our assumption. As $h(x,y-\lambda f(x))=p(x,y)$, we derive:

$$ $\displaystyle\begin{cases}\partial_{y}p(x,y)&=\partial_{v}h(x,y-\lambda f(x)), \\ \partial_{x}p(x,y)&=\partial_{x}h(x,y-\lambda f(x))-\lambda f^{\prime}(x) \partial_{v}h(x,y-\lambda f(x)).\end{cases}$ $$

Our first observation is that: $\mathbb{E}_{X,Y}\left(\partial_{y}\log p(X,Y)\right)^{3}$ is independent of $\lambda$ and $f$, as is shown below:

$$ $\displaystyle\mathbb{E}_{X,Y}\left(\partial_{y}\log p(X,Y)\right)^{3}$ $\displaystyle=\int_{\mathbb{R}^{2}}\left(\partial_{y}\log p(x,y)\right)^{3}p(x ,y)dxdy$ $\displaystyle=\int_{\mathbb{R}^{2}}\frac{\left(\partial_{y}p(x,y)\right)^{3}}{ p(x,y)^{2}}dxdy$ (By change of variables) $\displaystyle=\int_{\mathbb{R}^{2}}\frac{\left(\partial_{y}p(x,v+\lambda f(x)) \right)^{3}}{p(x,v+\lambda f(x))^{2}}dxdv$ $\displaystyle=\int_{\mathbb{R}^{2}}\frac{\left(\partial_{v}h(x,v)\right)^{3}}{ h(x,v)^{2}}dxdv.$ $$

On the other hand, by lemma 13 we get:

$$ $\displaystyle\left|\mathbb{E}_{X,Y}\left(\partial_{x}\log p(X,Y)\right)^{3} \right|=\left|A_{3}\lambda^{3}+A_{2}\lambda^{2}+A_{1}\lambda+A_{0}\right|.$ $$

Consequently, by assumption H $\left|\mathbb{E}_{X,Y}\left(\partial_{x}\log p(X,Y)\right)^{3}\right|
\rightarrow+\infty$ as $\lambda\rightarrow+\infty$, thus there exists a constant $\lambda^{*}$ s.t. for every $\lambda\geq\lambda^{*}$:

$$ $\left|\mathbb{E}_{X,Y}\left(\partial_{x}\log p(X,Y)\right)^{3}\right|>\left| \mathbb{E}_{X,Y}\left(\partial_{y}\log p(X,Y)\right)^{3}\right|.$ $$

$\blacksquare$

###### Remark 14 .

In particular, if $A_{3}\neq 0$, let $R=\left|\mathbb{E}_{X,Y}\left(\partial_{y}\log p(X,Y)\right)^{3}\right|$, we can set $\lambda^{*}=1+\max\left\{\frac{3\left|A_{2}\right|}{\left|A_{3}\right|},\sqrt{
\frac{3\left|A_{1}\right|}{\left|A_{3}\right|}},\sqrt[3]{\frac{3\left|A_{0}
\right|+R}{\left|A_{3}\right|}}\right\}$. It is easy to verify that

$$ $\left|A_{3}\right|\lambda^{3}-\left|A_{2}\right|\lambda^{2}-\left|A_{1}\right| \lambda-\left|A_{0}\right|>R,$ $$

when
$\lambda\geq\lambda^{*}$.

###### Lemma 13 .

###### Remark 14 .

## Appendix I Explicit computation of example 4

$p(x,y,z)\sim\exp\left(-\frac{z^{2}}{2}-\frac{(x-z)^{2}}{2}-\frac{(y-z-f(x)^{2}
)}{2}\right)$,

$$ $\displaystyle p(x,y)=\int_{\mathbb{R}}p(x,y,z)dz$ $\displaystyle=C\exp\left(-\frac{2}{3}\left(x^{2}+(y-f(x))^{2}-x(y-f(x))\right)\right)$ $\displaystyle=C\exp\left(-\frac{2}{3}\left(y-\frac{x}{2}-f(x)\right)^{2}-\frac {x^{2}}{2}\right),$ $$

$C=\frac{1}{\sqrt{3}\pi}$ is a normalization constant. Therefore,

$$ $\displaystyle\begin{cases}\partial_{y}\log p(x,y)&=-\frac{4}{3}\left(y-\frac{x }{2}-f(x)\right),\\ \partial_{x}\log p(x,y)&=\frac{4}{3}\left(y-\frac{x}{2}-f(x)\right)\left(\frac {1}{2}+f^{\prime}(x)\right)-x\\ &=\frac{2}{3}\left((1+2f^{\prime}(x))y-f(x)-2x-xf^{\prime}(x)-2f(x)f^{\prime}( x)\right)\\ &=\frac{2}{3}\left(\alpha(x)y-\beta(x)\right),\end{cases}$ $$

where $\alpha(x)=1+2f^{\prime}(x)$ and $\beta(x)=f(x)+2x+xf^{\prime}(x)+2f(x)f^{\prime}(x)$.

$$ $\displaystyle\textit{Skew}_{y}$ $\displaystyle=\frac{64C}{27}\left|\int_{\mathbb{R}}\int_{\mathbb{R}}(y-\frac{x }{2}-f(x))^{3}e^{-\frac{2}{3}\left(y-\frac{x}{2}-f(x)\right)^{2}}\mathrm{d}y~{ }e^{-\frac{x^{2}}{2}}\mathrm{d}x\right|$ $\displaystyle=\frac{64C}{27}\left|\int_{\mathbb{R}}\int_{\mathbb{R}}u^{3}e^{- \frac{2}{3}u^{2}}\mathrm{d}u~{}e^{-\frac{x^{2}}{2}}\mathrm{d}x\right|=0.$ $$

On the other hand,

$$ $\displaystyle\textit{Skew}_{x}$ $\displaystyle=\frac{8C}{27}\left|\int_{\mathbb{R}}\int_{\mathbb{R}}\left( \alpha(x)y-\beta(x)\right)^{3}e^{-\frac{2}{3}\left(y-\frac{x}{2}-f(x)\right)^{ 2}}\mathrm{d}y~{}e^{-\frac{x^{2}}{2}}\mathrm{d}x\right|.$ $$

Let $\gamma(x)=\frac{x}{2}+f(x)$. Notice that $\alpha(x)\gamma(x)-\beta(x)=-\frac{3x}{2}$

$$ $\displaystyle\int_{\mathbb{R}}\left(\alpha(x)y-\beta(x)\right)^{3}e^{-\frac{2} {3}\left(y-\frac{x}{2}-f(x)\right)^{2}}\mathrm{d}y=\int_{\mathbb{R}}\left( \alpha(x)(y+\gamma(x))-\beta(x)\right)^{3}e^{-\frac{2}{3}y^{2}}\mathrm{d}y$ $\displaystyle~{}~{}=3\alpha(x)^{2}(\alpha(x)\gamma(x)-\beta(x))\int_{\mathbb{R }}y^{2}e^{-\frac{2}{3}y^{2}}\mathrm{d}y+(\alpha(x)\gamma(x)-\beta(x))^{3}\int_ {\mathbb{R}}e^{-\frac{2}{3}y^{2}}\mathrm{d}y$ $\displaystyle~{}~{}=-\frac{9}{2}x\alpha(x)^{2}\int_{\mathbb{R}}y^{2}e^{-\frac{ 2}{3}y^{2}}\mathrm{d}y-\frac{27}{8}x^{3}\int_{\mathbb{R}}e^{-\frac{2}{3}y^{2}} \mathrm{d}y.$ $$

Let $\mu=\int_{\mathbb{R}}y^{2}e^{-\frac{2}{3}y^{2}}\mathrm{d}y=\frac{3\sqrt{6\pi}}
{8}$, $\int_{\mathbb{R}}e^{-\frac{2}{3}y^{2}}\mathrm{d}y$.

$$ $\displaystyle\textit{Skew}_{x}$ $\displaystyle=\frac{8C}{27}\left|\int_{\mathbb{R}}\left(-\frac{9\mu}{2}x\alpha (x)^{2}-\frac{27\mu_{0}}{8}x^{3}\right)e^{-\frac{x^{2}}{2}}\mathrm{d}x\right|$ $\displaystyle=\frac{4C\mu}{3}\left|\int_{\mathbb{R}}x\alpha(x)^{2}e^{-\frac{x^ {2}}{2}}\mathrm{d}x\right|$ $\displaystyle=\frac{4C\mu}{3}\left|\int_{\mathbb{R}}x(1+2f^{\prime}(x))^{2}e^{ -\frac{x^{2}}{2}}\mathrm{d}x\right|$ $\displaystyle=\frac{16C\mu}{3}\left|\int_{\mathbb{R}}x(f^{\prime}(x)+f^{\prime }(x)^{2})e^{-\frac{x^{2}}{2}}\mathrm{d}x\right|$ $\displaystyle=\frac{2\sqrt{2\pi}}{\pi}\left|\int_{\mathbb{R}}x(f^{\prime}(x)+f ^{\prime}(x)^{2})e^{-\frac{x^{2}}{2}}\mathrm{d}x\right|.$ $$

## Appendix J Additional experiments

Runtime
Table 1 presents the runtime, in seconds, for SkewScore using Sliced Score Matching (SSM) (Song et al.,, ), for various variable counts $d=\{10,20,50\}$, where the number of edges is twice the number of variables. All experiments are on 12 CPU cores with 24 GB RAM.

**Table 1: Wall-clock runtime (in seconds) of SkewScore on ER2 graph with $1000$ samples, comparing different gradient estimators.**
| Method | d=10 | d=20 | d=50 |
| --- | --- | --- | --- |
| SkewScore (SSM) | $310.24\pm 7.17$s | $704\pm 19.82$s | $1825.92\pm 39.23$s |

Effect of Non-Symmetric Noise We conducted experiments adhering to the same bivariate HSNM data generation process outlined in Section 6, except the noise now follows a Gumbel distribution. Our method remains relatively robust to Gumbel noise, a skewed distribution. The results in the table 2 demonstrate that the reduction in accuracy for SkewScore is less significant than for HOST, which operates under the assumption of Gaussian noise. LOCI and HECI are designed to handle various noise types, which unsurprisingly results in strong performance.

**Table 2: Accuracies of estimated causal relationships in the bivariate case with an non-symmetric Gumbel noise. The results are averaged over 100 independent runs.**
| Methods | HNMs (Gumbel) |
| --- | --- |
| SkewScore | 89 |
| LOCI | 100 |
| HECI | 99 |
| HOST | 33 |
| DiffAN | 92 |
| DAGMA | 27 |
| NoTears | 70 |

Performance on Additive Noise Model (ANM) Though SkewScore is designed for heteroscedastic noise, it is still valid on nonlinear models with homoscedastic noise. We explore the algorithm’s performance using data generated according to Additive Noise Models (ANM) with a GP as the nonlinear function and Gaussian noise.
After obtaining the topological order with SkewScore, CAM pruning (Bühlmann et al.,, ) is applied. We calculate the topological order divergence and the structural Hamming distance (SHD) across varying dimensions, focusing on graphs with $d$ and $2d$ edges, represented as ER1 and ER2 graphs, respectively. The results, presented in Tables 3 and 4, indicate that SkewScore maintains adequate performance in ANM settings. Results are averaged over 10 independent runs.

**Table 3: Topological order divergence across dimensionalities ($d\in\{10,20,50\}$) using data generated from ANMs.**
|  | ER1 | ER2 |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| Methods | d=10 | d=20 | d=50 | d=10 | d=20 | d=50 |
| SkewScore | $\textbf{1.00}\pm\textbf{1.00}$ | $\textbf{3.10}\pm\textbf{0.70}$ | $\textbf{8.50}\pm\textbf{3.17}$ | $3.80\pm 2.44$ | $11.00\pm 4.12$ | $\textbf{23.40}\pm\textbf{5.82}$ |
| HOST | $3.90\pm 1.45$ | $8.40\pm 2.11$ | $15.00\pm 3.74$ | $11.70\pm 2.45$ | $22.20\pm 3.79$ | $51.60\pm 3.77$ |
| DiffAN | $1.80\pm 0.98$ | $4.50\pm 3.07$ | $14.89\pm 2.18$ | $\textbf{2.70}\pm\textbf{1.19}$ | $\textbf{7.00}\pm\textbf{2.32}$ | $23.43\pm 2.50$ |
| DAGMAMLP | $3.70\pm 1.10$ | $7.70\pm 1.79$ | $15.70\pm 4.65$ | $10.20\pm 0.75$ | $17.10\pm 1.51$ | $35.00\pm 2.37$ |

**Table 4: SHD across different dimensionalities ($d\in\{10,20,50\}$) using data generated from ANMs.**
|  | ER1 | ER2 |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| Methods | d=10 | d=20 | d=50 | d=10 | d=20 | d=50 |
| SkewScore | $\mathbf{2.90\pm 1.47}$ | $\mathbf{7.08\pm 2.32}$ | $\mathbf{19.32\pm 7.26}$ | $8.26\pm 2.42$ | $30.50\pm 5.52$ | $\mathbf{76.50\pm 5.50}$ |
| HOST | $7.20\pm 2.04$ | $17.30\pm 3.95$ | $44.70\pm 6.56$ | $17.80\pm 2.44$ | $37.10\pm 4.04$ | $95.70\pm 6.42$ |
| DiffAN | $3.50\pm 1.40$ | $9.90\pm 1.87$ | $23.56\pm 1.27$ | $\mathbf{7.31\pm 1.42}$ | $\mathbf{30.10\pm 4.06}$ | $82.29\pm 5.95$ |
| DAGMAMLP | $6.90\pm 1.04$ | $15.60\pm 1.28$ | $37.60\pm 4.18$ | $12.30\pm 0.64$ | $37.50\pm 1.20$ | $92.40\pm 1.62$ |

Effect of Sample Size Table 5 illustrates the order divergence for various models on synthetic data generated, measured across increasing sample sizes $n=\{100,1000,10000\}$. This analysis was conducted with a fixed dimension $d=10$ with an equal number of edges, utilizing an ER graph model. The results indicate that our method, specifically SkewScore, demonstrates enhanced performance as the sample size increases. This trend suggests that SkewScore can efficiently leverage larger datasets to improve estimation accuracy. Notably, even with smaller sample sizes, SkewScore maintains commendable performance, underscoring its robustness and effectiveness in scenarios with limited data availability. This capability makes it particularly advantageous in practical applications where acquiring large volumes of data may be challenging.

**Table 5: Topological order divergence on synthetic data is analyzed as a function of sample size. We fix $d=10$ with the same number of edges, using an ER graph model, and vary the sample size across $n=\{100,1000,10000\}$. Results are provided for both the HSNM and the HSNM with latent confounders.**
|  | HSNMs | HSNMs with Confounders |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
|  | n=100 | n=1000 | n=10000 | n=100 | n=1000 | n=10000 |
| SkewScore | $\mathbf{1.90\pm 1.70}$ | $\mathbf{0.90\pm 0.75}$ | $\mathbf{0.60\pm 0.32}$ | $\mathbf{2.10\pm 0.83}$ | $\mathbf{1.30\pm 0.64}$ | $\mathbf{1.10\pm 0.58}$ |
| HOST | $4.40\pm 1.91$ | $4.40\pm 1.11$ | $3.80\pm 1.47$ | $4.56\pm 1.34$ | $4.20\pm 1.60$ | $3.90\pm 1.45$ |
| DiffAN | $4.70\pm 1.68$ | $2.30\pm 1.19$ | $2.40\pm 0.92$ | $3.60\pm 1.50$ | $2.80\pm 1.47$ | $2.75\pm 1.39$ |
| VarSort (Reisach et al.,, ) | $2.10\pm 1.45$ | $2.00\pm 1.38$ | $1.50\pm 1.81$ | $2.20\pm 1.54$ | $2.10\pm 1.51$ | $1.90\pm 1.37$ |

#### Real Data

Our method’s accuracy is 62.6% in the Tübingen cause-effect pairs dataset, excluding multivariates and binary pairs (47, 52-55, 70, 71, 105 and 107) as Immer et al., (). According to the reported results from Immer et al., (), our method performance is slightly better than LOCI. GRCI Strobl and Lasko, () performs the best, and QCCD Tagasovska et al., () also performs quite well. These methods are likely more robust to assumption violation in real data such as measurement error or low sample complexity. We will include this experiment and the advantages of these two methods in our revision.

Multivariate extension of the case study in Section 5 Finally, we conduct experiments to explore the multivariate extension of the case study in Section 5. The data generation process for experiments for multivariate HSNMs with latent confounders is based on . Let $Z\in\mathbb{R}^{d}$ represent the latent common cause. For each pair of distinct nodes $X_{i}$ and $X_{j}$, a Bernoulli random variable $C_{ij}\sim\text{Bernoulli}(\rho)$ is sampled, where $C_{ij}=1$ implies a confounding effect $Z$ between $X_{i}$ and $X_{j}$. The parameter $\rho$ determines the sparsity of confounded pairs in the graph. In our experiments, we choose $\rho=0.2$. Data generation follows a HSNM with latent confounders: $X_{i}=f_{i}(\text{pa}_{X}(X_{i}))+\phi_{i}(\text{pa}_{Z}(X_{i}))+\sigma_{i}(
\text{pa}_{X}(X_{i}))\cdot N_{i},$ for each variable $i=1,\ldots,d$. Here, $\text{pa}_{X}(X_{i})$ and $\text{pa}_{Z}(X_{i})$ represent the observed and hidden variable parents of $X_{i}$, respectively. The model incorporates a Gaussian processes with an RBF kernel of unit bandwidth to generate the function $f_{i}$, a linear function $\phi_{i}$, and a conditional standard deviation $\sigma_{i}$ modeled as a sigmoid function of the observed parents’ linear combination. For a given number of observed variables $d$, we adjust the sparsity of the sampled graph by setting the average number of edges to either $d$ (ER1) or $2d$ (ER2).
The results are shown in 6 and 7, where the lower topological order divergence means better performance. In general, SkewScore, Diffan, and NoTears have the best performances. In sparser graph (ER1), SkewScore performs the best when $d\in\{10,20,30,40\}$. In the ER2 graph, SkewScore and Diffan have similar good performances. The theoretical insights for this extension of the case study in Section 5 will be the focus of future work.

**Table 6: Topological order divergence: Synthetic data generated using multivariate HSNM with latent confounders with an ER1 graph.**
| Method | d=10 | d=20 | d=30 | d=40 | d=50 |
| --- | --- | --- | --- | --- | --- |
| SkewScore | $\mathbf{1.20\pm 0.98}$ | $\mathbf{2.40\pm 2.33}$ | $\mathbf{4.60\pm 1.80}$ | $\mathbf{6.60\pm 3.29}$ | $7.20\pm 2.64$ |
| HOST | $4.00\pm 1.48$ | $8.10\pm 1.64$ | $10.40\pm 2.29$ | $14.90\pm 2.66$ | $17.50\pm 2.73$ |
| DiffAN | $2.80\pm 1.17$ | $5.60\pm 2.62$ | $8.60\pm 2.91$ | $9.50\pm 2.46$ | $11.60\pm 3.67$ |
| DAGMAMLP | $6.00\pm 2.00$ | $14.50\pm 2.11$ | $19.30\pm 2.49$ | $27.10\pm 3.33$ | $34.50\pm 4.13$ |
| NoTears | $1.50\pm 1.02$ | $5.30\pm 2.45$ | $5.00\pm 2.57$ | $6.80\pm 3.97$ | $\mathbf{6.30\pm 2.45}$ |

**Table 7: Topological order divergence: Synthetic data generated using multivariate HSNM with latent confounders with an ER2 graph.**
| Method | d=10 | d=20 | d=30 | d=40 | d=50 |
| --- | --- | --- | --- | --- | --- |
| SkewScore | $\mathbf{3.00\pm 1.61}$ | $8.60\pm 2.11$ | $14.00\pm 4.73$ | $\mathbf{16.60\pm 2.42}$ | $25.80\pm 4.26$ |
| HOST | $11.80\pm 2.44$ | $22.80\pm 4.45$ | $32.90\pm 5.47$ | $41.90\pm 6.55$ | $51.80\pm 6.27$ |
| DiffAN | $3.30\pm 2.10$ | $\mathbf{7.20\pm 2.40}$ | $\mathbf{13.20\pm 3.49}$ | $18.10\pm 3.86$ | $\mathbf{22.10\pm 4.44}$ |
| DAGMAMLP | $17.60\pm 1.80$ | $34.90\pm 2.47$ | $54.20\pm 2.18$ | $74.40\pm 2.15$ | $91.20\pm 3.19$ |
| NoTears | $8.10\pm 2.81$ | $16.70\pm 6.81$ | $20.90\pm 5.73$ | $25.30\pm 7.06$ | $27.00\pm 4.45$ |

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: 1. Is the algorithm scalable?  It would be valuable to discuss if this approach can handle high-dimensional datasets effectively and whether its complexity scales polynomially or exponentially with the number of variables. 

2. Can you discuss the relation between the proposed method with other types of causal discovery algorithms, such as (1) constraints-based algorithms such as PC algorithm etc., and (2) continuous optimization based algorithms such as the one in https://arxiv.org/abs/1803.01422, which provides scalability? 

3. What are the effects of small sample sizes on SkewScore’s performance? It would be great if authors can discuss the convergence guarantee of the algorithm in terms of the size of samples. 

4. Can the proposed algorithm provide confidence intervals or uncertainty quantification in measuring the skewness or determining causal directionality?
soundness: 4
strengths: 1. An example in Section 5 provides a clear, practical illustration of the algorithm, connecting theoretical concepts with application. 

2. This paper provides a well-structured and self-contained summary of prior works in causal discovery algorithms based on functional causal models, allowing readers to clearly follow the evolution for learning causal graphs. As a reviewer, this summary enables me to track the advancements in the field and understand the recent progress. 

3. The core idea behind the proposed method is simple and intuitive, based on the asymmetry of score skewness to identify causal directions. This approach is conceptually accessible and computationally efficient solutions. 

4. The experiments provide strong empirical evidences of the paper’s contributions, with a wide range of scenarios that demonstrate the robustness of the proposed framework in heteroscedastic settings.
summary: This paper provides a novel causal discovery algorithm that leverages skewness to determine causal direction when heteroscedastic noise (HN) is present.
weaknesses: 1. The paper provides correctness guarantees only for the two-variable case, where skewness in the score function can reliably distinguish causal directions. However, for settings with more than two variables, there is no formal guarantee of correctness for Algorithm 1. 

2. Current empirical results seem insufficient in showcasing its performance in high-dimensional scenarios. In other words, it remains unclear whether SkewScore scales effectively as the number of variables grows significantly. 

3. While the authors claim robustness in the presence of latent confounders, this robustness appears to be limited to simple bivariate cases. In multivariate settings, SkewScore can still produce incorrect causal graphs when latent confounders are present. For instance, in the causal structure $X_1 \rightarrow X_2 \rightarrow X_3$ with a latent confounder $L$ influencing both $X_1$ and $X_3$ ($X_1 \leftarrow L \rightarrow X_3$), the method would infer the structure $X_1 \rightarrow X_2 \rightarrow X_3$ along with a spurious direct link $X_1 \rightarrow X_3$. This incorrect edge arises because the latent confounder $L$ induces a conditional dependence between $X_1$ and $X_3$, which SkewScore interprets as a direct causal link. As a result, while the method can handle some latent confounding, it does not offer sufficient mechanisms to detect and mitigate confounding in more complex structures, limiting its application in practical scenarios where latent variables are common.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 3
questions: Is it possible to provide identification results such as 'the skew score is larger under the wrong causal direction' without symmetric noise assumption?
soundness: 2
strengths: The notation is well-organized and the paper is easy to follow overall.
summary: This work proposes a causal discovery method for a class of restricted structural causal models with heteroscedastic symmetric noise. The identification of the causal graph is established using a score that quantifies the skewness of the score function. A causal discovery algorithm is then proposed based on the skewness score.
weaknesses: 1. Fact 2 is incorrect. Let X follow univariate exponential distribution with parameter $\lambda$, then $E[\frac{d \text{log} p(X)} { dx}] = -\lambda$. Its proof is also incorrect. 

2. 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. 

3. 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.

---

confidence: 2
contribution: 2
initial_rating: 3
presentation: 3
questions: 1. 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.
2. 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?
soundness: 3
strengths: 1. The introduction of the SkewScore, which is leveraged to identify causal structures without requiring the extraction of exogenous noise.
2. A theoretical extension of the criterion to multivariate cases, along with empirical validations showing its robustness, particularly in scenarios involving latent confounders.
summary: The paper presents a novel approach to uncovering causal relationships in data with heteroscedastic noise, which challenges traditional causal discovery models. The authors propose the use of a skewness-based criterion derived from the score of data distributions to identify causal directions, establishing that this measurement is zero in the causal direction but non-zero in the anti-causal direction.
weaknesses: See questions.

---

confidence: 4
contribution: 4
initial_rating: 10
presentation: 4
questions: ### On Proposition 5
Could you elaborate on why, “when $\sigma(x)$ is not constant (i.e., in the heteroscedastic model), $c \neq 0$ in most cases”? More broadly, could you comment on the condition $c=0$? This characterization could be crucial as it defines a class of distributions $p(u, x)$.

Proposition 5 considers $f$ in a finite-dimensional linear space, specifically polynomials. Could neural networks with hidden layers be considered also (since Example 7 doesn’t)?

It’s interesting that $f$ is required to lie in a “large finite-dimensional space.” I would assume that if $f$ were in an infinite-dimensional functional space, the space should instead be “small.” Could you share your thoughts on this apparent contrast?

### Additional Comments
Could similar theoretical results be proven for standard skewness?

Assumption 4: Typically, “regularity” refers to the properties of a well-behaved functional space. In Assumption 4, condition 1) is indeed a regularity condition, but it’s unclear if condition 2) qualifies as such. Additionally, could you explain why condition 2) is considered mild?

Minor: There is a duplicate bibliography entry for Mityagin (2015).
soundness: 3
strengths: The paper provides a neat and effective solution to a critical problem. The connection between score skewness and causal direction is both theoretically grounded and practically effective.

The extensions to multivariate cases and the potential applicability under hidden confounding are valuable contributions.

The theoretical analysis appears rigorous (though I haven’t verified every detail).

The paper is clearly written, with helpful elements like the two properties of the score function, the motivation of score skewness by analogy to standard skewness, and well-designed figures.
summary: ==Post-rebuttal update==

The authors' replies to my questions show their further insights into the problem and how the work could be deepened.

I also had a rough scan of other reviews and their respective replies, and I hold my original score.

This is an exceptionally well-written paper with valuable insights and high potential impact. Of course, no papers are perfect, but the current rating best describes my evaluation of the paper, i.e., it "should be highlighted at the conference".

==End update==

This paper proposes a method to identify causal directions in heteroscedastic symmetric noise models (HSNMs) using the skewness of the score function, which is the gradient of the log density. The skewness considered here is analogous to standard skewness but is distinct, offering computational advantages. Importantly, the paper proves that the score's skewness is zero only in the causal direction. This idea is later generalized to the multivariate case. The proposed method also shows robustness in the presence of hidden confounders, and a theoretical analysis is provided for an "additive" confounder in the bivariate case. The method demonstrates strong empirical performance.
weaknesses: The generality of Assumption 1, which defines the applicable HSNMs, could be explored further, although the paper presents sufficient results for a conference submission.

The method's applicability under hidden confounding might be limited to special confounding structures.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 4
questions: 1. Which ER graph model is used in the experiments, $G(d, p)$ or $G(d, m)$? 
2. Can you explain in more detail why $c\neq0$ in general in Remark 6?
3. What is the runtime for the order learning and independent tests separately?
soundness: 3
strengths: 1. The method is novel and can be potentially applied to causal discovery with latent confounders.
2. The theoretical justification of Assumption 1 seems sound, and the experiments showed that the proposed method works particularly well on sparse HSNMs.
3. The author extensively evaluated the method when the assumptions are violated, e.g. non-symmetric noise distribution, and latent confounder in the multivariate setting, and proved its effectiveness and robustness.
summary: In this paper, the author proposed a novel causal learning algorithm for heteroscedastic symmetric noise models (HSNMs), where the noise term of a variable $X_i$ can be parameterized by $\sigma(pa(X_i))N_{X_i}$ with symmetric noise $N_{X_i}$. Specifically, the author showed that the skewness of the score function (SkewScore) in the $X_i$-coordinate is zero if and only if $X_i$ is a leaf node under certain assumptions. Therefore, the topological order can be identified by iteratively removing the variables with the least SkewScore. Conditional independence tests are then applied to determine the edges. Its validity for the HSNM with latent confounder is also discussed. Extensive experiments showed that the proposed method is effective on HSNMs and is robust to assumption violations.
weaknesses: 1. The author mainly used ER1 to evaluate the methods on HSNMs without latent confounders, which is relatively sparse.
2. 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?

Meta-review: The paper tackles an important and challenging problem of causal discovery with heteroscedastic noise (HN), by utilising the third moment (skewness) of the score function (derivative of the log density). The idea with contrastive-based learning idea for the HN case is novel and substantial; the connection with causal discovery, score matching and skewness further strengthen the insights of the proposed approach. As the reviewers generally reach the consensus that the paper is nicely-organised and well-written; the contribution is significant for our community, the paper is recommended to appear and highlight in the ICLR conference.

Decision: Accept (Poster)