Distinguishing Cause from Effect with Causal Velocity Models
Authors: Johnny Xi
Hugh Dance
Peter Orbanz
Benjamin Bloem-Reddy
Venue: ICML 2025 poster (2025)
Keywords: causal discovery
functional causal models
dynamical systems
Primary area: general_machine_learning->causality
TL;DR: We use causal velocity functions, which specify distribution-free causal mechanisms, and which can be recovered from the score function of the data, for bivariate causal discovery.

Abstract: Bivariate structural causal models (SCM) are often used to infer causal direction by examining their goodness-of-fit under restricted model classes. In this paper, we describe a parametrization of bivariate SCMs in terms of a *causal velocity* by viewing the cause variable as time in a dynamical system. The velocity implicitly defines counterfactual curves via the solution of initial value problems where the observation specifies the initial condition. Using tools from measure transport, we obtain a unique correspondence between SCMs and the score function of the generated distribution via its causal velocity. Based on this, we derive an objective function that directly regresses the velocity against the score function, the latter of which can be estimated non-parametrically from observational data. We use this to develop a method for bivariate causal discovery that extends beyond known model classes such as additive or location-scale noise, and that requires no assumptions on the noise distributions. When the score is estimated well, the objective is also useful for detecting model non–identifiability and misspecification. We present positive results in simulation and benchmark experiments where many existing methods fail, and perform ablation studies to examine the method's sensitivity to accurate score estimation.

Distinguishing Cause from Effect with Causal Velocity Models
Johnny Xi Hugh Dance? 2
Peter Orbanz?
Benjamin Bloem-
Abstract
discovery possible, one common approach is make func- tional assumptions, for example, non-linear additive noise (Hoyer et al., 2008) or location-scale noise (Xu et al., 2022; Strobl & Lasko, 2023; Immer et al., 2023) on the underlying causal process, and decide on the causal direction based on model fit and/or complexity in both candidate directions.
Bivariate structural causal models (SCM) are often used to infer causal direction by examin- ing their goodness-of-ft fft under restricted model classes. In this paper, we describe a parametriza- tion of bivariate SCMs in terms of a causal ve- locity by viewing the cause variable as time in a dynamical system. The velocity implicitly defines counterfactual curves via the solution of initial value problems where the observation specifies the initial condition. Using tools from measure transport, we obtain a unique correspondence be- tween SCMs and the score function of the gener- ated distribution via its causal velocity. Based on this, we derive an objective function that directly regresses the velocity against the score function, the latter of which can be estimated nonparamet- rically from observational data. We use this to
Likelihood- -ased approaches seem natural here, as they can be used both for model estimation and evaluation. However, they require the full specification of both the mechanism and the noise distribution. This is undesirable due to the risks of model misspecification, particularly with respect to the noise distribution, which is not directly involved in causality, but can still lead to incorrect causal inference (Schultheiss & Bühlmann, 2023). To avoid that risk, the most promi- nent example of functional method is to fit an additive noise model (ANM) using nonparametric regression, which avoids needing to model the noise distribution. Instead of a likelihood- measure, an independence score (c.g.. HSIC, Gretton et al., 2005) between the estimated residuals and cause variable can be used to determine goodness-of of fit (Mooij et al., 2009; Peters et al., 2014). This approach has been extended to location-scale models (LSNM) by Strobl & Lasko (2023); Immer et al. (2023). Separately, Rolland et al. (2022); Montagna et al. (2023b) derived other goodness-of- fit criteria for ANMs that avoid structural model estimation altogether. Instead, those methods extract the signal directly from estimates of the score function V log p(r, y) under the assumption of Gaussian noise, Subsequent work used score estimation in conjunction with explicit model estimation in the case of non-Gaussian noise (Montagna et al., 2023a). However,
develop a method for bivariate causal discovery that extends beyond known model classes such as additive or location-scall noise, and that re- quires assumptions on the noise distributions. When the score is estimated well, the objective is also useful for detecting model non-identifiability and misspecification. We present positive results in simulation and benchmark experiments where many existing methods fail, and perform ablation studies to examine the method's sensitivity to ac- curate score estimation.
unaware of any work functional causal discovery methods beyond the ANM and LSNM classes, which may still be misspecified in many settings.
1. Introduction
Distinguishing cause from effect from purely observational data is challenging task. It is generally impossible to distinguish causal direction X ++ Y from an anti-causal direction Y + X without intervention, as the correspond- ing causal graphs are Markov equivalent. To make causal
In this work, we develop a new framework that treats bi- variate causal models as dynamical systems, and we use it to devise a new estimation procedure based on the score function of the data, applying it to causal discovery. Specifi- cally, we parametrize causal model via its causal velocity, which describes infinitesimal counterfactuals and can be directly recovered from the score in a simulation-free way analogous to the recent literature on flow- how-based generative modeling (Lipman et al., 2023; Albergo 2023), Assuming that the marginal and joint score functions can be estimated accurately, our proposed method combines
'Department of Statistics. University of British Columbia "Gatsby Unit, University College London. Correspondence to: Johnny Xi <johnny utat.ube.ca>. jtat.b..i@@stt..ubb.......
Vajden-Eijnden,
Proceedings of the ndnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 20255 by the author(s).

---

