Causal Discovery via Quantile Partial Effect
Authors: Yikang Chen
Xingzhe Sun
Dehui du
Venue: ICLR 2026 Poster (2026)
Keywords: causality
causal discovery
causal order
identifiability
normalizing flow
Primary area: causal reasoning
TL;DR: We propose a novel parametric assumption that affords cause-effect identifiability solely from the observational distribution, concomitantly generalizing and relaxing the Functional Causal Model assumption.

Abstract: Quantile Partial Effect (QPE) is a statistic associated with conditional quantile regression, measuring the effect of covariates at different levels. Our theory demonstrates that when the QPE of cause on effect is assumed to lie in a finite linear span, cause and effect are identifiable from their observational distribution. This generalizes previous identifiability results based on Functional Causal Models (FCMs) with additive, heteroscedastic noise, etc. Meanwhile, since QPE resides entirely at the observational level, this parametric assumption does not require considering mechanisms, noise, or even the Markov assumption, but rather directly utilizes the asymmetry of shape characteristics in the observational distribution. By performing basis function tests on the estimated QPE, causal directions can be distinguished, which is empirically shown to be effective in experiments on a large number of bivariate causal discovery datasets. For multivariate causal discovery, leveraging the close connection between QPE and score functions, we find that Fisher Information is sufficient as a statistical measure to determine causal order when assumptions are made about the second moment of QPE. We validate the feasibility of using Fisher Information to identify causal order on multiple synthetic and real-world multivariate causal discovery datasets.

