A Skewness-Based Criterion for Addressing Heteroscedastic Noise in Causal Discovery
Authors: Yingyu Lin
Yuxing Huang
Wenqin Liu
Haoran Deng
Ignavier Ng
Kun Zhang
Mingming Gong
Yian Ma
Biwei Huang
Venue: ICLR 2025 Poster (2025)
Keywords: Causal Discovery
Heteroscedastic Noise
Score Matching
Primary area: causal reasoning

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$ is modeled as $Y = f(X) + \sigma(X)N$, with $X$ as the cause and $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 \texttt{SkewScore}, an algorithm that handles heteroscedastic noise without requiring the extraction of exogenous noise. We also conduct a case study on the robustness of \texttt{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.

Published as a conference paper at ICLR 2025
A SKEWNESS-BASED CRITERION FOR ADDRESSING HETEROSCEDASTIC NOISE IN CAUSAL DISCOVERY
Yingyu Lin"!, Yuxing Huang Wenqin Liu Haoran Deng" Ignavier Ng²,
Kun Zhang Mingming Gong3. Yi-An Ma', Biwei Huang UCC San Diego, "New York University, Thhe University of Melbourne, *Zhejiang University, "Camegie Mellon University, Mohamed bin Zayed University of Artificial Intelligence {yi1208, yianma, bih007 J@ucsd. edu, edu, yh6004 @nyu.edu, (wengin. liu, mi mingming gong @un im b. edu au, denghaoranndzju. edu. cn, iigaviierng, kunz cmu. edu
edu,
ABSTRACT
Real-world data often violates the equal- 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 is modeled as Y f(X) o oX)),, with X the cause and N independent noise following a symmetric distribution. We introduce a novel crite- rion 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 computa- tionally tractable measurement that is zero in the causal direction but nonzero in the anticausal direction, enabling the causal direction discovery. We extend this skewness- criterion to the multivariate setting and propose SkewScore, algorithm that handles heteroscedastic noise without requiring the extraction of exogenous noise. We also conduct 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.
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., 2005) and Earth system science (Runge et al. 2019). The primary approaches to causal discovery include constraint-bar based methods (Spirtes & Glymour 199T, Spirtes al. (00T), which rely on conditional independence tests, and score-based methods (Heckerman et al.) 1995: Chickering 2002, Huang et al. 2018), which optimize a certain objective function. These methods often identify the causal structure only up Markov equivalence (Spirtes et al. 2001 Glymour al. 2019), since the causal directions in general cannot be determined without prior knowledge or additional assumptions (Pearl 2009).
Another line of approaches, based on functional causal models, impose restrictions the causal mechanisms to make the causal directions identifiable. The key idea is to show that, under these restrictions, 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 al.) 2006), nonlinear additive noise models (Hoyer et al., (008a), and post- nonlinear causal models (Zhang & Hyvärinen 2009). 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.) 2021) and robotics (Kersting et al. 2007).
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 & Lasko 2023 [mmer et al. 2023). This procedure needs be repeated for
Figual contribution.

---

Published as a conference paper at ICLR 2025
each variable, which can be computationally intensive. Alternative strategies that avoid extracting noise terms include the likelihood- distributions (Immer et al.] 2023; Khemakhem et al. 2021; Duong & Nguyen 2023), and the kernel- criterion proposed by Mitrovic et al. (2018). whose theoretical guarantees for handling heteroscedastic noise may not be clear.
approach. which often requires the noise to follow Gaussian
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 's t distributions, with applications in, e.g., finance (Anderson Arnold 1993), environmental science (Damsleth & El-Shaarawi 2018), and psychology (McGill 1962). Specifically, we leverage the skew- ness of the score (the gradient of the log density) of the data distribution to establish a mcasurement that is zero in the causal direction but positive in the anticausal direction, thereby uncovering the causal direction. We extend this skewness-basse 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 & Koller (2005). First, we estimate the topological order by iteratively identifying the sink node of the causal graph using the skewness-basee 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. (001)) 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 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 's impact within this triangular model.
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. (2022) 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. (020)); Strobl & Lasko (2023): y: [mmer et al. (2023) considered more general class of HNMs, where the conditional variance
deterministic function of its
parents. The works by Tagasovska et al. (2020) and [mmer et al. (2023) focus primarily on identifying pairwise causal relationships. Both (Strobl & Lasko (2023) and Duong & Nguyen (2023) include
noise-extraction step to estimate the sink node, with Strobl & Lasko (2023) using mutual information and Duong & Nguyen (2023) using normality as criteria. Yin et al. (2024) 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.)
2001] Colombo et al.] 2012
Akbari et al.] 2021), over-complete ICA-based techniques (Hoyer al.
2008b Salehkaleybar et al.
2020). Tetrad condition (Silva et al. 2006; Kummerfeld & Ramsey] 2016),
high-order moments
(Shimizu et al. 2009 Cai et |l.] 2019 Xie et al.] 20200] Adams al. [202]] Chen et al. 2022). matrix decomposition techniques (Anandkumar et al. 2013). mixture oracles (Kivva et al.) 2021) and rank
constraints (Huang et al.] 2022; Dong et al. 2023). 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 + Z + Y and X + Y, where Z is the latent variable.
Causal Discovery with Score Matching, A recent class of causal discovery algorithms employs constraints the score, i.e., the gradient of the log density function, to uniquely identify the
2

---