Distinguishing Cause from Effect with Causal Velocity Models
recover the correct causal direction.'
La Sa (v)) = (c) y(a')
The appealing properties of our method depend on accu- rate estimate of the score function at the data. The nonpara- metric score estimators that we use have known consistency properties and convergence rates (Zhou et al., 2020), which justifies the method asymptotically. However, some distri- butions, including those seen in common benchmarks, may not have favourable unite-sample behaviour even in the low- dimensional regime. We describe and probe failure modes of our method by way of poor score estimation.
fo (v)
(x,
v(, x)) = vz))) = (d/dx') (oddz')ssn fa
Outline In Section 2 we cover the necessary background on bivariate causal discovery and on dynamical systems. Section describes our framework of viewing SCMs as dynamical systems and the velocity parametrization. In Sec- tion 4, we use the continuity equation to derive an identity that we use to recover the causal velocity from the score function. Section 5 applies this to causal discovery and discusses consistency and identifiability. Finally, Section 6 and 7 describe related work and experimental evaluations.
Figure 1. Given a SCM X
Y and an observation (2, y). the causal curve represents the implied counterfactual outcomes had æ been x' (top). The derivative of the causal curve, which we call the causal velocity function z), generates all the causal curves and hence characterizes the counterfactuals of the SCM (bottom). In this paper, we show how the velocity can be estimated directly from data using the score function, without the need to evaluate the mechanism itself nor specify the underlying distributions.
2. Background
2.1. Structural Causal Models
A structural causal model (SCM) for X causing Y is typi- cally defined the tuple Px,)
the advantages of various state-of-the-arl methods in that it is agnostic to the noise distribution and avoids explicit esti- mation of the functional model, even for non-Gaussian and non-additive cases. Beyond this, the velocity parametriza- tion allows us to specify model classes that extend beyond ANM and LSNM in terms of dependence on the noise, but are still restrictive enough to be useful for causal discovery.
Y yxf(X,c))) X11lyy
(1)
with X Px and cy ~ Pa.. Throughout, we will assume that (=()) f(xr, +) is bijective for each ær and denote fo
its inverse. The pair (fz, P.,.) represents the conditional distribution of Y = z, and Px completes the specification of the joint distribution. The representation of a conditional distribution as a function of the conditioning variable and independent noise is always possible (Kallen- berg. 2021, Lem. 4.22), and it can be chosen to be bijective in the noise if, for each æ, the conditional distribution of Y|X is atomless.
We present new perspective on invertible SCMs dynamical systems through their implied velocity functions. In particular, counterfactual prediction becomes the solution to an initial value problem in which the factual observation is the initial condition. This allows us to specify novel causal models. It also allows us to use established ideas from measure transport to connect the causal velocity to the score functions of the data distribution. We use this to devise novel simulation and likelihood-free estimation procedure for bivariate SCMs and related goodness-of-ftt criterion for bivariate causal discovery.
Contributions
Bijectivity enables identification of counterfactuals. In particular, let (x, v(x)) = (I, fe(())) denote a realiza- tion of the SCM. Assuming no hidden variables, the coun- terfactual realization "had xr been xt is computed by the abduction-action prteriittttiiinnnriiiitionnn procedure (Pearl, 2009) as fa f= (v(2))) = (s', fa', = (x', v(x))) Without bijectivity, the conditional distribution of cylar, v(x) must be estimated. A bijective SCM point-identifies the counterfactual directly as a deterministic transformation:
To support the methodological contributions, we also ex- amine the causal direction identifiability problem from the dynamical perspective, where we obtain a general result that specializes those existing in the literature. We also establish statistical consistency of our estimation procedure, and show that the rate of convergence is dictated by the rate at which the score estimators converge. We show through synthetic and benchmark experiments that the added model flexibility is beneficial in situations where ANM and LSNM models can fail, but simple parametric velocity families
g()) - v(xx)=)=))))))) = Sar ST (y(x)).
(2)
et al. (2023) established the identifiability
in various settings of bijective SCMs, with y(x') given by
Codee supporting our experiments can be found Github https///ghhb..mm///jjnn////ull/
2

---

Distinguishing Cause from Effeet with Causal Velocity Models
Table Summary of existing and novel SCMs and their velocities. The velocity perspective allows for specification of novel model classes and can also be used as an alternative means of estimating existing model classes.
Model
SCM
Counterfactual fx.(fz (v))
Velocity
Parameters
f(X,cg)
(d/dr)) (5 (u))
ANM PNL
m(X) tey g(m(X)+ teg)
y m+(((=) ml(r)
riu(z) g(g (y))m(()
i R + R
g(g (y + m(xa)- m(x))
t. g : R + R
LSNM
m(X)= + c)i(x)) cy + Iz T ((((()) uidu
m(x)) m(x' + ch(=e)
m(xx) + (v) m(x)) a (
m. h IR -> R E RK
Basis
y+ l a (v)()) ujdu
Black box
Ey + 20 v(u(u), u) du
y + 12 v(u(u), u)du
v(y,
0: R² + IR
(2). Applying this transformation for all possible x' yields curve of counterfactual outcomes; see Figure 1.
Under regularity conditions, ODEs are known be in cor- respondence with two-parameter flows, defined as follows. Definition 2.1. A two-parameter flow (or flow) Pa, on Rd is a continuous mapping
2.2. Functional Bivariate Causal Discovery
Causal discovery aims to determine, from pairs of observed data (X, Y). whether X or Y is the cause. This problem is underdetermined when using observational data, due to the universality of the SCM representation of conditional distri- butions. One common approach, based on the argument that the data-genn process in the causal direction is simpler (Mooij et al., 2016), is to restrict the model class. It hoped that in doing so, the model will fit the causal direction but exclude the anti-causal direction (Goudet al., 2019).
(s, t, v) + (v,( that satisfies the flow properties
1. (r,(()) id(y) y. for all E R, y E Fd. Paa o PP,, 9s,u. for all 8, 8,1, u E R.
We note that yp,,1 = Pt.s. The ODE (3) is said to generate the flow Pa,z if
(or(y)) ( S. u(pp.u()), du
Because the distribution of the noise variables is not related to causality per se, most attention is paid to restricting the
E (4)
mechanism. A well-established framework that does not re- quire constraining or modelling the noise distribution is the regression and subsequent independence test (RESIT) ap- proach (Peters et al., 2014), which quantifies the goodness- of-fit of causal model by the independence of the cause variable to the residuals obtained from the model fit. Origi- nally proposed for nonparametric regression in ANMs. the RESIT approach has also been extended to LSNMs (Immer et al., 2023; Strobl & Lasko, 2023), PNL models by post- processing model estimates with non-lineat ICA (Zhang & Hyvärinen, 2009), and rank-based approaches (Keropyan et al., 2023). A key takeaway from this literature is that we should aim for methods that allow for though not unboundedly so-model classes that can be estimated without specifying the noise distribution.
If ya,s is differentiable in t and satisfies
d
dt Px,(()) = (U))).) and (x, = y
(5)
then Ss,t is said to be (classical) solution of (3). Con- versely, a flow can be used to define an ODE. If + (v) is differentiable in f. at I = for all y. then
d
u(y, 8) = dt G (v)
(6)
yields an ODE to which ss,, is a solution, and v is known as the velocity. Under regularity conditions (Appendix A.1), velocities and flows are in one-to-one correspondence.
Now suppose that an initial condition for the ODE (3) is chosen randomly as (s) ~ Pss where P, is the density of probability distribution on Rdd Letting Y(s) evolve according to the ODE yields (t)) Gx., (Y View- ing the flow measure transport, the density of Y()) is Pt (Pa,))=P. It can be shown that the family of densities (pt) solves the PDE (known as the continuity equation)
2.3. Flows and Differential Equations
The subsequent sections will expose connections between SCMs and dynamical systems generated by differential equations, so we briefly review some relevant results about the latter here. Our exposition largely follows Arnold (1998, Appendix B) and Santambrogio (2015, Ch. 4.1.2). Consider the ordinary differential equation (ODE) in Rd.
Oep: + V, v (pree) = 0,
(7)
with = v(.,t) and 0 = 0/Ot. It turns out that the solu- tion is unique up to (Lebesgue) almost everywhere equiv- alence. (See Appendix A.I for a formal statement.) Note that if P. is strictly positive then (7) is equivalent to
dy(l) = -(v(),, v(1) E Rd, LER. dt
0, log
-t - U2 V log.:
(8)
3

---

Distinguishing Cause from Effeet with Causal Velocity Models
3. SCMs as Dynamical Systems
point wo. which plays the role of the noise distribution as follows
Generally, flows describe the time evolution of dynamical systems: the usual statement of Definition 2.1 yields the evolution of y between time points, and causality is associ- ated with the forward flow of time. Here, we take different flow
(x) = (Cy): by ~ p(y)zo)
Under regularity conditions, the SCM is unique and holds over all of I Technically, an SCM in the sense of Section also requires Px to fully specify the joint distribution but it plays no role in the causal mechanism or counterfactuals. The following theorem formalizes the relationship between bijective SCMs and dynamical systems. See Appendix A. for technical details on the regularity conditions.
perspective, viewing counterfactual transformation in which the role of time is played by cause variable xr Whilst a similar connection has been made in previous work (Dance & Bloem- Reddy, 2024), by specializing to the bi- variate real valued case and viewing the cause explicitly as time, we are able to make a direct connection to ODEs. This enables us to derive novel velocity- based parametrization of SCMs and simulation- free learning procedure. For the remainder of the paper, we work in the setting of X E IR and Y E R.
Theorem 3.2. Let X cause Y, with X,Y E R. Under regularity conditions, a bijective SCM uniquely determines a velocity-densigg pair (v. )o)) Conversely, (u. (()))) determines a SCM uniquely up to changes in Px-
Definition 3.1. Let X cause Y. with bijective SCM (I). The flow generated by the SCM, or SCM flow, is defined as
The above equivalence shows that we can view bijective SCMs in terms of their underlying velocity functions with- out loss of generality. See Table 1 for the velocity functions associated with familiar classes of SCMs.
(z, H for (f (y)) Fx, (v)
is easy to see that Equation (9) satisfies the axioms of
flow Definition 2. In particular, for an observation (ar, y),
3.1. Velocity Parametrization of SCMs
the causal curve (Figure 1) is the function of x' defined by
The dynamical perspective allows causal model to be specified by the velocity and its implied counterfactuals instead of in terms of first (ANM) or second (LSNM) condi- tional moments, which can be more interpretable. Viewing p(y zo) as sampling individual baseline measurements from some population, we see that ANMs indicate that the relative effect of an intervention is the same for all individu- als, regardless of their measurement value y. This manifests graphically parallel curves (Figure 2). LSNMs relax this, but implies that the difference between individual trajecto- ries are proportional. Distributionally, ANMs and LSNMs imply that x controls at most the second moment of p(y x). In particular, if p(y xo) were Gaussian for any zo. then it is also Gaussian for all x. On the other hand, the velocity func- tion ex implics a more complicated mechanism that can model insensitivity in the tails of y. no matter the value of x. Distributionally, this also modifies higher order moments of the conditional distribution; see the rightmost panel of Fig- ure 2n practice, the velocity model can be specified in an interpretable way via domain knowledge of the underlying process, or with basis functions or neural networks for more flexible models. The velocity parametrization automatically specify bijective SCMs, and as such their counterfactuals are identifiable (Nasr-Esfahan et al., 2023).
x' + y(z') Par. (3)
(10)
Causal curves represent the counterfactual outcomes under the SCM from factual observation (I: y). Generally, each observation specifies a different causal curve. Based on the dynamical systems perspective reviewed in Section 2.3, (x:, y) can be viewed as initial conditions.
In analogy to (6), if (or,a'.)) + (a.zr (y) is continuous and x + fz (sy) is differentiable for all ey then the causal velocity of the SCM flow is (see Figure 1, bottom)
u(y, x) = di for (t
(11)
This holds for all (x, in the support of Px.y. The resulting v describes the local behaviour of the flow at (xr, v): x() is the effect of an infinitesimal intervention x' ='=r+t + ö as 8- + O, which, based on the algebraic proper- ties of the flow, generates the counterfactual curves. Thus, we see that the mechanism of bijective SCM uniquely determines a flow and its associated velocity.
Conversely, the causal mechanism of bijective SCM can be parametrized by a velocity via its corresponding flow from an arbitrary To R acting the noise state,
4. Score Functions and SCM Flows
a (v) = y I" v(Paro. (V)) w) du.
(12)
In the previous sections, we established an equivalence be- tween SCMs and an associated pair (v. p(ylz))) Here, we show that when a joint distribution is generated by an SCM, the velocity leaves an explicit signature on the derivative of
To complete the specification of the SCM, it remains to spec- ify the conditional distribution p(y | mo) at some arbitrary
4

---

Distinguishing Cause from Effeet with Causal Velocity Models
ANM
LSNM
v(y,z) = e²²
Figure Distributions of causal curves obtained by sampling 100 initial conditions from ply | == ro) (o, 1). Note ANM and LSNM cross-sections ply æ)) are necessarily Gaussian also, no matter how flexible the underlying function classes are. On the other hand, the exponential velocity yields non- con-Gaussian conditionals at any x o 0. exhibiting right and left skew as æ < and æ respectively.
the log-density, i.e., the score function. We leverage this observation to derive a goodness-of of- criterion based on the score, which can then be used for model-fitting and checking directly at the level of the causal velocity. Re- markably, this allows to check whether data can possibly be generated from an SCM with causal velocity without ever having evaluate the mechanism, making any assumptions about the underlying distributions beyond dif- ferentiable score functions.
'True Score
True Mechanism
Let v be the velocity of a SCM that generates the conditional density p(y]æ)), and assume that the joint distribution has full support, with differentiable marginal and joint log-densities. Let (y) == O₂ log (z))2)) and (y(zz) = O, log p(y (z) denote the partial derivatives of the log conditional density, and similarly 8 (2, y) and 85 (xr) the partial derivatives of the log joint and marginal densities, respectively. Viewing the cause variable x as time, the log version of the continuity equation (8) yields
Figure 3. If the model is well-specified and the true score is given. velacity functions as obtained by minimizing Equation (15) (up Monte Carlo approximation error) can be integrated numerically (left) recover the true mechanism (right).
sx 2 r) -0, x) vyy. ())((()))))
(13)
Observe that (14) is entirely in terms of observable scores, and thus is agnostic to the unobserved noise distribution Pey. Indeed, SCMs with the same mechanism but different noise distributions will still satisfy (14) with their respec- tive scores. This makes it an especially suitable objective for functional causal discovery, which asks whether distri- butions can be generated by a restricted mechanism class, while leaving the noise distributions unspecified.
By known results (see Appendix A. every solution p(yl arises from some initial condition density p(y) trans- ported by the SCM flow, for fixed but arbitrary zo E TR. Moreover, that solution is unique. Therefore, (13) can be used to characterize when a conditional density p(ylz) could possibly have been generated by a SCM with velocity v. For practical purposes, (13) can also be stated in terms of marginal and joint scores, which can be estimated from data.
5. Velocity-Ban -Based Causal Discovery
Theorem 4.1. Let X cause Y, with X, Y R. Assume that Px, has full support, with differentiable joint, conditional, and marginal log densities. Then
For methods that resolve the direction of causality by fitting model to the data in both candidate directions, some re- striction to model complexity must be made. As discussed in Section 2.2, many existing methods assume some variant of ANM in order to avoid making specific distributional as- sumptions about the unobserved noise. Our velocity-
(x -0,u(y, v(y, x: (g( (x(x) (14)
for all (x,y) E ? if and only if Рx,Y can be represented by a SCM with velocity v and marginal density p(z:).
5

---

Distinguishing Cause from Effeet with Causal Velocity Models
method, described in this section, also makes no assump- tions on the noise distribution by capturing its influence via score estimation, and allows all modeling and complexity control to be imposed via the parametrization of the velocity. As established in previous sections, this allows the model to be substantially more flexible than ANMs or LSNMs.
Causal
Anticausal
GoF (True Score)
GoF (Estimated Score)
100
9
0.5
Empirical GoF
10-1
0.3-
We define a goodness-of-fi (GoF) statistic of a velocity
10-2
function v as,
8
8
0.1-
10-1
L(u)= == (X) - (x, s. (X, r))
(15)
100
500
1000
2500
5000
10000
100
500
1000
2500
5000
10000
Sample Size
which is derived from (14) and uses the notation
Figure 4. If the ground truth score known, the anticausal GoF orders of magnitude larger than the causal GoF. In this case, cause can he distinguished from effect with certainty even with as few as 100 samples. Right: As the score estimate quality improves with larger sample sizes, the causal GoF decreases (Theorem 5.1), and the gap between the causal and anticausal GoF widens.
$,(2, y) == y) + (y,,r)s (2,
(16)
This notation reflects the fact that s, the directional deriva- tive of p(a, v) along the causal curve, i.e., in the direc- tion O₂ (x,, (1, "(v(2)), ):) (see Appendix A.2).
We propose to minimize the GoF directly (15), which is zero if and only if (14) holds Pa.y-almost everywhere, to estimate the velocity given the scores. See Figure 3 for an example; when the true score is given, minimizing (15) recovers the true mechanism accurately. In practice, the score function is unknown and two-step approach is required. First, we estimate the score nonparametrically (Zhou et al., 2020); second, we estimate the velocity by minimizing (15), (17) where the scores are replaced by estimators and the expectation is estimated empirically (18).
and the model is well-specified, (17) is orders of magni- tude larger than (15), which identifies the causal direction with certainty regardless of the sample size. In practice, the causal discovery performance improves as the score estima- tion improves in sample size (see Section 7.3 for a more in depth evaluation). In the following section, we show that the empirical estimator of (15) converges rate controlled by the rate which the score estimator converges.
5.1. Consistency of Empirical GoF Estimation
For causal discovery, we also estimate the model the Y X direction by minimizing
Let (aaa, ya y:) denote the available data. The two-step pro- cedure described above yields the following objective fune- tion for estimating the velocity
( = (Y) 0x6(X, - so (X,x))) (17)
In principle, an estimated velocity could be integrated to a base point in order to estimate the noise variables, i.e., & = xx.. (Y), which could then be used as residuals in independence tests for causal discovery (as in RESIT). However, for the remainder of this paper, we will use the more direct approach of using the value of the objective itself to determine causal direction. The resulting two-step method is as follows:
(v)(v) = Σ ((s((1))) - Oyy vs) - 8 (0 (1))) 1m
(18)
Recall that $, v.) o(ys, 22), 91) + S y.). and here šn, by are the score estimators, for which we use either the Stein estimator (Li & Turner, 2018) or one based on kernel density estimator (Wibisono et al., 2024). These estimators can both be viewed as doing regularized, vector- valued reproducing kernel Hilbert space (RKHS) regression on the true scores (Zhou al., 2020). Based on this, the fol- lowing result shows that our estimator for the loss function converges at a rate no worse than the convergence rate of the score function estimators. In the following theorem we denote $..y := sy (2, y). s2.9 == Sy (z. y) for convenience. and similarly for the corresponding estimators.
1. Estimate joint and marginal scores from data.
2. Estimate causal velocity in both directions by mini- mizing (15), (17), with expectation taken over the data. Choose as causal whichever direction has a lower value of objectives (15), (17) evaluated on the data.
As might be expected, this method requires accurate esti- mates of the scores. However, since (14) expresses a fixed relationship between the scores and the velocity, an inaccu- rate score estimate will produce an inaccurate velocity that nonetheless could still yield a small value of (15), thereby potentially introducing error into the subsequent causal dis- covery. Figure+ 4 shows that when the true score is given
(v)
Theorem 5.1. Let 8 E Hx and s SEn E Hxx where H.x, Hx are RKHSs induced by bounded kernels be < B,k = k2 @ bi < B. Additionally, let E CR (R2,R) R). If the score estimators have bounded RKHS norm that con- (a) |1 80. l
verges, such that each lf - sxl1pxx
6

---

Distinguishing Cause from Effeet with Causal Velocity Models
22.9 Hxx is Op(n for some a > 0, then [v) L(v)= = Op(n" }). 3 max(a, 23.
|
Expressions for the relevant partial derivatives of in terms of the forward model parameters 4p u are given in Ap- pendix C. We show in Appendix C.1 how the known ANM non non-identifiability differential equation follows naturally from this characterization. In Appendix C.2, we use it to derive a criterion for non-identifiabiiity that holds uniformly over the class of LSNMs, albeit at a loss of interpretability. Thus, we expect that general analysis of identifiability- one is possible may require new techniques.
Zhou al. (2020) showed that under smoothness assump- tions on each of log y), kzz ky, then convergence at rate E [3, 4] can be achieved. See Figure 4 (right) for an empirical demonstration of the behaviour of L, (v) in well-specified model, and Section 7.3 for the implications on causal discovery performance.
5.2. Identifiability
6. Related Work
In the limit of infinite data, the causal direction of an ANM is known to be identifiable under certain conditions, formu- lated as the non-existence of solution to a differential equa- tion for log p(z) (Hoyer et al., 2008: Zhang& & Hyvärinen, 2009). For fixed ANM parameters m, Pe, the set of marginal densities that satisfy the differential equation (and hence are non-identifiable) is three-dimensional space contained in an nntinite-iimensionall space, and therefore it is believed that ANMs can identify causal direction in "most cases" (Peters et al., 2014). Identifiability in ANMs is the strongest known result, as the characterization does not depend simul- taneously on the parameters of both model directions, and therefore holds uniformly over the model space.
Score Matching ANMs A recent line of work also use sig- natures in the score function for causal discovery ANMs. (Rolland et al., 2022; Montagna et al., 2023b) use properties of Gaussian noise to derive conditions on the score that are satisfied by non-cause variables. (Montagna et al., 2023a) show that even under arbitrary distributions the score is de- terministic function of the noise variables, and fit the model using nonparametric regression to evaluate this condition. These methods are restricted to the ANM case, but are fo- cused on sink node (non-cause) identification in multivariate SCMs, while our focus is on more general model classes in the bivariate setting. The continuity equation (13) simplifies for ANMs, and we show how to recover the identities used in this line of work in Appendix A.3.
As argued by Tagasovska et al. (2020); Xu et al. (2022); Strobl & Lasko (2023); Immer et al. (2023) and others, ad- ditive noise can be an overly restrictive assumption: even basic heteroscedastic noise (i.e., a LSNM) can lead ANM- based causal discovery procedures to mis-identify the cor- rect causal direction. While LSNMs extend ANMs, even that straightforward generalization significantly complicates the identifiability analysis. In analogy to ANM identifiabll- ity, Strobl & Lasko (2023); Immer et al. (2023) derived par- tial differential equation (PDE) whose solutions characterize non-identifiability. In contrast to ANMs, the LSNM PDE depends simultaneously on the parameters of both model di- rections. Therefore, is useful on a case-by-case basis, but does not yield a characterization of non-identifiabilityy uni- formly over the model class. As we show in Appendix C. similar characterization for velocity models can be obtained by analyzing the continuity equation (14) in both directions.
Cocycles in Causal Modeling Connections between coun- terfactual causal models and dynamical systems have been previously established by Dance & Bloem Rrodm- (2024) in more general setting, using class of maps called cocycles. A woo-parameter flow is example of cocycle. Dance & Bloem-Redds (2024) focused on distribution-rooust in- ference with a known causal ordering with multiple cause variables, rather than bivariate discovery. Furthermore, they do not analyze the cocycle as dynamical system nor make any connection to the score function.
Continuous Normalizing Flows A popular class of velocity-based probabilistic model are continuous normal- izing flows (CNF) (Chen et al., 2018). There, instead of conditioning, time is an auxiliary variable, and the genera- tive model is for a marginal density (30), where yo is drawn from an arbitrary base distribution. Our learning ob- jective (18), which targets conditional distributions instead, is similar to recently proposed simulation-free objectives, which attempt target the velocity directly in CNFs, e.g., Flow Matching (Lipman et al., 2023). Normalizing flows have been used for causal inference (Khemakhem et al., 2021; Javaloy et al., 2024) with a known causal graph. The causal augo-regressive flow model (Khemakhem et al., 2021) when restricted to the bivariate case represents a flexible LSNM, and was used for likelihood-based causal discovery. Tu et al. (2022) propose to analyze the CNF velocity as an alternative decision rule for causal discovery in ANMs.
Proposition 5.2. Assume that the joint distribution of the observed data can be expressed as causal velocity model in both directions,
Y 2 ( X
x 2 ((x), Yilez,
with corresponding velocities v and i, respectively. Then, as long as π ( y) = log p(a, y) and the required derivatives exist, the velocities satisfy on the support of p(a, y).
03- r + 8y(v- x) = -8.0 + 2z(0 x).
7

---

Distinguishing Cause from Effect with Causal Velocity Models
Table 2. Accuracy (and AUDRC) of velocity models for deter- mining the correct causal direction, evaluated on a collection of synthetic datasets (columns). Score is estimated using the Stein method with the Gaussian kernel. The results shown for LOCI (NN/Spline), IGCI, and CGCI represent the best performance over different variants. The best performing velocity model, and, if ap- plicable, the overall best model, are indicated in bold. KDE results and a subsampled a = 1000 setting can be found in Appendix E.
Score Estimation MSE
Causal Discovery
1.0
0.40
AJDRC
AUDRC/SucceS, %
Success
Score MSE
0.30
0.8
0.20
0.10
0.6
Model
Velocity
Sigmoid
ANM
LSNM
0.01
100
500
1000
2500
5000
10000
100
500
1000
2500
5000
10000
B-LIN
91 (97)
(95)
66 (70)
(77)
B-QUAD
97 (98)
80 (94)
63 (70)
73 (79)
Sample Size
V-ANM
86 (97)
36 (41)
92 (95)
59 (67)
V-1.SMM
89 (99)
81 (95)
86 (92)
82 (92)
Figure 5. Score estimation and causal discovery performance as sample size increases in well-specified ANM benchmark. Upper and lower limits indicate 1st-3rd quartile over 100 datasets. The larger error at 10000 is likely due to an increased chance of encountering outliers; note the median continues to decrease.
V-NN
96 (99)
56 (69)
65 (72)
57 (64)
LOCI (HSIC)
40 (25)
70 (86)
97 (99)
83 (89)
LOCT (I.ik)
43 (66)
43 (63)
49 (56)
40 (31)
CDS
56 (47)
12 (5)
92 (99)
59 (68)
IGCI
66 (77)
73 (83)
29 (21)
37 (31)
which we found to result in the best performance.
RECI
33 (26)
13 (4)
34 (36)
38 (55)
CGCI
52
41 (33)
85 (96)
(89)
Computation Experiments are written in JAX (Bradbury et al., 2018) and carried out either on a M1 Mac or NVIDIA V100 GPU. The most costly operation is the Stein score estimation, which is O(z²) due to the inversion of an n x matrix (Li & Turner, 2018). Note this only has to be performed once, prior to training the velocity model. With n = 1000 and n = 5000, the median time is 0.02s and 0.3s respectively on GPU, and 0.9s and 15s on MI.
7. Experiments
We evaluate our method empirically on new synthetic datasets as well as on existing benchmarks from the lit- erature (Mooij et al., 2016; Immer et al., 2023). We study the relative performance of different velocity model classes, including velocity- based parametrizations of known classes. In Section 7.3, we show empirically that the performance of our method hinges on the accuracy of the score estimators. When the true score is given, our method achieves perfect causal discovery with as few as n = 100 samples (Figure 4).
Evaluation Causal direction is determined by comparing the empirical loss (18) in both directions and selecting the smaller value. Either the squared or absolute values can be used; here we use the squared value. We follow the literature and report the raw causal discovery accuracy as well as the associated area under the decision rate curve (AUDRC).
Velocity Estimation We use three novel parameterizations of the velocity: a 2-layer MLP (V-NN), and linear/quadratic basis families (B-LIN/P CEIIN/--QUAD We also experiment by fitting ANM and LSNMs via their velocity parameteriza- tions (V-ANM and V-LSNM). See Table 1 and Appendix D for details. Given an estimate of the score, we minimize (18) using the Adam optimizer (Kingma & Ba, 2015) to estimate the velocity. The term 0, x²) is obtained by au- tomatic differentiation. The datasets we consider are small (n < 10²), and thus take full batch gradient steps. To mitigate potential non-identifiability and regularize estima- tion, we penalize the complexity of the model via higher order derivative terms d² y(x) dx (Kelly et al., 2020) in the case where both directions are specified to encourage selecting the simpler model. In practice we use k: = 2.
Other Methods We compare our method to LOCI, an LSNM-based model which represents the state-of-the- thef-tth-artttt in functional causal discovery based on the results reported by Immer et al. (2023). It has two major variants, one based on the Gaussian conditional likelihood (Lik) which makes strong parametric assumptions on the conditional distribu- tion, and one based on the RESIT approach (HSIC), which is robust to the noise distribution and thus more comparable to our method. We also compare our method to con-funciional methods that do not make parametric assumptions such as CDCI (Duong & Nguyen, 2022) and other methods that are implemented in the Causal Discovery Toolbox (Kalainathan et al., 2020), which include CDS (Fonollosa, 2019), IGCI (Daniušis et al., 2010), and RECI (Blöbaum et al., 2018).
Score Estimation For estimation of the score, we use the Stein gradient estimator of Li & Turner (2018) (STEIN), with the Gaussian kernel and bandwidth selected by the median heuristic. We also consider the derivative of the log of a standard kernel density estimator with empirical Bayes smoothing (Wibisono et al., 2024) (XDE). For the KDE- based estimator, we use the exponential (Laplace) kernel,
7.1. Datasets
Synthetic We design two synthetic settings for which ex- isting functional methods (e.g., LOCI) are misspecified and are expected to fail. All benchmarks consist of 100 datasets of size n = 5000 with randomly sampled mechanisms. The
8

---

Distinguishing Cause from Effect with Causal Velocity Models
Table Accuracy (and AUDRC) of velocity models on benchmark data. LOCI (NN/Spline). IGCI, and CGCI represent best results over all variants. The best performing velocity model, and, if applicable, the overall best model, are indicated in bold. Additional comparisons on these benchmarks can be found in Mooij et al. (2016); Immer et al. (2023).
Model
SIM
SIM-C
SIM-G
Tü
Tü (cts.)
KDE
STEIN
KDE
STEIN
KDE
STEIN
KDE
STEIN
KDE
STEIN
B-LIN
40 (39)
58 (62)
56 (50)
55 (54)
83 (93)
71 (85)
71 (76)
54 (59)
78 (81)
55 (68)
B-QUAD
39 (37)
60 (64)
43 (37)
62 (54)
94 (99)
71 (82)
52 (65)
56 (58)
58 (68)
61 (68)
V-ANM
36 (39)
63 (77)
52 (49)
76 (71)
58 (77)
77 (93)
68 (78)
60 (64)
74 (82)
69 (72)
V-LSNM
39 (39)
62 (67)
52 (49)
76 (71)
74 (92)
72 (88)
58 (70)
62 (65)
57 (73)
68 (73)
V-NN
37 (38)
57 (69)
48 (44)
60 (54)
84 (95)
68 (78)
58 (70)
59 (65)
65 (72)
64 (68)
LOCI (HSIC)
79 (89)
83 (93)
81 (91)
60 (56)
57 (58)
LOCI (Lik)
52 (68)
50 (63)
78 (89)
57 (66)
60 (66)
CDS
68 (85)
78 (89)
73 (82)
66 (62)
60 (55)
IGCI
37 (33)
(35)
86 (96)
60 (71)
53 (64)
RECI
45 (48)
(58)
42 (41)
72 (88)
71 (88)
CGCI
68 (84)
76 (91)
75 (88)
67 (64)
61 (57)
Velocity benchmark is generated by numerically integrating periodic velocity functions with initial values given by the noise variable, as in (12). The Sigmoid benchmark can be seen as a variation on LSNMs with additional post-nonlinear and affine transformations.
well-estimated, and is competitive with other methods over- all (Table 3). Notably, taking only the continuous subset of the Tübingen dataset improves performance in our method but not existing methods, which emphasizes the importance of the log-density assumption. Interestingly, using the KDE estimate improves causal discovery performance in cases where the noise distributions are Gaussian (SIM-G), but also the real data benchmark (Tü).
In all cases, noise variables are sampled as a randomly sampled invertible transformation (i.e., a random CDF trans- form) of Gaussian noise. As such, methods based Gaus- sianity are also expected to fail. To illustrate our method as an alternative to HSIC in existing functional classes, we also design synthetic ANM and LSNM benchmarks for which V-ANM and V-ISNM are well-specified. Sample plots of datasets are given in Appendix D.1.
7.3. Sample Size and Score Estimation
To study the effect of score estimation on velocity estima- tion and downstream causal discovery, also designed synthetic ANM and LSNM datasets with Gaussian noise, for which the true marginal cause and joint scores can be obtained analytically. For the effect variable, we compute the marginal density by Monte Carlo integration of the con- ditional and differentiate to obtain the marginal score.
Following previous work (Blöbaum et al.. 2018; Tagasovska al., 2020; Immer et al., 2023), we also evaluate our method on the SIM-series of simulated benchmarks and the Tübingen Cause-Efféee pairs of Mooij et al. (2016). Since our method requires the existence of log- densities, we also test removing instances from the Tübingen collection with inteeer-valued observations, which we refer to as Tü (cts.) (details in Appendix D.4).
Benchmark Data
In Figure 4, we saw that using the ground truth scores as input to estimate the velocity determined the causal di- rection with certainty, regardless of the sample size used in minimizing (18). On the other hand, when the score is esti- mated, causal discovery performance appears to improve as sample sizes grows. Here, we provide additional evidence that the improvement causal discovery performance is due to an improving score estimate. We use the ground truth scores to evaluate the quality of our estimated scores by eval- uating the mean squared estimation error. Figure 5 shows that the score estimation largely improves in sample size, and this corresponds to an improvement in both measures of causal discovery performance. Note the figure shows estima- tion of the marginal score of the effect variable-joint and marginal cause score estimation follows similar patterns. See Appendix E.3 for full tables of results. We believe improvements to score estimation will directly benefit the applicability of our methodology in future practice.
7.2. Results
The experiments show that our two-step method, in par- ticular with B-LIN and B-QUAD velocity models which represents novel model classes, are able distinguish cause from effect when existing methods fail (Table 2). For non- Gaussian ANM and LSNM data, velocity parametrizations achieve performance that is competitive with the HSIC vari- ant of LOCI. As expected, the (Gaussian) likelihood vari- ant of LOCI performs poorly in our non-Gaussian simu- lations. Our method also achieves state-of-the-art perfor- mance certain benchmarks, in particular the Gaussian SIM-G dataset which we suspect is due to the score being

---

Distinguishing Cause from Effeec with Causal Velocity Models
Acknowledgments
Goudet, O., Kalainathan, D., Sebag, M., and Guyon, I. Learning Bivariate Functional Causal Models, pp. 101- 153, Springer International Publishing, 2019.
The authors are grateful to the anonymous reviewers, whose comments helped improve the paper. This research was supported in part through computational resources and ser- vices provided by Advanced Research Computing at the University of British Columbia. JX is supported by Nat- ural Sciences and Engineering Research Council of Canada (NSERC) Canada Graduate Scholarship. BBR acknowl- edges the support of NSERC: RGPIN2020-0499 RGPAS- 2020-00095. HD and PO are supported by the Gatsby Char- itable Foundation.
Gretton, A., Bousqeet, O., Smola, A., and Schölkopf, B, Measuring statistical dependence with Hilbert-Schmie norms. In ALT, 2005.
Hodgkinson, L., van der Heide, C., Roosta, F,, and Mahoney, M. W. Stochastic continuous normalizing flows: training SDEs ODEs. In UAI, 2021.
Hoyer, P., Janzing, D., Mooij, J. M.. Peters, J., and Schölkopf, B. Nonlinear causal discovery with additive noise models. In NeurIPS, 2008
Impact Statement
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.
Immer, A., Schultheiss, C., Vogt, E., Schölkopf, B., Bühlmann, P., and Marx, A. On the identifiability and es- timation of causal location-scale noise models. In ICMI., 2023.
Javaloy, A., Sánchez-Martin, P., and Valera, I. Causal nor- malizing flows: from theory to practice. In NeurIPS, 2024.
References
Albergo, M. and Vanden Vanden-Eijnden, E. Building normalizing flows with stochastic interpolants. In ICLR, 2023.
Kalainathan, D., Goudet, O., and Dutta, R. Causal Discov- ery Toolbox: uncovering causal relationships in Python. Journal of Machine Learning Research, 21(37), 2020.
Aroold, L. Random Dynamical Systems. Springer Berlin Heidelberg, 1998.
Bartlett, P. L. and Mendelson, S. Rademacher and Gaussian complexities: Risk bounds and structural results. Journal of Machine Learning Research, 3(11), 2002.
Kallenberg, O.
Foundations of Modern Probability. Springer International Publishing, 2021.
Kelly, J., Bettencourt, J., Johnson, M. J., and Duvenaud, D. K. Learning differential equations that are easy to solve. In NeurlPS, 2020.
Blöbaum, P., Janzing, D., Washio, T., Shimizu, S., and Schölkopf, B. Cause-effect inference by comparing re- gression errors. In AISTATS, 2018.
Keropyan, G., Strieder, D., and Drton, M. Rank-based causal discovery for post-nonlinear models. In AISTATS, 2023.
Bradbury, J., Frostig, R., Hawkins, P., Johnson, M. J., Leary, C., Maclaurin, D., Necula, G., Paszke, A., VanderPlas, J., Wanderman-Milne, S., and Zhang, Q. JAX: composable transformations of Python+NumPy programs, 2018. URL http: /ttp://githbb com/ jax-ma/jaxa.....
Khemakhem, I., Monti, R., Leech, R., and Hyvarinen, A. Causal autoregressive flows. In AISTATS, 2021.
Kingma, D. and Ba, Adam: A method for stochastic optimization. ICLR, 2015.
Chen, R. T. Rubanova Y., Bettencourt, J., and Duvenaud, D. K. Neural ordinary differential equations. In NeurlPS, 2018.
Li, Y. and Turner, R. E. Gradient estimators for implicit models. In ICIR, 2018.
Dance, H. and Bloem-Reddy, B. Causal inference with cocycles. arXiv preprini arXiv:2405.1 13844, 2024.
Lipman, Y., Chen, R. T. Q., Ben-Hamu, H., Nickel, M., and Le, M. Flow matching for generative modeling. In ICLR, 2023.
Daniušis, P., Janzing, D., Mooij, J., Zscheischler, J., Steudel, B., Zhang, K., and Schölkopf, B. Inferring deterministic causal relations. In UAI, 2010.
Montagna, F., Noceti, N., Rosasco, I.., Zhang, K., and Lo- catello, F. Causal discovery with score matching on addi- tive models with arbitrary noise. In CLeaR, 2023a.
Duong, B. and Nguyen, T. Bivariate causal discovery via conditional divergence. In CLeaR, 2022.
Fonollosa, J. A. Conditional distribution variability mea- sures for causality detection. In Cause Effect Pairs in Machine Learning. Springer, 2019.
Montagna, F., Noceti, N., Rosasco, L.,, Zhang, K., and Lo- catello, F. Scalable causal discovery with score matching. In CLeaR, 2023b.
10

---

Distinguishing Cause from Effeet with Causal Velocity Models
Mooij, J., Janzing, D., Peters, J., and Schölkopf, B. Re- gression by dependence minimization and its application to causal inference in additive noise models. In ICML, 2009.
Zhang, K. and Hyvärinen, A. On the identifiability of the post-nonlinear causal model. In UAI, 2009.
Zhou, Y., Shi, J., and Zhu, J. Nonparametric score estima- tors. In ICML., 2020.
Mooij. J. M., Peters, J., Janzing, D.. Zscheischler, J., and Schölkopf, B. Distinguishing cause from effect using observational data: methods and benchmarks. Journal of Machine Learning Research, 17(32), 2016.
Nasr-Esfahany, A., Alizadeh, M., and Shah, D. Counterfac- tual identifiability of bijective causal models. ICML, 2023.
Pearl, J. Causality. Cambridge University Press, 2009.
Peters, and Bühlmann, P. Identifiability of Gaussian structural equation models with equal error variances. Biometrika, 101(1), 2014.
Peters, J., Mooij, J. M., Janzing, D., and Schölkopf, B. Causal discovery with continuous additive noise models. Journal of Machine Learning Research, 15(58), 2014.
Rolland, P., Cevher, V., Kleindessner. M.. Russell, C., Janz- ing, D., Schölkopf, B., and Locatello, F. Score matching enables causal discovery of nonlinear additive noise mod- els. ICMI,, 2022.
Santambrogio, F. Optimal Transport for Applied Mathe- maticians. Birkhääserr Cham, 2015.
Schultheiss, C. and Bühhlann, P. On the pitfalls of Gaussian likelihood scoring for causal discovery. Journal of Causal Inference, 11(1), 2023.
Silverman, W. Density estimation for statistics and data analysis. Routledge, 2018.
Steinwart, I. Support Vector Machines. Springer, 2008.
Strobl, V. and Lasko, T. A. Identifying sptient-specific root causes with the heteroscedastic noise model. Journal of Computational Science, 72, 2023.
Tagasovska, N., Chavez-Demoulin, V., and Vatter, T. Dis- tinguishing cause from effect using quantiles: Bivariate quantile causal discovery. In ICML, 2020.
Tu, R., Zhang, K., Kjellström, H., and Zhang, C. Optimal transport for causal discovery. In ICLR, 2022.
Wibisono, A., Wu, Y., and Yang, K. Y. Optimal score estimation via empirical Bayes smoothing. In COLT, 2024.
Xu, S., Mian, O. A., Marx, A., and Vreeken, J. Inferring cause and effect in the presence heteroscedastic noise. In ICML., 2022.

---

Distinguishing Cause from Effeet with Causal Velocity Models
A. Additional Details
A.1. Background on Flows
The key results we use from dynamical systems are as the following two theorems. We refer the reader to Arnold (1998); Santambrogio (2015), and references therein, for more details.
Theorem A.1 (Amold, 1998, Thm. B.3.1, B.3.5). Suppose that (t. v(y, 1) is locally Lipschitz continuous in y, and satisfies the local linear growth condition,
v(yyt))= & att)lll + b(t)
where and b are non-negative locally integrable functions. Then the maximal solution to (3) generates unique solution flow Pa.t in Equations (4) and (5). Ifv(+,() is k-times continuously differentiable in for each i, then is (P,,( ..). Conversely, if a flow (y) is differentiable in at ott= = for every then the ODE defined by (6) both generates and is solved by Ps,t-
Flows generated by velocities also characterize the possible solutions to the continuity equation.
Theorem A.2 (Santambrogio, 2015, Thm. 4.4). Fix u and assume the conditions of Theorem A.I. Then for any Pss the family of densities pt (Ps,t)-ps uniquely solves the continuity equation (7) with initial condition density ps- Moreover, every solution (Pr)LER to for fixed v is obtained from the corresponding flow of some initial condition density Pa-
A.2. Density Along Causal Curve
The continuity equation in (14) can be viewed an Eulerian perspective on the problem, tracking velocity and density points (z. v). Alternatvely, may consider the Lagrangian perspective, viewing the problem along trajectories, which in this case are causal curves. The result is the same, but we include the analysis here to shed additional light on the problem.
Let y( be a causal curve with velocity v. We evaluate the log joint density along the curve (x, (())):
log ((),((())) log p(x)
(19)
@e (v(x)) 0y(xx)
llgg(y((o))) (ro) + log
where wo is an arbitrary origin point. This follows from the usual change-of-variables formula that results from transporting along the causal curve from (z, (z)) to (Io, v(zo)). Compared the usual SCM construction, log p(y(zo) ) xо) the
log-density of the noise, and (Pz,zo represents the residual map. Notice this "noise" term does not depend on æ. This is because the noise realization is the same no matter where the counterfactual prediction is made.
Now, we take (total) derivative in and use the instantaneous change-of-variables formula commonly used in the neural ODEs literature (Chen et al., 2018; Hodgkinson et al., 2021) to obtain
d dloog p(x, v())) dx
d logp(r) dx
dn (v(((t),
dy(x)
= (2(()) @v(u(x),z)
(20)
Oy(xx)
Notice the first term is precisely the marginal score sz (x). We take the same total derivative using the chain rule
d do pl,, (x))
dx a log p(z, v()) dx =s, y(x)) + ())(),, vz))
dy(x) Olg p(a, v(z))
dx
Oy(I)
(21)
In fact, this is a directional derivative of the joint log-density along the causal curve, i.e., in the direction O, (z, y(x))
(1, x)), Hence, denoting
(x, = sx (x, vx)) (()())) (x,xz)))
(22)
12

---

Distinguishing Cause from Effeet with Causal Velocity Models
and equating (20) and (21), obtain an expression equivalent (14),
z (z) u = v(x).
(23)
This is relation between the marginal and joint log-density functions that is satisfied when the conditional distribution arises from an SCM with velocity v. It is intuitive that this is in terms of u and the score, which both characterize local change.
A.3. Score-based ANM Methods
Rolland et al. (2022): Montagna et al. (2023a;b) also use the score function for functional causal discovery. There the focus on finding leaf nodes in a multivariate causal graph, that is, non-cause nodes. This is equivalent to finding the effect variable the bivariate setting. Here, we show how the conditions derived there can be interpreted via the continuity equation. Throughout, let y be the effect variable so that X + Y is the causal graph.
The continuity equation in setting states (7)
8 (v|x) (,((x)),,,(((rr)
(24)
where note for an ANM, o(y,x) = mxx and hence 0y x) = 0. Recall that for an ANM we also have
ply x) = p. (y m())),
(25)
where pr is the density of the noise variable ey. Let 8 = 0 log pe (t), then the continuity equation becomes m (( = - z),
(26)
where we wrote y m(x) the LHS. Now, noting that sy .r) sy (2, v), the component of the joint score, have
S.(E) = -y(( v),
(27)
which holds for all (x, y) where ri(x)) # 0 (thus, for all x if m is injective, corresponding to Condition (2d) in Appendix B.I). The above expression is precisely the general one used by Montagna et al. (2023a) for general non-Gaussian noise. There, they use score estimation to estimate sy (a.,)) y), then fit non-parametric regression model to estimate m and obtain the residuals ē. Then, the equation above says that for the effect variable, it is possible perfectly predict s, y) from é.
Under the Gaussian noise assumption cy NN(,, Rolland et al. (2022); Montagna et al. (2023b) derive more specific equation that is easily derived independently of the continuity equation. Consider the component of the joint score,
(,, v) = sy(y x) = s. - -m()))
(28)
Under the Gaussian noise assumption, se(e) = -e/o². Thus, we have
m(x) x²
sy v)
(29)
This indicates that (x, v) = -2 which is constant over (x, y). Rolland al. (2022) hence devise an algorithm to estimate the Hessian of the log-likelihood (i.e., Jacobian of the score) and select the node with minimum empirical variance as the effect.
B. Proofs
B. Regularity Conditions of SCM-Velocity Correspondence
The correspondence in Theorem 3.2 between SCMs and velocity-initial condition density pairs (v. p(y(2o)) essentially follows from the known results on dynamical systems reviewed in Appendix A.I. The only subtleties are under what conditions the relationships hold over all of I rather than on subset. We discuss them here before proving Theorem 3.2.
1. Px, has full support on 2 and is absolutely continuous with respect Lebesgue measure.
13

---

Distinguishing Cause from Effeec with Causal Velocity Models
2. For any SCM mechanism f:
(a) For each x E I., fæ, )) is bijection I- R.
(b) The mapping (a, 1) + fzr (v)) is continuous for each æ, x', y E R
(c) The mapping x' + for (fz (v)) is differentiable at x' = x for all x, y E R
(d) There is a unique zo E IR such that for all Ey E R, f(zo, Ey) Ey:
3. For any velocity field e:
(a) For each x E IR, v( x) is locally Lipschitz continuous or k-times continuously differentiable.
(b) Satisfies
ee x)1 & a()))vl + b(x).
where and are non non-negative locally integrable (integrable on every compact subset of I) functions.
Within the confines standard practice, these assumptions are not restrictive. Assumption 2(d) may require some care, but easy to achieve in practice. For example, with LSNMs, it becomes
m(xr)) + Mo
which requires each of m, to have unique zero at vn. Assuming that each of m, h have at least one zero, if they are injective then the zero is unique.
In general, assuming that least one such zu exists, a sufficient condition for uniqueness is that x * xf implies that (x,,5y) # [(x', syy for some by IR.
Assumption 3(b) the velocity ensures that a local solution to the ODE dy/dr v(y, x) extends to unique global solution
(v) = y + S u(Peoo.(U)) u) du
whereas assumption 3(a) ensures that
d dx Fz.. (V) Pzo.zo (v) = y
holds for all in the solution domain. See, for example, Arnold (1998, Appendix B) for details.
Proof of Theorem 3.2. First, assume that Px, has full support on R² and is specified by bijective SCM
Y f(X,
Xiley
with noise density Ey Po- Then
(r) for (am (1))
defines a continuous flow. Along with the assumed differentiability in 2(c), by (Arnold, 1998, Thm. B.3.5),
d vly.x) dx far (fz
r'=rr
defines the velocity that generates the flow. As long as the mechanism is bijective on all of IR for each x then this relationship holds over all of ? By the uniqueness zo. p(y(zo) (r)) = (y) is uniquely specified by the SCM.
Conversely, fix pair (v, p(y (o))) such that p(y (2)) has full support on y E R. As long as e is sufficiently regular so as to yield flow (xr v) - P (v) over all of R22 then the flow is unique and therefore
f(X, cg) := Pxx.
with cy ~ po)) uniquely specifies bijective SCM, except for the marginal distribution Px. Assumption 3(b) above guarantees the global uniqueness of the flow- mechanism.
14