Published as a conference paper at ICLR 2026
CAUSAL DISCOVERY VIA QUANTILE PARTIAL EFFECT
Yikang Chen, Xingzhe Sun, Dehui Du
Shanghai Key Laboratory of Trustworthy Computing, East China Normal University
ABSTRACT
Quantile Partial Effect (QPE) statistic associated with conditional quantile re- gression, measuring the effect of covariates at different levels. Our theory demon- strates that when the QPE of cause on effect is assumed to lie in finite linear span, cause and effect are identifiable from their observational distribution. This gener- alizes previous identifiability results based on Functional Causal Models (FCMs) with additive, heteroscedastic noise, etc. Meanwhile, since QPE resides entirely at the observational level, this parametric assumption does not require considering mechanisms, noise, the asymmetry of shape characteristics
the Markov assumption, but rather directly utilizes
the observational distribution. By per- forming basis function tests on the estimated QPE, causal directions can be dis- tinguished, which is empirically shown to be effective in experiments on a large number of bivariate causal discovery datasets. For multivariate causal discovery, leveraging the close connection between QPE and score functions, we find that Fisher Information is sufficient
statistical measure to determine causal order when assumptions are made about the second moment of QPE. We validate the feasibility of using Fisher Information to identify causal order on multiple syn- thetic and real-world multivariate causal discovery datasets.
INTRODUCTION
Multivariate causal discovery aims to elucidate inter-variable structures, yielding causal graphs or causal orders crucial for non-parametric identification and effect estimation in causal inference (Pearl) 2009). The sparsity of these graph structures aids downstream tasks such as feature se- lection (Guyon al., (007) and disentangled representation learning (Yang al. 2021). To identify hidden Directed Acyclic Graphs (DAGs), various causal discovery methods have been developed, including constraint-based (Spirtes & Glymour [991] Spirtes 1995) and score-based (Cooper & Herskovits, 1992 Chickering, 2002) approaches. However, these methods typically identify only up to an equivalence class, failing distinguish causes from effects without additional assumptions.
Functional Causal Models (FCMs) are a class of Structural Causal Models (SCMs) that impose con- straints on causal mechanisms. These include Additive Noise Models (ANM) (Hoyer et al. 2008), Heteroscedastic Noise Models (HNM) (Tagasovska et al. 2020), and Post-Nonlinear (PNL) models (Zhang & Hyvärinen 2009). The inherent asymmetries in FCMs ensure that models satisfying these assumptions can distinguish cause from effect, except in some marginal cases. However, FCMs re- quire strong assumptions about the underlying mechanism, noise, and Markov property, which may not hold in real-world scenarios. Recently, causal velocity (Xi et al., 2025) has generalized FCMs without requiring assumptions on functional form and noise. Nevertheless. reliees on counterfactual concepts, making it challenging test the validity of the underlying counterfactual assumptions.
Inspired by causal velocity, this work discovers a more fundamental concept called Quantile Par- tial Effect (QPE), which reflects the shape characteristics of the observational distribution. This is a statistic purely at the observational level, defined by the conditional quantile function (Koenker) 2005) (Definition 3. It includes two equivalent definitions (Propositions 3.2) and (.3). We also find cause-effect identifiability presents if QPE is assumed to lie in finite linear span when basis functions are given (Assumption 3.5). This assumption generalizes past FCM assumptions (detailed in Table Tp. In contrast to previous methods, however, identifiability via QPE relies purely on the observational distribution, independent of the underlying mechanism, noise, or Markov properties
"Corresponding author; dhdu @ sei.ccnu.edu.cn

---

Published as a conference paper at ICLR 2026
Joint of X,
Conditional of X
Conditinal of X
QF X
23
hb
y F-
X
Figure 1: Distributions and their QPEs for ANM Y (in + U and HNM Y =x= + tanh((X - 1)²)) (a) Joint distribution (heatmap) and samples (scatterplot); (b) Conditional density of Y I |_X (heatmap), conditional quantiles (white curves), and their gradients (white arrows); (c) QPE of | (3D surface) and its intersection with the Y-Z. plane (white curves); (d) Conditional density and conditional quantiles of X Y: (e) QPE of X Y (3D surface) and its intersection with the X-Z plane. ANM guarantees that the intersection
the QPE of Y X in the Y-Z plane is a constant function, while HNM guarantees is affine function, due to restrictions the QPE form (see Table I; the converse generally does not hold Section 3.2).
Two algorithms for bivariate causal discovery based on QPE estimation and ba- sis testing were subsequently developed (Sections 4.2 and 4.3) Experiments on bivariate causal discovery datasets demonstrate the effectiveness of these algorithms (Section 6.1).
Corollary 3.8
When discussing multivariate causal discovery, we turn to indirect statistical measures to determine causal order, estimating QPE becomes increasingly difficult in high dimensions. Given the close relationship between QPE and the score function (Lemma 3.4), we find that the second moment of QPE influences Fisher information (Theorem 5.2 Under certain assumptions regarding only QPE, Fisher information can sufficiently distinguish between cause and effect (Corollary 5.3). Based on this theory, we develop simple and efficient non-pan non-parametric algorithm to identify causal orders, provided that /Assumption 5.4 is satisfied. Finally, We validate the feasibility of this method on multiple synthetic and real-world multivariate causal discovery datasets (Section 6.2).
PRELIMINARIES
Notation and General Assumptions Generally, this paper discusses random variables RRd, We denote one-dimensional random variable and its realization as X and respectively. For multi- dimensional random variables, we use X = (X1..... X²) and æ = (e1....,aa) In some cases, lowercase boldfacing also represent vector-valued functions, while uppercase boldfacing may denote matrices. Set operations hetween random variables refer to operations on their index sets and will also be used. For all probability density functions px appearing in this paper, we assume they are always existent, strictly positive, and at least C -f-ncfuttoos when involving k-th derivatives.
Structural Causal Models For set of variables X, SCM is triplet (f.X,U) U), where X and are endogenous and exogenous variables, respectively. For each X; E X, we have X; = (P,, U,). where P. X (Xx, are the parent variables and f, is the causal mechanism. An SCM is recursive if the graph G with nodes X and directed edges (X,, X,) for every X, E P, is DAG. A causal order any total order consistent with the topological sort of G. A recursive SCM satisfies the Markov assumption if the exogenous variables are jointly independent (i.e., causal sufficiency or no latent confounders (Pearl 2009)). which implies P. U, for each X₁.
Score Function and Stein's Identity This paper refers to the gradient of the logarithmic density, V z log px, the score function. When mild boundary condition, lim h(x:) px (æ:) 0.
2

---

Published as a conference paper at ICLR 2026
General Assumptions
ASSUIIIOO 3
Assunntion 3.7
Prymsstren Fod
Rentee 3.
Conntar
Menoomicity Mankowanity
Asumattin
Progosition 3.
Thorrem S.2
Conolary
Figure 2: For clarity, we present a dependency graph of the assumptions and theorems in the main text, where "A-+B" indicates that depends on A: (a) Theorems in cooresppond /Section 3 (b)
Theorems in correspond to Section 4 (c) Theorems in correspond to [Section 5
holds for class of functions h, we have Stein's Identity (Stein 1972p:
E [, log + v h(x]] =.
direct consequence is that the score function has zero expectation. Its variance, which is equiva- lent its second moment, constitutes the Fisher Information. For single variable X,, the Fisher Information is the variance of the partial score 0z, Pxg... equals the negative expected second derivative: E[(O,, log px) ) = -E [o? log px] where the second term is the negative expectation of diagonal entry of the Hessian matrix of the log-density.
Px. Under the same regularity conditions,
Causal Normalizing Flows Javaloy et al. (2023) summarized the connection between autoregres- sive normalizing flows and SCMs, proposing to use flows model SCMs and fit observational distributions, and showing their capacity for counterfactual inference. Taking single causal mech- anism X, (P, U,) as an example, assuming it is strictly monotonic w.r.t. U, and the SCM satisfies the Markov assumption, then according to the change-of-varian formula:
logPx, log
lug PU. log Og,uul,
where the causal flow is modeled as (p,, r:). which maps the endogenous variable X, to latent variable Ue. and PU, is called the latent distribution. This formula also allows us use Maximum Likelihood Estimation (MLE) to estimate parameters 0 by maximizing Eloog Pu, + log |0,,4|] Mol] Javaloy et al. (2023) also discussed how to use autoregressive flows to fit multivariate SCMs and proved that under the above assumptions, such modeled causal flows ultimately lead to latent vari- ables identifiable up to element-wise invertible transformations. Chen & Du (2025) then proved that causal flows and the true SCM are counterfactually consistent under these assumptions.
3 QPE FOR CAUSAL DISCOVERY
3.1 DEFINITION OF QPE AND EQUIVALENT CONCEPTS
We use QPE, a statistical object induced from conditional quantile regression, for causal discovery. QPE describes the sensitivity of quantiles to covariates, quantifying covariate effects at different levels. Its visualization, as shown in Figure I] reflects the shape characteristies of the observational distribution. Let Frix be the conditional cumulative distribution function (CDF) corresponding to the conditional distribution pr|x, and Qylx be the conditional quantile function. Then:
Definition 3.1 (Quantile Partial Effect). The Quantile Partial Effect (QPE) of random variable Y given X is x (3 x) = V Q (T xx), where the quantile 7 = Qyi Frix.
Based on the equality Fyix = Qr|x we can immediately find an equivalent description of QPE: Proposition 3.2 (QPE from CDF). d V FxFFixx0
-VaFyx//PyIx
Xt al. (2025) introduced the concept of causal velocity to generalize FCM-based bivariate causal discovery. Consider a Markovian causal mechanism Y s(X, U), where / is strictly monotonic
3

---

Published as a conference paper at ICLR 2026
Table 1: The functional forms of FCMs and their corresponding QPE forms. For HNM, > and for PNL, g' o'(g The QPEs of these FCMs can all be expressed in finite-rank forms.
FCM
Functional Form
QPE Form
QPE Basis Functions
LiNGAM
Y = ctæ +
ANM
Y= a(x) +
Va (Va Vb)
1
HNM
== =(x) + b(x)
+
Vb
0 y
1,
PNL-ANM
Y = =(((()) +
Vag
9
PNL-HNM
Y g(o(x) + b(x) u)
Vaa 5 V6) g
4b
b g g
9.9 g
Assumption 3.5
Perhaps no closed-form
Σ eg (1)
Φ
w.r.t. U. The term f(x', u(x,, y)) is referred to the counterfactual outcome (or SCM flow). where 21 (f -))= Causal velocity is then defined as fz,(xx u). According to (Nasr- Esfahanny et al, 2023 under these assumptions, the counterfactual outcome is identifiable from the observational distribution, meaning causal velocity is independent of the exogenous distribution. In fact, we find
Proposition 3.3 (Causal Velocity is QPE). 4b1x = vaf(= u) u. 0/ a.e.
According to the Causal Hierarchy Theory (Bareinboim al. 2022), this implies that causal ve- locity, a counterfactual quantity, can "collapse" entirely into an observational quantity, despite its definition through counterfactual concepts. In contrast causal velocity, QPE does not require monotonic mechanisms or the Markov assumption. Its definition (Definition 3.1) and cause-effect identifiability (Section 3.2) depend solely on the observational distribution.
Xi et al. (2025) also found PDE relationship between causal velocity and the score function, based on the SCM flow and the continuity equation. This PDE naturally applies to QPE due to their equivalence. In fact, by directly taking the partial derivative of the implicit function defined in Proposition 3.2, we can derive this equation:
Lemma 3.4. V Py + & a, Pog + 0,5 = 0 and lim Eprix = -ffff = d/y|x-
Now consider each covariate X,, and denote the corresponding component in QPE as dyix 0,,Qy Qy Since 0z, log independent Y, according to log Px,y log px + log PYJx. we can obtain the equality 0y log x = 0g log Px.y between the conditional log density and the joint log-density. Therefore, second-order mixed PDE can be derived from Lemma
from 3.4
oz 0y logpx.y + 0y (vvlxx; 0g logpx,y + 0gVrix.)
which only involves the score function of the joint distribution and QPE.
3.2 CAUSE-EFFECT IDENTIFIABILITY BY QPE IN FINITE LINEAR SPAN
For a function f R + IR, if it can be represented as linear combination of a set of basis functions & (di...., 2x) such that there exist coefficients , satisfying f'xx) Σラ= ej (x)),
then f is said to be in the finite linear space spanned by Φb denoted as f E sp)n(
Assumption 3.5. For each cause variable X; and any x, t/y x,(( known set of basis functions that depend only on the effect variable Y.
(: ax) span(()) where & is
In other words, this assumption requires that each component thy x- x) can be represented in finite-rank form as a sum of finite c. (v)). where a.j are coefficient functions w.r.t.
We observe that previous FCMs with constrained functional forms can in fact be expressed a finite-rank form, detailed in Table 1 However, note that Assummtion 3.5 requires known Φ; thus, although it generalizes ANM and HNM, does not directly apply to PNL where g is unknown. Furthermore, this assumption is entirely independent of monotonicity or the Markov property, dis- tinguishing our approach from the generalization path taken by Xi et al, (2025). Specific examples that violate these properties yet still satisfy Assumption 3.5 are provided in Appendix A.3
4

---

Published as a conference paper at ICLR 2026
Next, consider constructing equations solely dependent on the joint distribution using linear
lationships. For a set of multivariate functions fI..... : fk. where each fa Rd + R, the determinant
f1
f2
fk
Jzh1
0z12
Jzf
Wx(f1.... (a) = det
0 f
0 f2
0 Sk
is called the Wronskian determinant w.r.t. the variable X. y 02, 0y log Pogx and |||x 0,( 02 logg.... Y Φ¹). Here, sx,,Y Y is an off-diagonal element in the Hessian matrix, while fy is the partial derivative of ø's Stein operator w.r.t. y. Then, according to the second-order mixed PDE (Equation 2 the following theorem holds:
Theorem 3.6 (Identifiability of QPE in Finite Linear Span). For each variable X; and any æ, x x)) E span(() implies (sx,,, ( = 0 for any y. If: (i) The components in xT(x))
Tyix are linearly independent; (ii) There exists y such that Wy(ny + 0; (iil) For each basis function dz. limy-tto 93 Py|x 0; then the converse is also true.
Theorem 3.6 a necessary (or sufficient, assuming well-behaved conditions) condition for the validity of Assumption 3.5 Crucially, Theorem 3.6 relies solely on the Wronskian determinant, which is function only of the joint distribution and the basis set Φ (as well as their higher- derivatives). By leveraging this Wronskian determinant, we can formally characterize the asymme- try of shape characteristics in the observational distribution given Φ. and further achieve cause-effect identifiability that relates only to the observational distribution:
Assumption 3.7. The set of basis functions Φ is known, and for any Z X and any X; X we have W 2 (sx., m2 # O.
Corollary 3.8 (Cause-Effect Identifiability by QPE in Finite Linear Span). If Assumption 3.5 and Assumption 3.7 hold simultaneously, then Y is the effect variable.
Proof Theorem 3.6 establishes a necessary condition for Assumption 3.5
(dyixx; E span ())
to hold is Wz(sx, M = 0, so Assumption 3.7 implies 12 span (). Consequently, Assumption 3.5 does not hold for any Z YY, Y as the only possible effeet variable.
It is worth noting that both Assumption 3.5 and Assumption 3.7, utilized in Corolllly 3.8 depend solely on the observational distribution (assuming Φ is pre-specified). The former is independent of monotonic mechanisms and the Markov assumption, as illustrated in [Appendix A.3] the latter describes a PDE system over the observational distribution, reflecting its asymmetry. Consequently, via QPE, we achieve cause-effect identifiability based strictly on the observational distribution.
QPE FOR BIVARIATE CAUSAL DISCOVERY
4.1 QUANTIFY CAUSE-EFFECT IDENTIFIARILIT THROUGH THE LENS OF FCM
While Corollary 3.8 establishes cause-effect identifiability through Assumption 3.5 and Assump
tion 3.7 weaknesses remain. Most notably, characterizing the probability that Assumpisnmiiinn 3.7 holds chal challenging. To resolve this, we revert to the bivariate FCM framework, replacing Assumption 3.7 with the assumptions of monotonic mechanisms and the Markov property. This enables quantita- tive description of cause-effect unidentifiability grounded in manifold theory.
Assumption 4.1. A bivariate SCM (f, X,,XUU U)) such that Y = [XX,U, U), where the causal mechanism and exogenous distribution py are given. Furthermore, f is strictly monotonic w.r.t. U. and the Markov property or X 1 U holds.
Assumption 4.1 describes a family of FCMs where only Px of X is unspecified. All possible px distributions form an infinite-dimensionl manifold, O denoted as dim(@) = oc, representing the scale of this FCM family. We then introduce forward and backward versions of Asssuppiion 3.5 on QPE: "for any x, d (- [z) E span(p) (en"" "for any y, (- [1) E span( These conditions induce two submanifolds, O x and θ x db. respectively. These submanifolds contain all Px that satisfy the corresponding assumption and must least satisfy Theorem 3.6
5

---

Published as a conference paper at ICLR 2026
Corollary 4.2. Assume that there are k basis functions in Φb and certain regularity conditions hold. Then: (i) dim x = =o; (ii) dim (om(@xxxx) xx-Y; 2
In other words, under the premise of an FCM satisfying Assumption 4.1 if the forward version of Assumption als holds, the choice for px remains arbitrary. However, if the backward version of Assumption holds, Px is constrained to lie within submanifold of at most & + dimensions. It is obvious then that the unidentifiable scenario, where both the forward and backward versions the assumption hold simultancously, also falls within such at most k + 2-dimensional submanifold. Since this manifold is miiee-immnsiooal, is typically of measure zero for non-degenerate mea- sures anfnite-dimensionall space. Therefore, generally, except for extremely marginal cases, Assumption 3.5 universally provides identifiability of causal direction.
While Corollary 4.2 retains the need for monotonic mechanisms and the Markov property, it offers guarantees that serves as a generalization of quantifiable cause-effect identifiability under the classic FCM framework. In particular, this corollary generalizes the conclusion by Hoyer et al. (2008) that unidentifiable ANMs lie in a three- ftrce-dimensional affine space (since ANMs correspond to k = where the Wronskian simplifies linear ODE) any finite-linearly spanned QPE (k: oo).
4.2 KERNEL-BASED QPE WITH LEAST SQUARE BASIS TEST
With the identifiability guarantees established above, we can determine the causal direction by es- timating the QPE and checking whether lies within the assumed span. The original definition of QPE suggests calculation via conditional quantile regression, but its accuracy and efficiency depend on the underlying quantile regression. Below an efficient and entirely non-paramen non-parametric method for bivariate causal discovery, independent of quantile estimation.
By Proposition 3.2 we can use non- non-paramente con-parametric methods to estimate V FxFyyx and 0, Fx. thereby indirectly computing dbyix x- Given N samples (zr,, drawn from px y, the conditional CDF can be estimated using kernel estimator as:
Kernel-based QPE
EN K (xx, x) S(y.. y)
Frix x))
EN K(x;,
where is smoothed version of the indicator function (wx E d(x)) and S smoothed version of the 1 13y S y). We choose K as Gaussian kernel and S as a sigmoid function. Based on this equation, since the kernel functions are known and smooth, we can derive closed-form expressions for V, Fy x (yl a)) and F,Fr X m)
Least Square Basis Test After obtaining the estimated QPE vhyix. we can test whether the es- timate satisfies Assumption 3.5 at fixed (æ, y) pairs. Specifically, we can pre-select {vi..... yM } locations for and construct basis matrix B such that Bm.j - (y Concurrently, select (æ1..... ) as test samples for æ. For the i-th component, based on these fixed locations, we construct response matrix L such that V.,, = / : (N/m xx)) for each ym and 2₁. the true QPE x (y x) E span(Ø). then for a coefficient matrix C, such that Cr,,, (.)() it must hold that V = CTBT, Since C, is unknown, this problem can be modeled as an Ordinary Least Squares (OLS) problem where the noisy L; is treated as linear response w.r.t. B. Then,
V Č,B = Ć. : B(BBBB)
arg min C:
where, because the S are not necessarily linearly independent across the M test locations, which means that the inverse of BT may not exist, we use the pseudoinverse BB+ instead to ensure numerical stability. Finally, we average the residuals for each component:
1
Σ
L 10---BB VBB(BBB)+Br
This value reflects the degree to which Assumption 3.5 satisfied: Ex-Y 0 only dyyx
x) E sp)n(Φ) holds for all covariates X,. In bivariate causal discovery, identify the causal direction by comparing Ex and Ey- We infer Y as the effect if Ex > E and X as the effect otherwise.

---

Published as a conference paper at ICLR 2026
True QPE of Y
OPE-K
Causal Velocity
QPET
(0
Figure 3: True and estimated QPEs of Y X at samples from HNM Y X + (1 + tanh((XX
))²)) U. From left to right: (a) True QPE: (b) QPE-k (Section 4.2): (c) Causal velocity model (Xi
al. 2025) (V-NN); (d) QPE-f (Section 4.3). The black lines represent the intersection of the QPE
surface with the Y-Z plane. Only QPE-f's trend tends to match the true QPE in high-density areas.
4.3 FLOWDBASED QPE WITH NEURAL BASIS TEST
While kernel methods combined with the least squares basis test offer fast and effective non- parametric approach, their performance is limited by the choice of kernel bandwidth and sample size. QPE estimation tends to be less accurate and prone to overly smooth predictions. We will now introduce neural network-based parametric method, which, despite requiring neural network training, often yields more accurate QPE estimates, as shown in Figure 3) In experiments, this method outperforms most ofate-of-the-arr bivariate causal discovery methods.
Flow-based QPE Propossiinio 3.3 for any SCM that satisfies the assumption, we can calculate the QPE using the function that maps the observational distribution to an exogenous distribution. According to [Section 2| any causal flow modeled by Equation I is an SCM satisfying Proposi- fion 3.3, and these SCMs are counterfactually identifiable. Thus, we can estimate the QPE through any flow parameterized as in Equation and optimized via MLE. We use the standard normal distribution as the latent distribution and parameterize it with a neural network. After obtaining the ug using automatic differentiation techniques, thereby x V uo (0y2. Note that Proposition 3.3 implicitly guarantees that the QPE obtained this way is always consistent, even if different latent distribution is chosen.
causal flow ug. can compute ux and ug
calculating the QPE y
Neural Basis Test Although the previously described OLS- OLS-based test still applies to the estimated d here, requires fixed y values and cannot directly test based on irregularly distributed sam- ples. Furthermore, OLS corresponds to MLE under the assumption that errors follow Gaussian distribution. To address these issues, we consider using a parameterized neural network to directly model the coefficient function o (3)), which corresponds to optimizing the objective:
arg minex-yX + A||l|| where EX-Y = E d - Co ΦT 0
where À is regularization hyperparameter, and C. the matrix formed by C.j.a- In contrast to the OLS test, here the samples 2t, 2/ can be irregularly distributed. The optimal E will serve as a measure how well Assumption 3.5 satisfied, and its use in determining the causal direction is the same as described previously.
This testing method can be simplified, as proposed by (Xi et al. 2025). to directly fit s/yix using unconstrained neural networks, called V-NN models. Specihcally, they consider directly optimizing PDE in Lemma 3.4, which is guaranteed to converge under the vanishing assumption. However, Lemma 3.4) involves first-order gradients, and the performance heavily relies on the accuracy of score function estimation algorithm.
5
FISHER INFORMATION FROM THE QPE PERSPECTIVE FOR MULTIVARIATE CAUSAL DISCOVERY
5.1 THE CONNECTION BETWEEN FISHER INFORMATION AND QPE
Due to the curse of dimensionality, QPE estimation in high-dimensional settings performs poorly, limiting its direct application in multivariate causal discovery. Since Lemma 3.4 describes de-

---

Published as a conference paper at ICLR 2026
terministic relationship between the score function and QPE, we consider constructing an equation between the Fisher information and QPE. We then indirectly infer the Fisher information by assum- ing QPE, enabling its for causal discovery. For variables X, Y, denote the score functions of the joint and marginal distributions, 0z, log Px, and 0z log px, as 8x, and rx; respectively. Assumption 5.1. Both vy x Px, and vhy x (,,0 xx, px,Y vanish at to0.
Theorem 5.2. For each covariate x, Assumption 5.7 holds then
E ((Vr)x,)) (sy) E (sx,)) E [(xx,) )² + + (0, / 83
Corollary 5.3. For each covariate X;,
Assumption 5.7 holds then
E [1, ) + 2 p 0, vyxxi < [((ev)x,;) 1) (sx) + (rxx))
if and only E [(xx, | > E [9.
This implies that if the marginal Fisher Information is sufficiently large, or the second moment of yyjx, & sufficiently large and the moments of the higher-order partial derivatives of r w.r.t. y are sufficiently small, we can directly distinguish the effect variable by the Fisher Information.
Qualitative analysis under the heteroscedastic Gaussian assumption By introducing additional assumptions on the second-order partial derivative (y (i.c., HNM), we can isolate its second moment from the inequality in Corollary 5.3
E [ E [sy )2]
E ((,,
+ Var [(brxx))) Var [(sy) }²)] E [(x)]
(+'yix,)))
Then, for Corollary 5.3 to be as valid as possible, we need to make assumptions about QPE: (( [J must be sufficiently large, and one of the following must hold: (i) Varl(wyixx " sufficiently small when vyyxx (e.g., ANM); (ii) Var((eyix,))
sufficiently small and (I,,,)rix))) is sufficiently small when de = 0 (c.g.. HNM). Qualitatively, when assum- ing [[[0,,vyxx, is relatively stable, this trend requires the coefficient of variation (CV) of the squared QPE, [ xx a) [E(//((())) ?]. to be sufficiently small.
Furthermore, if we assume heteroscedastic Gaussian conditional distribution, i.e., Y X ~ ((((z), (a(z)) we can derive the exact value of the CV of the QPE under linear μt, σ. or its upper bound in the general case (see Appendix B.2 for the complete derivation):
Var [(by]x))
/4 1442 2x² 1 x²
Var [(
/1 + 6-2 33+
v E (
[n],
E [( <
E [(])
1 12
E [(x')]]
where k = 0//p²' and < |x| Suu
Based on the above analysis, qualitatively, under the heteroscedastic Gaussian assumption, the CV of QPE is sufficiently small when: (i) |k] = |oo//u| is sufficiently small, meaning rather than o dominates the shape the distribution; and (ii) μ and are close linear functions. To further con- fimm the correctness of this qualitative analysis, we conducted experiments on datasets with widely varying average |x] and linearity. The results are presented in Appendix D.4
5.2 FISHER INFORMATION CAUSAL ORDERING
Assumption 5.4. For each variable X,, any of its parents P, E P: satisfy Assumption 5.1 and
E (0,0x, p..)) + +2 (2xxp/1..) j..j 8 (.
E [Vx. 1)² (sy) )2| ++++ (rx,))
8

---

Published as a conference paper at ICLR 2026
According to Corollary 5.3 if both Assump-
Algorithm Fisher Information Causal Ordering
tion 5.1 and Assumption 5.4 hold for a given set of variables X, then the variable X; with the minimal Fisher information [x(xx,)) must have no child variables, called the leaf variable. Therefore, peereforr.Assumpttion fication, implying that once the score function is estimated, there is a possibility of finding the leaf node. After selecting leaf node X, a sub- problem on X is formed, which allows for recursion until all variables are removed. The resulting sequence of variables is the reverse of the causal order. This process is called FICO (Fisher Information Causal Ordering), and the detail is described in Algorthhh I
Input: a set of random variables X. Output: causal order n
X <- X, nT -
for j = to d do
[Assumption 5.4 is crucial simpli-
+ arg min; E (0,, log
x x 1, + [1
end for return π
Next, given the causal order, for d variables, we only need to perform most d(d-1)/2 conditional independence tests to prune the full DAG induced by this topological order, thereby retaining only the necessary edges. Under the assumptions of causal sufficiency and faithfulness, independence tests ensure that the DAG lies within the Markov equivalence class, and the direction of each edge is determined by CCrolllry 5.3 Thus, the resulting DAG is fully identifiable.
Notably, CaPS (Xu et a 2024), another score function based causal discovery method, is algo- rithmically fully equivalent to FICO. However, distinct differences exist between them: (i) CaPS is specifically tailored for the ANM, and its theoretical derivation relies on the ANM assumption. In contrast, FICO applies to any model, provided that the vanishing assumption (Assumption 5.1) holds and the QPE satisfies Assumption 5.4 (ii) CaPS utilizes -E 8 logPx| whereas FICO employs E [(.,, log )²]. Although both are equivalent, FICO involves lower-order derivatives, resulting in higher computational efficiency (see Section 6.2 for details).
It is important to emphasize that although FICO's assumptions and theory have been extended to general cases, and despite its simplicity and computational efficiency, it still faces challenges worth addressing: (i) The assumption relied upon by FICO (Assumption 5.4) has only been analyzed qualitatively under the heteroscedastic Gaussian setting. We still lack a precise understanding of what the validity of [ssuupption 5.4 implies in broader contexts. (ii) Calculating the QPE in high- dimensional settings difficult, which renders the testability of Ass/sssutioon 5.4fproblematic. Given these limitations, caution still advised when applying FICO in practice.
6 EXPERIMENTS
6.1 BIVARIATE CAUSAL DISCOVERY EXPERIMENTS
Datasets We conducted experiments on 24 synthetic and real-world benchmarks, including: (i) AN, AN-c, LS. LS-c, MNU from (Tagasovska et al.) 2020); (ii) SIM, SIM-c, SIM-g, SIM-In from (Mooij et al. 2016); (iii) Cha, Multi, Net from (Guyon et al., 2019); (iv) Per, Sig, Vex generated by random SCM flows from (Xi et al. 2025); (v) Qd-V, Sig-V, Rbf-V, NN-V generated by constrained QPEs (ppenddxx D.2): (vi) Tübingen cause-effect pairs challenge (Mooij et al.) 2016);; (vii) Gene network reverse engineering challenge from (Marbach al. 2009). Datasets (vi) and (vii) are real- world datasets, while the underlying SCMs of (iv) and (v) are not necessarily ANM, HNM or PNL.
Baselines We consider the following open source methods as baselines: (i) ANM and PNL with independence test, implemented by (Zheng et al., 2024); (ii) ANM-based RECI (Bloebaum et al.) 2018) and CDS (Fonollosa 2019). with implementations from (Kalainathan et al. 2020);; (iii) HNM- based CDCI (Duong & Nguyen 2022), НЕCI (Xu et al. 2022), and LOCI (Immer et al. 30 (iv) The causal velocity model CVEL (Xi et al.] 2025). In Appendix D.2, we will cover additional methods (21 baselines in total based on different theories).
Results Table 2 shows the accuracy of QPE-k (Section 4.2) and QPE-f Section 4.3 on 12 datasets,
with more detailed and complete results available in Appendix D.2 For QPE-f, we present the results of causal flows under their best configurations, as detailed in Appendix D.I As shown in Table 2. QPE-f performs best these benchmarks due to its stronger expressive power and more

---

Published as a conference paper at ICLR 2026
Table 2: Accuracy of QPE-k, QPE-I, and baselines on 12 bivariate datasets. The best is bolded.
Method
AN
LS
SIM
SIM-c
Cha
Net
Per
Sig
Qd-V
NN-V
Tue
D4-s1
Time (s)
ANM
0.43
0.46
0.45
0.49
0.41
0.47
0.49
0.44
0.49
0.48
0.65
0.50
0.250
PNI
0.30
0.33
0.46
0.54
0.45
0.51
0.42
0.43
0.46
0.41
0.51
0.33
37.770
RECI
0.18
0.22
0.44
0.53
0.56
0.60
0.00
0.63
0.49
0.64
0.58
0.002
CDS
0.99
0.76
0.71
0.76
0.71
0.78
0.18
0.08
0.78
0.52
0.67
0.58
0.017
CDCI
1.00
I.00
0.84
0.76
0.67
0.81
0.48
0.42
0.74
0.72
0.68
0.67
0.001
HECI
0.98
0.92
0.49
0.55
0.57
0.72
0.01
0.13
0.59
0.45
0.61
0.42
0.026
LOCI
1.00
1.00
0.78
0.81
0.73
0.87
0.96
0.70
0.71
0.78
0.61
0.58
14.981
CVEL
1.00
0.98
0.63
0.72
0.68
0.62
1.00
0.84
0.91
0.87
0.64
0.67
1.597
QPE-K
0.99
1.00
0.83
0.79
0.60
0.89
0.77
0.89
0.42
0.53
0.54
0.58
0.009
QPE-f
1.00
1.00
0.88
0.88
0.85
0.86
1.00
0.90
0.91
0.90
0.70
0.79
7.804
accurate fitting. but it is relatively time-consuming due the need to train causal flows. QPE-k performs similarly to QPE-f on various benchmarks and runs very fast, but its identification capacity is limited. Additionally, methods based on ANM HNM do not perform well as CVEL and QPE- f on causal flow and constrained QPE datasets. This suggests that our identifiability theories based on QPE have a broader scope of applicability beyond just common FCMs.
6.2 MULTIVARIATE CAUSAL ORDERING EXPERIMENTS
Table 3: Runtime efficiency of FICO and CaPS, in seconds per sub-test.
Method
d
d= 10
20
d= 50
100
CaPS
0.455 + 0.037 0.425 + 0.322
1.074 + 0.056
2.761
0.285
10.822 = 1.037
33.794 = 3.501
FICO
0.797 + 0.364
1.727 = 0.523
5.550 + 0.943
13.538 = 1.248
Results Given that CaPS and FICO are algorithmically equivalent, their performance is nearly identical, with minor differences attributable only to numerical computation errors. For simplicity, therefore report only their computational efficiency in the main text.
Performance comparisons between CaPS, FICO, and other causal ordering baselines are detailed in the appendices: Appendix C.2 describes the datasets and baselines used; Appendix D.3 reports the relationship between performance and sample size: Aizeeddppeiiix D.5]dttails perlormance on real-world datasets; and Appendix D.6 covers performance all score function based causal ordering algorithms perform very similarly across various settings, demonstrating the robustness reported by Montagna et al. (023). This empirically suggests that the underlying assumptions or characteristics of these methods are implicitly satisfied.
synthetic datasets. In summary, find that
Table 3 reports the average time consumed by both methods across all tests under different dimen- sionalities, showing that FICO significantly outperforms CaPS. For a comparison of computational efficiency between FICO and all other causal ordering methods, see Appendix D.7
CONCLUSION
In this work, building upon research into quantile partial effects, propose novel parametric assumption (Assumption 3.5) that enables cause-effect identifiability (Theorem 3.6) solely from the observational distribution. This assumption simultaneously generalizes and relaxes the Func- tional Causal Model assumption. Consequently, we develop two algorithms, QPE-k (Section 4.2) and QPE-f (Section 4.3). for effective bivariate causal discovery, and evaluate their performance in numerous experiments (Section 6.1). For multivariate causal discovery, we investigate indirect sta- tistical criteria related to QPE for efficient causal ordering. We propose FICO (Section 5.2). which performs causal ordering efficiently when AAssumptinn 5.4) holds, validated on both synthetic and real datasets (Section 6.2). For causal discovery based on quantile partial effects, future work will first investigate cause-effect identifiability under more general conditions, such as relaxing the fixed basis function assumption. In multivariate causal discovery, high-dimensiona quantile partial effect estimation is challenging; while this paper explores an indirect alternative, its underlying assump- tions lack interpretability and practical verifiability. Future work should develop more intuitive and inherently suitable information measures for multivariate causal discovery.
10

---

Published as a conference paper at ICLR 2026
ACKNOWLEDGMENTS
This work was sponsored by the Open Project of Shanghai Driverless Train Control System Engi- neering Technology ResarrhhenterroosSUTC--22200
REFERENCES
Elias Bareinboim. Juan D. Correa, Duligur lbeling, and Thomas Icard. On Pearl's Hierarchy and the Foundations of Causal Inference, pp. 507-556, Association for Computing Machinery, New York, NY, USA, 1 edition, 2022. ISBN 978 9781450395861.
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-8239 2022.
Patrick Blocbaum, Dominik Janzing, Takashi Washio, Shohei Shimizu, and Bernhard Schoelkopf. Cause-effect inference by comparing regression errors. In Amos Storkey and Fernando Perez- Cruz (eds.), Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics, volume 84 of Proceedings of Machine Learning Research, pp. 900-909. PMLR, 09-11 Apr 2018.
Maxime Bocher. Certain cases in which the vanishing of the wronskian is sufficient condition for linear dependence. Trans. Am. Math. Soc., 2(2):139, April 1901.
Peter Bühlmann, Jonas Peters, and Jan Emest. Cam: Causal additive models, high- -igh-dimensional order search and penalized regression. The Annals of Statistics, 42(6)) 2014. ISSN 00905364.
Yikang Chen and Dehui Du. Exogenous isomorphism for counterfactual identifiability, In Forty- second International Conference on Machine Learning, 2025.
David Maxwell Chickering. Optimal structure identification with greedy search. Journal of machine learning research, 3(Nov):507-554, 2002.
Gregory F Cooper and Edward Herskovits. A bayesian method for the induction of probabilistic networks from data. Mach. Learn., 344:30993347, October 1992.
Povilas Daniušis, Dominik Janzing, Joris Mooij. Jakob Zscheischler, Bastian Steudel, Kun Zhang, and Bernhard Schölkopf. Inferring deterministic causal relations. In Proceedings of the Twenty- Sixth Conference on Uncertainty in Artificial Intelligence, UAI' pp. 143-150, Arlington, Vir- ginia, USA, 2010. AUAI Press. ISBN 9780974903965.
Laurent Dinh, Jaseha Sohl-Dickstein, and Samy Bengio. Density estimation using real NVP. In International Conference Learning Representations, 2017.
Bao Duong and Thin Nguyen. Bivariate causal discovery via conditional divergence. In Bernhard Schölkopf, Caroline Uhler, and Kun Zhang (eds.), Proceedings of the First Conference on Causal Learning and Reasoning, volume 177 of Proceedings of Machine Learning Research, pp. 236- 252. PMLR, 11-13 Apr 2022.
Bao Duong and Thin Nguyen. Heteroscedastic causal structure learning. In Frontiers in Artifi- cial Intelligence and Applications, Frontiers in artificial intelligence and applications. IOS Press, September 2023,
Conor Durkan, Artur Bekasov, lain Murray, and George Papamakarios. Neural spline flows. In Hanna M. Wallach, Hugo Larochelle, Alina Beygelzimer, Florence d' A'Ahh-BBuu Buc, Emily B. Fox, and Roman Garnett (eds.), Advances in Neural Information Processing Systems 32: Annual Con- ference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada, pp. 7509-7520, 2019.
Josè A. R. Fonollosa. Conditional Distribution Variability Measures for Causality Detection. pp. 339-347. Springer International Publishing, Cham, 2019. ISBN 978-3-030-21810-2 doi: 10. 12.
11

---

Published as a conference paper at ICLR 2026
Olivier Goudet, Diviyan Kalainathan, Philippe Caillou, Isabelle Guyon, David Lopez-Paz, and Michele Sebag. Learning functional causal models with generative neural networks. In Explain- able and interpretable models in computer vision and machine learning, pp. 39-80. Springer, 2018.
Siyuan Guo. Viktor Toth, Bernhard Schölkopf, and Ferenc Huszar. Causal de finetti: On the identi- fication of invariant causal structure in exchangeable data. In A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine (eds.), Advances Neural Information Processing Systems, Curran Associates, Inc., 2023.
volume 36, pp.
Isabelle Guyon, Constantin Aliferis, et al. Causal feature selection. In Computational methods of feature selection, pp. 79-102. Chapman and Hall/CRC, 2007.
Isabelle Guyon, Alexander Statnikov, and Berna Bakir Batu. Cause effect pairs in machine learning. Springer Nature, Cham, Switzerland, October 2019.
Patrik Hoyer, Dominik Janzing, Joris M Mooij, Jonas Peters, and Bernhard Schölkopf. Nonlinear causal discovery with additive noise models. In D. Koller, D. Schuurmans, Y. Bengio, and Bot- tou (eds.), Advances in Neural Information Processing Systems, volume 21. Curran Associates, Inc., 2008
Chin-Was Huang. David Krueger, Alexandre Lacoste, and Aaron Courville. Neural autoregressive flows. In Jennifer Dy and Andreas Krause (eds.), Proceedings of the 35th International Con- ference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pp- 2078-2087. PMLR, 10-15 Jul 2018.
Aapo Hyvärinen, Kun Zhang, Shohei Shimizu, and Patrik O. Hoyer. Estimation of structural vector autoregression model using non-gaussianity. Journal of Machine Learning Research, 11 1599:177- 2010.
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 Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett (eds.), Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pp. 14316-14332. PMLR, 23-29 Jul 2023.
Kasra Jalaldoust, Saber Salehkaleybar, and Negar Kiyavash. Multi-domain causal discovery in bijective causal models. In Fourth Conference Causal Learning and Reasoning, 2025.
Adríán Javaloy, Pablo Sinchez-Martin and Isabel Valera. Causal normalizing flows: from theory to practice. Advances in Neural Information Processing Systems, 36:58833-58864 2023.
Diviyan Kalainathan, Olivier Goudet, and Ritik Dutta. Causal discovery toolbox: Uncovering causal relationships in python. Journal of Machine Learning Research, 21(37):1-5, 2020.
Ilyes Khemakhem, Ricardo Monti, Robert Leech, and Aapo Hyvarinen. Causal autoregressive flows. In Arindam Banerjee and Kenji Fukumizu (eds.), Proceedings of The 24th International Confer- ence Artificial Intelligence and Statistics,, volume 130 of Proceedings of Machine Learning Research, pp. 3520-3528. PMLR, 13-15 Apr 2021.
Roger Koenker. Quantile Regression. Econometric Society Monographs. Cambridge University Press, 2005.
Yingyu Lin, Yuxing Huang, Wenqin Liu, Haoran Deng, Ignavier Ng. Kun Zhang, Mingming Gong, Yian Ma, and Bíwei Huang. A skewness-based criterion for addressing heteroscedastic noise in causal discovery. The Thirleenth International Conference on Learning Representations, 2025.
Daniel Marbach, Thomas Schaffter, Claudio Mattiussi, and Dario Floreano. Generating realistic in silico gene networks for performance assessment of reverse engineering methods. Journal of Computational Biology, (6(2):229-239 2009. doi: 10. 10.10089/mbbbb www.:/cmm...0....mm 2008. 0TTT. PMID: 19183003.
12

---

Published as a conference paper at ICLR 2026
Alexander Marx and Jilles Vreeken. Identifiability of cause and effect using regularized regression. In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, KDD '19, pp. 852-861, New York. NY. USA, 2019a. Association for Computing Machinery. ISBN 9781450362016. doi: 10. 1411/3292500033336533
Alexander Marx and Jilles Vreeken. Telling cause from effect by local and global regression. Knowl. Inf. Syst, (003):1277-1130000 September 2019b.
Francesco Montagna, Nicoletta Noceti, Lorenzo Rosasco, Kun Zhang, and Francesco Locatello. Causal discovery with score matching on additive models with arbitrary noise. In Mihaela van der Schaar, Cheng Zhang, and Dominik Janzing (eds.), Proceedings of the Second Conference on Causal Learning and Reasoning, volume 213 of Proceedings of Machine Learning Research, pp- 726-751. PMLR, 11-14 Apr 2023.
Joris M. Mooij, Jonas Peters, Dominik Janzing, Jakob Zscheischler, and Bernhard Schölkopf. Dis- tinguishing cause from effect using observational data: Methods and benchmarks. Journal of Machine Learning Research, 17(32):1-102, 2016.
Arash Nasr-Esfahany, Mohammad Alizadeh, and Devavrat Shah. Counterfactual identifiability of bi- jective causal models. In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett (eds.). Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pp. 25733-25754. PMLR, 23-29 Jul 2023.
Weronika Ormaniec, Scott Sussex, Lars Lorch, Bernhard Schölkopf, and Andreas Krause. Stan- dardizing structural causal models. In The Thirteenth International Conference on Learning Rep- resentations, 2025.
Judea Pearl. Causality: Models, Reasoning and Inference. Cambridge University Press, USA, 2nd edition, 2009, ISBN 052189560X.
Jonas Peters, Joris M. Mooij, Dominik Janzing, and Bernhard Schölkopf. Causal discovery with continuous additive noise models. Journal of Machine Learning Research, 15(58):2009-2053 2014.
Alexander Reisach, Christof Seiler, and Sebastian Weichwald. Beware of the simulated dag! causal discovery benchmarks may be easy to game. In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan (eds.), Advances in Neural Information Processing Systems, volume 34, pp. 27772-27784 Curran Associates, Inc., 2021.
Alexander Reisach, Myriam Tami, Christof Seiler, Antoine Chambaz, and Sebastian Weichwald. A scale-invariant sorting criterion to find causal order in additive noise models. In A. Oh, T. Naumann, A. Gloherson, K. Saenko, M. Hardt, and S. Levine (eds.), Advances in Neural Information Processing Systems, volume 36, pp. 785-807. Curran Associates, Inc., 2023.
Patrik Reizinger, Yash Sharma, Matthias Bethge, Bernhard Schölkopf, Ferene Huszár, and Wieland Brendel. Jacobian based causal discovery with nonlinear ICA. Transactions on Machine Learning Research, 2023. ISSN 2835-8856.
Paul Rolland, Volkan Cevher, Matthaus Kleindessner, Chris Russell, Dominik Janzing, Bernhard Schölkopf, and Francesco Locatello. Score matching enables causal discovery of nonlinear additive noise models. In Kamalika Chaudhuri, Stefanic Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato (eds.), Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pp. 18741-18753. PMLR, 17-23 Jul 2022.
Pedro Sanchez, Xiao Liu, Alison Q O'Nei, and Sotirios A. Tsaftaris. Diffusion models for causal discovery via topological ordering. In The Eleventh International Conference on Learning Rep- resentations, 2023.
Shohei Shimizu, Patrik O. Hoyer, Aapo Hyvärinen, and Antti Kerminen. A linear non non-gaussian acyclic model for causal discovery. Journal of Machine Learning Research, 7(72):20332033 2006.
13

---

Published as a conference paper at ICLR 2026
Shohei Shimizu, Takanori Inazumi, Yasuhiro Sogawa, Aapo Hyvärinen, Yoshinobu Kawahara, Takashi Washio, Patrik O. Hoyer, and Kenneth Bollen. Directlingam: direct method for learn- ing linear non-gaussian structural equation model. Journal of Machine Learning Research, 12 (33):1225-1242)) 2011.
Peter Spirtes. Directed cyclic graphical representations of feedback models. In Proceedings of the Eleventh Conference Uncertainty in Artificial Intelligence, UAI pp. 491-498, San Fran- cisco, CA, USA, 1995. Morgan Kaufmann Publishers Inc. ISBN 1558603859.
Peter Spirtes and Clark Glymour. An algorithm for fast recovery sparse causal graphs. Social science computer review, 9(1):62-72, 1991.
Charles Stein. A bound for the error in the normal approximation to the distribution of sum of dependent random variables. Proceedings of the sixth Berkeley symposium on mathematical statistics and probability, volume 2: Probability theory, volume 6, pp. 583-603. University of California Press, 1972.
Eric V. Strobl and Thomas A. Lasko. Identifying patient-specific root causes with the heteroscedastic noise model. Journal of Computational Science, 72:102099, 2023. ISSN 1877-7503.
Xiangyu Sun and Oliver Schulte. Cause-effect inference in location-scale noise models: Maximum likelihood vs. independence testing. Advances in Neural Information Processing Systems, 36: 5447-5483, 2023.
Natasa Tagasovska, Valérie Chavez-Demoulin, and Thibault Vatter. Distinguishing cause from effect using quantiles: Bivariate quantile causal discovery. In Hal Daumé III and Aarti Singh (eds.), Pro- ceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pp. 9311 13323. PMLR, 13-18 Jul 2020.
Ruibo Tu, Kun Zhang, Hedvig Kjellstrom, and Cheng Zhang. Optimal transport for causal discovery. In International Conference on Learning Representations, 2022.
Antoine Wehenkel and Gilles Louppe. Unconstrained monotonic neural networks. In Hanna M. Wallach, Hugo Larochelle, Alina Beygelzimer, Florence d' Alché-Bue, Emily B. Fox, and Roman Garnett (eds.), Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada, pp. 1543-1553, 2019.
Johnny Xi, Hugh Dance, Peter Orbanz, and Benjamin Bloem-Reddy. Distinguishing cause from ef- fect with causal velocity models. In Forty-second International Conference Machine Learning, 2025.
Sascha Xu, Osman A Mian, Alexander Marx, and Jilles Vreeken. Inferring cause and effect in the presence of heteroscedastic noise. In Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato (eds.), Proceedings of the 39th International Conference Machine Learning, volume 162 of Proceedings of Machine Learning Research, pp. 24615- 24630, PMLR, 17-23 Jul 2022.
Sascha Xu, Sarah Mameche, and Jilles Vreeken. Information-tthorrticc causal discovery in topolog- ical order. In The 28th International Conference Artificial Intelligence and Statistics, 2025.
Zhuopeng Xu, Yujie Li, Cheng Liu, and Ning Gui. Ordering- based causal discovery for linear and nonlinear relations. In A. Globerson, L. Mackey, D. Belgrave, A. Fan, U. Paquet, J. Tomezak, and C. Zhang (eds.), Advances Advances in Neural Information Processing Systems, volume 37, pp. 4315 4340. Curran Associates, Inc., 2024.
Mengyue Yang, Furui Liu, Zhitang Chen, Xinwei Shen, Jianye Hao, and Jun Wang. Causalvae: Disentangled representation learning via neural structural causal models. In Proceedings of the JEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pp. 9593-9602, June 2021.
14

---

Published as a conference paper at ICLR 2026
Yue Yu, Jie Chen, Tian Gao, and Mo Yu. DAG-GNN: DAG structure learning with graph neural networks. In Kamalika Chaudhuri and Ruslan Salakhutdinov (eds.), Proceedings of the 36th International Conference Machine Learning, volume 97 of Proceedings of Machine Learning Research, pp. 7154-7163. PMLR, 09-15 Jun 2019.
Keli Zhang, Shengyu Zhu, Marcus Kalander, Ignavier Ng, Junjian Ye, Zhitang Chen, and Lujia Pan. geastle: python toolbox for causal discovery. arXiv preprint arXiv:2111.15155 2021.
Kun Zhang and Aapo Hyvärinen. On the identifiability of the post-nonlinear causal model. In Proceedings of the Twenty-Fifth Conference Uncertainty in Artificial Intelligence, UAI '09, pp. 647-655, Arlington, Virginia, USA, 2009. AUAT Press. ISBN 9780974903958.
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.
Yujia Zheng, Biwei Huang, Wei Chen, Joseph Ramsey, Mingming Gong, Ruichu Cai, Shohei Shimizu, Peter Spirtes, and Kun Zhang. Causal-learn Causal discovery in python. Journal of Machine Learning Research, 25(60):1-8, 2024.
15

---

Published as a conference paper at ICLR 2026
A DISCUSSION ON QPE AND IDENTIHABILITY
A. PROOFS FOR SECTION 3.2
Proposition 3.2 (QPE from CDF), 4 V /
Proof. For any function y f(x, u). if f(x,-) is strictly monotonic for any æ, then its inverse is (f-((x;)) (x.-)) -) Fixing and taking the total derivative w.r.t. of the implicit function y = f(x, u(ary)) yields
dy = Oz, f(x u) + 0af( u) 0s, u(x, v) = O,
dx.
Thus, - of Given the inverse relationship between the conditional quantile function and the conditional CDF, Qrix and F x are inverses of each other. Therefore, V= = - seelvix/Og / Furthermore, from the relationship between the conditional CDF and condi- tional PDF, 0g Fvix = Prjx, which implies V y) = Fr /aFFFx//PPxxxX
Proposition 3.3 (Causal Velocity is QPE). 1/yxx = a u) = u/,00 a.e.
Proof. From the previous derivation, we have V = - -V v //// Fixing æ and taking the total derivative w.r.t y of the implicit function y = =0 f(x, u( u)) yields
dy = 0., u) 0y (x, y) 1,
dy
Thus, = Therefore, // v. which proves the latter half of the identity.
Consider another SCM Y = g(X, W) where is uniformly distributed on [0, 1], such that g is strictly monotonic w.r.t. W and X I1 W. Now, let any be given. According the proba- bility integral transform, the conditional CDF Fyix transforms the conditional distribution py|x into pw. The inverse function w = ()(=,)))) also transforms the conditional distribution Py(x into pw. Since both Frix and w are strictly monotonic for continuous variables, they are both Knothe transports from Pylx to pw. Due to the a.c. uniqueness of Knothe transport, Fy w a.e Furthermore, as ENasr-Esfahany assumption, for any SCM that generates the observation distribution Px,., the exogenous vari- ables are identifiable up to an invertible transformation. Let this invertible transformation be- tween exogenous variables U and W be h such that U h(W).. Then V.u. h V..w and 0,u = S,w, which implies V / Vaw/O,w. Since Fy = w a.e., by Proposition 3.2 / F. F
al. 2023), under strict monotonicity and the Markov
By x a.e., which proves the first half of the identity.
Lemma 3.4. V I + E O log + 0,E 0 and tim EPYx 0 iff 1/1|x-
V Farri + Φr Pr = O. Under assumptions of smooth- ness and strict positivity, taking the partial derivative w.r.t. // gives
Proof. According to Proposition 3.2
V Vaprix + ( VaPy PrIx d + dYIi OxPYIx
PY|x
Pr
V log pr|p by x 0, logPy|x + 0, 0,
(3)
which means the equation holds when & = bby x- Additionally, according to Proposition 3.2 lagain, byyix x Pr|x V. Fy
lim v/ Prix +0o
lim V VaFy = V lim FYix 4-o
lim 1 = 0.
(4)
y-04000
40000000
The case for + - xo is analogous. From Equation 3 and Equation 4 the necessity of the lemma established.
Now, for any & satisfying the conditions, subtract the equation it satisfies from the equation that d satisfies obtain
( x/rxx) O log 1ogr 00, x/yx) 0.
16

---

Published as a conference paper at ICLR 2026
Then, multiply both sides of the equation by pr|x and expand 0 log pr|x, which yields
(5 - x--by1x) QgPy x + PrIX 0,(5 x = -,5((5 - 4blxx) 0.
Integrating this w.r.t y from some yo gives
(E (yixx)piXX + C
where C(x)) is a function independent of Taking - too for this equation on both sides, since E S tim
and lim b/ X PY|X = 0, it follows that limy-++o C = C 0. Thus,
(E x py|x 0 holds for all æ, y. Due to the strict positivity assumption, Py > 0, so must be that 5 x = O, which establishes the sufficiency of the lemma.
Theorem 3.6 (Identifiability of QPE in Finite Linear Span), For each variable X; and any æ, Vy|x,(11 æ) E xp implies (xxxx, ( 0 for any y. If: (i) The components in Trix are linearly independent; (ûi There exists y such that (y # 0; (iii) For each basis function 0j. limy-+ioe Q3 PrIx 0; then the converse also true.
Proof. Let arbitrary be given in this proof. Substituting 1 xxxxyyla x)) Eド= (15)) (5() into the second-order cross PDE (Equation 2 and rearranging terms yields
k
2z 0g logXx, + Σ
ca. 0z(0; 0, logpx,y + t,) = -X,,Y + Σ Ca j=1
Tr, = 0, (5)
j=1
where 7/y denotes the j-th component of n/Ιx. This implies that sx.. and fy are linearly dependent. Taking the m-th partial derivative w.r.t. for m = O,..., k gives
0g SX.Y
+
Cij Og TriX =0
for all O...., k. Note that these ls equations (including the zero-th order) share the same coefficients. Therefore,
SX..Y 0
1 0
M 8, X
Ga We =
:
01
Cis
for any Since e is not identically zero, det W = 0, which means Wy Wy 7xxx) 0.
Assume (sx,, 7)(xx) =0 for any y. By the Bacoer Theorem (Bocher 1901), the two additional conditions (i) and (ii) guarantee that sx..Y and 77 linearly dependent. Thus, there exist 1, G.I..... Cas such that Equation 2 holds. We only consider the zero-th order equation, i.e., Equation holds. Integrating it w.r.t. y from some yo yields
2z log
Σ (os By logxxx + 13) + + 0,
(6)
where C(x) is a function independent of Now, taking the expectation of both sides w.r.t. Py|X, notice that
[ 0x.Px,y dy = =d,, log Px.
E 0r log X 2]
px
Furthermore, Φ Og logpx.y+ is the Stein operator acting on Φ w.r.t. x (since oy log Pr|x 0, logPx.y) Thus, by Stein's 's identity (under the vanishing condition (iii)),
o dy log + o |X = x = 0,
for each j. Also, E (x)) X m| C(x)). For Equation 2 to hold, C = -d.; logpx. Rearrang- ing Equation further, we get
dz logprix +
Σ ai os 0y logpy py x + 0,
Σ cis 23 0.
17

---

Published as a conference paper at ICLR 2026
Given the boundary condition tim Og pr|x Σnm Gij Øj. Then
0, let each component of be &.
k
lim 4ー0
PX
Caj
lim oy = 0. y
j=1
Therefore, by the equivalence stated in Lemma 3.4 E = Vylx,, i.e, v/1xx.(112)) E span(().
A.2 PROOFS FOR SECTION 4.1
Corollary 4.2, Assume that there are basis functions Φ and certain regularity conditions hold. Then: (i) dim(Ox oo;
oo; (ii) dim(Oxxr < k
Proof. Here, we can consider the log-likelihood logpx: given the assumption of absolute positiv- ity. First, by Markovianity and strict monotonicity, there exists equation by change-of-variabl formula:
log PY log Pa log T,f log log 0yu,
where U can be expressed by X, Y. and f. Thus, log Prjx only depends on f and pv, which are assumed to given. So log py
known function. Since the joint log-density is log px,y = log px + log Pyix. for any given l/o. any PDE on logpx,y can be simplified to ODE solely on logpx. Furthermore, any partial or mixed derivative w.r.t. y will only involve the given log py x-
Next, when the forward version of Assumption 3.5) holds, the Wronskian determinant described in Theorem 3.6 is an identity involving only partial or mixed derivatives w.r.t. y. This is completely independent of log px. In other words, the forward version of Asssmmttoo 3.5 not impose new constraints, which implies 0x-Y; = @, and thus dim (xx-Y: X0.
Conversely, when the backward version of Assumption 3.5) holds, the Wronskian determinant in Theorem requires an ODE for log to be satisfied. This ODE can be rewritten the form
(log px) (k+2) = G(x,, log (logPxY,,.. (x) (
whose highest order is k 2++2 This high-order nonlinear ODE can be reduced k 2-dimensional first-order ODE system. According to the Picard-L.nneeie theorem, under certain regularity condi- tions (e.g., global Lipschitz continuity), the solution to this ODE exists and is unique, determined by k + initial conditions (log px (2o), (log (xz)), .... (log px ) (K++1) (xx))) Therefore, log P.x be described by these k 2 parameters. If degenerate cases exist, fewer parameters may be required.
Hence, dim e x
< 2
Discussion on Corollary 4.2
As stated in the main text, this corollary generalizes the conclusion by (Hoyer et al. 20008) that unidentifiable ANMs lie in three-dimensional affine space (since ANMs linear ODE). It also fills the gap left by
correspond to and the Wronskian simplifies previous identifiability works (Immer et al., 2023 Strobl & Lasko 2023) regarding quantifying unidentifiable HNMs (HNMs correspond to derived identifiability ODEs for ANMs and HNMs using causal velocity. but they did not further quantify the solutions for HNMs. In contrast to these works, [orollary 4.2 provides quantified identifiability for HNMs and extends this conclusion to any finite-linearly spanned QPE (k < 00).
2). Similar
work, Xi et al. (2025) also
A.3
QPE IS IRRELEVANT TO MONOTONIC MECHANISM AND MARKOV PROPERTY
Xi et al., 2025) demonstrated that causal velocity, specifically QPE as in Proposition 3.3] inde- pendent of the precise form of noise and latent mechanisms. This relaxes the strength of parametric assumptions on FCMs. Given that QPE is an observational quantity, inferable solely from obser- vational distributions, we have relaxed the counterfactual assumptions typically required for causal velocity, including strict monotonic mechanisms and the Markov assumption, which are crucial for counterfactual identification (Nasr-Esfahany al. 2023 Chen & Du 2025). Therefore, as claimed in the main text, QPE fully relaxes assumptions regarding the latent causal mechanisms, shifting fo- entirely to the observational distribution's 's shape. This implies that "observational causal discov- ery," even requiring identifiability of causal direction, can operate entirely within the observational layer of the causal hierarchy.
18

---

Published as a conference paper at ICLR 2026
Next, we provide two simple examples demonstrating that even if the underlying causal mechanisms do not satisfy strict monotonicity or the Markov property, their observational distributions' QPEs are identical and still satisfy Assumption 3.5]
Example A.1 (Strictly Monotonic Mechanism but Semi-Markovian). Let X = Z + W and Y = exp(X X Z + U), where the confounding variable and the exogenous variables U, W independent standard normal. Then the QPE d æ) - (2x + 0.5)y.
Proof. Utilizing the property of linear Gaussians, consider Z X W. Then (Z X z) (0..5. 0.5), and simultaneously (Z + U X = x) ~ (0.5a, 1.5). Let T ~ (0, 1.5), so is equivalent T + 0.5x. Thus, the conditional event Y
7 + U X
X holds if and
only if x² + + 0.5z & log(y). Therefore, Fr (ylx) = P(Y S X = x) = P(T < log(y) - x² -²2 0.5x) = (logl)( x² - -0.5z) = h(z y). We can calculate O = hzz, ((z,))(22r 0.5) and 0yFi = h(,, v) which are Fyix (u/z) and Fr(x(ylz) respectively. By Proposition 3.2 obtain dbyxx (y) x) = (2r 0.5)y-
(~2r
Example A.2 (Non-Monotonic Mechanism and Semi-Markovian), Let X Z + W and Y exp(X + Z + + where the confounding variable Z and the exogenous variables U, W are inde- pendent standard normal. Then the QPE x x) 1.5y..
Proof. Utilizing the property of linear Gaussians, consider Z X W. Then (
(0.55a 0.5). Since [ ~ xi and is independent of Z and X, (2+ = x) ~ NO 0.5)+
xi. Let T ~ (0, 0.5) + xí. then Z + |
is equivalent to T 0.5x.. Thus, the conditional
event Y < y X = holds if and only if x + T + 0.5x < log(v). Therefore. Fr|x(ylz)= P(Y S
y X = x) = P(T < log(y) 1.5z) = (v) - 1.5x) = h(r,y) y). We can calculate Fr =
12, y) and 0, = h(,, (1/b). which are 0, x (r) and 1, (2/z) respectively. By
Proposition 3.2, obtain d (y 2) 1.5y.
In Example A.1 Z x correlates with x, implying that the non-endogenous term exp(Z + U) and the endogenous term exp(X) are not independent, thus violating the Markov property. However, since it can be equivalently written some random variable T + 0.5z. where T is conveniently eliminated when computing QPE. it "coincidentally" still satisfies [Assumption 3.5) This means its observational distribution can even be equivalent to an HNM (the basis is y). In Example A.2 the principle for violating the Markov property is similar to Example A. with the difference being that the non-endogenous term exp(2 Z + 2²2) is non- non-monotonic w.r.t the exogenous variable. Yet, it can still be equivalently written as some random variable T + 0.5z.. The existence of these two examples reveals series of special cases where, even the underlying SCM does not satisfy strict monotonicity or the Markov property, its observational distribution still exhibits par- ticular characteristics in QPE. In such cases, assumptions on FCMs fail, but assumptions on QPE remain valid.
is crucial to emphasize that despite QPE being independent of the underlying causal mechanisms, it still represents parametric assumption. This parametric assumption is empirically verifiable when the analytical form of the observational distribution is fully known. For instance, we can compute the Wronskian determinant given in Theorem 3o empirically determine if this assump- tion holds. However, in reality, one typically only has samples from the observational distribution. Therefore, whether the QPE parametric assumption holds still requires careful consideration in real- world problems, similar to FCMs, to mitigate potential risks arising from assumption violations.
19

---

Published as a conference paper at ICLR 2026
B DISCUSSION ON FICO
B.1 PROOFS FOR SECTION 5.1
Theorem 5.2. For each covariate X,, Assumption 5.1 holds then
E (vy)xx)) (sr) = E (sx,)² E ((xx,)) +
2 2 0.yrI,,)
Proof. Rewrite the equation described in Lemma 3.4 as:
8x = rx, - v sy 0y1,,
where #fyjx, is abbreviated as vh. Squaring both sides and taking the expectation, the following equation still holds:
E [(xx = -
[( By + -
(rx, -2 +0,4)]
where E |xx, (V/ sy + 0,,,) = |x E ( sy + + X = æ1]. And E ( + 0, X x m]
0 due to Stein's 's identity, assuming y 16 pxx, 0. Thus, the entire equation simplifies to:
E [(sx:)]] = I3. [(xx,)
E [(y sy + 0,4))²]
Now consider E [(V) sy + 2)
= E (sy )2| + [0, + +2 &l 2,00. Also, according
to Stein's 's identity, if x 3 0, pxp = 0, then
18 [sy a 0, + ( 0,,ر =0
Therefore, E (ay & [,yو = - [ [(0 m5. Equation 7 yields the theorem.
Substituting these simplified terms back into
Corollary 5.3. For each covariate X3,
Assumption 5.7 holds then
E [(00,vypx,)) 2)² + x (y < <E [((vy]x,)) (sx) if and only F [(xx)]] E. [(x )3].
9 + [r]
Proof. According to Theorem 5.2 E (sx,) ) - E [(s)) equals
E [(xxx,) )²] + E [((v)Ix,;) ) 1) (sv) E [(0 (0, v 0, x
Therefore, this inequality holds iff E [(xx)]]] - E [(>v )²] > 0, which proves the corollary.
B.2 FICO UNDER HETEROSCEDASTL GAUSSIAN ASSUMPTION
First, we provide the formal statement and proof of the corollary to Corollary 5.3 regarding the HNM assumption, discussed in the main text:
Corollary B.1. For each covariate X;, given Assumption 5.1 holds. If OyyIx, = 0 and
E (ry)
E [(0,Ayixx))) Var [(VYIx))) Var ((x)))
E [(WV)Ix))'] >
E [sy )
E (sy )
then E [ssx,)=]= > E [s.
Proof. When 836 [sxx, )
[[sy]] equals
E (rxx.)
+ [( (sx) - [(a,))}]
where
E
2 (sy )21 = E 2 E [sy) + Cov (v² (sy)
And according to the Cauchy- Schwarz inequality.
Cov (sv))) 2 Var 2 Var (sy }²)]
20

---

Published as a conference paper at ICLR 2026
Assuming E (sy }2] 0, then follow the condition in this corollary,
E [(xx,))] + E 2 E (sy ) E [(0,))]] Var [+ Var [(s)] ) B[(sy)
Since the inequality for Cov (v/2, (sy holds, it implies
E [(xx,)]²] + [ E [(sY) E [(),b)]]] + Cov (+i (sr)3)) [[(sy) Rearranging this gives
[(xx.)*] + [(V2 1) (sy) [{0,,)] [(sx,) E
[[sy )"] > 0,
which completes the proof.
Corollary B.I qualitatively requires the squared coefficient of variation (CV) of the QPE to be as small possible. Next, we show how to derive the simplified explanation presented in the main text under the linear heteroscedastic Gaussian assumption, and how to extend this to the nonlinear case to obtain qualitative explanation. Without loss of generality, for simplicity, we only consider the bivariate case below. The multivariate case only requires replacing the univariate functions or their derivatives in this section with multivariate functions their partial derivatives.
Assumption B.2. The conditional distribution PYxx is Gaussian for any x. i.c. I -x=-- (u((z), (o()))
Using the expression for the Gaussian distribution and Proposition 3.2) is straightforward derive the QPE dy x =p² (y fo)). For any H. 7 = (V -M)/a ~ (N0 1), so v/yx -k¹ - Zoo. Thus, it can be derived that
E [wby x) 2 X x]
EE [(p' + Zo')² (1/)² (d')² E [(p' + Z (u')* 6(')) +3(0))
E [(by x) x =
where x: is given, so yl and a' are treated as constants. For standard normal variable Z, its moments are E [Z] the squared QPE as:
[z]] = E, [ 0, [ = 3. This further leads to the mean and variance of
E [(] E E [(Wr]x 3 =
=]] =x]
E [(W]] +
Var [Wyxx))²
= E Var (1xy|ix)
+ Var E [ x xx=+]] =x]]
E E [(Vr)x)* |X= ---E- E² =x [ =2]1 + (4')² + (a')²
E [((+')* + (())) ((')))+ + Var [u') + Var (('') (0')3)
2 (2)))))(()2++"+)))))))
by the law of total expectation and law of total variance. Let = 0' //'.. Then the mean and variance can be written as:
E [Vy|x E (1+))
Var [( 2E [m'y'* x² (2+ Var (4d) x²)]
When μa and a are linear functions, 2e' and o' are constants w.r.t. X, and thus k is also a constant. In this case, the CV simplifies
Var [(Vy[))]]
/ 2x
|x]]
(8)
[Wyix)²
1 1 x²
Note that this function strictly increases as [x| + oo and its limit is V2. Therefore, to minimize
the CV under the linear heteroscedastic Gaussian assumption, |k], i.e., /0//1'), should be as small as possible. This means either o' is as small as possible (e.g.. degenerating to ANM when o' = 0) or pl is sufficiently large relative a' such that y dominates the distribution shape.
21