Published as a conference paper at ICLR 2025
causal graph in additive noise models (ANMs). This line of thought determines the topological order by iteratively identifying sink nodes (i.c.. leaf nodes) using the score (Rolland et al. 2022
Montagna et al. 2023c; Sanchez et al. 2023 Montagna et al. 2023b). Interestingly, MMontagna et al. (2023a) empirically indicates that the presence of latent confounders may not significantly disrupt the inference of the topological ordering; score-matching-ban
approaches still provide reliable ordering. Effective gradient estimation techniques are crucial in connecting the gradient of the log- oggllikelhhoo to successful causal discovery. Various methods have been proposed estimate the score from samples generated by an unknown distribution, including score matching (Hyvärinen & Dayan 2005 Vincent Song & Ermon 2019 [Song al., 2020), and kernel score estimators based on Stein's 's methods (Li & Turnrr] 2017: Shi et al., 2018: Zhou et al. 2020). These methods have been successfully applied to generative modeling, representation learning, and addressing intractability approximate inference algorithms.
PROBLEM SETUP AND PRELIMINARIES
Symbols and notations. Let R. = (0, oo). We denote by C (x*) the set of all continuously differentiable functions on Rd and by C* (R* I++ the subset of (1R) consisting of functions that
take only positive values. We denote by - (IR)) the set of all essentially bounded functions on IR
Furthermore, let - = ( and C = - (R*)).. For brevity, we write ūz; as Ōz,- 2
represent (XI,.... Xi-1, X.+1,.... X Xa) as X-i, and similarly, (x1...... 2;-1 I4+1, T.) as -- We use px (r), or simply p(x), to denote the probability density of X, and px, xx_ (0)(2=+) as p(=)(=_;) for conditional densities. We denote by pa( X,) ) the set of parents (direct causes) of X; in a directed acyclic graph (DAG), and by child(X,) the set of children (direct effects) of X,, A node with no children, i.c., no outgoing edges, is called a sink or leaf node.
For random variable X (X1, X2, ................................ Xa) in Rd,
3.1
HRTEROSCEDASTIC SYMMETRIC NOISE MODELS (HSNM)
We consider structural causal model over d random variables {X,] where each variable X; corresponds to node in the directed acyclic graph (DAG) G. The Heteroscedastic Symmetric Noise Model (HSNM) is defined
X, = (pa( + o ==
(1)
where N; an exogenous noise variable that follows a symmetric distribution with mean zero, i.e.,
( (-n)) PN, (n,). The noise terms {N,) are jointly independent, and each N, is independent of pa(X.) when pa(X,) * 0. We assume that px CRR. (RR, R4), and for each E 21,2.... .d. PN. C* (F, R.). with f., of E C¹ and a; E L Additionally, all partial derivatives Oz of L We also assume there exists a constant > 0 such that inf inf (1......d) inl a;(x) 2 r. Furthermore, we assume that px. PN. fi, and o; satisfy weak regularity conditions to ensure that the derivatives and expectations involved in our analysis are well-defined.
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 distributions, with applications in, e.g., finance (Anderson & Arnold 1993). environmental science (Damsleth & El-Shaarawi 2018). and psychology (McGill 1962).
3.2 SCORE FUNCTION, SCORE MATCHING, AND SKEWNESS
In this paper, the score function refers to the gradient of the log density, Vp(z), which can be estimated using various score matching methods Hyvärinen & Dayan 2005 Vincent (201 1); Song & Ermon (2019): Song et al. 2020); Li & Turner (2017); ;); Shi et al. 2018); Zhou et al. (2020). We introduce two important properties of the score function and define the skewness of the score.
.a)) E CRd R2) ) be a density function that is strictly postive on Rd. The Facts Let p( 22,
following two facts hold for every i = 1,2...,,0.
Fact Oz logp(1,(<1, I1-1, #141, ... ,T)) Or log Log,((11, x2. ,.)).
Fact 2. Ex-p 10n, (ogg(X)) 0.
The proofs of Fact T and Fact are provided in Appendix D These two facts are crucial for deriving
Theorem [] skewness-ban based criterion for HSNMs. Fact 1 aallww us work with derivatives of the
3

---

Published as a conference paper at ICLR 2025
log joint density rather than directly estimating conditional densities, while Fact shows that the expectation the score is zero vector, when E C² CR RR).
Skewness measures the asymmetry of 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 Skew(W) E (W - E[W])" Inspired by
this, we introduce the notion of skewness for the score function.
Definition 1. Given a probability density p for random variable X = 12. ;a)) Rd with p(x) E C* (R4, R+), we define its skewness of the score the x,-coordinate as:
SkewScore (p) = E xx (0z, log - Exap (Ox, log px))
(2)
= Exap (Oz, log x())
where the second equality follows from Fact
This measure captures the asymmetry within each coordinate of the score. If the conditional density p(ax x--)) is symmetric function, then SkewScore,, (p) = 0, as discussed in Appendix
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. ( 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 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. We then 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, shown in Corollary P
Consider the following bivariate heteroscedastic symmetric noise model (HSNM):
Y = f(x) - (
(3)
where N symmetric exogenous noise independent of with mean zero, i.e., Pr (-n) = Pn (n). and X IL N. We assume the same regularity conditions stated in Section 3.1 specified Assumption
and in Appendix
The key insight for the identifiability of the heteroscedastic symmetric noise model (HSNM), given in Eq. 3 is that the conditional distribution P(Y X) is symmetric, while P(X Y is generally asymmetric under mild assumptions f and a. This contrast is illustrated in Figure I 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 V log p( X, Y), where (X, Y) follows the data distribution Specifically, Fact 1 allows us to convert the partial derívatives of the log conditional density into the gradient of the log joint density, which 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 with the model restriction given in Assumption I
Assumption 1 (Model Restriction). Denote A(x)
(w) d'(x) B(x) - f')) a(z) and C(x)
P. (x) a(x)
d'(x) ( where p'g. o', and f' represent the derivatives of px, o, and f respectively. Assume the following inequality holds:
I (A + C(x)
PLtu)u Pala)
+ B(xx (1(aa ? '()))(()))) (A(x)= + C(x)
Pl (L(())) (w) Pn. (1))) o
4

---

Published as a conference paper at ICLR 2025
ply]))
ththe aftar
p(xly)
X: the cause
Figure Identifiability of heteroscedastic symmetric noise model (HSNM). causal direction X (left), the conditional distribution p(yla) is symmetric any æ (top right), and variance of p(y(z) varies with different values of due to heteroscedasticity. In contrast, p( y) is asymmetric for some y (bottom right).
I 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 [ Appendix B we demonstrate that Assumption [ hholds generically by showing that the set of models violating this assumption has property analogous to having zero Lebesgue measure. More discussion is provided in Appendix B
[
We now present 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. 33 Let p(z, v) 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:
SkewScore,, Ex, (00, logg(( Y>> 0.
2. Moreover, under Assumption 7 the score function's 's component in the X-coordinate (the cause direction) follows skewed distribution:
SkewScore, = Exxy (0= logp(XX #
(5)
The proof is provided in Appendix B The asymmetry in the above theorem can used to deter- mine the causal direction between X and Y. Specifically, one could estimate the score function V log p(x, y) and determine the causal direction based on the magnitude of the skewness of its x- and y-coordinate. This result naturally generalizes the multivariate case for d > stated in the following corollary: the score's 's component only along the sink-node coordinate is unskewed, under weak assumptions.
We first clarify the notations: Let X (X1, X2. E Rd follow the multivariate heteroscedastic symmetric noise model given by Eq. ① and p(.) = p(1 1,51,,22 U2. ... .aa denote the joint distribution of X. For simplicity, denote dæ = dx dxy dxa. fo = f((pa(())) a; =,,(p((z)), Pt
Ti-fi and pl = P'v, (
PN.
Corollary 2. Assume that the following condition holds for every node with direct causes (i.e..
child(k) + +0):
pl.
o,0m, fo (z.
f1)02.
Ony σi
]
3
II
Pj dz 0.
(6)
TEPk
Pi
of
sEchild(k)
Then, for every node j. it is sink if and only if the score's component in the co-coordinate follows an unskewed distribution, ie..
SkewScore x, (p) 0,
where SkewScorez, (p) is defined in Definition
5

---

Published as a conference paper at ICLR 2025
The proof is provided in Appendix This corollary enables to find the sink nodes of 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 provide practical algorithm that leverages this Corollary for identifying the causal order.
4.2 ALGORITHM
We apply Corollary 2nnd and propose an algorithm determine the causal ordering, named Skewscore (Algorithm D Our approach follows the topological order search framework of Rolland et al. (2022), 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. 2023 Duong & Nguyen 2023 Montagna al. (023b) or to estimate the diagonal of the Hessian of the log density (Rolland et al. 2022 Sanchez et al. 2023). Additionally, after finding the minimum SkewScorez, in Line 6 of Algorithm I can check whether SkewScore is close to zero, where a value not close to zero indicates a violation of the symmetric noise assumption or other model assumptions. We note that the algorithm's numerical stability and sample complexity can depend on the test odd function chosen in the skewness definition, see Appendix
Once the topological order is obtained using Algorithm T the complete DAG be constructed by adding edges with conditional independence tests (Zhang et al. (001T) associated with the topological order. In this order-based search, (dd)) CI tests are performed. compared exponential complexity in classical methods (Spirtes & Glymour 1991 in worst-case scenarios.
Algorithm 1 SkewScore
Input: Data matrix X E Raxd test odd function (+()) = 8², conditional independence test oracle CI (returns p-value), and significance level (0, Initialization: nodes = {l..... .d......, d}, topological order = [ adjacency matrix A E IRdxd (all
zeros).
for k = L....,.,dd d do
Estimate the score function (XX
Vlog Px.... (X).
Estimate SkewScorez, [Ex.. [s) (X))]l > Skewness of the score's projection on zy D> Find the sink node
+ arg gnndes SkewScore,;
< [e,
> Update the topological order
nodes f- nodes- - {E} Remove (-th column of X
10. end for
for =
1
d do
12:
for j = i +,,1..... 1,....de do
if CI( CICIIxR. ( Xx X (1)))))))))) <a then
13:
14:
An,.,
D. Add edge n; + T
end 16: end for
15:
17: end for
18. Output: Topological order स, and adjacency matrix A.
5 CASE STUDY: ROBUSTNESS TO THE LATENT CONFOUNDER IN BIVARIATE SETTING
It has been shown that various types of assumption violations may affect the performance of causal discovery (Reisach et al. 2021 Montagna et al. 2023a Ng et al. 2024). In this section, we conduct a case study to analyze the robustness of Algorithm to the latent confounder in a pairwise causal model given by Eq. We provide theoretical insights into this robustness in Proposition 3 with experimental results presented Section 6
We focus on a three-variable model consisting of two observed variables, X and Y, and one latent variable, Z, depicted in Figure
2
The model is formally expressed as:
= Φo(Z) No.
Y = xf(x) + Φ1 (2) + (X)NN1

---

Published as a conference paper at ICLR 2025
where N, ds with (n) c[(-n) and qi E C' CR,, for 0, and Z ~ p=()). In this model, the effect of the latent confounder has an additive structure, and the conditional standard deviation or (>X) only depends on the observed variable X, without interaction with the latent variable. The function is replaced with Xf, where À 0 scales the effect of X on Y. To analyze the confounding effects of Z, we vary À while keeping a, ts Us and (2(2) fixed. As in Section 3.1
we assume f. a E C² (R), o, o' E 1 (TR), and that there exists a constant > such that a (x) 2 r.
Y
xf
Figure 2: Three-variable HSNM model with latent confounder given by Eq.
The criterion in Theorem D Iidettffies the causal order without requiring exogenous-nois extraction. This simplicity allows us to quantify the impact of latent confounders. Notably, to recover the causal direction in the infinite sample limit, Line in Algorithm[]] is is equivalent to observing that the score's s component in the effect coordinate is less skewed than in the cause. The following proposition shows that in Model 7 with À large enough, the causal direction X Y can be identified by Algorithm in the infinite sample limit, with only mild regularity assumptions on the function f, o.. qi for i E {0, 1]. and pz. These conditions are discussed in Appendix H aloogg with the proof.
Proposition 3. Suppose two variables X and Y follow the causal models in Eq. 7 Given f. a. di
q. for i {0, 1], and pz, let px (x( y) denote the joint distribution of (X, Y). Under assumption 4 there exists X E R,,, such that for every X 2 X", we have SkewScore (Px) > SkewScore (Px);
where SkewScore, SkewScore,,() is defined in Definition 7 identified Line in Algorithm iin the infinite sample limit.
This indicates the effect Y can be
Roughly speaking, when À is sufficiently large, similar to high signal-to-noise ratio (SNR), the effect of the latent variable Z can mostly be treated as noise. In this setting, serves as an indicator of how much the observed variables dominate over the latent confounder. The threshold X 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 + Y consistently surpasses that of the reverse direction Y - X. For À < A", this snewness--assee measure may not clearly distinguish the causal direction. However, once À exceeds X", 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.
For the multivariate extension, notable limitation is that the DAG recovery method in Lines 11-17 of Algorithm[] | may lead to spurious direct links when there are latent confounders. A potential way mitigate this issue is to leverage the idea from FCI (Spirtes al. 2001), which will involve more sophisticated ways of performing conditional independence tests (after obtaining the topological ordering) and output the Partial Ancestral Graph (PAG) instead of just a DAG.
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 Consider the following model:
x = = 2+ No, {x
= Z Z f(X)+Ni
where Z, N: {0, 1) are independent variables, all distributed as (0, 1). Then, we have:
Skew
2v2n
(x = =z
(8)
Skewy = O.
7

---

Published as a conference paper at ICLR 2025
The explicit calculation Eq. (8) can be found in Appendix f f is linear function, then Skew 0 and the model is unidentifiable by Algorithm [ However, in most cases, Skew, > 0. making the model identifiable by Algorithm|] [ For example, for a quadartic function f(x) ax² + be + with a + O, Skew, > 0 if and only if + 2- SSiilarly, for almost all p polynomials with a degree greater than? 2, we have 5 > O, indicating that Algorithm ean effectively determine the causal direction in the infinite sample limit.
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) coeet--coffouded triangular HSNM setting (bivariate model with latent confounder, as discussed in Section and (3) multivariate HSNMs. For each scenario, noise is sampled from two types of symmetric distributions: Gaussian distribution and Student's distribution. The degree of freedom of the Student's 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, a, 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 a set to the absolute function (denoted as "Sig-abs"). 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 2 2. The sample size for all experiments is 5000. The score estimator employs sliced score matching (SSM) with 3-layer MLP. We use a batch size of 128 and configure the hidden dimensions to 128 for d < 10 and 512 for d 2 10. Optimization is performed using the Adam optimizer with learning rate of 10-3, and subsample 1000 points for conditional independence testing after obtaining the topological order.
process-sigmoie
Comparisons. We empirically evaluate SkewScore e against state-of -he--hhhrr bivariate causal discov- ery methods designed for heteroscedastic noise models: HECI (Xu al.) [022]; y; LOCI (Immer et al., 2023) and HOST (Duong & Nguyen 2023). Given that Since HOST can also handle multivariate settings, we include in those comparisons. We also test DiffAN (Sanchez al., 2022). which uses the second-order derivative of the log- ligllikliihood, estimated through diffusion models. Additionally, we evaluate NOTEARS (Zheng et al.) 20188), which uses continuous optimization for linear DAG learning, and the MLP version of DAGMA (Bello et al., 2022), 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 2.
Metrics. For the pairwise causal discovery settings, report the accuracy of the estimated causal directions. For multivariate cases, we compute two quantities: the topological order divergence (Rol- land et al. 2022), 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 π and binary ground-trui -ruth adjacency matrix A. the topological order divergence is defined as Drop (T, A) = Ei Ey.,,>,, Aij. Lower values of Drop (T, A) indicate higher accuracy in the estimated topological order, reflecting more precise causal discovery. The structure Hamming distance measures the difference between two 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, antent-coffounde---rinnggii and multivariate scenarios are shown in Figure Bal Figure 3b 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 z distribution. In contrast, methods like NoTEARS, DAGMA, and DiffAN- are not specifically designed for HNMs-fae performance limitations in both pairwise and multivariate scenarios. SkewScore demonstrates consistently strong performance in the ctennt-oofouudee triangular setting, supported by the theoretical insights discussed in Section
8

---

Published as a conference paper at ICLR 2025
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 r-distributed noise effectively, and HOST performs well for the Gaussian noise confounded data is generated using the GP-sig formulation, these HNM baselines do not match the performance of SkewScore. In more challenging multivariate scenarios, method outperforms the others, especially as the dimensionality increases. HOST depends on noise extraction and normality tests, which can be difficult with complex data. DifTAN performs well with small number of variables, but its performance drops dimensionality grows. This is likely due to its non-retraining approach and the score derived from ANM, which may lead error accumulation over iterations. Additionally, all experiments are 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
However, when the latent-
Accurocy
60
Accurscy
60-
SkewScore
SvewSrore
GP-4ig
GP-sig
Sigeats
Signats
GP-sig
GP-sig
Sig-abe
Sig-abs
(0) Hivarinte beteroscedastic noise model.
(b) Latett-coufouuded trinngular heseroscedastie noise model.
Figure 3: Accuracy causal direction estimation across different data generation processes for (a) bivariate heteroscedastic noise models, and (b) latent confouneddd triangular heteroscedastic noise models, Results are averaged over 100 independent runs.
Mehret
Vatoen
Skow
SlowSoore
HOBI
HOST
Topologica divergencce
DIPAN
DAGMAA
Number Vanablee
Numbor Variabies
Figure 4: Topological order divergence and structural Hamming distance (SHD) across different number variables (the dimension d). Lower values indicate better performance for both metrics.
CONCLUSION
In this study, we introduced a novel criterion-the skewness-of of-sssor----roreo 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 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.

---

Published as a conference paper at ICLR 2025
ACKNOWLEDGEMENTS
The research is partially supported by the NSF awards: SCALE MoDI-2134209, CCF-21 16665 (TILOS). is also supported by the U.S. Department Energy, the Office of Science, the Facebook Research Award, as well as CDC-RFA-FT-23-ON Outbreak Analytics. BH is supported by NSF DDMS-242800 We would also like to acknowledge the support from NSF Award No. 2229881, 1, AI Institute for Societal Decision Making (AI-SDM), the National Institutes of Health (NIH) under Contract ROIHL159805, and grants from Quris AI, Florin Court Capital, and MBZUAI-WIS Joint Program. YL thanks Sumanth Varambally for the discussion. YH thanks Xuyang Lin for the discussion.
from the CDC's Center for Forecasting and
REFERENCES
Jeffrey Adams, Niels Hansen, and Kun Zhang. Identification of partially observed linear causal models: Graphical conditions for the non-gaussian and heterogeneous cases. Advances in Neural Information Processing Systems, 34:22822-22833, 2021.
Sina Akbari, Ehsan Mokhtarian, AmirEmad Ghassami, and Negar Kiyavash. Recursive causal structure learning in the presence of latent variables and selection bias. Advances in Neural Information Processing Systems, 94:10119-10133 2021.
Animashree Anandkumar, Daniel Hsu, Adel Javanmard, and Sham Kakade. Learning linear bayesian networks with latent variables. In International Conference on Machine Learning, pp. 249-257.
PMLR, 2013.
Dale N. Anderson and Barry C. Arnold Linnik distributions and processes. Journal of Applied Probability, 30(2):330-344 1993.
Kevin Bello, Bryon Aragam, and Pradeep Ravikumar. Dagma: Learning dags via m-matrices and log-determinant acyclicity characterization, Advances in Neural Information Processing Systems, 35:8226-8233 2022.
Peter Bühlmann, Jonas Peters, and Jan Ernest. Cam: Causal additive models, nigh-dimensional order search and penalized regression. 2014.
Ruichu Cai, Feng Xie, Clark Glymour, Zhifeng Hao, and Kun Zhang. Triad constraints for learning causal structure of latent variables. Advances in neural information processing systems, 32, 2019.
Zhengming Chen, Feng Xie, Jie Qiao, Zhifeng Hao, Kun Zhang, and Ruichu Cai. Identification of linear latent variable model with arbitrary distribution. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 36, pp. 6350-6357,
David M. Chickering. Optimal structure identification with greedy search. Journal of Machine Learning Research, 3(Nov):507-555 2002.
Diego Colombo, Marloes Maathuis, Markus Kalisch, and Thomas S Richardson. Learning high- dimensional directed acyclic graphs with latent and selection variables. The Annals of Statistics, pp. 294-321, 2012.
Eivind Damsleth and A. H. El-Shaarawi. Arma models with double-expooeentiall distributed noise. Journal of the Royal Statistical Society: Series B (tttauddooloica)))((()))) 12 2018.
Xinshuai Dong, Biwei Huang, Ignavier Ng, Xiangchen Song, Yujia Zheng, Songyao Jin, Roberto Legaspi. Peter Spirtes, and Kun Zhang. A versatile causal discovery framework to allow causally- related hidden variables. arXiv preprint arXiv:2312.1100. 2023.
Bao Duong and Thin Nguyen. Heteroscedastic causal structure learning. arXiv preprint arXiv:2307.07973, 2023.
Clark Glymour, Kun Zhang, and Peter Spirtes. Review of causal discovery methods based on graphical models. Frontiers in Genetics, 10, 2019.
David Heckerman, Dan Geiger, and David Chickering, Learning Bayesian networks: The combination of knowledge and statistical data. Machine Learning, 20:197-243, 1995.
10

---

Published as a conference paper at ICLR 2025
Patrik Hoyer, Dominik Janzing, Joris M Mooij. Jonas Peters, and Bernhard Schölkopf. Nonlinear causal discovery with additive noise models. Advances in neural information processing systems, 21, 2008a.
Patrik O Hoyer, Shohei Shimizu, Antti J Kerminen, and Markus Palviainen. Estimation of causal effects using linear non-gaussian causal models with hidden variables. International Journal of Approximate Reasoning, 49(2):362 -378 2008b.
Biwei Huang, Kun Zhang. Yizhu Lin, Bernhard Schölkopf, and Clark Glymour. Generalized score functions for causal discovery. In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, 2018.
Biwei Huang, Charles Jia Han Low, Feng Xie, Clark Glymour, and Kun Zhang. Latent hierarchical causal structure discovery with rank constraints. Advances in Neural Information Processing Systems, 35:5549-555 5561,,55122222 2022.
Aapo Hyvärinen and Peter Dayan. Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research, 6(4), 2005.
Alexander Immer, Christoph Schultheiss, Julia E Vogt, Bernhard Schölkopf, Peter Bühlmann, and Alexander Marx. On the identifiability and estimation of causal location-scale noise models. In International Conference on Machine Learning, pp. 14316-14336 PMLR, 2023.
Kristian Kersting, Christian Plagemann, Patrick Pfaff, and Wolfram Burgard. Most likely het- croscedastic Gaussian process regression. In International Conference Machine Learning, 2007.
llees Khemakhem, Ricardo Monti, Robert Leech, and Aapo Hyvarinen. Causal autoregressive flows. In International conference on artificial intelligence and statistics, pp. 3520-3528. PMLR, 2021.
Bohdan Kivva, Goutham Rajendran, Pradeep Ravikumar, and Bryon Aragam. Learning latent causal graphs via mixture oracles. Advances in Neural Information Processing Systems, 34:1801 18:1178811110, 2021.
Erich Kummerfeld and Joseph Ramsey, Causal clustering for 1-factor measurement models. In Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining, pp. 1655-1664, 2016.
Yingzhen Li and Richard Turner. Gradient estimators for implicit models. 17Xx007705.07107 2017.
arXiv preprint
William McGill. Random fluctuations of response rate. Psychometrika, 27(1):3-17, 1962.
Bruno Merz, Günter Blöschl, Sergiy Vorogushyn, Francesco Dottori, Jeroen C.J.H. Aerts, Paul D. Bates, Miriam Bertola, Matthias Kemter, Heidi Kreibich, Upmanu Lall, and Elena Macdonald. Causes, impacts and patterns of disastrous river floods. Nature Reviews Earth & Environment, 2021.
Jovana Mitrovic, Dino Sejdinovic, and Yee Whye Teh. Causal inference via kernel deviance measures. Advances in neural information processing systems, 31, 2018.
Boris Mityagin. The zero set of real analytic function, 2015. URLBEttps:/arxcvv.orog///
1512.07273
Francesco Montagna, Atalanti Mastakouri, Elias Eulig, Nicoletta Noceti, Lorenzo Rosasco, Dominik Janzing, Bryon Aragam, and Francesco Locatello. Assumption violations in causal discovery and the robustness of score matching. In Advances Neural Information Processing Systems, 2023a.
Francesco Montagna, Nicoletta Noceti, Lorenzo Rosasco, Kun Zhang, and Francesco Locatello. Causal discovery with score matching additive models with arbitrary noise. In Conference Causal Learning and Reasoning, pp. 726-751. PMLR, 2023b.
Francesco Montagna, Nicoletta Noceti, Lorenzo Rosasco, Kun Zhang, and Francesco Locatello. Scalable causal discovery with score matching. arXiv preprint arXiv:2304.00 03382, 2023c.
11

---

Published as a conference paper at ICLR 2025
Ignavier Ng. Biwei Huang. and Kun Zhang. Structure learning with continuous optimization: A sober look and beyond. In Proceedings of the Third Conference on Causal Learning and Reasoning. 2024.
Judea Pearl. Causality: Models, Reasoning and Inference. Cambridge University Press, 2009.
Alexander Reisach, Christof Seiler, and Sebastian Weichwald. Beware of the simulated dag! causal discovery benchmarks may be easy game. Advances in Neural Information Processing Systems, 1-227722227844 2021.
Paul Rolland, Volkan Cevher, Matthäus Kleindessner, Chris Russell, Dominik Janzing, Bernhard Schölkopf, and Francesco Locatello. Score matching enables causal discovery of nonlinear additive noise models. In International Conference on Machine Learning, pp. 18741-18753, PMLR, 2022.
Runge, S. Bathiany, E. Bollt, G. Camps-Valls, D. Coumou, E. Deyle, C. Glymour, M. Kretschmer, M. Mahecha, J. Muñioz-Mari, E. V. van Nes, J. Peters, R. Quax, M. Reichstein, M. Scheffer, B. Schölkopf, P. Spirtes, G. Sugihara, Jie Sun, Kun Zhang, and J. Zscheischler. Inferring causation from time series in earth system sciences. Nature Communicationss 10, 2019.
Karen Sachs, Omar Perez, Dana Pe'er, er, Douglas Lauffenburger, and Garry P Nolan. Causal protein- signaling networks derived from multiparameter single-cell data. Science, 308(5721):553--29, 2005.
Saber Salehkaleybar, AmirEmad Ghassami, Negar Kiyavash, and Kun Zhang. Learning linear non- gaussian causal models in the presence of latent variables. Journal of Machine Learning Research, 21(39):1-24, 2020.
Pedro Sanchez, Xiao Liu, Alison Q O'Neil, and Sotirios A Tsaftaris. Diffusion models for causal discovery via topological ordering, arXiv preprini (...iv:22100062001 2022.
Pedro Sanchez, Xiao Liu, Alison Q O'Neil and Sotirios A. Tsaftaris. Diffusson models for causal discovery via topological ordering. In The Eleventh International Conference on Learning Repre- sentations, 2023. URL httpa//sseepaaeeeeeeeeeeeeeeeeea
wet.eewwwwwwiieeeeeeeoo.e.eeeeee n www.:
Jiaxin Shi, Shengyang Sun, and Jun Zhu. spectral approach to gradient estimation for implicit distributions. In International Conference on Machine Learning, pp.
PMLR, 2018.
Shohei Shimizu, Patrik O Hoyer, Aapo Hyvärinen, and Antti Kerminen. A linear non-Gaussian acyclic model for causal discovery. Journal of Machine Learning Research, 7(0e:)2003,,0000, 2006.
Shohei Shimizu, Patrik O Hoyer, and Aapo Hyvärinen. Estimation linear non-gaussian acyclic models for latent factors. Neurocomputting, 72(7-9):2024-2000 2027, 2009.
Ricardo Silva, Richard Scheines, Clark Glymour, Peter Spirtes, and David Maxwell Chickering. Learning the structure of linear latent variable models. Journal of Machine Learning Research, 7 (2), 2006.
Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. Advances in neural information processing systems, 32, 2019.
Yang Song, Sahaj Garg, Jiaxin Shi, and Stefano Ermon. Sliced score matching: A scalable approach to density and score estimation. In Uncertainty in Artificial Intelligence, pp. 574-584. PMLR, 2020.
Peter Spirtes and Clark Glymour. An algorithm for fast recovery of sparse causal graphs. Social Science Computer Review, 9:62-72, 1991.
Peter Spirtes, Clark Glymour, and Richard Scheines. Causation, prediction, and search. MIT press, 2001.
Eric Strobl and Thomas A Lasko. Identifying satient-specific root causes with the heteroscedastic noise model. Journal of Computational Science, 72: 102099, 2023.
12

---

Published as a conference paper at ICLR 2025
Natasa Tagasovska, Valérie Chavez-Demoulin, and Thibault Vatter. Distinguishing cause from effect using quantiles: Bivariate quantile causal discovery. In International Conference on Machine Learning, pp. 9311-9322 PMLR, 2020.
Marc Teyssier and Daphne Koller. Ordering- based search: simple and effective algorithm for learning bayesian networks. In Proceedings of the Twenty-First Conference on Uncertainty in Artificial Intelligence, pp. 584-590, 2005.
Paseal Vincent. A connection between score matching and denoising autoencoders. Neural computa- lion, 23(7):1611-1674, 2011.
Feng Xie, Ruichu Cai, Biwei Huang, Clark Glymour, Zhifeng Hao, and Kun Zhang. Generalized independent noise condition for estimating latent variable causal graphs. Advances in neural information processing systems, 15:14891-44902 2020.
Sascha Xu, Alexander Marx, Osman Mian, and Jilles Vreeken. Causal inference with heteroscedastic noise models. In Proceedings of the AAAI Workshop on Information Theoretic Causal Inference and Discovery, 2022.
Naiyu Yin, Tian Gao, Yue Yu, and Qiang Effective causal discovery under identifiable heteroscedas- tic noise model. In Proceedings of the AAAI Conference Artificial Intelligence, volume 38, pp. 16486-16494, 2024.
Zhang, Peters, D Janzing, and B Schölkopf. Kernel-based conditional independence test and application in causal discovery. In 27th Conference on Uncertainty in Artificial Intelligence (UAI 2011), pp. 804-813. AUAI Press, 2011.
Kun Zhang and Aapo Hyvärimen. On the identifiability of the post-nonlinear causal model. In Conference on Uncertainty Artificial Intelligence, 2009.
Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric Xing. Dags with no tears: Continuous optimization for structure learning. Advances in neural information processing systems, 31, 2018.
Yuhao Zhou, Jiaxin Shi, and Jun Zhu. Nonparametric score estimators. In International Conference on Machine Learning, pp. 11113-11522 PMLR, 2020.
13

---

Published as a conference paper at ICLR 2025
A ADDITIONAL DISCUSSION ON THE SkewScore MEAUSURE
This section discusses SkewScore (Definition including intuition and motivation for considering this measure.
the conditional density p(æ; x==s) is a symmetric function, then we have SkewScore,, (p) = 0.
The proof of this result is deferred the end of this section. For asymmetric distributions, we provide two one-dimensional examples with explicit computations at the end of this section:
1. For the Gumbel distribution (2)
exp (- "
(-5" c x
with B > we have
StewScore(PGummee) gs.
2. For the Gamma distribution [[k, 0), with k > 3, and > 0, i.c.
TTKjOr -
120 0
PGamma (2))
x <<0 0,k>3),
we have SkewScore(ppammm) 0-3
We provide additional discussion the motivation for considering the skewness of the score, i.e., the SkewScore measure. Consider the bivariate HSNM model where the variable X causes Y, i.e. X + Y. Figure shows that to identify the causal direction, rather than examining the skewness of the residuals, we can focus on the skewness of the conditional distribution. However, working directly with the conditional random variables X and Y|X is often challenging. Therefore, introduce an auxiliary random variable. V log p(X, Y), where (X, Y) follows the data distribution p. This wwo-dimensional random variable, induced by the score function on the data distribution, is represented as the black directed line segments in Figure 5 This auxiliary random vector's s first element, 0₂ log p(X, Y). is visualized as the orange directed line segment, while the second element, 0, log p(X, Y). ). is shown in blue. Fact I connects the conditional distribution p(y(ær and p(zlyy with this auxiliary random vector, allowing us to transform the conditional- cintonalddisstiuttonessstte to joint-distribution one. [
problem that under Assumption IT the orange directed segment (representing O2 log p(X, Y)) follows a skewed distribution, whereas the blue directed segment (representing 0, log pl X, Y)) follows symmetric distribution. This means that SkewScore, # 0 while ScowScore, (p) = O. We summarize the motivation in Figure
points (X,)
Score function: V Voog Vlog)(t)))
projection log plt,) Score's projection on d, log p(x.)
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 Texamines the skewness of these projections across the data distribution.
Skewness of the Con- ditional Distributions
Skewness of the score V log p(X, Y): SkewScore, (p) and SkewScore, (P) (Our approach)
Skewness of the Exogenous Noise N (Difficult to extract N)
p(zly), p (r) (Challenging to pirically estimate)
Figure 6: Conceptual flow from skewness of exogenous noise skewness of the score.
We present the deferred proof and derivation as follows.
14

---

Published as a conference paper at ICLR 2025
Proof of: SkewScorez, (p) = 0 for symmetric p(z;; x=-)). By Definition 1 SkewScores..p) [x- (00,
we have
(ogg(())))
1 / JR p(11, Id)
{ [rr,
ITa) ,3a)) dr,dr-a
The assumption that p(z; xi-)) is symmetric implies: for every E R
p( - Fia"" :* ,d) - p(x1, ... Fz T.).
Therefore,
O2, p(21, 2i, rd) 0(11( 2 Fd)
Now we fix Ø- E 12-1, and let h(z;)) play, ,a ma) derive that h(x;)) -h(-.;).. So is odd function and thus
(0,p( sa From the previous analysis, we
f h(z;)dz; O.
Consequently, for every x-; E TR
10,,(111 ra)] dx; = 0. p(e1 ,z))
0
log 1og,(X1,Xp) X2, xa) = Jau
8m, p(11, X2, : (d)
E
p(11.122 02. " ,ra)de_
Ox,
p( I2, "2."" Id]
Ra-1 [
p(x1 22, rg)dd,d ",a)d,d
J Oxi
Thus we prove SkewScorea, (p) O.
Derivation of SkewScore(pGGmtell Recall that p(x) exp
3
-
Then logp(z)
log }
and d log p(z) = = de
(-b)
(-1+e 30"
3
SdewScorrePGaabee
/
1 +
3 )
3 exp
1
1
c
-t
- dr
3
" 1
(-1+e (=tee=))
dz
g3
R
-=
1
(-1 y) ddd..
03
Finally, since So y = r(k + = kl, for l E N, we conclude that SkewScoreep)) = g
Derivation of SkewScore (teeScoorePcamma.) (PGamma). We first note that k > 3 implies that p(z) E C² and ensures the integrals involved are all well-defined. We have log p(x) = (k- - 1) In x $ + Im [(k)a= and thus
de logp(z) 1-1 J. Then we have
(
d log p(x)))
SkewScore(p) E
dx
1
[( -l
e dx
J0
1)
1
1
(k: 1)30 -3 - -3(k- 1)22 0
3(k 1)
dxr
[(k)ØR
82
03
0
(k
1) 3) - -3k 1)2 2) + 3(k - 1)[ 1) T(k)
[(
15

---

Published as a conference paper at ICLR 2025
By the property of Gamma function [(2 + 1) = =2] we have
E
logp(x)) dr
r(k)Ok 1(k))k-3
1
(k: 1)3
3(k: 1)²
+
3(k 1)
(ke - 1)(k - 2)( 3)
(k: 1) - 2)
k
0 4
(k - (): - 2)
B DISCUSSION ON MODEL ASSUMPTION 1
The following proposition demonstrates that Assumption is not overly restrictive and can generally be satisfied; specifically, the set of models that violate this assumption relatively small. To establish this, we fix Pa. Pri and (zz), and consider f within finite-dimensional space with sufficiently large dimension to capture the features of the data. We then demonstrate that the set of functions f() that do not satisfy the assumption has zero measure.
Proposition 5. Fix Pa. Pre and (x). Consider f within an M- -d-dmensiioal linear space V. this space V. define the non-compliant
Suc(V) = (SEV: f does not satisfy Assumption 7
With the notation in Assumption define
p.(())² du) A(z)p=(x) + (JR Ea du) h(x): := =3 (JR p.(()) a())²
C(r)pa((t) a(x))
c: (A(x) C(x) (
(u (w)u)
Pn (())p(())ddds
Pa(u))
which do not involve f.
Since dimV = M, V is isomorphic to R and thus naturally inherits the Lebesgue measure on RR This induced measure does not depend on the choice of basis and is denoted by Assume there exists 2o E I such that h(xr)) * O, then the following holds.
1. If + O, then the non-compliant set has zero measure, i.e.,
m(Suc(V)) =0
2. If O, assume V is the space of polynomial functions of degree most M 1, iee E {R.. Assume JR (p=(x))² + + (x)22 exo² dx < oc. for some μ E (0, +00). Then there exists My E N, such that for every M 2 Mo, the non-compliant set has zero measure, i.e..
V = E M-1 axxt: 00,
m(Suc
This proposition indicates that Assumption 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
We note that when (x)) is not constant (Le., in the heteroscedastic model), e + 0 in most cases. However, if (2)) is constant and ps is symmetric, then In this scenario, if X is Gaussian variable, the assumption in the second case is always satisfied; that is, there exists some
Remark 6.
(0, +00) such that Sg (x)) + (2:)²
eux² da < cc.
For the multivariate case, intuitively, the general validity of the assumption in Corollary i similar to Proposition 5] where the set of functions (fi, o;) that violate Eq. (6) exhibits property analogous to having Lebesgue zero measure. Specifically, we parameterize the causal model by parameterizing functions fe/ and of. where 0 belongs to the parameter space θ C R (e.g.. neural network; see Example 7). The left-hand side of Eq. 6 defines a function of the parameter 6, and we write it as F(0). With this parametrization, the model that violates Eq. (6) corresponds to the solution of the equation F(0)= = O. Under regular assumptions on ft. di Px. and pN,, F is a analytiayytic function.
Then by applying (Mityagin 2015 Proposition 0), the violated set (Q) (0 E O F(0) = 0)
16

---

Published as a conference paper at ICLR 2025
has zero measure as long the function F is not constant function. is worth noting, however, that in the homoscedastic Gaussian linear model, is identically zero, and thus the proposition from Mityagin 2015) no longer applies.
We provide the following examples to give more intuitive understanding of Assumption T
Example 7. In this example we consider Homoscedastic model, where (x) is constant. Let f be a two-layer Neural Network with Softplus activation function:
f(z) w2 log(1 + c(wratbi) + b2, W1, W2 b1, b2 E IR
Assuming Pa is an unimodal distribution symmetric with respect to O, then Assumption 7 is violated if and only w, or W2 0, i.e, when is constant function.
Proof The LHS of Assumption T eqqual to
Swfwz
JR (1 + ty (o)de /
(F'x(u))² PN(u) du.
Notice that (
Ja (Ite Con -- Pa (r) # when wi # O. (xx(t)) du > 0. We only need to prove pv(a)
(a) is always non-negative, so S
Notice that (z( is
function thus (x) is an odd function, thus
P da
(1 + (wee-wez P/ (x) + S (1 + -w(w
(x)
p'((x)
(1 + e (w)
Ps (x)de
(1+ (w
(1 (1 -1
00
(2 + e (wix++) + e (-ww-)) (e - W1()
(1 - + (-
Since * 0, and P. is an unimodal distribution symmetric with respect to 0. The integrand is always of the same sign (always strictly greater than strictly less than 0) over the integration interval (0, co). Therefore, the above integral is non-zero.
Example 8. In this example we consider the Heteroscedastic case with axx) = (=== where 1 + O, X > Let f(x)) ax b, X \( 1), N ~ (O, Then f violates Assumption 7 and only if
where 7 = 31d
and c is defined in Proposition 5 Especially, when
(4X-1)v44+1 (6A+1)/6X41
0, always satisfy Assumption
evident that the non-compliant set Snc(V) = - x R always has
measure (it is an empty set when < 0), providing an example for proposition 5 Additionally, due to the simple structure, we already have m(SuccV)) 0 when = 0. Therefore, the discussion of the second case Proposition iss unnecessary. The proof of this example is provided as follows.
Proof Recall in the proof of proposition to
(Eq.
16
we proved that f violates Assumption I Jeuuls
(())) -c.
As f(x) a. is linear, it equals
2 M()dd -c &
17

---

Published as a conference paper at ICLR 2025
On the other hand, we can simplify h(z) as follows:
(x)
9px(x)o(x)
3x
4Mc 2
18A(x 11)
GACe 3
(xx) = xx)² a
V
2
Thus we can compute that
J h(x)d==
4e 104
6e 200
3rd =3p
(4h+ 1)
(6A+ + 1V66++1
Then we know that Ja h(x)dr # O. Therefore, we conclude that / violates Assumption if and only if
a² = Ju h(x)dx
which equals to
+
In the following example, we discuss a case where f does not satisfy Assumption T Example 9. Let (xx) be a constant and be linear, then I does satisfy Assumption if and only if:
J (1(()) P, dx 0.
When X follows symmetric distribution, Eq. (9) holds. Specifically, when X follows a Gaussian distribution, Eq. (9) holds, implying that f fails to satisfy Assumption [ This indicates that the Gaussian Linear model is unidentifiable by our algorithm, consistent with previous research such as (Hoyer et al., 2008a). However, as illustrated in example if (xx) 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
| Notice that when o(z) =
a constant, (x) aæ + bs linear, C(x) = O, B(x) -.. A(x)
(x(()) So f does satisfy
P=(z)*
Assumption equals
[
(x)3 dr + (x)
3a² /
_ / pl((x)dx== = =0 0.
o² J Pu(
Notice that P/y (x) = 0, it then equals to
(x)33 J
J p=)
de = 0,
thus completing the proof.
C DISCUSSION ON THE SKEWNESS DEFINITION
To enhance the flexibility of our algorithm, we redefine the skewness of a random variable W as Skeww(W) (W - [, employing nonlinear odd test function that satisfies (0(s) -s)s). Specifically, in Theorem the chosen function is v0(s) = 3 Note that for any symmetric random variable W. Skew
0 for any odd function dh, thus upholding the validity Eq. 44)
under this generalized definition. To ensure that gi-analog of Eq. (5) also holds in the anti-causal direction, it is essential to verify the vi-analog of Assumption [ A necessary condition for & is nonlinearity due to Fact 2 The choice of the odd test function affects the robustness of the algorithm: an effective odd test function captures the asymmetry of the random variable while requiring less sample complexity.
18

---

Published as a conference paper at ICLR 2025
D PROOFS OF SUPPORTING LEMMAS AND FACTS
Proof of Fact 1 Without loss of generality, it suffices to provide the proof when p E C² Notice that the derivative of the log conditional density translates to the derivative of the log joint distribution, then
0 dx:
1(211,22,.-- Idd
log p())
log
0z,
p(x-s)
0
,dd) ((,,
log p(x
0z,
0
log 1o "","",,
0x,
Proof of Fact 2
8 Or,
Du. 1 :. (4)
E
logp(Xn, X2, Xx)
Jzu p(1 12,-1 add
p( 22, sa)dd,d
Jz J Oz.
Since p is continuously differentiable, we apply the Newton-Leibniz formula to Ty. For any (x1, Ti- 1,41, ... Ta) Rd-1, we have
S JD p( T za)dæ; p(x1, (
dx,
where limz;-+too p(11,". .. .a) lim p(x ... Te ....... T.) 0 because p is probability density. Therefore, have shown that E
at, log pl (X,, _2, Xa)
We also note that the assumption (x(x)= e C ("", can relaxed. For example, Fact holds for uniform distribution 1-1, despite having singularities at -1 and 1, the densities at these two endpoints cancel out after the Newton-Leibniz step. More generally, Fact A holds when p is piecewise continuously differentiable with singularity points w1 102
Wm, and
>n [p(w;+) - p( O.
Attentive readers may observe, however, that for the exponential distribution (Pexp(x) = de- Ar for 2 and 0 otherwise). Fact 2 does not hold. This is because of its singularity at x. O, where Pexp (0-) = O and Pexp (0+)*
E PROOF OF THEOREM
Assumption (Regularity assumption for Pay Pr.) Pro Pe C.
1. Jn. [rr, (1))* (u)
du < oo,
2. Jx
(2)3 p=( dr co. J&
[( va(a))
du < co.
Assumption (Regularity assumption for f). EEC!
J [1())) < xc.
Remark 10. Assumption and ensures all the integrals we discuss here are well-defned.
Proof of Theorem We first prove Eq.
4 As Y = f(X) + o(X)N the density of (X, Y) should
be
(x(x)
y yy-f(a)
plx, v) = Pa (x)p((()))
o(x)
(xx)
19

---

Published as a conference paper at ICLR 2025
Phn (リ-((a)
Ex, (0, logp(XXX)))
/ (x(x)
o(x)
drdy
Jree (x))
ph (-f(())
(x(x)
(p/ (u))
drdu
R2 (x)3
Ph (w)
(x(
(F/(
a { dx
du
(u(
For Eq. 55 we first expand it explicitly:
{
(z)
d'
(u(()
f((z)
d'
Exy (0, logp(X
Pr (a
R2 (x) (xx)
Pa (n(u)
oxx)
(x))
(10)
f po(x)) / A(x) + B(x) (u)
P' (u)u
+ ((x)
Pn (u)duudx.
(u(
pn
(11)
As Pn (u) (n), derive f& (u(2²² o(u) and C(x) o (x) oxz)² as is denoted Assumption
where A(x) (x) d'(x) B(x)= -(z)=-f(a
(w) o(w)
Pu(u)zdu Sg P/a (u)du 0,
expand Eq. II and
get
Ex,y (0, log = Pa( / [[((( [((x) + C(x) (n)
(A(o)
dæ
(12)
+3 3 B(x) Khlw))
AA(x) + C(x) 1' Ph(uuj (n(
pa pu(uudd..
(13)
Pm Pn(u))
By assumption, Exy (0, loggll,,,,.... Y * O, thus completing the proof.
F PROOF OF COROLLARY
Proof of corollary 2 As in prior work, assume all integrals below well-defined. We first compute that: for every k E (1,2,-.. ......... d}.
P's
Pl.
0,00 fo + (x: fi)O2a
Oey
2n log p(x)
Σ
(14)
O&Pk
Pi
o7
igchitd(k)
It suffices to notice that the joint distribution of X can be written
d
p(z)= = II
(z. (x),)))
1=1
lox)(x)
log (z; pxa()))
i=1
a fi (pa, (z) σ; (p),(x))
log
PN.
log o. (pa,x)))
We then differentiate with respect to 2y. and directly derive Eq. (14)
We now prepared prove the first assertion. Considering sink Xx. know that child(k) 0. Therefore, by the formula in the first line of the proof:
P1
O, log p(x)
OPPk
20

---

Published as a conference paper at ICLR 2025
Then compute:
P'v
Mh fe-Fe (2))
3
Ex (dng log p(X))
ox(ppa(0))
p (z))pa, dx
@o (2)) P.ve
Dy 34 (pa=f((pa(2))) (21) ox(paa(())
i=1
(P'v
/ ox (pax(()))
(x)a
(pa
II p (ar pa, (x)) de iffh
(a-Salpa, (x(pp,(()))
=f p(pa (x))
(pw, (u))
3
Je 0% (pa (x))3
du
II (x.laaa) dx,)
Piv. (u)
(x))
II (o (z//pa,(x)) dx;) / (P'v, (u)) Phx (u) du
Rd-1 (x(pa (x))3
17k
0.
(Fyy, (u))² px,(" (w) du 0. More explicitly, due to the assumption that
We obtain the last equality because Jn P'N, (u) is symmetric,
(N'v (w)) PR (M)
is a odd function and thus its integral is O.
We then prove the second assertion. Recall that
Ozy log p(x) = P/s
plc o,02, fo + (x, - f.)0
Σ
oā
T&Pk
Pi
σi
(Echild(R)
Therefore,
Ex (0,, logp(X))²
Pk
Pl o, (z+ f.)d2. 01
Ou
13
Pj dx.
Rz
TkPk
Pi
o²
(Echild(k)
Thanks to our assumption co. we derive that:
Ex (0,, log)(X))² 40,
thus completing the proof.
G PROOF OF PROPOSITION 5
Lemma 11. J (Pa(1)² + +z()) gx dr < 00 for some E +0o), recalling
M(x): =3 (1 ph (u)2 (u)
A(x)pa(T)
+3 ( // (u) du
((x)px())
du
a(r)²
Pn(u)² o())²
then have
[ h(x)2e dx <
Proof As A(x)
(x) d'(z) (z))) x(x) and C(z) = d'(x)
a(a) and or are bounded and (xx) 2 r > 0,
there exists a constant s.t. AA()) S
(x)
+ b. [(x)] < b. Notice that fR (l(x))
du and Pa(a]
palx)
fx pl, (u)"M du are all constant, then there exists constant e > 0 s.t.
[Mxx) (px (x) (z))
Therefore
[Mxx))2 S 20 (xx(2)² (x)2 +1
+ f h(x))e e dx < / (po(x)²+1'((1)) (x)))) e0s dx < 00.
21

---

Published as a conference paper at ICLR 2025
Defining
et = he LR( ( ((x) do < co).
lemma 11 equals to saying that h E (2 cn=²). The following lemma is interesting as the integral is taken over I instead of finite interval, and thus the Stone Weierstrass theorem cannot be applied.
Lemma 12. For h E (R, el if x²x(d) = 0 for every n E N. then h 0 almost everywhere.
Proof Without loss of generality, we assume μ = well known that L²RR, e is Hilbert space. Then by knowledge of Hermite polynomials, there exist series of polynomials {H,),
which is a orthogonal basis L (R, e!Te) ). fa =2))ddd = 0 for every n E N implies that S& Hn. (z))(zz)) = 0 for every n E N, thus proving = 0.
Proof of Proposition For every finite dimensional subspace V Vo. Snc equals to (()))( ( pa(a)
0 0=
(((z) + EL(N 3 Pa(M)) (A)
+ B((r) (2(w) ) (A(x) + C(x)
Pn Pn(u)po (x))ddu.
(15)
In the above equation, only (xx) involves (x)) So we can simplified Eq.
15 into
J (x)) (
(16)
where
Fl,()) A(x)pa(s) h(x) =3 (Jn
Po(u) du) du)
(' du) (x) (Ja
o(2))
Pataj a(x)²
3
c= ( (A(x)) + C(x)
Pa(u)
Pm
(n(p (.)))((()dxd))
Let fa (1 sis d) be basis of V. = E, a: then Eq. T6 equals to
Σ aayy ISi
/ [(o0)((M))) +c=
Let Jy(a1, a,,) = Eisissa 0,a, (()) + c. J is polynomial and thus a real analytic function. Define KerJy = ((11, ... ,a.) Ju a2, .2, ,2,",,,,n) = 0}, then
f (vc = (a1, a2, "............ E KerJy.
We first consider the second case where fa (px (x))) eura dx for some μ E (0, +00). By lemma 11 Ja 2el) < oc. If f x" ((x)dr= for every n N. then by lemma 12
almost everywhere. Notice that also continuous, implying that h(x) contradicts with our assumption on h. So we know that there exist Mo. s.t. fa x. (oMoh(x)rr * 0. Therefore, for M 2 Mo + 1, if we take V to be the space of polynomials with degree smaller than M and fi(z) x² then Jv is a non-zero function. By basic knowledge of real analysis Mityagin Propossionn 0j, KerJy has zero measure. So (n(S = m m(KerJy) =0
for every
R. This
For the first case, if / 0, then (0, 0) / 0. So Jv is a non-zero function for every finite dimensional space V C (RR) Similar to the analysis above, we derive m(SnccV)) = and thus completing the proof.
H PROOF OF PROPOSITION 3
Let
1(x, u) = / -(())o((=)0(()))) (x)(=))((x))-)) ("-
dz.
a(x)
22

---

Published as a conference paper at ICLR 2025
Under assumption in Formally, we denote
5
it is easy to verify that h is well-defined and belong C (R?) L
A3 = - Sza BuM,)e dxdv,
h(a.u)²
A2 = 3 sRa ['()) (e drdv,
k(a,1)
A -3 Sez
drdv,
h2,0)
Ao SE O.klo,)) dxde. M(z0))²
The well-defineness of A, (0 I i S 3) will be verified in the following lemma.
Assumption (Restriction on Model with Latent Confounder).
1. (Regularity assumption.) It exists an open subset U C +oo) s.t. for i E {x,y) and every
AE Ex,y (0, Po Y))* < 0c.
2. A7 + 43 + A3 #
Lemma 13. Under assumption 4
A; (i 0, 1,2,3) all well-defined integral and Ex,y (0, log Y))
AgA* + A2X² + A1 + Ag.
Proof
/ 22
E (02 (X,
(0, (ogaat)))) Px(E, vddrdy
(x)))
P)(x,U)2 dady
(0,h(a, - (f(())0,((z, v))² u)²
dado
Consequently we can write the integrand as:
(0.h(a, -)(()),,(,2)))))
Σ ()(
h(I,0)2
where
(2 v) ((
h(2,0)2
(2 v) Pu.hiei
A(e,0)z
91 (x)v) -33 M(z.v)
h(a,)))
(0 = 2ch(a.0)
b(a,o))
As h(, v) > 0, are continuous. Taking d E C d + X₂ for + j. By assumption 4 for < & < Szz 9f u) dxdv < oo, i.. Σ X E L Notice that the matrix A = (A.))ss.ke
such that A₁ X is Vandermonde matrix and thus invertible. Then since L (IR2) is a linear space, we conclude that gi E L' (R2) for < < 3. Consequently A; Jg (xx dxdd is well-defined and
Ex, (8, log A2X3 + A2X² + AiA An-
Proof of Proposition
p(z, y, 2) p())p(xl=))oull,, =x(=)go(x (co))(=)))))) (リ- (リ~Aʃ(x) (1(2))
Af(x))
o(x)
23

---

Published as a conference paper at ICLR 2025
We write (1 v) p(a, y) for short, therefore,
p(x. v) p(=)go(x do(2))o(z)) 21 y (V-A((x)= - 2x(33) dz.
o(x)
Let (x, u) = p(=))(()(())) - ( 'n1 (=)
dz. Notice that h is independent of X, f
and h E (RR)( under our assumption. As h(T, - {f(x)) p(a,, v), we derive:
@,p(a.y) a, Xf(z)),
0.p(r.y) 0,h(z, xf())) (
Our first observation is that: Ex, (0g log Y))² is independent of A and f. is shown below:
Ex,y (0, logglpp (X,)) 8))²
(0, log u))² p(x, v)drdy
R²
(Oypla, v\) drdy
J/2
p(a,
(0,plz, xf(=))))) dxdv
(By change of variables)
p(a,u (x)))
R2
(0.1(z, v))
drdv. 1)z,0)2
On the other hand, by lemma 13 we get:
Ex,y (0, log )))² = AgA* + A22² + A, Aol-
Consequently, by assumption
Ex,y (0g log constant s.t. for every 2
+ + as A + too, thus there exists
Ex.y (0, (ogg)))(X
ExY (0,, log PXx
Remark 14. In particular, if As + 0, let R = Exy (a,, log p( (XX, X))
we can set X = 1
max
3|A21 Agl
3|4 TAs
3.Agl +R Aa
J. 11 is easy verify that
[A X |2 X2 - |ool > R,
when X 2 A
I EXPLICIT COMPUTATION OF EXAMPLE
p(a, v.2) ~ exp (- x² (a-2)² ly 1y-z-f(a)3)
2
p(x, = p(I, y. p(I, 2)dz Cexp
(x² + (y f(x))² -x(y - f(x)))
JR
3
2
( T. - 1(x))
x²
= Ceep
3
2
2
C' I is normalization constant. Therefore,
v3x
Oy log p(æ, y)
3 (y = - f(x))) (y 5- f())) ( + fx))) -
O2 log (gg((r.))
(1)
- -f(x) - 2x - xxf((z) 2f(x)f(x)))
(a(x) -(x)),
24

---

Published as a conference paper at ICLR 2025
where a<x) 1 + (z) and 3(x) f(x)
+ 2x +
64C
/
/ (y
((z)))c-
dy e dz dx dy
Skew,
27
R
64C
uee fu² du e & dx =.
27
R
On the other hand,
Skew = 27 [
8C||
[ (a(x))-
(((--))(() + dy -
Let (2x) = + + f(x). Notice that a(x)(()) - 3(x) =
3 2
/ R (((x)y-B(x)) ( dy / (a(z)(++ B(x))* e +tddy
8e--3v' by - 8(0))3 / c -33
= (()) -
9 -=aa(s))
/ 32 fv dy
27 x3 [ E 3v dy.
Let μ fu e2= 32 dy 3vddy
3/65 SR e- 322² dy.
8C
/
27p0 nnr)
8 e 4 do
Skew,
na(x))
27
4Cp
xa(x) 2dx
3
4Cp
x(1 2x)) (
3
16Cp
x((((x)) + (())))e dr
2v
/ x(f((x) + (x) 2 dr
25

---

Published as a conference paper at ICLR 2025
J ADDITIONAL EXPERIMENTS
Runtime Table 1 presents the runtime, in seconds, for SkewScore re using Sliced Score Matching (SSM) (Song al. 2020), for various variable counts d = (10, 20, 50], where the number of edges twice the number of variables. All experiments are on 12 CPU cores with 24 GB RAM.
Table I: Wall-clock runtime (in seconds) of SkewScore ER2 graph with 1000 samples.
Method
d = 10
d = 20
= 50
SkewScore (SSM) 310.24 + 7.17s 704 = 19.82s 1825.92 + 39.23s
Effect of Non-Symmetric Noise We conducted experiments adhering to the same bivariate HSNM data generation process outlined in Section except the noise now follows Gumbel distribution. Our method remains relatively robust to Gumbel noise, skewed distribution. The results in the table 2 demonstrate that the reduction in accuracy for SkewScore 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. However, these two methods are designed for bivariate models.
Table Accuracies of estimated causal relationships the bivariate with non-symmetric Gumbel noise. The results are averaged over 100 independent runs.
Methods
HNMs (Gumbel)
SkewScore
89
LOCI
100
HECI
HOST
33
DiffAN
DAGMA
27
NoTears
Performance on Additive Noise Model (ANM) Though SkewScore is designed for heteroscedas- tic noise, it is still valid on nonlinear models with homoscedastic noise. We explore the algorithm's '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 pruming (Bühlmann et al. 2014) is applied. We calculate the topological order divergence and the structural Hamming distance (SHD) across varying dimensions, focusing graphs with d and 2d edges, represented as ERI and ER2 graphs, respectively. The results, presented Tables 3 and 4 indicate that SkewScore maintains adequate performance ANM settings. Results are averaged over 10 independent runs.
Table 3: Topological order divergence across dimensionalities (d E {10, 20, 50}) using data generated from ANMs.
ER1
FR2
Methods
d=10
d=20
d=50
d=10
d=20
d=50
SkewScore
1.00 + 1.00
3.10 + 0.70
8.50 + 3.17
3.80 + 2.44
+1.00
23,40 + 5.82
HOST
3.90 + 1.45
6.40
15.00 = 3.74
11.70 2.45
51.60 2..70
DiffAN
1.80 0.98
4.50 3.07
14.89 2.18
2.70 + 1.19
7.00
2.32
23.43 1: 2.50
DAGMAMLP
3.70 1.10
7.70 + 1.79
15.70
4.65
10.20 + 0.75
17.10
35.00 2.37
Effeet 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 ER graph model. The results indicate that our method, specifically SkewScore, demonstrates enhanced performance the sample size increases. This trend suggests that SkewScore efficiently leverage larger
26

---

Published as a conference paper at ICLR 2025
Table 4: SHD across different dimensionallitis (d {10, 20, 50}) using data generated from ANMs.
ERI
ER2
Methods
d=10
d=20
d=50
d=10
d=20
d=50
SkewScore
2.90 1.47
7.08 +: 2.32
19.32
7.26
8.26 2.42
30.50 = 5.52
76.50
=::00:::0 5.50
HOST
7.20 +: 2.04
17.30
3.95
44.70
6.56
17.80
2.44
37.10
4.04
95.70
+: 6.42
DiffAN
3.50
1.40
9.90 + 1.87
23.56 1.27
7.31 + 1.42
30.10 + 4.06
82.29
5.95
DAGMAMIPP
6.90 1.04
15.60
1.28
37.60 4.18
12.30 0.64
37.50 1.20
92.40 1.62
datasets to improve estimation aceuracy. 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 ER graph model, and vary the sample size across {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
1.90 1.70
0.90 + 0.75
0.60 = 0.32
2.10 1: 0.83
1.30 1 0.64
1.10 1 0.58
HOST
4.40 = 1.91
4.40
1.11
3.50 = 1.47
4.56 + 1.34
4.20
+ 1.00
3.90
+ 1.45
DiffAN
4.70 1.68
2.30
+ 1.19
2.40 - 0.92
3.60
+ 1.50
2.80
+ 1.47
2.75
+ 1.39
VarSort (Reisach et
2021
2.10 1.45
2.00
1.38
1.50
1.81
2.20 + 1.54
2.10 + 1.51
1.90 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. (2023). According to the reported results from Immer et al. (003), our method performance is slightly better than LOCII GRCI [Strbb & Lasko (203) performs the best, and QCCD Tagasovska et al. (2020) also performs quite well. These methods are likely more robust to assumption violation in real data such as low sample complexity. The Tübingen cause-effect pairs dataset contains discrete datasets, where the score function is not well-defined, and datasets that do not satisfy the symmetric heteroscedastic assumption, which is challenging for our method.
Multivariate extension of the case study in Section Finally, we conduct experiments to explore the multivariate extension of the case study in Section | The data generation process for experiments for multivariate HSNMs with latent confounders is based on Montagna et al. (2023a). Let E R represent the latent common cause. For each pair of distinct nodes X; and X,, a Bernoulli random variable Ci3 ~ Bernoulli(p) is sampled, where Cij = implies confounding effect Z between X; and Xj. The parameter o determines the sparsity of confounded pairs in the graph. our experiments, we choose p 0.2. Data generation follows HSNM with latent confounders: de (pa)((.)) + (pp (X.)) NN for each variable = -,...., d. Here,
X; = ff (x.) pax (X,) and paz (X:) represent the observed and hidden variable parents of X,, respectively. The model incorporates the Gaussian processes with an RBF kernel of unit bandwidth to generate the function fi, a linear function d,, and a conditional standard deviation Ø; modeled as 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 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 NoTaa have the best performances. In sparser graph (ER1), SkewScore performs the best when d E (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.
The data used for the barplots Figures Ba and 3b are provided in Table SkewScore maintains at least 95% accuracy across all eight settings in Table g while other baselines fall below 85% accuracy at least one setting.
27

---

Published as a conference paper at ICLR 2025
Table 6: Topological order divergence: Synthetic data generated using multivariate HSNM latent con- founders with an ERI graph.
Method
d=10
d=20
d=30
d=40
d=50
SkewScore
1.20 + 0.98
2.40 2.33
4.60
1.80
6.60
3.29
7.20 +: 2.64
HOST
4.00 1.48
8.10
10.40
14.90
2.66
77.00
DiffAN
2.80 1.17
2.60
8.60
2.91
9.50 2.46
1 3.67
DAGMAMLP
6.00 2.00
14.50 2.11
27.10 3.33
34.50
4.13
NoTears
1.50 + 1.02
5.30
5.00 2.57
6.80 3.97
6.30 2.45
Table 7: Topological order divergence: Synthetic data generated using multivariate HSNM with latent con- founders with ER2 graph.
Method
d=10
d=20
d=30
d=40
d=50
SkewScore
3.00 + 1.61
8.60 + 2.11
14.00 + 4.73
16.60 + 2.42
25.80 + 4.26
HOST
1..80
22.80 4.44
32.90 5.47
41.90 + 6.55
51.80
6.27
DiffAN
3.30 2.10
7.20 = 2.40
13.20
1.10
22.10
DAGMAMLP
17.60 1.80
34.90 2.47
54.20 2.18
74.40
2.15
91.20 + 3.19
NoTears
8.10
16.70 6.81
2.79
25.30 - 7.06
27.00 4.45
Table Accuracy (%) causal direction estimation across different data generation processes bivariate and lotent-confounded triangular heteroscedastic noise models. These data were used for the barplots in Figures and 3b Additional notes: Dark green bold indicates the hest performance, and dark green regular indicates the second best performance.
Bivariate (no latent confounders)
Latent-confounded Triangular
GP-Sig (Gauss)
GP-Sig
Sig-abs (Gauss)
Sig-abs (t)
GP-Sig (Gauss)
GP-Sig
Sig-abs (Gauss)
Sig-abs (t)
(t)
SkewScore
99
100
98
95
97
96
100
LOCI
99
99
98
96
68
83
99
100
HECI
100
95
100
64
54
64
99
83
HOST
100
81
100
89
89
91
100
91
DiffAN
92
89
92
48
76
90
80
35
DAGMA
72
58
58
26
45
61
40
NoTears
85
90
100
48
55
55
13
52
28


Reviews:
[{"review_id": "fq0kj2POjX", "reviewer": "Reviewer_h4sn", "summary": "This paper provides a novel causal discovery algorithm that leverages skewness to determine causal direction when heteroscedastic noise (HN) is present.", "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. \n\n2. 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? \n\n3. 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. \n\n4. Can the proposed algorithm provide confidence intervals or uncertainty quantification in measuring the skewness or determining causal directionality?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "1. An example in Section 5 provides a clear, practical illustration of the algorithm, connecting theoretical concepts with application. \n\n2. 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. \n\n3. 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. \n\n4. 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.", "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. \n\n2. 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. \n\n3. 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Ii1OnvfUOh", "reviewer": "Reviewer_b2UL", "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.", "questions": "Is it possible to provide identification results such as 'the skew score is larger under the wrong causal direction' without symmetric noise assumption?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "The notation is well-organized and the paper is easy to follow overall.", "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. \n\n2. 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. \n\n3. 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "0DU81KXvQh", "reviewer": "Reviewer_8P9r", "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.", "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.\n2. 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?", "limitations": "", "rating": 5, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. The introduction of the SkewScore, which is leveraged to identify causal structures without requiring the extraction of exogenous noise.\n2. A theoretical extension of the criterion to multivariate cases, along with empirical validations showing its robustness, particularly in scenarios involving latent confounders.", "weaknesses": "See questions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "jDyhcGiqMT", "reviewer": "Reviewer_8vJA", "summary": "==Post-rebuttal update==\n\nThe authors' replies to my questions show their further insights into the problem and how the work could be deepened.\n\nI also had a rough scan of other reviews and their respective replies, and I hold my original score.\n\nThis 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\".\n\n==End update==\n\nThis 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.", "questions": "### On Proposition 5\nCould 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)$.\n\nProposition 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)?\n\nIt’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?\n\n### Additional Comments\nCould similar theoretical results be proven for standard skewness?\n\nAssumption 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?\n\nMinor: There is a duplicate bibliography entry for Mityagin (2015).", "limitations": "", "rating": 10, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 4, "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.\n\nThe extensions to multivariate cases and the potential applicability under hidden confounding are valuable contributions.\n\nThe theoretical analysis appears rigorous (though I haven’t verified every detail).\n\nThe 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.", "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.\n\nThe method's applicability under hidden confounding might be limited to special confounding structures.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "g7FZntaOhT", "reviewer": "Reviewer_iLTQ", "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.", "questions": "1. Which ER graph model is used in the experiments, $G(d, p)$ or $G(d, m)$? \n2. Can you explain in more detail why $c\\neq0$ in general in Remark 6?\n3. What is the runtime for the order learning and independent tests separately?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "1. The method is novel and can be potentially applied to causal discovery with latent confounders.\n2. The theoretical justification of Assumption 1 seems sound, and the experiments showed that the proposed method works particularly well on sparse HSNMs.\n3. 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.", "weaknesses": "1. The author mainly used ER1 to evaluate the methods on HSNMs without latent confounders, which is relatively sparse.\n2. 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?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)