---

Distinguishing Cause from Effeect with Causal Velocity Models
B.2. Proof of Theorem 4.1
Proof of Theorem 4.1. By Theorem A.2, the continuity equation is uniquely solved by densities generated by the flow associated with and some initial condition density p(y]z)). (o)) Hence, if Px,y can be represented by SCM with velocity then its conditional p(y|z) will satisfy the continuity equation
d.plylæ) -B, (p(vlx) e(y,x)).
Since Y (and hence Py|x) (x) is assumed to have full support, its density will be strictly positive and the continuity equation can be written
Daplylz) -,, x) v(y, x)) Oyplylz) p(yl+)
p(y(x) x)
02 log p()]2r) = (v,( x) - cly, x )d,, pggp(yl
Noting that 0yp(ulæ) = 0yp(2, y) and adding 0₂ log to both sides, we have
a₂ log v) = -0, x) - x), log p(z, v) + 02 log p(z),
which is (14).
Conversely, if Px,y satisfies (14) then p(y|z) must satisfy the continuity equation and therefore the associated velocity can be used to represent p(ylæ) with the SCM constructed as Theorem 3.2.
B.3. Proofs for Section 5.1
Theorem B.1. Let 82 E Hxx and 82.91 Sa. E Hx.y where Hx, Hx,y are RKHSS induced by bounded kernels liz < B, ka.y = ka @ ky < B. Additionally, let v E C. (R2 ; R). Then, if the score estimators have bounded norm and converge
(x) [v]
in RKHS norm, such that each of || sz s
sl1) s.. Hxxy
((r) $2.9 Hxx
is Opp(n" for some a > 0, then
(v) L(u) On for 3 maxfa,
Proof. For convenience in this proof, we will use the notation 81 = 82, 82 == -8, 83 == -2.9. h1 =
(x)
(3)
8v ho
By
and Z - -X, Pn.y. Note that by assumption, we have that s (81, 82, 83) H. where H Hi H2 X
Hs is
tensor product reproducing kernel Hilbert space with components H₁ Hxx H2 = H3 = H.x. Using this notation, can express Equation (15) and its estimator Equation (18) as
L(v) Bf Z)
(30)
= 1 Σ
L.(())
[(8, 2,)
(31)
n 1=1
where [(s, 2) dataset, and å = (31, $2, $1) are the estimated scores using D,,. We can expand the deviation between the estimated and true loss into two bounding terms
h2 (2)) (2))²
and D, == = id px, is an i.i.d.
[in(()) - L(v)l 1 n
Σ
f(3, Z.) Ezf(s 2)
+ (3)
(32)
n Σ i=1
sup
f(s, 2,) - Ef((s 2) + |E,f(,, Ef(s, Z)1
(33)
where for clarity B(H, = B(H,, M) U B(H2, M) U B(H,, C
H. f(s
(z(())+h(()) + 52(2) + (=))))(())) and Z is an iid copy.
To analyze the LHS term, we first show that the function class {f(s.)): E B(H, } is a separable and complete Carathéodory family. This requires that (i) B(H, M) is a separable, complete metric space, and (ii) 8 + f(8 2)
15

---

Distinguishing Cause from Effeet with Causal Velocity Models
continuous for every 2 E 2 (Steinwart, 2008). (i) Follows immediately from the properties of RKHS's s and the fact that B(H, M) is closed ball in this space. (ii) Follows from the fact that 81, 82, 82, hi hi are all bounded and continuous (note the boundedness of 82, $2, 83 follows from the boundedness of the kernels kay ky). This property also means that F := {f(s..) : B(H,M) E (2) and (B,,, (l >l As result, by Proposition 7.10 in Steinwart (2008) we have
1 " Ef(s, Z)
E sup ME
-
Σ i=1
[8, Z.)
< 2ERadp, (F,
(34)
Where RadD, (F, n) E. (H. Rad(1/2)) Now, note that we can expand
1 - ff Z,)l is the empirical Rademacher average (i.e., {o,) idd
n
[(s, 7) +2 (Z) (Z) + 25 28 (2))2 + (0)( + (2)())++())))))))))))))) (35)
+ (Z) + (2) + 83
+ (Z)))
If we substitute in this definition of into (34), we get the inequality.
D. (..,,,FF,n)) (F.)) & Radp, (Bxq. n)) + Rado, (Bxg. n) Rado, (fh2] @ Bug. n) + 2Rad (BH, @ &BH,, n)
+ 2Radp, (BH, @ BHg @ {12] n) + 2Rad (BHr, @ BHg, @ {h2], n) + 2Radp, (BH, @ {hi).n) + 2Rady, (Bu, @ {hr], n) + 2Radp, (BH, @ [h1] @ {h2), n) + Rado, ((h1) @ (h1).n)
(36)
Where in the above we use the shorthand B+- B(H, M). The tensor product spaces with singletons satisfy the property that if Φ E g @ {h} where G is space of functions and {4} is a singleton, then there is some E G such that 6(z) = g())((>) Vz E Z. Note that the above inequality follows from the fact that, for any two function classes Gr G2. we have by the triangle inequality.
sup (01 (01.92)GG
Σ (o + (2
S sup gn E9
1 7 t=1
ig
(7:)
+ sup go E9
n
0: (4.)
(37)
Now, note that every term in (36) is the Rademacher average of a closed ball in an RKHS of bounded functions. For the terms involving tensor- tensor-produse product spaces without singletons like [hz] } this is by true by definition (since we work with tensor products of bounded RKHS's). To briefly show that this holds for the terms involving tensor product spaces with singletons
(e.g. H/1 @ (hr}), note in general that if H4 is an RKHS with bounded kernel k, then so is H4 @ {h2], if h2 E (Z,RR, because the evaluation functional 6= h & hn + h(2)/1 (=) remains continuous: 82 8-hh) )
(h, k(z, MH hn (1 <
||/a ((, =)]l/|||| Ar. I B B||A||HH Ht. B > O. Similarly, if || I M, then |/h/ Nt. S M also. This is because one can simply choose the norm of H4 @ {h1] as the norm of He. due to the boundedness of h1-
Now, by Lemma 22 in Bartlett & Mendelson (2002), we have Rado, (B(He, M), < C for some C > 0, where B(Ha, M) (h- E Hd || S M) and He is an RKHS with bounded kernel k. This means that
n Σ f(s, n
sup BE
Bf(s, 2)
s 2ERad,, (F, S /2/V2
(38)
for an appropriate constant D > 0, which gives us the desired result for the first term of Equation (33).
Now for the second term. For this section, all expectations are taken over Z only (i.e. not 8). To start, we can expand out f.
Es(s, 2) = (Z)² + s2 (2) + (2) sa 2s (Z) + + 2(sss (Z) + 82 (7) + 83
+ (82(2)ss(2)hall
(39)
Taking differences with the same quantity & gives
Ef Z)- - Ef Z) = - +111::/1000000001100000000000010 (2())) (Z)
), (40)
+2E[(si(Z) (Z()) (|000)/)--BBB((000----0011111111 (z())) (4))
16

---

Distinguishing Cause from Effeet with Causal Velocity Models
Since v E (IRI,,
(IR2,)), R), we know 1 < A1, b2 A2 are bounded, where A1, A2 > 0 are constants. Using Jensen's inequality
therefore have
+ 2E[ss
82(2)8,(2))
(s))) - (2)) + 2A, - 22 83(
(42)
Which can be simplified using a change of notation to
Ef(s, 2) EE, |) - || silLe, + ||s] 8211e, + ||$ -53lle,
(43)
+ 2|| 81š211e, + 2A2l|s153 3183 L. +2A2||s253 $2331121
+2A | 31lle, + 2A1lls32 8allz, + 2A11|s3 $21le,
(44)
where we define sa : 2 s;(2)² and (s,8;) : (z))) There are two kinds of summands above: (i) ||s; 8,1lz, and (ii) ||s;sy - , ll,, for i.j E {1,2,3] 3}. All that remains is to bound each of these terms by sums of terms like ||s; - 8,112,
This is immediate for (i) by the properties of lpp norms. To show (ii) we can simply use the triangle inequality and Cauchy Schwartz:
||s88 8,8,1la, = [|s:(s) - 8,) + -
(45)
I E[s
+ ( I || Syllz, + +15s|leaal1s 8ell
- 5,(2))
(2)) (())
(46)
(47)
s As(Ils) - 8,lle, + ||s: 8:1lea)
(48)
Now, note that for any bounded positive-definite kernel h z² + R with associated RKHS Ha. if l: < B we have ll, I IIf Nu: because [f(=)] = (f(2, I |k(2, 2)ll S B. This means that
[)( -) BJ(, = O 8: 8:11 1, On 2).
(49)
Where the last equality follows from the convergence assumption of the score estimators in the theorem. Combining this result with the convergence result in Equation (38) for the first term of Equation (33) completes the proof.
C. Identifiability of Causal Direction with Bivariate Velocity Models
In order for the causal direction to be identifiable, it cannot be the case that p(æ, y) can be expressed in terms of SCM from the same class both causal directions. For the model in the Y - X direction, let ö(x. y) denote the velocity. Theorem 4.1
can be used to determine conditions under which the causal direction cannot be identified.
For convenience, let π denote log p. for example (Ix y) log p(z, 3), ((y(=) taking a partial derivative in y yields
log ((g))))), and so on. Starting with (14),
0 0,0,T(I,V) v) - 850 -0 ((,
(50)
Similarly, in other model direction, (14) becomes (,T(x, = -0,0(x,00 - derivative with respect to x, we have
+ 0,)( so taking
0,0y5(r, u) = - Pi6(, v) - 02ū(2, y) . (xx v) ū( (f,,))0
(51)
Equating them, we find that the direction is not identifiable if and only if
Zo(y,z) + 0, (),,(,,)0,mtt x)0, v) + 1(,,))037(2, y) = )v((,,) v) + O2 0( (0,0,,)),,711 v) + 6(x T (52)
This proves Proposition 5.2.
17

---

Distinguishing Cause from Effeec with Causal Velocity Models
To this, write the various partial derivatives of ((; y) in terms of the forward model. For SCM with velocity and flow (e, so that the SCM is y Pzo.= (fy), the log joint density can be written
π(x,V) (x) + (o( (v)) + log 1 (v)) = (x) + mu (v)) log 10y.P (1))
and therefore,
(z y) = #o(p2.2. (v)) (v) By 10 0 (2, y) = + to (1. 0n log O 0y (y)
(53)
(54)
- =(
(,
(55)
(z y) =
8e log 2y (v)
(56)
(2 y) #(z) + tv (z0 y) = #o(P2.2. (0)
to 8z (1) 0,
(57)
(u 0z.0p (v) + 0,00 log (v)
(58)
= (1) ())
(v)) - (())
(59)
We have used the identities Pa (3) u0, Pa, (v) and 0g log p2 (v) ,y x) to simplify some of the expressions. Substituting these into (52), if we view Tu. s,, and (x) as given (specified by nature), the result is a PDE for the reverse model velocity, 0. Alternatively, as is common in the literature (Peters & Bühlmann, 2014), we allow x(x)
to vary then we might manipulate some combination of Equations (14) and (50) to (52) to obtain a differential equation for (1) in terms of only the fixed forward model.
more thorough general analysis of identifiability is beyond the scope of this work, though we analyze below the special cases of ANMs and LSNMs. In doing so, obtain new characterization of non-identifiability that holds uniformly over the model class. However, even in that somewhat simple extension of ANMs, the characterizing equation is much more complicated and does not yield an easy interpretation.
C.1. Additive noise models
For ANMs, write the models Y mXX) + Ey with Ey ~ po. and X = mmY) + ez with ez ~po. As shown in the main text, v(y, x² = mz) Putting this into (52) yields
m(x ( v) = m
Hence, the structural functions are mutually constrained, with mnyy) satisfying
m = mn(x)) 8, (x,
(60)
0 (I,
We observe that in the special case of linear ANMs, m(x) = and myyy b. so that (60) implies that each of ( v) (2(xx and TZ v) a (rly) must be constant and therefore Gaussian, i.e., b/a / We also see that if (ar, y) is jointly Gaussian then the ANM not identifiable if and only mz))/m(y) is constant, i.e., each of the model
directions is linear.
Continuing from (60) by taking another partial derivative in x, find that
a m )2)00p(< =.
Moreover,
<(x = x(z) + - (
so that
0g7(2,8) = #o - m(x))
(z,(,)) = #o(y m(x))
= (x)
#oly m(x))
m(=)) = #o -o((v-(<z)m(e)
18

---

Distinguishing Cause from Effeect with Causal Velocity Models
Carrying out the algebra, find
#(z) #===#(=GGGaa H
(61)
where
G(z.v) = m(xx) m(xx) H o - m(x))
r))x)
Ry m(x))
(zx = (y (z )in + fo
- m(x))m(())
to m())) # (ly m(z)) #oly - (m()=)())((())))
#oly m (x))
ri(x))
This is the same differential equation obtained by Hoyer et al. (2008) in their analysis of identifiability using ANM models.
Altemaatvely, Eq. (6) from (Hoyer al., 2008), which also leads the identifying differential equation (61), can be obtained directly from the continuity equation (14) for the backward model Y X. In that case, (51) yields
0, y) m v).
Hence, solving for 1/m(y) and differentiating with respect x,
0, O O,OyT(Z,M) 0, (L
which is Eq. (6) in (Hoyer et al., 2008), from which the rest of their identifiability results on ANMs follow.
C.2. Location-Scale Noise Models
For LSNMs with y = m(x) + ch(m) and = mmyy) + M we have that v(y,x)) = ml(z) + h(x)(v) -m(x)). and similarly ö(x, v) = m)(y) + h(yy) (x= m(v))) Therefore, (51) yields
020 y) = (())a=((x v) + (m(r) + M(v)( ((x)))27(E11))
Dividing by TZTT,, y) and differentiating with respect to x,
(000T(2,))) n(w) (1 (1+ On
(0,T(Z.)) v) BZ m)
-0,
Therefore,
o 0 (,0)
(1+ 0,
O v) BT y)
)
0.
Carrying out the differentiation and simplifying yields (suppressing the (z, y) arguments of m))
[3(0 200 + ((0,T)(0 + [2(0]T)) -00
-
= 0.
Since by using the forward model we have
(1, v) T(x) + To (y- h(x),
in principle this can be solved to characterize, for fixed To, h the set of x(x) for which the causal direction is not identifiable. Alternatvely, substituting
0,0,T(Z, 3) h(x)0, 9) + (m/()) + h(x) - m(x)))af=(xx y)
yields different but equivalent equation.
19

---

Distinguishing Cause from Effeec with Causal Velocity Models
Table 4. Parameter settings for synthetic benchmarks.
Benchmark
fo(2, (v)
Velocity
0 IR
Cy + I 0 (v((),, uddu,
mg
sin (u)
sin(y)
(yy,u)==
cos(u)
cos(y)
sin(x + y)
Sigmoid ANM
Par 0b 0c: 0g: 2x64 MLP weights
c(x) a(x)) Φ (sigmoid(a(z)) + t b(x) (,)
0.2
6n 3x64 MLP weights
m(x) + Ey m(x) + h(=) + 0.2)cy
0.2
0.2
LSNM
mm 02 2x64 MLP weights
0.2
0.2
D. Experiment and Simulation Details
We describe extra details of the experiment and simulation here. All experiments are conducted using the JAX library (Bradbury et al., 2018).
D.1. Synthetic Data Generation
All synthetic benchmarks can be described as generating from an SCM Y = ff((X,cy). In each case 100 samples are drawn from NO,, oß) to generate the 100 datasets. The noise distributions X, Cy are generated as follows.
Ex, Ey ~ (0, I2),
(62)
x= To, = To. (Ey)
(63)
where By is a noise scaling parameter that represents the signal-to-noise ratio (larger oy indicates noisier data). See Table for specifics in each benchmark. To, are randomly sampled (for each dataset the benchmark) triangular monotonic increasing (TMI) maps parametrized as follows:
To, (z)= 1// (fftplus(ff,(z))dd
(64)
and 0, (.(00 are 3x64 MLP parameters.
Analytic Score Calculation For Section 7.3, additionally designed ANM and LSNM datasets with Gaussian noise variables, X ~ (0, 1), ~N,, Other settings are in Table 4. The Gaussian noise variables allow for an analytic score calculation of the marginal and joint score functions, as follows:
Sq (x)
ggussion (x) = -T.
(65)
Denote the mechanism by (ey).. We have
(x, = 02 lolo bogp(zz.) = 02 logp(x)) + Og log x)
(66)
= + 0. log Pe,, (f. (v)) + 02 (log (g))
(67)
- -
(68)
0yfr (v)
and
(, 1) = 0, log logp(1, v) = 0, log (x)) + 0, logp(y x)
(69)
logpe, (f (v)) + 0 (log 0,J 'y)=)(y))
(70)
for (1)0,f (v) + 0gufz (v)
(71)
0, (v)
20

---

Distinguishing Cause from Effeet with Causal Velocity Models
Velocity
Sigmoid
y->x
x->y
y->x
x->y
X->y
y-2x
X->y
y-x
y-2x
X->y
y->x
x>y
ANM
LSNM
y->x
X->y
y->x
X->Y
X->y
y->x
X->y
y->X
y->x
X->y
ソ->X
X->y
Figure 6. Example datasets from each synthetic benchmark.
21

---

Distinguishing Cause from Effeet with Causal Velocity Models
where the derivatives calculated with automatic differentiation. Finally, the marginal density is given
p(u) S. ply x)p(z)de,
(72)
and the score is approximated by Monte carlo approximation of the integral and automatic differentiation of the log-density.
D.2. Velocity Parametrization and Training
We distinguish two types of models: those where the velocity function is a linear combination of basis functions, and those where the velocity function are neural network based.
For the basis models, we use the Adam optimizer with a base learning rate of 0.1, scaled by factor of 1/ log(# of parameters).
The basis functions we use as follows:
P (v, x)
v). T (quad(y
8² xy
(73)
For each of these, our experiments also include appending the following exponential terms:
Peep x²) [e-" e- E E-(22+0?)]
(74)
The linear basis and quadratic basis models are hence parametrized by K 3, 6 real-valued parameters, respectively, When the exponential terms are added, the parameter count is increased further by 3.
All neural networks involved are 3 layer fully connected MLPs with hidden size of 64, and tanh activation functions. For the ANM, we directly parametrize the velocity (i.e., ri) as an MLP, while for the LSNM, we parametrize the functions m and h, and evaluate their derivatives using automatic differentiation to obtain the LSNM velocity (see Table for the specific form).
D.3. Score Estimation
For the Stein score estimate, we use the suggested regularization parameter À = 0.1 (Li & Turner, 2018). For the KDE, we use a regularization parameter of = -2 as suggested (Wibisono al., 2024). For the KDE estimator, we found that the standard Silverman rule of thumb (Silverman, 2018) for obtaining the bandwidth works well, as the optimal bandwidths as proposed in (Wibisono et al., 2024) include unknown constants.
D.4. Benchmarks
For the benchmarks in Table 3, we also remove the data points corresponding to the most extreme 5% of marginal values as an automated procedure, following score estimation. We found that this improved causal discovery performance especially for the real data in the Tübingen dataset. Furthermore, all other benchmarks have sample size of 1000 besides the Tübingen dataset. To minimize hyperparameter search, all datasets are sub-sampled, or re-sampled if necessary, to uniform size of n = 1000.
Continuous Tübingen In addition the datasets with binary variables (pairs 47, 70, 107) and multivariate settings ("pairs" 52, 53, 54, 55, 71, and 105), which are filtered by default, we also identified 29 additional datasets with discrete variables that expected score estimation to fail on, for example, the cause variable in pair 5 is integer-valued (# of tree rings). The full set of removed pairs as follows: =55.8.84..3411102221133333663311
E. Additional Figures and Results
E. Synthetic Results
Here we present additional results on our synthetic dataset: Table is the KDE counterpart to Table the main text, and Table 6 contains results when subsampled 1000.
22

---

Distinguishing Cause from Effeec with Causal Velocity Models
Table 5. KDE results on synthetic data (n 5000).
Model
Velocity
Sigmoid
ANM
LSNM
LIN
68 (59)
30 (27)
45 (29)
40 (32)
B QUAD
65 (61)
33 (29)
30 (19)
15 (5)
V ANM
46 (33)
25 (15)
31(18)
39 (26)
V-I.SMM
65 (78)
47 (54)
(33)
(38))
V-NN
69 (65)
56 (48)
30 (24)
36 (18)
Table Results on synthetic data (subsampled
1000).
Model
Velocity
Sigmoid
ANM
LSNM
KDE
STEIN
KDE
STETN
KDE
STEIN
KDE
STEIN
B-LIN
68 (59)
87 (96)
30 (27)
(88)
35 (22)
39 (30)
40 (32)
49 (60)
B-QUAD
65 (57)
88 (97)
25 (28)
(84)
30 (19)
49 (53)
(5)
51 (58)
V-ANM
46 (33)
81 (96)
25 (15)
40 (23)
31 (18)
77 (85)
39 (26)
52 (54)
V VSLSM
65 (78)
87 (96)
47 (54)
69 (86)
37 (32)
73 (82)
38 (38)
66 (64)
V-NN
69 (65)
90 (98)
56 (49)
58 (69)
30 (24)
58 (66)
36 (18)
48 (57)
LOCI (HSIC)
40 (30)
70 (83)
91 (98)
69 (86)
LOCI (Lik)
46 (66)
49 (61)
44 (52)
31 (31)
CDS
46 (42)
28 (13)
90 (98)
48 (48)
IGCI
66 (76)
53 (67)
34 (22)
32 (20)
RECI
36 (26)
18 (8)
36 (34)
43 (56)
CGCI
48 (42)
66 (71)
71 (88)
55 (69)
E.2. Additional Benchmarks
In Table 7, we present results on the benchmarks of Tagasovska et al. (2020). These are variants of ANM/LSNMs with Gaussian noise variables; the methods that assume Gaussian noise perform well in these settings.
Table 7. Results on benchmark data Tagasovska (2020).
Model
AN
AN-s
LS
LS-s
MNU
KDE
STEIN
KDE
STEIN
KDE
STEIN
KDE
STEIN
KDE
STEIN
B-LIN
76 (87)
14 (5)
97 (100)
15)
63 (79)
19 (25)
84 (96)
52 (53)
68 (79)
52 (51)
B-QUAD
98 (100)
21 (8)
98 (100)
(4)
89 (99)
22 (25)
92 (99)
48 (53)
96 (100)
57 (71)
V ANM
93 (99)
68 (85)
73 (91)
59 (67)
77 (93)
63 (71)
57 (75)
(72)
55 (59)
10
V-LSNM
(98)
(40)
89 (98)
31 (19)
90 (99)
69 (86)
81 (92)
62 (82)
61 (68)
41 (40)
V-NN
97 (100)
22 (6)
94 (100)
33 (17)
84 (97)
19 (20)
74 (92)
63 (79)
55 (53)
53 (64)
LOCI (HSIC)
100 (100)
100 (100)
95 (99)
89 (97)
100(100)
LOCI (Lik)
100 (100)
100 (100)
100 (100)
100 (100)
100 (100)
CDS
99 (100)
92 (99)
77 (87)
60)
66 (70)
IGCI
91 (98)
95 (100)
90 (98)
92 (98)
82 (93)
RECI
19 (6)
32 (18)
28 (14)
43 (43)
23 (7)
CGCI
100 (100)
95 (99)
100 (100)
85 (94)
(99)
E.3. Sample Size Experiments with Known Score
Here, we report full tables of results corresponding to Figure and Figure 5 in the main text, which Table 8. We also repeated the same experiment for a well-specified LSNM, which is reported in Table Finally, Figure 7 shows the visual effects of increasing the sample size, Figure 3 for its counterpart when using the ground truth score. Note score estimation here refers to using the Stein score estimator with the Gaussian kernel.
23

---

Distinguishing Cause from Effeet with Causal Velocity Models
Table Known score sample size experiment results for ANM data. For MSE/GoF, values shown following the convention computed over the 100 datasets in the benchmark: Median (QI, Q3).
n
Success (AUDRC)
MSE (Cause)
MSE (Effect)
MSE (Joint)
GoF (Causal)
GoF (Anticausal)
100
55 (58)
0.26 (0.17, 0.36)
0.36 (0.18, 0.39)
2.66 (1.46, 4,49)
0.39 (0.35, 0.46)
0.42 (0.36, 0,51)
500
61 (75)
0.08 (0.06, 0.12)
0.09 (0.06, 0.14)
2.27 (1.05,4.1) 4.))
0.29 (0.24, 0.32)
0.31 (0.27, 0.36)
1000
74 (88)
0.06 (0.04, 0.08)
0.05 (0.04, 0.08)
2.17 (0.99, 4.26)
0.23 (0.19, 0.27)
0.26 (0.23, 0.30)
2500
(96)
0.03 (0.01, 0.04)
0.03 (0.02, 0.04)
2.07 (0.94, 4.20)
0.18 (0.15, 0.20)
0.22 (0.19, 0.28)
5000
86 (97)
0.02 (0.01,
0.02 (0.01, 0.03)
2.06 (0.92, 4.13)
0.14 (0.12, 0.16)
0.19 (0.16, 0.23)
10000
91 (96)
0.02 0.04)
0.02 (0.01,
2.03 (0.90, 4.21)
0.12 (0.11,
0.19 (0.16, 0.23)
Table Sample size experiment results for LSNM data. For MSE/GoF, values shown following the convention computed the 100 datasets in the benchmark: Median (Q1, Q3).
Success (AUDRC)
MSE (Causc)
MSE (Effect)
MSE (Joint)
GoF (Causal)
GoF (Anticausal)
100
59 (51)
0.25 (0.16, 0.43)
0.29 (0.16, 0.43)
1.28 (0.85, 1.92)
0.65 (0.46, 0.87)
0.67 (0.50, 0.83)
500
53 (65)
0.09 (0.05, 0313
0.09 (0.06, 0.14)
0.85 (0.65, 1.17)
0.31 (0.27, 0.37)
0.33 (0.28, 0.40)
1000
(78)
0.05 (0.03, 0.07)
0.05 (0.03, 0.08)
0.72 (0.59, 1.13)
0.25 (0.22, 0.29)
0.27 (0.24, 0.31)
2500
67 (81)
0.03 (0.02, 0.03)
0.03 (0.02, 0.04)
0.64 (0.49, 1.05)
0.20 (0.17, 0.22)
0.21 (0.18, 0.25)
5000
(87)
0.02 (0.01, 0.02)
0.02 (0.01, 0.03)
0.61 (0.47, 1.08)
0.16 (0.14, 0.18)
0.18 (0.16, 0.21)
10000
(91)
0.02 (0.01,
0.02 (0.01, 0.04)
0.58 (0.45, 1.01)
0.14 (0.13, 0.17)
0.17 (0.14, 0.21)
Sample Size: 100
Sample Size: 500
Sample Size: 1000
Sample Size: 5000
Sample Size: 10000
Figure 7. The estimated causal curves better resemble the ground truth as the score estimation improves in a well-specified LSNM.
24


Reviews:
[{"review_id": "7SwogTyAg4", "reviewer": "Reviewer_xfsW", "summary": "This work studied a class of bijective structural causal models from the perspective of dynamical systems. The identifiablity of the causal model was shown through velocity functions. A loss function is proposed to solve for the velocity function, as well as being used to quantify how well a bijective causal model can be fitted.\n\nOverall, this work is well-written and theoretically grounded.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their comments. \n\n**Extension to multivariate settings**: The multivariate case is not a straightforward extension of the bivariate setting. To extend the velocity interpretation we would need to carefully define a multivariate time. We agree that this extension is interesting, but we have not identified a best approach yet and consider this future work.\n\n**Scalability**: The entire procedure, even for $n=5000$, takes less than 3 seconds for KDE estimators, and less than 20 seconds for the Stein estimator on a 2020 M1 Macbook air. The bottleneck is the score estimation, which is asymptotically $O(n^3)$. However, KDE only requires a single matrix multiplication, while Stein only requires several matrix multiplications and inversions, ahead of optimization. One can also use nearest-neighbour based subsampling to reduce the complexity to $O(n \\log n)$, which enables scalability to datasets with $n \\gg 10000$. Once the score is estimated, optimization of Eq. 18 is a cheap $O(n)$ regression task given the score estimate.\n\nTo give this discussion some detail, we report the average computation times/steps for the experiment reported in Table 9rbC.2 of our [response to reviewer 9rbC](https://openreview.net/forum?id=gV01DWTFTc&noteId=yMKZKiiUzH).\n\n**Table xfsW.1**: Computational performance of the proposed method (averaged over 100 replications; velocity optimization convergence tolerance is 1e-6)\n\n|                                                 | n = 10 | n = 100 | n = 500 | n = 1000 | n = 2500 | n = 4000 |\n|-------------------------------------------------|--------|---------|---------|----------|----------|----------|\n| Score estimation time  (Stein w/ Gaussian kernel)                        | 0.027s | 0.034s  | 0.297s  | 0.865s   | 4.600s   | 11.56s   |\n| Velocity optim time (avg)                       | 0.87s  | 0.75s   | 0.63s   | 0.63s    | 0.75s    | 0.79s    |\n| Velocity optimization steps until convergence (causal) | 150.73 | 133.85  | 78.28   | 71.95    | 73.23    | 71.69    |\n| Velocity optimization steps until convergence (anticausal)            | 165.37 | 130.57  | 82.07   | 76.69    | 71.86    | 73.13    |", "claims_and_evidence": "The claims seem well-supported by theoretical evidence.", "methods_and_evaluation_criteria": "The designs of the loss and causal discovery procedure are sound.", "theoretical_claims": "The theoretical results look solid to me, while I did not check correctness.", "experimental_designs_or_analyses": "The experiment design seems reasonable to me.", "relation_to_broader_scientific_literature": "The work established identifiablity for a wider class of bivariate SCMs.", "essential_references_not_discussed": "No", "other_strengths_and_weaknesses": "No", "other_comments_or_suggestions": "No", "questions_for_authors": "I have no major concerns about this work. I am curious about the extension of the method to the multivariate case. The current causal discovery procedure does not seem scalable. What are potential ways to improve it in the multivariate setting (if identifiablity can be shown)?", "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": "YdIFa18cyt", "reviewer": "Reviewer_9rbC", "summary": "This paper delves into how to distinguish between causes and effects in causal relationships through Causal Velocity Models. It proposes a novel framework that treats bivariate Structural Causal Models (SCMs) as dynamical systems and parameterizes these models using causal velocity. The core idea of this method is to infer the causal direction by estimating the score function, without making any assumptions about the noise distribution.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their careful reading and analysis of our paper. The reviewer seems to have appreciated the main benefits of the proposed method (bivariate causal discovery with minimal model assumptions), with some reservations about its reliance on nonparametric score estimation. Our rationale for the reliance on a nonparametric estimate is as follows:\n\n- We view this as an instance of the fundamental statistical trade-off between strength of assumptions and required data/computation. In causality, ANMs sit near the strong assumptions/small sample end of this spectrum. LSNMs already require more data/computation (i.e., two-stage regression) in general. Our method sits further into the weak assumption/more data regime. So do others, such as doubly-robust estimation methods. We would argue that all of these methods are interesting; each is the right tool for some class of problems, with weak assumptions being preferred in settings with limited a priori evidence for favoring specific models.\n- Our approach is to reduce causal discovery to score estimation. We do not claim to introduce a tool that works well on every data set. Our goal is to design a tool that turns reasonable score estimates, where available, into conclusions about cause and effect. This is conceptually similar to methods that reduce causal discovery to independence tests.\n- We find this conceptually appealing since it cleanly separates the causal aspect of the problem (computed once the score estimate is obtained) from the statistical ones. All statistical aspects are, in a sense, encapsulated in the score estimate—the dependence on sample size, for example, only enters through the score estimate.\n- Since score estimation is itself an active research area, we also point out that our method is agnostic to how the score estimate is obtained; it only depends on the quality of the estimate. As new or improved score estimators for a given problem become available, these can simply be plugged into our method.\n- We also note that recent work on causal discovery in ANMs (e.g., NoGAM [Montagna et al., 2023]) also relies on score estimation as a first stage of the discovery method.\n- For computational complexity, please see our [response to reviewer xfsW](https://openreview.net/forum?id=gV01DWTFTc&noteId=pI8MqzAGXf). \n\n### Discovery with known score\n\nTo highlight the separation of causal/statistical components, we performed the following experiment. Using synthetic data generated from an ANM (Gaussian noise; mean function is a 3-layer MLP with random weights and tanh activation) such that the score could be computed analytically/numerically in both directions, we compute the ground truth scores as input to our method and find that sample size plays no role (beyond $n>10$, which is required to evaluate the criterion at a sufficient number of approximation points), and causal direction can be inferred with certainty. We also report the GoF statistic from Eq. 18, showing that it is an order of magnitude lower on average in the causal direction.\n\n**Table 9rbC.1**: Velocity-based causal discovery with known score (results averaged over 100 replications)\n\n|                     | n = 10 | n = 100 | n = 500 | n = 1000 | n = 2500 | n = 4000 |\n|---------------------|--------|---------|---------|----------|----------|----------|\n| GoF Stat Causal     | 0.0459 | 0.0051  | 0.0042  | 0.0041   | 0.0039   | 0.0023   |\n| GoF Stat Anticausal | 0.1067 | 0.0390  | 0.0365  | 0.0365   | 0.0362   | 0.0283   |\n| Success Rate        | 89%    | 100%    | 100%    | 100%     | 100%     | 100%     |\n| AUDRC               | 0.90   | 1.0     | 1.0     | 1.0      | 1.0      | 1.0      |\n\n### Discovery with estimated score\n\nWe repeated the experiment using the same synthetic data as above, using varying sample sizes to estimate the score and subsequently perform velocity-based causal discovery (as in the experiments reported in the paper). We find that as score estimation improves, the success rate of causal discovery also improves.\n\n**Table 9rbC.2**:  Velocity-based causal discovery with score estimated by Stein score/Gaussian kernel (results averaged over 100 replications)\n|                     | n = 10 | n = 100 | n = 500 | n = 1000 | n = 2500 | n = 4000 |\n|---------------------|--------|---------|---------|----------|----------|----------|\n| Score MSE Cause     | 0.0575 | 0.0432  | 0.0169  | 0.0083   | 0.0043   | 0.0033   |\n| Score MSE Effect    | 0.0638 | 0.0417  | 0.0154  | 0.0099   | 0.0049   | 0.0036   |\n| Score MSE Joint     | 0.1436 | 0.1000  | 0.0386  | 0.0219   | 0.0111   | 0.0081   |\n| GoF Stat Causal     | 0.2279 | 0.1165  | 0.0749  | 0.0568   | 0.0429   | 0.0366   |\n| GoF Stat Anticausal | 0.2413 | 0.1235  | 0.0780  | 0.0627   | 0.0495   | 0.0457   |\n| Success Rate        | 49%    | 55%     | 54%     | 58%      | 68%      | 76%      |\n| AUDRC               | 0.54   | 0.59    | 0.59    | 0.69     | 0.77     | 0.87     |", "claims_and_evidence": "The claims made in the submission are supported by clear and convincing evidence.", "methods_and_evaluation_criteria": "The proposed methods and evaluation criteria make good sense for the problem or application at hand.", "theoretical_claims": "I have examined all the theoretical proofs to ensure their accuracy.", "experimental_designs_or_analyses": "The theoretical analysis of this paper primarily focuses on bivariate causal models. Although the authors mention the potential to extend this method to the multivariate case, specific extension methods and theoretical guarantees have not been discussed in detail.", "relation_to_broader_scientific_literature": "The key contributions of this paper are related to traditional causal discovery methods, such as the Additive Noise Model (ANM) or the Location-Scale Noise Model (LSNM). Unlike traditional causal discovery methods, the method proposed in this paper does not require any assumptions about the noise distribution. This makes the method applicable even when the noise distribution is unknown or complex.", "essential_references_not_discussed": "Most of the essential references have already been cited within this paper.", "other_strengths_and_weaknesses": "Strengths:\n1. Unlike traditional causal discovery methods, such as the Additive Noise Model (ANM) or the Location-Scale Noise Model (LSNM), the method proposed in this paper does not require any assumptions about the noise distribution. This makes the method applicable even when the noise distribution is unknown or complex.\n2. The method proposed in this paper is not only applicable to traditional ANM and LSNM models but can also be extended to a wider range of model categories. By introducing the concept of Causal Velocity, the authors demonstrate how to parameterize more complex causal mechanisms using basis functions or neural networks.\n3. The validity of the method is verified through extensive simulation experiments and benchmark datasets in this paper. The experimental results show that the method can infer causal directions well under various complex data generation mechanisms, especially when existing methods (such as ANM and LSNM) fail, the method proposed in this paper still performs excellently.\n\nWeaknesses:\n1. The method proposed in this paper heavily relies on the accurate estimation of the score function. Although the authors have utilized non-parametric estimation methods (such as KDE and Stein's estimator), the accuracy of score estimation may be affected in finite sample scenarios, especially in the tail regions of data distributions.\n2. Due to the need for non-parametric estimation of the score function and the use of automatic differentiation to compute derivatives of causal velocity during the optimization process, the computational complexity of this method is relatively high, particularly for high-dimensional data or large-scale datasets.", "other_comments_or_suggestions": "I do not have any other comments or suggestions.", "questions_for_authors": "See Weakness.", "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": "z83U7kq7l3", "reviewer": "Reviewer_sNSP", "summary": "The paper proposes a bivariate causal discovery algorithm utilizing velocity models, viewing structural causal models as dynamical systems where the cause variable acts like time. The approach establishes a relationship between causal velocity and score functions of data distributions, which is exploited for distinguishing cause from effect. The method generalizes beyond traditional additive and location-scale noise models, with promising performance on synthetic and real-world datasets.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their careful reading and analysis of our paper. The reviewer lists a few weaknesses and questions that we would like to address. \n\n**Extension to multivariate settings**: Please see our [response to reviewer xsfW](https://openreview.net/forum?id=gV01DWTFTc&noteId=pI8MqzAGXf).\n\n**Performance on cause-effect pair dataset collection**: Indeed, there is room for improvement. There are a number of possible explanations.\n\n- If datasets are generated by a process for which a specific model class (e.g., ANM or LSMN) is well specified, we would expect the specific model-based methods to perform better. Whether that is the case for datasets in the collection is not something we can test.\n- Some of the constituent datasets in the Tuebingen collection have discrete or ordinal data, while our framework requires densities that are continuous on $\\mathbb{R}$. We re-ran the experiments after discarding 29 total datasets that have discrete or ordinal data (e.g., “rings” in pair 5 is integer-valued). On this data, our best performing methods (B-QUAD and B-LIN with KDE), obtained 83% (B-QUAD) and 79% (B-LIN) accuracy (up from 59% and 69%, respectively) which is more competitive with other methods. We also re-ran experiments with competing methods in this setting and found that their performance generally did not improve after removal of discrete/ordinal datasets. Please see our [response to reviewer H4Vs, specifically Table H4Vs.1](https://openreview.net/forum?id=gV01DWTFTc&noteId=V3mM0TgYlj) for full results.\n\n**Violation of assumptions**: \n- Non-invertible mechanisms are not an issue for our method. They would be an issue for making unit-level counterfactual inference, but our method only relies on the distributions and their *representations* via functional models (i.e., SCMs), so there is no loss of generality in assuming that the conditional distribution can be represented by a bijective SCM. This is noted after Eq. (1) in the paper, but we will be sure to emphasize the point more clearly in revisions. \n- Confounders: We have not studied this setting and therefore cannot say anything specific, though we consider it important future work. We note that all methods that assume causal sufficiency would likely have issues in the presence of confounders.\n\n**Computational complexity**: We thank the reviewer for pointing this out. We will add a discussion to the appendix. Please see our [response to reviewer xsfW](https://openreview.net/forum?id=gV01DWTFTc&noteId=pI8MqzAGXf) about computational considerations.\n\n**Performance of B-QUAD**: We believe that the B-QUAD model is generally flexible enough to fit most data, while still being simple enough (it only has 9 scalar parameters, see Eqn 67 + 68) to distinguish the causal direction. See Figure 6 of the Appendix for an example of fitting B-QUAD causal curves to real data. Developing a better understanding of this trade-off is an interesting direction for future work.\n\n**Using the velocity parameterization for SCMs**: A great point–we will add a discussion in revision. In short: In the bivariate scalar setting, the velocity parameterization could be used to define novel classes of functional causal models for causal modelling. Since they generate bijective SCMs, they would inherit the counterfactual identifiability results of BCMs [Nasr-Esfahany et al., 2023]. In particular, they can be used to compute queries at all levels of Pearl’s ladder of causation. We can use numerical integration to evaluate the causal curve (Figure 1). Note also that the abduction and prediction steps are especially simple in BCMs, only involving inverting and forward evaluation of the model. In particular for velocity models, this just corresponds to integrating the velocity from an observed condition $x$ to a query condition $x'$, see the discussion above Eqn 2.", "claims_and_evidence": "- Theoretical claims connecting SCMs are well-supported by the theoretical statements\n- Empirical results are convincing\n- Proofs (although I have only skimmed over them) seem to be sound\n- Authors acknowledge limitations", "methods_and_evaluation_criteria": "- Standard data sets and evaluation criteria\n- Careful evaluation, although not all baseline values are directly reported (only referenced in the original papers)", "theoretical_claims": "- Correspondence between SCMs and velocity-density pairs, and the connection to score functions, is sound\n- The work has an identifiability result (core piece for any causal discovery approach)", "experimental_designs_or_analyses": "Good and common selection of benchmark datasets (also see Methods And Evaluation Criteria)", "relation_to_broader_scientific_literature": "While there are other relevant works, the most important ones were discussed. Overall, a fair selection.", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "Strengths:\n- Novel approach with an interesting perspective using velocity models\n- Generalization of existing models\n- Great theoretical justification\n- Good figures to follow the ideas\n\nWeaknesses:\n- Unclear how to extend it to non-bivariate settings (minor)\n- The cause-effect pair performance is not that impressive (related methods are partially better)\n- Unclear how well the method performs when assumptions are violated", "other_comments_or_suggestions": "- A brief discussion on computational complexity using the empirical approaches (perhaps in an appendix) would be helpful.\n- It is unclear whether the approach could extend beyond the bivariate setting; a small remark on this would help.", "questions_for_authors": "Overall a great paper and well justified approach. Only have a few minor questions/remarks:\n\n- The performance in the cause-effect pair dataset is not that impressive. Are there any insights on the reason why?\n- With respect to the previous point, more insights toward violations of assumptions (e.g., non-invertible mechanisms, violation of causal sufficiency etc.) in a more systematic way could be helpful. E.g., the SIM-C dataset contains confounders.\n- It seems that the B-QUAD model performs the best on average, is there any insight regarding why?\n- An outline of how the proposed approach could be used as functional causal models when modeling a SCM to, e.g., compute Rung 2 and Rung 3 queries in Pearl's ladder of causation can be very insightful (e.g. to reconstruct the noise given a (x,y) pair).", "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": "ZbAmzoMH3w", "reviewer": "Reviewer_H4Vs", "summary": "The authors proposed a novel solution to the bivariate causal discovery problem. The key idea is to view the SCM as a flow. The flow model is learned by posing the continuity constraints (minimizing an objective that forces the continuity equation). The value of this objective is further used to decide the causal direction (a smaller violation of the continuity suggests the causal direction).", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their critical reading of our paper. The reviewer raises a few concerns and points of confusion that we believe can be addressed.\n\n**Reliance on score/density estimate**: Please see [response to reviewer 9rbC](https://openreview.net/forum?id=gV01DWTFTc&noteId=yMKZKiiUzH).\n\n**Questions for Authors**:\n1. **Do flows represent ANMs?** Yes. We can regard the velocity/flow parameterization as a new way to express SCMs. ANMs, PNLs, and LSNMs are special cases (see Table 1). __The results in Tables 2-3 refer to the velocity parametrization/estimation of ANM/LSNM__. We can also express new SCMs easily via their velocities (see Table 1, Fig. 2, and Sec. D.2 for examples). In fact, Thm. 3.2 tells us that bijective SCMs and velocity models are in one-to-one correspondence. One of the advantages of our method is that entirely new classes of SCMs based on velocities can be fit and used for causal inference/discovery. ANMs are known or suspected to be mis-specified in many settings, in which case they may not be reliable for causal inference/discovery. \n1. **Performance of method with varying sample sizes**: The method is sensitive to sample size mainly through the accuracy of score estimation. Given the exact score, Eq. 14 holds pointwise, so there is no statistical element remaining in minimizing Eq. 18. We designed an experiment to evaluate this empirically. Please see our [response to reviewer 9rbC, Tables 9rbC.1 and 9rbC.2](https://openreview.net/forum?id=gV01DWTFTc&noteId=yMKZKiiUzH) for details.\n1. **Comparisons to additional methods**: We have evaluated ANM + independence tests and other methods from the `cdt` package, as well as CGCI, on the Tuebingen dataset (after filtering out discrete/ordinal datatsets, for which our model is mis-specified; see our response to reviewer sNSP) as well as the three simulated datasets in the paper (Table H4Vs.1 below). We find that our method performs best overall, although IGCI works well due to the Gaussian noise when generating the data.\n1. **Difficulty of optimizing loss**: Please see our [response to reviewer xsfW](https://openreview.net/forum?id=gV01DWTFTc&noteId=pI8MqzAGXf) about computational considerations and optimization.\n\n**Table H4Vs.1**: Performance of our method compared to other methods\n| Method | Tuebingen (continuous distns. only) | Periodic | Sigmoid | Velocity |\n|---|---|---|---|---|\n| B-QUAD + KDE | 83% (0.89) | 99% (1.0) | 72% (0.87) | 95% (0.99) |\n| B-LIN + KDE | 79% (0.89)  | 98% (1.0) | 71% (0.82) | 95% (1.0) |\n| LOCI (best setting) | 63% (0.66)  | 86% (0.95) | 50% (0.72) | 61% (0.87) |\n| ANM | 60% (0.59) | 17% (0.10)  | 37% (0.25) | 12% (0.02) |\n| CDS | 61% (0.56) | 22% (0.09)  | 20% (0.08) | 12% (0.02) |\n| IGCI (Gaussian) | 53% (0.65) | 96% (0.99)  | 65% (0.81) | 87% (0.97) |\n| IGCI (Uniform) | 63% (0.67) | 0% (0)  | 16% (0.08) | 94% (0.99) |\n| RECI | 71% (0.88) | 0% (0) | 11% (0.03) | 92% (0.99) |\n| CGCI (best setting) | 61% (0.57) | 47% (0.37) | 65% (0.66) | 72% (0.74) |\n\n**Questions on Experimental Design/Analysis**\n- **LOCI not only SOTA method.** See answer 3 above.\n- **Noise level:**  To our knowledge, methods such as NOTEARS have been shown to be sensitive to marginal variances, in particular where the variance in the effect is larger than that of the cause. Following recent convention in the field first established by [Reisach et al., 2021], we standardize the data prior to causal discovery. Note that we use a noise scale of 3 when generating data to reduce the signal-to-noise ratio and increase the difficulty of causal discovery, but this affects all methods equally after standardization.\n- **Sample size.** See answer 2 above.\n- **Optimization of loss.** See answer 4 above.\n\n\n**Questions on Claims and Evidence**: **Optimized loss in causal direction as sample size diverges.** Please see [Table 9rbC.1](https://openreview.net/forum?id=gV01DWTFTc&noteId=yMKZKiiUzH), which shows that when the exact score is known and used, the optimized loss is small even for relatively small $n$, and appears to tend towards 0 for the causal direction, while it stays an order of magnitude larger in the anti-causal direction. For estimated scores, [Table 9rbC.2](https://openreview.net/forum?id=gV01DWTFTc&noteId=yMKZKiiUzH) shows a similar trend towards increasing discovery accuracy as score estimation improves with increasing $n$, which is supported theoretically by Thm 6.1.\n\n**On Other Comments/Suggestions/References:** Thank you for the detailed suggestions. Due to space restrictions we cannot respond to each in detail. We will add the suggested clarifications and references in revisions. We would like to clarify that in the context of Sec. 2.3, $s$ can be viewed as a “starting time” and $t$ as a time variable at which we evaluate the flow that was started at time $s$. Sec. 2.3 is background on flows and ODEs. Starting with Sec. 3, we substitute a cause variable $x$ to play the mathematical role of time, as stated at the start of Sec. 3.", "claims_and_evidence": "Experimental evidence is not sufficient to support the claims. For example, I'd like to see if the optimized loss (in the causal direction) will be close to zero when the sample size goes to infinity.", "methods_and_evaluation_criteria": "The methods and evaluation criteria are acceptable.", "theoretical_claims": "I only read through the claims but didn't check the proofs.", "experimental_designs_or_analyses": "I have several suggestions for improving the experiment presentation.\n\n1. I don't think LOCI is the only one \"SOTA\" method. The comparison with other methods is necessary rather than referring to other papers. You may put the results in the supplementary.\n\n2. On synthetic data. I don't think the authors should fix the noise level. Some methods may be sensitive to the noise level, say NOTEARS.\n\n3. I'd like to see how the sample size will affect the accuracy of your methods of causal discovery.\n\n4. Will the optimization of the loss be an issue? I'm not sure how fast the algorithm is.", "relation_to_broader_scientific_literature": "TBD", "essential_references_not_discussed": "I believe the following paper should be discussed.\n\nTu, Ruibo, et al. \"Optimal transport for causal discovery.\" International Conference on Learning Representations, ICLR, 2022.\n\nAdditionally, there is an extensive body of research on bivariate causal discovery, including methods such as Conditional Divergence-based Causal Inference (CDCI) and Maximal Correlation-based PNL. The author may wish to engage in a comparative analysis with these approaches.", "other_strengths_and_weaknesses": "Strengths:\n- The flow view of SCM is comprehensive.\n\n\n\nWeaknesses:\n-  The method relies on the accuracy of existing density estimation methods.", "other_comments_or_suggestions": "What makes me confused:\n1. Notations in sec.2.3 are not clearly explained. $s$ for start time? $t$ for termination? I often got confused with $x$ and $t$ in this paper.\n\n2. \"Scores\" in Section 4. What do the \"scores\" refer to? It seems the \"scores\" are term from the \"score function\" in Statistics rather than the causal score for distinguish the cause and effect. I suggest that the authors point out the causal score and its physical meaning explicitly in this section. To me, the physical meaning is simply the violation magnitudes of the continuity equation. \n\n3. Notations like $\\dot{m}$ are not defined. In convention, it refers to the first-order derivative w.r.t time. Is that true here?\n\nIn Figure 2, please consider showing the axis names and the SCM expression for each subfigure. BTW, I'm not sure what the purpose of showing quadratic and quadratic LSNM is.", "questions_for_authors": "1. Take ANM for example, can we regard the flow model as another representation? If so, the difference of the proposed method lies in how to learn the underlying functions and the decision criterion (proposed loss / HSIC test). Is it possible to furthur analysis the advantage of the proposed method compared to ANM? More effective learning or the decision rule is better? \n2. Experiments when increasing/decreasing sample size.\n3. Comparison with independent test-based methods and other recent methods.\n4. Is the objective easy to optimize? What about the learning time? Hope to see some discussion.\n\n\nI've updated my score accordingly.", "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:

Meta-review: Reviewers overall agree this submission is a solid theoretical contribution supported by experiment. I therefore recommend acceptance.
Decision: Accept (poster)