---

Published as a conference paper at ICLR 2026
When j and QPE, yielding
general, assuming / S |x|
We can still extract [x] from the moments of
E [((vvix)
2 (1+ [(x')]]
Var [+ - 2 (1')* x² (2 + x²)] + E (1+ {
$ (1 + 6u + 3x²) E [('''1]
Substituting this into the CV expression, we obtain upper bound for the CV:
Var [(vy x² E [(by)x)
vI + 6u 3+ 1+12
T (1'']] E ((1))
To minimize the CV, the upper bound must be minimized. This requires two conditions: (i) The upper bound of |x] should be as small as possible, and the lower bound as large as possible. This means |x] should not vary significantly (it is constant under the linear assumption). Additionally,
because the growth rate of the numerator's 's upper bound is greater than that of the denominator's s
lower bound, [t| should generally be small as possible. (ii) v ((x")*)/ [(=')²] should be sufficiently small. According to the Cauchy- Schwarz inequality, its lower bound is 1 (with equality when y' is constant, i.e., under the linear assumption).
22

---

Published as a conference paper at ICLR 2026
C EXPERIMENT DETAILS AND RELATED WORKS
C.I BIVARIATE CAUSAL DISCOVERY
Datasets Section 6. 6. enumerates the bivariate causal discovery datasets used in our experiments, along with their sources. Most of these datasets follow conventions from previous work (Immer) et al. 2023 Xi et al. 2025) and are standardized. For the Tübingen cause-effect pairs challenge, we manually removed discrete pairs with IDs "47, 52, 53, 54, 55, 70, 71, 105, 107". These pairs visibly and severely violate the continuity assumption of distributions, which most bivariate causal discovery algorithms rely on. This selection may differ from previous literature, potentially leading to lower baseline results reported in this paper.
Additionally, we created a set of constrained QPE synthetic datasets, including Qd-V, Sig-V V, Rbf-V, and NN-V. These datasets are generated from the random optimization of the following objective:
Ux
arg min 0,0
VY +
+ Var [vy]xx,]]
0yus
where no is randomly initialized causal flow (see Equation T and choose a random Gaus- sian mixture (3 components, with variance 0.25 S S. 2.0 for each component) as the exoge- distribution py for the causal flow. x is constrained QPE satisfying [ssumption 3.5 (i.e., given basis functions), randomly initialized with parameters & (which control the coelfccent functions). Qd-V, Sig-V, Rbf-V, and NN-V correspond quadratic basis (1, y, y²), sigmoid basis (1, syymoid(y), sigmoid(2y)), RBF basis (1, (1,xxp(- exp((y"), cxpp=""p)xxppp-pp-ep + 0.5)²), exp(- - )), and neural network basis (with tanh activation function), respectively. À regularization parameter, set to À 0.1, to prevent the random optimization from making the QPE too extreme, which could lead to excessively irregular observational distributions.
Baselines In the main text, we primarily use baseline methods as examples. In Appendix D.2 we provide a comprehensive evaluation using 21 open-source baselines based on various theories to fully and transparently demonstrate the superiority of our method. These baselines include:
Linear Non-Gaussian Acyelic Model: ICA-LiNGAM (Shimizu et al. 2006).
-AR-LINGAM
(Hyvärinen et al. 2010), Direct-LiNGAM (Shimizu et al., 2011).
Additive Noise Model: ANM (Hoyer et al. 2008). CAM (Bühlmann et 2014), RESIT (Peters
et al. 2014), RECI (Bloebaum et al.) 2018). CGNN (Goudet et al. 2018)) CDS (Fonollosa 2019).
Post Non-Linear Model: PNL (Zhang & Hyvärinen
2009).
Heteroscedastic Noise Model: QCCD
Tagasovska et al
2020), CAREFL Khemakhem et al.
202I). HECI (Xu et al. 2022). GRCI Strobb & Lasko
2023). CDCI (Duong & Nguyen 2022)
LOCI (Immer et al.) (023)..
Minimum Description Length: IGCI Daniušis et al.
2010), SLOPE Marx & Vreeken 2019a)
SLOPPY (Marx & Vreeken 2019b).
Optimal Transport: DIVOT (Tu et al. 2022
Causal Velocity Model: CVEL (Xi et al. 2025
Among these methods, SLOPPY includes variants using different information criteria (AIC and BIC). CAREFL and LOCI include variants using likelihood and HSIC (Hilbert-Schmidt Indepen- dence Criterion) (Immer et al. 2023; Sun Schulte 2023). CDCI and CVEL involve multiple configurations or combinations. For all these methods, we tuned their hyperparameters to select the best configuration for fair comparison with our method. However, we only tested the ANM variant of DIVOT due to its deprecated dependencies for PNL implementation.
While several methods share our goal of generalizing the functional forms of LINGAM, ANM, or HNM, they achieve identifiability via assumptions that differ from our setting, such as require- ments on the data generation process (Guo et al. 2023). contrastive learning frameworks (Reizinger et al. 2023), access to multi- data (Jalaldoust et al. 2025). In contrast, method, QPE, identifies the causal direction using only the observational distribution. Accordingly, our main com- parisons are with methods designed for this purely observational context.
23

---

Published as a conference paper at ICLR 2026
Runtime Regarding hardware, PNL, CGNN, CAREFL, LOCI, CVEL, and QPE-f utilize GPU acceleration for neural networks, while others run on CPU. As for environments, implementations of SLOPE, SLOPPY, QCCD, and RECI are based on R, while others based on Python. Each experiment is conducted using the same hardware, system, and default runtime configurations, The average time taken for cause-effect pair under these conditions has been reported in Table 2
Complexity analysis For QPE-k, can analyze its complexity. Let N be the number of samples, M be the number of test locations, and T be the number of test samples used in the OLS test. The most computationally intensive part is the estimation of the response matrix in OLS, which has complexity of (TT. In our experiments, we set M = T = 20, and test locations are uniformly distributed within [-2.5, 2.5] since the datasets are standardized.
For QPE-f, which is based on neural networks and random optimization, its training process is simpler compared CVEL involving higher order terms. QPB-f only needs to be trained according to Equation | Although calculating log [0z, ug| also involves gradients, it can be expressed in closed form using pre-defined discrete transformations of the causal flow, as detailed in/pppeniix D.I
C.2 MULTIVARIATE CAUSAL ORDERING
Datasets
The configurations for the synthetic datasets are detailed as follows:
ANM-GP and HNM-GP are generated by the following processes, respectively:
X, = GP(PP;; + GMM(U, o.), X; = =GP(PP;)..) + GP(P;; GMM o.),
where P; denotes the parent variables in the DAG, GP the random fourier features Gaussian process, and GMM is a Gaussian mixture model. All parameters, including ,,, 0,,,, 1, 0r,2, M4g i are randomly selected according to a prior.
The gCastle (Zhang et al. 2021) synthetic process can be represented as X; = NN(P;: 0,) + (0, where NN Is a neural network mechanism, Ao 1) is the standard normal distribution, and the neural network parameters are randomly chosen.
The LiNGAM (Shimizu al., 2006 synthetic process can be represented as X, a. + 6; + Gumbel 1), where Gumbel 0,, 1) is the standard Gumbel distribution, and the coefficients and the bias of the linear equation are randomly selected.
To demonstrate robustness under FCM assumption violations, we also include confounded ver- sions of ANM-GP and HNM-GP, denoted ANM-GP-c and HNM-GP-c, respectively:
X; GP(Z, P:; 0.) + GMM( o:), X: GP(P:; + GP(Z, P.p.,.) GMM(u; o:),
where Z ~ GMM do) is common confounding variable.
Except for gCastle, other datasets are internally standardized to iSCM (Ormaniec et al. 2025) to eliminate potential sortability. Specifically, during data synthesis, we are required to add extra step such that X, + (X, E(X,))/ Var[X.] before processing its children.
Baselines In the main text, we list the open-source baselines involved in our experiments. Below are the assumptions underlying these methods, providing context for interpreting the experimental results:
as synthetic SCMs may contain unintentional "fingerprint" information. Ormaniec al. (2025) proposed a simple method to eliminate such artifacts.
Sortability: Var-Sort (Reisach et al. 2021 and R2-Sort (Reisach et al. 2023).
2006) and Direct- Linear Non-Gaussian Acyelic Model: ICA-LiNGAM (Shimizu al.
LiNGAM (Shimizu et al.)
(01T) for linear mechanism and non-Gaussian noise.
Additive Noise Model: SCORE (Rolland et al. 2022), based score function with Gaussian 2024), based on score function with strong parent influence (and its criterion is equivalent FICO according to Stein's identity). RESIT (Peters et al., 2014), based on HSIC test of residuals. NoGAM (Montagna et al. 2023), based on score function of residuals, supporting arbitrary noise.
noise. CaPS (Xu et al.
24

---

Published as a conference paper at ICLR 2026
Heteroscedastic Noise Model: HOST (Duong & Nguyen 2023).
based on Shapiro-Win Wilk test for gaussian residuals. SKEW (Lin et al. 2023). based on score function with symmetric noise.
Minimum Description Length: Topic (Xu et al. 2025 based on cntemattiongeommetne theory.
Compared to the above baselines, FICO is characterized by fewer assumptions, being solely related to QPE rather than based a specific FCM. This broadens its theoretical applicability. How- ever, while FICO has a wider theoretical scope, it does not guarantee superior performance over these models. This is because, ent characteristics within the baselines may still enable them to perform effectively, For instance, Montagna al. 2023) showed that score function based methods exhibit robustness assumption violations in a wide range of experiments.
when their specific theoretical assumptions are violated, inher-
We do include other algorithms that primarily identify causal graphs, such NOTEARS (Zheng itt al. 2018), Grad-GNN (Yu et al. 2019), and DAGMA Bello et al. (2022), despite their ability identify unique causal graphs (due an implicit ANM assumption). These methods mainly output causal graphs, not causal orderings, which preclude accurate calculation of the OD and ODR metries. Furthermore, this paper focuses exclusively on the causal ordering task. The graph pruning process is omitted because, under consistent pruning strategies, OD and ODR sufficiently reflect the discrepancy between the final DAG and the ground truth.
Metrics For causal ordering tasks, we primarily use Order Divergence (Rolland et al. 2022 to measure how well the output causal order aligns with the underlying true DAG. Specifically, for an output causal order ıT and adjacency matrix A of the underlying true DAG, define
d ODR(T, A) = OD x A)
OD( R A)
Σ Σ A.jj i=1 n.>my
A..--
In other words, OD reflects how many edges in the DAG violate the causal order and ODR scales OD to the unit interval. When they are both zero, is exactly a topological order of the true DAG.
Runtime Regarding hardware, all methods can be run on CPU. Additionally, all score functions are estimated using a kernel-based score matching algorithm from (Rolland et al.) (022). This eliminates the need for extensive GPU training, providing fast and fair baseline, albeit with some sacrifice in accuracy. For the environment, all methods are based on Python. We run each experiment using the default configurations of these methods. The average time taken for causal ordering is detailed in TTble 11
Complexity analysis In our implementation, the primary bottleneck for FICO comes from the score function estimation in each iteration. According to the implementation by (Rolland et al.) 2022 assuming N the sample size and d is the dimension, each estimation step requires O + 22d). Since FICO requires d estimations in total, the overall complexity is O(N' N2d).. This implies that when >> d, the sample size becomes the bottleneck, leading to severe inefficiency. For large samples, denoising score matching (Sanchez et al. (023) can be applied, which amortize the complexity over multiple steps, resulting in a single-step optimization complexity of only O(N max(d, c)), where is the number of model parameters.
25

---

Published as a conference paper at ICLR 2026
D ADDITIONAL RESULTS
D.1 PYPER-PARAMETER TUNING FOR QPE-F
Hyperparameters for QPE-f A primary source of hyperparameters for QPE-f arises from the necessity to initially fit a causal flow (Equation[] D]. Specifically, we employ discrete causal flow, where the causal flow ug is parameterized as composition of several invertible transformations:
Ti T2,0
T2,0,
where each invertible transformation o (a. y) is a monotonic function w.r.t. y. Candidate choices for these transformations include affine transformations (Dinh et al.) 2017). Rational-Quadi Spline (RQS) transformations (Durkan et al. 2019), Monotonic Neural Networks (MNN) (Huang
et al. 2018 and Unconstrained Monotonic Neural Networks (UMNN) (Wehenkel & Louppe, 2019 These methods exhibit progressively increasing representational capacity. In our experi- mental setup. affine transformations were excluded they restrict the model to representing only
HNM (Khemakhem et al. 202T), thereby failing to generalize to broader model classes
The second source of hyperparameters for QPE-f pertains to the selection of the hypothesized QPE model employed for the linear span test. We consider the following choices for the hypothesized model: (i) Hypernetwork: an unconstrained neural network that implicitly satisfies the hypothesis; (ii) Low-rank network: ΣK (x( 0.) (y:: Φ:). where both a, and b, are unconstrained neural net- works, explicitly representing finite number of basis functions in the hypothesis; (iii) Polynomial network: Σ (z: 0.) y" explicitly assuming polynomial basis functions.
For training both the flow and the hypothesized QPE models, the Adam optimizer was employed with a learning rate of 0.01 and a weight decay of 0.001. To ensure continuity, networks utilized the Sil.U activation function and comprised two hidden layers, each with a width of 100. Training proceeded for total of 1000 epochs, and the model corresponding to the epoch with the minimum loss function value was selected for basis test. The training and validation sets were not partitioned.
Hyperparameter tuning We conducted hyperparameter tuning across several dimensions: the choice of transformation type (RQS, MNN, UMNN), the number of composite transformations (t E {1, 5)), the hypothesized QPE model (Hypernetwork, Low-rank network, Polynomial net- work), and the specific hyperparameters for the hypothesized QPE models (e.g., rank for Low-rank networks, degree for Polynomial networks). Due to space constraints, the hyperparameter tuning results for the SIM and Tue datasets are presented in Tables 4 and respectively.
Table 4: Hyperparameter tuning on the SIM dataset. Best per flow configuration is bolded.
ROS
MNN
UMNN
Test
Hyper
0.75 (0.76)
0.67 (0.60)
0.57 (0.52)
0.80 (0.83)
0.75 (0.77)
0.78 (0.82)
0.79 (0.78)
0.76 (0.75)
0.67 (0.68)
2
0.76 (0.77)
0.68 (0.60)
0.57 (0.53)
0.87 (0.87)
0.75 (0.
0.75 (0.76)
0.79 (0.80)
0.78 (0.81)
0.73 (0. 73)
Lowrank
0.76 (0.77)
0.64 (0.56)
0.54 (0.47)
0.50 (0.81)
0.73 (0.74)
(0.72)
0.83 (0.82)
(0.74)
10
0.75 (0.76)
0.61 (0.56)
0.51 (0.18)
0.51 (0.81)
0.79 (0.79)
0.73 (0.75)
0.77 (0.73)
0.76 (0.76)
0.72 (0.72)
0
0.87 (0.88)
0.82 (0.76)
0.76 (0.74)
0.83 (0.83)
0.85 (0.85)
0.85 (0.85)
0.88 (0.86)
0.87 (0.87)
0.73 (0.77)
Poly
1
0.84 (0.s3)
0.78 (0.75)
0.69 (0.08)
0.83 (2)
0.79 (0.80)
0.79 (0.84)
0.85 (0.83)
0.86 (0.85)
0.75 (0:77)
0.81 (0.79)
0.73 (0.70)
0.09 (0.09)
0.81 (0.81)
(0.76)
0.78 (0.63)
(0.83)
0.81 (0.82)
(0.73)
Table 5: Hyperparameter tuning the Tue dataset. Best per flow configuration is bolded.
ROS
MNN
UMNN
Test
5
5
Hyper
0.46 (0.41)
0.47 (0.49)
0.56 (0.48)
0.59 (0.61)
0.59 (0.53)
0.49 (0.51)
0.58 (0.56)
0.48 (0.45)
0.51 (0.45)
2 0.43 (0:36)
0.52 (0.52)
0.52 (0.47)
0.81 (0.61)
0.56 (0.54)
0.51 (0.46)
0.53 (0.4)
0.43 (0.46)
0.52 (0.45)
Lowrank
0.47 (0.37)
0.53 (0.65)
(0.46)
0.66 (0.67)
0.49 (0.47)
(0.57)
0.42 (0:45)
0.51 (0.45)
0.47 (0.37)
0.51 (0.52)
0.56 (0.46)
0.57 (0.66)
0.63 (0.65)
0.48 (0.44)
(0.66)
0.42 (0.46)
0.48 (0.43)
0
0.52 (0.45)
0.56 (0.51)
0.60 (0.49)
0.56 (0.57)
0.56 (0.53)
0.58 (0.62)
(0.4))
0.45 (0.35)
0.47 (0.46)
Poly
0.48 (0.40)
0.56 (0.53)
0.60 (0.53)
0.62 (0.60)
0.59 (0.57)
0.58 (0.19)
0.57 (0.52)
0.56 (0:48)
0.54 (0.42)
0.55 (0.522)
(0.51)
(0.52)
0.70 (0.79)
(0.61)
0.49 (0.45)
(0.57)
0.58 (0.80)
0.58 (0.52)
26

---

Published as a conference paper at ICLR 2026
Optimal hyperparameter configurations
Table 6 provides the empirically determined optimal
hyperparameter configurations for each dataset, thereby ensuring experimental reproducibility.
Table 6: Empirically optimal hyperparameter configurations for each dataset.
Dataset
Transform
Test
K
Dataset
Transform
Test
K
Dataset
Transform
Test
K
AN
UMNN
Poly
0
SIM-In
MNN
Poly
D4-s2c
MNN
Lowrank
ANe
MNN
Poly
0
Cha
MNN
Lowrank
10
Per
RQS
Hyper
LS
UMNN
Poly
0
Net
ROS
Poly
Sig
RQS
Poly
LS-s
UMNN
Poly
Multi
UMNN
Poly
Vex
UMNN
Lowrank
10
MNU
RQS
Poly
Tue
MNN
Poly
2
Qd-V
MINN
Poly
SIM
UMNN
Poly
D4-s1
UMNN
Lowrank
Sig-V
UMNN
Poly
SIM-C
M.N
Poly
0
D4-s2a
MNN
Poly
Rbf-V
UMNN
2
Poly
SIM-B
Poly
Poly
NN-V
Lowrank
D.2 COMPARATIVE EVALUATION QPE-K, QPE-F, AND BASELINES
Tables to present the comprehensive results of 20 distinct methods across 24 datasets. In this extensive comparison, QPE-f consistently achieves state-of-the theth--rtt performance across the vast ma- jority of datasets. with only minor exceptions where is marginally outperformed by certain meth- ods. direct comparison between QPE-f and CVEL,, as well as QPE-k, empirically validates the assertion made in (Section 4.3) regarding the superior accuracy of flow-based methods in QPE esti- mation. Furthermore, CVEL and QPE-f consistently outperform other methods in both non-ANM and non-HNM datasets (specifically, flow synthesized Per, Sig, Vex; and constraint QPE synthesized Qd-v, Sig-V, Rbf-V, NN-V). This empirically corroborates the theoretical insight that the identifia- bility conditions for causal velocity or QPE are more broadly applicable than those only for ANM and HNM.
D.3 CONVERGENCE OF FICO
Figure 4 fillustrates the convergence behavior of score-based methods (SKEW, SCORE, and FICO) w.r.t. sample size HNM-GP datasets. All methods utilize the same score function estimation algorithm, which is known to asymptotically approach the true score function with increasing sample size. Figure 4 specifically demonstrates how the accuracy of these methods in the causal ordering evolves as the score function estimates become increasingly precise.
HNM-GP (SF. 100)
HNM GP
200
100
0
100
2000
500
800
1000
3000
3000
100
300
500
600
1000
2000
3000
Mullivaritts Classal Ordering Methods
SKEW
SCORE
Figure 4: Convergence behavior of SKEW, SCORE, and FICO HNM-GP datasets.
Notably, the SKEW method exhibits degradation in performance the score function estimate improves, primarily because underlying assumptions (e.g., symmetric noise) are violated this context. In contrast, SCORE and FICO maintain consistent convergence profiles, likely attributable to the satisfaction their intrinsic assumptions or properties. Furthermore, FICO consistently out- performs SCORE the sample size gradually increases.
27

---

Published as a conference paper at ICLR 2026
Table 7: Accuracy (and AUDRC) of QPE-K, QPE-f, and 21 baselines on 9 bivariate datasets.
Method
AN-1
1.8
1.S+
MNU
SIM
SIM-c
SIM-g
SIM-In
ICA-LINGAM
0.63 (0:81)
6:69 (0.61)
0:63 (0.61)
0.63 (0.61)
0:63 (0.61)
0.63 (0.61)
0.63 (0.51)
0.63 (0.61)
0.63 (0.61)
VAR-LINGAM
0.66 (0.07)
0.00 (0.00))
0.11 (0.11)
0.001 (0.00)
0.00 (0.00)
0.42 (0.36)
0.47 (0.25)
0.27 (0125)
0:23 (0.14)
Divees LINGAM
0.06 (0.07)
0.00 (0.00)
0.11 (0.11)
0.001 (0.00)
0.00 (0.00)
0.42 (0.36)
0.47 (0.25)
0.20 (0.25)
0.23 (0.14)
ANM
0.45 (0:35)
0.47 (0.42)
0.46 (0.50)
0.45 (0.45)
0.40 (0.37)
0.45 (0.56)
0.49 (0:46)
0.41 (0.42)
0.46 (0.45)
CAM
1.00 [1.00]
1,00 (1.00)
1.00 (1.00)
0.51 (0.58)
0.91 (0.58)
0.50 (0.58
0.50 (0.5
0.50 (0.79)
0.85 (0.90)
RESIT
0.59 [1 00)
1.00 [1:00)
0.72 (0:70)
0.08 (0.07)
0.01 (0.00)
0.78 (0.78
0.52 (0.55
0.76 (0.721
0.67 (0.00)
RBC1
0.18 (0.27)
0.35 (0.34)
0.22 (0.16)
0.44 (0.88)
0.23 (0.36)
0.44 (0.40
0.53 (0.533
11.39 (1.330
0.64 (0.40)
CGNN
0.96 (0.94)
0.57 (0.57)
0.92 (0.88)
0.64 (0.64)
0.54 (0.55)
0.75 (0.76
0.76 (0.7%
0.72 (0.590
0.75 (0.07)
CDS
0.89 (0.58)
0.99 (1 00)
0.76 (0.79)
0.05 (0.06)
0.70 (0.74)
0.71 (0.55
0.76 (0.55
0.75 (11776
0.05 (1.63)
PNL
0.30 (0.31)
0:40 10.46)
0.30 (0.35)
0.40 (0.53)
0.58 (0.58)
0.46 (6.
500
0.63 (0.40
0.42 (0.57)
QCCD
1,000 [0:
0.65 (0.75)
1.00 00)
1.00 (1.00)
1.00 (1.00)
0.68 (0.60)
0.77 75)
0.67 (0/63)
0.57 (0.87)
CAREFL
1:00
1.00
(1.00)
1.00
1.00 (1.00)
0.82
JICCI
0.98 (0.97)
0.55 (0.60)
(0.86)
0.35 (0.00)
0.53
(0.54)
0.56 (0.52)
0.55
GRCI
(0.61)
(0.65)
0.47
(0.87)
0.52
1.00
1.00
(0.57)
(0.97)
(0.53)
(0.90)
(0.56)
(0.91)
0.56
(0.50 (0.37)
0.96
(0.54)
0.14
(0:21)
(0.11)
(0.20)
0.00
0.65
SLOPPY
1,00 (1.00)
1.00 (3.00)
1.00 (3.00)
0.56 (0.50)
0.96 (0.97)
(0.70)
0.822 (0.62)
0,82 (U.80)
0 (1.50)
DIVOT
0.62 (0:58)
0.69 (0:73)
0:45 (0.47)
0.69 (0.68)
1.00
0.68 (0.72)
0.47 (0
0.80 (0.61)
0.63 (0.68)
CVEL.
1,00 (1.00)
0.98 [1 00)
0:08 (0.97)
0.93 (0.91)
0.94 (0.98)
0.63 (0.66)
0.72 (0.77)
0.90 (0. (05
0.76 (0.73)
QPE-
0.99 (0:06)
0.88 (0.75)
1.00 (1.00)
0.78 (0.50)
1.00 (1 11.00)
0.53 (0.85)
0.70 (0-53)
0.83 (0.81)
0.68 (0.69)
O
1.00 (1.00)
1,00 (1.00)
1.00 (1.00)
0.99 (1.00)
1.00 (1.00)
0.59 (0.56)
0.88 (0.52)
0.96 (0.84)
0.92 (0.91)
Table 8: Accuracy (and AUDRC) of QPE-K, QPE-f, and 21 baselines on 8 bivariate datasets.
Method
Cha
Multi
Tue
D4-11
D4-42h
T)4-s2e
ICA- ICA-LINGAM
0.52 (0.57)
1152
1.52
0.64
0.67 (0.63)
0.50
0.52 (0.06)
0.51 (0.96)
VAR-LINGAM
0.53
(0.62)
Direct LINGAM
(0.52)
0.35
(0.72)
(0.52)
(0.55)
ANM
(0.37)
0.47 (0.66)
(0.67)
(0.70)
0.35 (0.52)
(0.47)
CAM
(8.48)
0.78 (1.83)
0.345 (0.36)
0.55 (1.54)
0.42 (0.57)
0.35 (0.31)
0.44 (0:29)
(0.49)
RESIT
0.74 (1.79)
0.76
0.37 (0.4.)
(0.61
(0.72)
0.63 (0.00)
0.55 (0.58)
RECI
0.56
0.90
0.85 (0.85)
0.64 (0.55)
0.58 (0.64)
0.58 (0:09)
0.50 (0.64)
0.52 (0.55)
CGNN
0.61 10.50
0.75 (0. 761
0.84 (0.83)
0.60 00.68)
0.50 (0.65)
0.50 (0.50)
0.47 (0.57)
0.50 (0.54)
CDS
0.71 10. 77)
0.78 (0.90)
0.54 00.000 (0.65)
0.67 0..00
0.58 (0.20)
0.50 (0. 60:
0.54 (0.53)
0.58 (0.57)
PNL
0.45 (0.36)
0.50 (0.44)
0.05 (0.43)
0.51 (0.54)
0.28 (0.19)
0.45 (0:112
0.59 (0.62)
0.29 (0.40
QCCD
0.55 (0.54)
0.81 (0.455
0-09 (0.90)
0.66 (0.73)
0:35 (0.52)
0.55 00:53:
0.53 (0:01)
0.53 (0.58)
CAREPL
0.72 (0.
0.85 (0.983)
0.76 (0.79)
0.463 (0.60)
0.08 (0.64)
0,09 150
0,63 (0.68)
0.56 (0.
EIECT
0.57 00.
0.72 (0.711
0.01 (0.90)
0.61 (0.57
0.02 (002)
0.56
(0 00:
0:50 10 05)
0.47 10.50)
GRCI
0.53 (0.5 551
0.58 00.53)
0.61 (0.55)
0.55 (0.56)
0.07 (0.70)
0.51 (0.57)
0.1 10.70)
0.50 10.5
LOCE
0.73 (0.77)
11ST (11888)
0.79 (0.30)
[L:11 (0.67)
0.58 (0.55)
0.899 (0.74)
0.61 (0.06)
0.54 (0.47)
CDCT
0.67 (0.78)
0.84 (0.50)
0.92 (0.055
0.68 (0.78)
0.67 (0.09)
0.08 (0.65)
0.60 (0.59)
0.61 (0.55)
IGCI
0.55 (0.54)
0.57 (01.55)
0.68 (0.05)
[LG2 (0.67)
0.42 (0.27)
0.44 (0.46)
0.43 (0.42)
0.40 (0.4-)
SLOPE
0.56 (0.39)
0.61 (0.00
0.88 (0.85)
0.57 (0.55)
0.35 (0.52)
0.51 (0.50)
0.40 (0:52)
0:30 (0. (5)
SLOPPY
0.48 (0.40)
(0.89)
0.06 (0.47)
0.64 (0.65)
0.35 (0.16)
0.35 (0.25)
0.44 (0:40)
0-18
DIVOT
0.49 (0.38)
0.34 (0.58)
0,38 (0.42)
(0.54)
0.57
0.55 (0.49)
0.56 (0.51)
CVEL
(0.77)
0.02 (0.60)
0.97 (0.94)
(0.50)
(0.70)
(0.63)
(0.57)
QPE
(0.53)
0.99 (0.90)
0,88
0.54 (0.47)
0.58 (0.75)
0.67 (0.69)
0.61 (0.64)
0.64 (0.61)
OPE-T
(0.97)
0.96
(0.75)
(0.72)
(0.60)
(0.18)
Table 9: Accuracy (and AUDRC) of QPE-K, QPE-f, and 21 baselines on 7 bivariate datasets.
Method
Per
Sig
Vex
Qd-V
Sig-V
Rbf-V
NN-V
ICA LINGAM
0.63 (0.61)
0.63 (0.01)
0.63 (0.61)
0.63 (0.61)
0.63 (0.51)
0.63 (0.51)
0.63 (1.81)
VAR LINGAM
0.6 (0.63)
0.36 (0.36)
0.48 (0.40)
0.57 (0.55)
0.39 (0.54)
0.61 (0.35)
0.85 (0.89)
DINGSLINGAM
0.6 (0.63)
0.37 (0.36)
0.47 (0.47)
0.58 (0.50)
0.50 (0.56)
0.500 (0.54)
0.56 (11.27)
ANM
0:40 (0.51)
0.41 (0.50)
0:30 (0:38)
0.49 (0.57)
0.50 (0.36)
0.43 (0.56)
0.15 (0.50)
CAM
0.00 (0.00)
0.08 (0.06)
0.24 (0.11)
0.12 (0 16)
0.47 (0.46)
0.30 (0.22 22)
0.23 (0.22)
RESIT
0.70 (0.72)
(0.22)
0.03 (0.00)
0.50 (0.52)
(0.45)
RDCI
6.00 (0.00)
0.07 (0.06)
0.94 (0.94)
0.63 (0.50)
(0.52)
20)
0.49 (0.52)
CGNN
(0.80)
(0.66)
0.71
(0.87)
0.554
CDS
(0.00)
(0:43)
(0.46)
(0.53)
QCCD
6.02 (0.01)
(2)
(0.00)
0.34 (0)
(0.56)
26)
0.53
CARHEL
0.85 (0.97)
(0.65)
(0.922
0.72 (0.78)
0191 (0.88)
(0.54)
(0.82)
HBCT
0.01 (0.00)
(0)
0.94 (0.94)
0.20 (0 2011
0150 (0.48)
0.31 (0.20)
1145 (0.51)
GRCI
0.56 (0.48)
(0.54)
0.54 (0.55)
0-47 (0
0.51 (0.32)
(0-16
100 (0.55)
LOCI
0.56 (0.97)
0.70 (0.08)
0.87 (0.56)
0.71 (0-72)
0.57 (0.58)
0.61 (6)
0.75 (0. (0.79)
CDCI
0.48 (0.47)
0.42 (0.45)
0.40 (0.52)
0.74 (0:73)
0.50 (0.77)
0.57 (0 00:58
0.72 (0.72)
10C1
1,00 (1.00)
0.77 (0.s1)
0.57 (0u80)
0.47 (0:24)
11.99 (11222
0.22% (1:23)
11:00 (0.54)
SLOPE
0.000 (0.00))
0.06 (0.0)
0.59 (0.50)
0.61 (0.37)
(1522 (1.46)
0.18 (0.19)
11:4 (0.45)
SLOPPY
0.62 (0.02)
0.11 (0.06)
0.10 (0007)
0.17 (0.14)
0.48 (0.19)
0.32 (0.26)
1135 (0.54)
DIVDT
0.87 (0.06)
0.52 (0.64)
0.05 (0.04)
0.52 (0:35)
0.44 (0.19)
0.63 (0. (6.55)
1147 (0.44)
CVEL
1.00 (1:00)
0.54 (0.81)
0.96 (0.97)
0.91 (0.52)
0.94 (0
0.52 (0.34)
0.57 (0.90)
QPC-k k
6.77 (0.70)
0.50 (0.84)
0.63 (0.60)
0.42 (0.44)
0.67 (0.73)
0.68 (0.96)
0.53 (0.49)
QPE-1
1.00 (1.00)
0.90 (0.
0.91 (0.92)
0.91 (0.91)
0.91 (0.90)
0.94 (0.94)
0.90 (0.90)
28

---

Published as a conference paper at ICLR 2026
D.4 FICO'S PERFORMANCE UNDER HETEROSCEDASTER GAUSSIAN ASSUMPTION
We analyze the interpretability of FICO under the heteroscedastic Gaussian assumption concerning the mean and variance functions in Appendix B.2 Figure 5 details this relationship. The synthetic (
dataset with heteroscedastic Gaussian noise is generated by X; (P,(.) + os (P)) 1), where N4 h(:; a, b a, 1) and 0; (og(exp(h(:: ci d,, a, B)) + 1). Here, h(; a, bb, a, B) is an affine function with sin perturbations:
E PJE
C.y Cas P +a Σ
sin sin(Pr.))
PJEP.
where coefficients a.j,, Ci.j ~ (O, 1) and biases b;, d; ~ NO,, 2). Hyperparameters a and / control the magnitude of the sin perturbation and the gradient of the affine function. respectively. Indirectly, 3 controls the magnitude of |^| = = (since / for 14 is always 1), while a controls how closely 14 and Ø; approximate linear functions. Each cell in Figure 5] shows the average 100 sub-tests, where corresponding sub-tests share the same coefficients and biases, varying only hyperparameters and
11:55
0:44
3:43
0.43
5:43
09
18.0
0.9
644
3.4
5:41
08
10.00
0140
8:47
1.47
2.235
06
039
0:33
9.71
10.9
4:41
0.7
8:
4.55
0.49
45
11:55
0:46
1.49
07
0.3
6:1
0:44
04
0.44
0.44
6.42
4:45
0.45
06:2
0:43
2.43
04-7
0.50
2.03
9.475
06
037
3.35
AUR
4:43
0:44
3.40
0.5
0:47
1.40
4.
4:40
MAAZ
0:44
2:43
MARY
0:45
2.02
0.5
125
9.77
0:44
3.40
044
0:30
5.50
0:42
9:45
0.51
006
3.
4.455
0:37
9.99
0.6
0.26
011
023
085
112
06
69
1.2
1.8
(a) ER,
(b) ER, d 20
Figure 5: Relationship between FICO's ODR and hyperparameters a and 3 under the heteroscedas- tic Gaussian assumption.
(a) 10-variable ER graph: (b) 20-vari ER graph. The expected num- bers of edges in these graphs is 4 times their dimensions.
The results indicate that FICO's performance gradually degrades as a and 3 increase. This reflects that the magnitude of [s] and the linearity of μu; and a; indeed affect the validity of Assumptioo 5.4 FICO performs best when |x| sufficiently small and both /4 and a are linear. As the assumption is progressively violated, the performance weakens, which empirically consistent with our analysis.
D.S FICO ON WORLD DATASETS
Table 10: ODR on real-world datasets for different methods. The best is bolded.
As shown in Table TO) on real-world datasets, score function based methods generally perform poorly on Sachs, likely due to the assumption not holding. However, CaPS and FICO achieve the best perfor- mance Syntren.
Method
Sachs
Syntren
R2-Sort
0.29
0.89 = 0.07
Var-Sort
0.82
0.50 0.17
ICA-LINGAM
0.59
0.42 0.19
Direct-LiNGAM
0.47
0.62 1 0.10
HOST
0.18
0.42 0.10
D.6 FICO SYNTHETIC DATASETS
RESIT
0.47
0.74
TOPIC
0.59
1.183
NoGAM
0.65
0.39 0.09
presents the experimental results across 8 synthetic datasets under types random graphs (SF and ER). It is important to note that ODR exhibits hidden baseline at 0.5, which corresponds to random ordering (since the probability of n, My is 0.5 in
Figure 6
SKEW
0.71
0.49
SCORE
0.71
0.38 + 0.10
CaPS
0.71
0.33 0.08
FICO
0.71
0.33 - 0.08
29

---

Published as a conference paper at ICLR 2026
&Custile (SF)
@astile (ER)
ODK
LINGAM (SF)
LINGAM (ER)
1.0
0.8
0DR
0.5
0.4
0.2
0.0
ANM-GP (SP)
ANM-GP (ER)
ODR
0.2
0.0
HNM-GP (SF)
HNM-GP (PR)
1.0
0.5
OBR
0.5
ANM-GP-c (SE)
GP (ER)
0.5
0.6
0.4
0.2
0.0
HNM-GP- (SF)
HNM-GP-c (ER)
ODK
Mulrowiitte Oaikming Methods.
BKD
STORE NoGAM
SKEW
Dive-UINAM ICN
VanSon R2-Ses
CUPS
TOPIC
HOST
Figure 6: ODRs of FICO and baselines on 12 multivariate causal discovery datasets. Lower is better.
random ordering). An ODR value less than 0.5 indicates that the method effectively retains more than half of the edges, implying that the underlying assumptions characteristics met within the dataset. Conversely, an ODR greater than 0.5 suggests that these assumptions or characteristics likely to be violated. When ODR approaches 0.5, the method essentially performs level equivalent to random ordering.
The results presented in Figure reveal several key observations:
30

---

Published as a conference paper at ICLR 2026
The elimination of "fingerprints" by iSCM leads ODR values exceeding 0.5 for both Var-Sort and R2-Sort, indicating performance worse than random ordering. This suggests low sortability in these datasets, thereby mitigating the possibility of "hacking" through these specific metrics. ICA-LiNGAM and Direct-LiNGAM perform well exclusively in LiNGAM datasets, approximat- ing random ordering in other contexts.
RESIT exhibits unexpectedly strong performance GP datasets. This could be due to GMM noise, particularly in high gigh-dimnnsional settings. This phenomenon warrants further investigation. Conversely, its performance in gCastle and LiNGAM datasets is moderate, nearing random order- ing, potentially due to the Gaussian and Gumbel distributions.
HOST demonstrates significantly superior performance in low-dimensional gCastle scenarios, which may also explain best performance on Sachs. However, in other contexts, its ODR consistently remains above 0.5, violating assumptions, largely because its normality check on noise becomes progressively challenging in higher dimensions.
TOPIC exhibits variability. It outperforms other methods in gCastle. Yet, its performance is suboptimal elsewhere. For instance, it notably above 0.5 in LiNGAM and approximates random ordering in other datasets.
SKEW's performance is generally slightly lower than other score function based methods due to the violation of its symmetry noise assumption in most datasets, with the exception of the Gaussian noise in gCastle.
The remaining
function based methods exhibit comparable performance. In high- dimensional settings, CaPS and FICO are marginally superior to SCORE and NoGAM (with a potential difference of s 0.05). Given their formal equivalence, any performance disparities be- tween CaPS and FICO are solely attributable to precision differences.
Overall, score function based methods demonstrate the robustness reported in Montagna et al. (2023). Despite potentially being slightly outperformed by other methods on specific datasets, their ODR values consistently remain below 0.5 across all experimental settings. This empirically sug- gests that the underlying assumptions or characteristics these methods are implicitly satisfied.
D.7
FICO'S RUNTIME EFFICIENCY
Table 11 presents FICO'S runtime performance in comparison other baseline methods.
Table 11: Runtime efficiency of FICO and baselines, in seconds per sub-test.
Method
d= 5
d= = 10
d= 20
d = 50
d= 100
R2-Sort
0.000 0.000
0.000
0.000
0.000
0.000
0.000 + 0.000
0.005 0.007
Var-Sort
0.000 0.000
0.000 0.000
0.000 0.000
0.000:00.000
0.000 0.000
IIN.M
0.010
0.043 0.035
0.128 0.073
0.638 0.155
2.200 0.636
Direct-
LINGAAL
0.010
0.050
0.338 0.000
5.006
39.942 0.507
HOST
0.625:0.000
1.967 0.106
4.159 0.207
7.177 0.330
17.332 - 0.640
TOPIC
0.318
1.984 0.541
9.381 2.604
60.869 13.916
250.343 = 45.157
RESIT
0.268 4
1.091 0.051
4.772 0.268
00.000220022
227.662 9.178
NoGAM
1.243 0.299
4.251 + 0.593
14.864 1.028
92.578 6.251
370.316 + 16.281
SKEW
0.543 + 0.362
1.201 0.624
2.548
0.876
8,022 + 1.248
20.043 + 1.136
SCORE
0.831 0.462
1.883 0.737
4.306
0.982
14.180 1.914
37.609 2.007
CaPS
0.455 + 0.087
1.074 + 0.056
2.761 + 0.285
10.822
1.037
33.794 + 3.501
FICO
0.425 1 0.322
0.797 + 0.364
1.727 1 0.523
5.550 1 0.943
13.538 1 1.248
Among the methods evaluated, excluding those based on sortability and LiNGAM, FICO demon- strates the fastest execution, particularly in high-dimensional settings.
31

---

Published as a conference paper at ICLR 2026
THE USE OF LARGE LANGUAGE MODELS
Paper Writing: LLM was used to polish the text of this paper. We ensure that all modifications were manually verified to be accurate and consistent with the authors' intended meaning.
Programming Assistance: L.M assisted in implementing the QPE-k algorithm (Section 4.2p and generating code for some figures (Figures 1 3 and 6 All code developed with LLM assis- tance underwent unit testing and manual debugging. All other experimental procedures and data processing were performed manually.
Runtime Environment Fixes: LLM provided suggestions for resolving issues with older baseline runtime environments in newer versions and assisted in connecting and Python, allowing all experiments to be conducted within Python.
32


Reviews:
[{"review_id": "J0UUB035d9", "reviewer": "Reviewer_mFRa", "summary": "The paper studies causal discovery from observational data using distribution shape constraints. Authors motivate non-graphical assumptions about the gradient of the joint distribution that emphasize the asymmetry between the cause and effect in a bivariate case, and extend the results to multivariate causal discovery by iteratively finding the last variable in the causal order. Through examples, they show how the shape constraints introduced in their work is distinct from the graphical constraints commonly employed in this problem, and experiments validate usefulness of their approach.", "questions": "Please address my concerns above.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "1. Clear statement of the assumptions underlying the theoretical results.\n2. Motivating examples throughout the paper allow the more curious reader to understand the motivation.\n3. Algorithm 1 is impressively simple, and the connection between fisher information and the shape constraints are quite insightful.", "weaknesses": "1. Appendix A3 contains very important remarks that deserve to be presented in the main text. The paper argues that the shape assumptions are more *relaxed* than the existing assumptions employed in causal discovery, and this subtle points needs to be established carefully, or else the reader may find the contributions not well-motivated.\n2. Sections 3.3 and 3.4 are not contributing to what comes later in the manuscript; they are useful extensions of 3.2, though, I recommend a restructuring as they hurt the coherency.\n3. Missing discussion about related work [1,2,3,4] that would enrich the work.\n\n[1] Hyvarinen, A., Sasaki, H., & Turner, R. (2019, April). Nonlinear ICA using auxiliary variables and generalized contrastive learning. In The 22nd international conference on artificial intelligence and statistics (pp. 859-868). PMLR.\n\n[2] Jalaldoust, K., Salehkaleybar, S., & Kiyavash, N. (2025). Multi-domain causal discovery in bijective causal models. In Proceedings of the Fourth Conference on Causal Learning and Reasoning.\n\n\n[3] Guo, S., Tóth, V., Schölkopf, B., & Huszár, F. (2023, December). Causal de finetti: On the identification of invariant causal structure in exchangeable data. In Thirty-seventh Conference on Neural Information Processing Systems.\n\n[4] Reizinger, P., Sharma, Y., Bethge, M., Schölkopf, B., Huszár, F., & Brendel, W. (2023). Jacobian-based causal discovery with nonlinear ICA. Transactions on Machine Learning Research.", "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": "lrhexnM618", "reviewer": "Reviewer_W1w7", "summary": "The paper proposes an approach to causal discovery from observational data, based on the notion of quantile partial effect (QPE). In particular, three distinct approaches utilizing QPE are proposed:\n\n(i) A kernel-based, non-parametric approach for bivariate discovery which assumes that the QPE can be represented as a linear combination of known basis functions,\n\n(ii) A parametric approach for bivariate discovery using the idea of normalizing flow, where each flow is parameterized by a neural network,\n\n(iii) A parametric approach, which recovers a causal ordering among a set of variables, based on the Fisher information matrix.", "questions": "(Q1) Is the representation in Table 1 actually true?: specifically, for PNL-ANM and PNL-AHM models, the basis functions $\\bar g, g^{-1}\\bar g$, etc. are needed. Isn’t it generally the case that the function $g(\\cdot)$ is not known in these settings? How does this tie with the fact that the basis functions in $\\Phi$ need to be known?\n\n(Q2) Relating to the previous point, is there any nuance w.r.t. choosing the basis? The choice of the kernel bandwidth parameter is mentioned as a limitation, but the choice of basis (which is even more fundamental?) is not really mentioned.\n\n(Q3) The performance of QPE-f seems to be very strong empirically. I wonder, for settings where the data-generating model is actually an ANM, does the ANM method perform better? (similar question for LiNGAM & others). If I understand, current experiments would seem to imply the answer is no? If so, it is curious.\n\n(Q4) Relating to previous point — given the strong performance of QPE-f, is there any way to move beyond bivariate discovery with such an approach?", "limitations": "", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "(S1) Causal discovery is an important topic in causal inference, and further advances in this context are needed. The paper deals with an important topic.\n\n(S2) The framing of previous methods (LiNGAM, ANM, HNM) as specific cases of the general setting in which QPE lies in the linear span of fixed basis functions, seems interesting and novel. Getting a generalization of current approaches assuming different functional constraints is a good direction to explore.\n\n(S3) The performance of the flow-based (neural) approach for QPE testing seems to perform quite well empirically.", "weaknesses": "(W1) The exposition of the paper could be substantially improved. A range of assumptions are discussed, and in the current version, it is non-trivial to follow for which of the proposed approaches exactly which assumptions are used. This could be streamlined through an overview table.\n\n(W2) Causal Velocity comparison is misleading: Some of the language around causal velocity / bijective causal models seems concerning (abstract, after Proposition 3.3). \n\nWhen a bijective causal model assumption is invoked, some collapse of the layers does occur. This seems to be an important point for causal velocity.\n\nQPE is defined purely at the observational level, and does not pertain to any structural assumptions — it is simply a functional of the underlying observational distribution. Therefore, comparing these two in such a way (saying that “QPE does not need a monotonicity assumption”) seems quite confusing. It seems more likely that (i) the definition of QPE does not require monotonicity; but (ii) for the proposed discovery methods to work, monotonicity is assumed?\n\n(W3) Removal of the Markovian assumption: the claims on the removal of the Markov assumption are very curious. Most work in causal discovery using functional assumptions considers the Markovian case. If the authors are truly relaxing this assumptions, this seems like a big deal, and should be emphasized strongly. If not, explaining exactly how the Markov assumption is “not needed” is very important — is the case again that QPE can be defined without invoking the Markov assumption, but to apply QPE-based discovery methods, the Markov assumption is needed? If the latter is true, these claims need to be reframed properly.\n\n(W4) It is not uncommon for causal discovery methods to provide some form of formal guarantees that the method will work with high probability etc. Are the formal guarantees needed / missing from the paper?\n\n(W5) Line 289 states that “OLS assumes a Gaussian distribution”. OLS is a method, which does not assume any distribution implicitly. Of course, with Gaussian noise, OLS may exhibit specific properties (become the MLE etc.). Clarifying this would be important.\nFurthermore, in this setting arbitrary distribution of $x_t, y_t$ can now be handled — is this previously not the case? This is definitely worth clarifying.\nThis lack of clarity on the exact assumptions connects back to point (W1).\n\n(W6) In Lines 338, it is said that $\\partial_y \\psi_{Y | X,i} = 0$ is “e.g. ANM”. To be clear, in this context, additive noise is assumed in addition to the fully parametric model used in this setting? \n\nANM usually allows arbitrary functions $a(X)$; which is not the case here? It would be more precise to say “additive noise” than ANM in this case?\n\nSame comment for the “e.g. HNM” model.\n\n(W7) Assumption 4.5 seems to be important for the FisherInfo-based approach. I currently have no idea how one might judge the validity of such an assumption, or even start thinking about it? It seems that this is insufficiently discussed, and raises the concern whether the exact condition needed for the method is assumed, regardless of its plausibility.\n\nI note that I would be willing to reconsider my grade in light of the authors' responses to weaknesses / 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": "yhnrW0ZE2N", "reviewer": "Reviewer_vus6", "summary": "The authors propose estimating the quantile partial effect (QPE) as a means of functional causal discovery. The QPE generalizes causal velocity to the distributional level, which avoids making SCM assumptions. They derive novel conditions for both identifying causal direction in the bivariate case as well as sink nodes in a multivariate setting. Based on these conditions, they propose two methods for bivariate causal discovery (QPE-f, QPE-k) and an order-based approach for the multivariate case (FICO).", "questions": "- Question on your Appendix A.2: You show that, given a joint distribution (up to $p(x)$) the class of forward models $p(y|x)$ are non-parametric, while the class of backward models $p(x|y)$ are parametric, based on an extension of the ANM argument in Hoyer et al., 2008. My question is this: does the (k+2)-th order ODE for $\\log p(x)$ somehow depend on the anti-causal model, either via $\\eta_{X \\mid Y}$ or via the basis functions $\\phi(x)$ in Assumption 3.5? And if it does, will changing these unknown properties also change the ODE? \n\n- Question on FICO: can you clarify the relationship between FICO and SCORE? You still have to estimate the score function for FICO, correct? SCORE looks for the smallest variance in score across nodes, while FICO looks for the smallest FI, or squared score. But asymptotically under mild conditions the FI and variance are the same (the experimental results also look very similar). Is Corollary 4.3 a condition on when we expect SCORE to work for non-linear, non-additive noise models?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 2, "strengths": "- The authors propose a Wronskian criterion for bivariate causal discovery that is not only theoretically interesting but also directly informs novel methodology based on a practical test, which shows strong performance in the common benchmarks. \n- Mathematical erivation and discussions (especially in the Appendix) are complete and rigourous\n- Experiments are exceptionally comprehensive and easy to follow \n- Overall writing is clear, the paper is well-structured and logical", "weaknesses": "- The theoretical generalization beyond causal velocity is limited in practice, since the discussion in Appendix A.2 on identifiability still requires the SCM assumption which equates QPE and causal velocity. \n- Assumption 4.5 for FICO is specifically designed for Corollary 4.3 and is difficult to interpret in practice (although I appreciate the authors acknowledging this).\n- The resulting algorithm, FICO, seems very similar to the SCORE algorithm (Rolland et al., 2022).", "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": "vzaUDfuQOX", "reviewer": "Reviewer_w7Hq", "summary": "The paper proposes Quantile Partial Effect (QPE) as an observational-level object for causal discovery. QPE is defined via conditional quantiles and admits multiple equivalent forms, enabling estimators that do not assume a specific functional causal model (FCM). The core theory shows identifiability of causal direction when each component of the QPE lies in a finite linear span of basis functions, formalized by Theorem 3.6. Building on this, the authors give two bivariate procedures: (i) a kernel QPE with an OLS basis test (QPE-k) and (ii) a flow-based QPE with a neural basis test (QPE-f)—and report strong empirical results. For multivariate discovery, they leverage a PDE link between QPE and the score function to derive Fisher-Information Causal Ordering (FICO), which greedily removes leaf nodes using marginal score variances under an additional assumption (Assumption 4.5). Experiments show QPE-f competitive or superior to many bivariate baselines, and FICO performs competitively on synthetic and real benchmarks.", "questions": "a. How is the assumption on lines 96–97 used downstream? e.g., in which lemmas/theorems or in the estimator design. and what would break without it?\n\nb. The paper claims that QPE does not require the (causal) Markov assumption. However, some steps appear to rely on independence of noise terms (exogeneity)?\n\nc. How should practitioners select basis families for Assumption 3.5 in new domains? Could the authors provide data-driven diagnostics or model selection among candidate bases?\n\nd. How sensitive are ODRs to the estimator class and its hyperparameters in FICO?\n\ne. Why Sachs is hard? Can the authors analyze where Assumption 4.5 breaks on Sachs?", "limitations": "", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "a. The equivalence of QPE and the Wronskian-based identifiability criteria are clear and self-sufficient at the observational level.\n\nb. Unifies causal-velocity style reasoning with quantile/CDF views and score-based analysis.\n\nc. Two practical bivariate pipelines: kernel, OLS basis test (fast), and flow, neural basis test (accurate). And, for multivariate cases, FICO is simple to implement given score estimators.", "weaknesses": "a. The central finite-span assumption (Assumption 3.5) is elegant but potentially restrictive; practical guidance on choosing/validating basis sets is limited. How sensitive are conclusions to basis mis-specification?\n\nb. Assumption 4.5 for FICO appears to lack easy a-priori diagnostics.  The performance on Sachs suggests it fails in practice, which narrows applicability.\n\nc. Multivariate ranking depends on accurate scoring function estimation; the paper uses first-order scores for speed but does not deeply analyze robustness to estimator bias or noise.", "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)