Year: 2021

Abstract:
Approaches based on Functional Causal Models (FCMs) have been proposed to determine causal direction between two variables, by properly restricting model classes; however, their performance is sensitive to the model assumptions, which makes it difficult for practitioners to use. In this paper, we provide a novel dynamical-system view of FCMs and propose a new framework for identifying causal direction in the bivariate case. We first show the connection between FCMs and optimal transport, and then study optimal transport under the constraints of FCMs. Furthermore, by exploiting the dynamical interpretation of optimal transport under the FCM constraints, we determine the corresponding underlying dynamical process of the static cause-effect pair data under the least action principle. It provides a new dimension for describing static causal discovery tasks, while enjoying more freedom for modeling the quantitative causal influences. In particular, we show that Additive Noise Models (ANMs) correspond to volume-preserving pressureless flows. Consequently, based on their velocity field divergence, we introduce a criterion to determine causal direction. With this criterion, we propose a novel optimal transport-based algorithm for ANMs which is robust to the choice of models and extend it to post-noninear models. Our method demonstrated stateof-the-art results on both synthetic and causal discovery benchmark datasets.

heading: 1 INTRODUCTION
text: Determining causal direction between two variables is a fundamental and challenging causal discovery task (Janzing et al., 2012). Conventional constraint-based and score-based causal discovery methods identify causal structures only up to Markov equivalent classes (Spirtes et al., 2001), in which some causal directions are undetermined. To address this challenge, properly constrained functional causal models (FCMs) have been proposed. FCMs represent the effect as a function of its cause and independent noise, and can help identify the causal direction between two variables by imposing substantial structural constraints on model classes, such as additive noise models (ANMs) (Shimizu et al., 2006; Hoyer et al., 2008) and post-nonlinear models (PNLs) (Zhang and Hyvärinen, 2009b). While some of the models, such as PNLs, are highly flexible, the constraints are still restrictive and difficult to interpret and relax. Inevitably, the performance of these methods is sensitive to model assumptions and optimization algorithms, especially in real-world applications.
To handle the mentioned issues, we consider FCMs from a dynamical-system view. By augmenting a time dimension for the static causal discovery task, we interpret FCMs with dynamical causal processes under the least action principle (Arnol’d, 2013). The new interpretation connects FCMs with a large class of models in dynamical systems. It then provides more freedom to model causal influences, possibilities to derive new causal discovery criteria, and a potential direction to generalize causal models with identifiable causal direction.
In particular, we exploit the above idea by leveraging the intrinsic connection between FCMs and optimal transport. Optimal transport is originally introduced by Monge (1781), which has been ap-
FIX (GA77, pBTi) plied in a large range of applications, not only because it is a natural way to describe moving particles (Ambrosio et al., 2012) but also because of its recent improvement of the computational methods (Cuturi, 2013; Kolouri et al., 2019). Recently, it has also been largely applied to generative models for measuring the distance of probability distributions (Arjovsky et al., 2017; Kolouri et al., 2018; Genevay et al., 2018). Among different optimal transport definitions, the L2 Wasserstein distance got extensive applications in statistics (Rachev and Rüschendorf, 1998), functional analysis (Barthe, 1998), and multiple disciplines of science (McCann, 1997; Otto, 1997). The dynamical formulation
of the L2 Wasserstein distance is introduced by Benamou and Brenier (2000) for relaxing the computational costs. We find that in the context of the dynamical formulation, FCMs can be connected with optimal transport. Furthermore, with the dynamical interpretation of optimal transport, one can naturally understand FCMs from a dynamical-system view, which makes it possible to derive new criteria to identify causal direction. Moreover, it also enables us to develop practical algorithms with optimal transport for static causal discovery tasks in the bivariate case without learning a regression model. Accordingly, our main contributions are:
1. Dynamical interpretation of FCMs in the bivariate case. We provide dynamical interpretations of optimal transport under the constraints of FCMs. Furthermore, we introduce a time variable, determine the underlying dynamical process under the least action principle (Arnol’d, 2013) for the static bivariate causal discovery task, and characterize properties of the corresponding dynamical systems (Sec. 3.1 and Sec. 3.2).
2. A criterion of determining causal direction between two variables for ANMs. We study the corresponding dynamical systems of FCMs, and prove that ANMs correspond to volume-preserving pressureless flows. Moreover, based on the divergence of their velocity fields, we propose a criterion to determine causal direction, and show that under the identifiability conditions of ANMs, it is a valid criterion for ANMs, which can be extended to PNLs directly (Sec. 3.2).

heading: 3. An optimal transport based approach (DIVOT) for distinguishing cause from effect between two
text: variables. DIVOT inherits the advantages of one-dimensional optimal transport. It can be computed efficiently, and does not require independence tests, learning a regression model, or deriving likelihood functions for complicated distributions. Experimental results show that our method is robust to the choice of models, and has a promising performance compared with the state-of-the-art methods on both synthetic and real cause-effect pair datasets (Sec. 4 and Sec. 6).

heading: 2 PRELIMINARIES

heading: 2.1 OPTIMAL TRANSPORT
text: In this work we mainly use the L2 Wasserstein distance and we introduce its framework following the notations and definitions of (Benamou and Brenier, 2000) in this section.
FIX (pBTi)The underdetermined Jacobian problem. Suppose that two density functions, ρ0(x) and ρT (x) where x ∈ Rd, are non-negative and bounded with total mass one. The transfer of ρ0(x) to ρT (x) is realized with a smooth one-to-one map M : Rd → Rd. The Jacobian problem is to find M that satisfies the Jacobian equation:
ρ0(x0) = ρT (M(x0))|det(∇M(x0))|, xT =M(x0), (1) where det(·) denotes determinant. This is an underdetermined problem as many maps can be the solutions. Thus, a natural way is to choose the optimal one in the sense of lowest cost depending on applications. A common cost function is the Lp Wasserstein distance (Kantorovich, 1948).
Lp Wasserstein distance and its one-dimensional closed-form solution. The Lp Wasserstein distance between p0 and pT , denoted by Wp(p0, pT ), is defined by Wp(p0, pT )p = infM ∫ |M(x0) − x0|pp0(x0)dx0, where p ≥ 1. We say that the optimal transport M∗ achieving the infimum solves the optimal transport (or Lp Monge-Kantorovich) problem. In this work, we mainly use the square of the L2 Wasserstein distance, denoted by W 22 . Moreover, following the approach of Kantorovich (1948), the one-dimensional (1D) Lp Wasserstein distance has a closed-form solution, e.g., the 1D optimal solution of W 22 is M∗ = P−1T ◦ P0, (2) where P0 and PT are the cumulative distribution functions for p0 and pT , and “◦” represents the function composition. In practice, the 1D optimal solution in Eqn. (2) can be computed with the average square distance between the sorted samples from p0 and pT (Kolouri et al., 2019).

heading: 2.2 FUNCTIONAL CAUSAL MODELS
text: A FCM represents the effect Y as a function f(·) of the direct cause X and independent noise Ey , where function f describes the causal influence of X on Y , and Ey is the exogenous vari-
able/noise. Without any additional assumption on the functional classes, the causal direction is not identifiable (Hyvärinen and Pajunen, 1999; Zhang et al., 2015a). Roughly speaking, this is because given variable pair (X,Y ), one can always construct Y = f(X,Ey) and another different FCM, X = f̃(Y,Ex), such that both of them are consistent with the data distribution and have independent “noise” (Hyvärinen and Pajunen, 1999; Zhang et al., 2015a). Several works further introduce proper assumptions on the functional form of the FCMs and the distributions of the involved variables, and then identify causal direction in the bivariate case, i.e., the independence of cause and noise only holds in the causal direction. Examples include the ANM and PNL models:
ANM: Y = g(X) + Ey; (3) PNL: Y = h(g(X) + Ey), (4)
where g and h are nonlinear functions and h is invertible.

heading: 3 DYNAMICAL INTERPRETATION OF FUNCTIONAL CAUSAL MODELS
text: We first show the connection between FCMs and optimal transport in Sec. 3.1. In Sec. 3.2, we further elaborate the analogy between the optimal transport problem and the causal direction determination problem. We then study the optimal transport under the constraints of FCMs, show the corresponding dynamical systems of FCMs, and characterize the properties of such systems.

heading: 3.1 CAUSALITY AND OPTIMALITY OF THE MAPS IN THE JACOBIAN PROBLEM
text: As introduced in Sec. 2.2, FCMs are used to approximate the true data generation process. Given the FCM, Y = f(X,Ey), we rewrite it in a Jacobian equation way,
xT = [ X Y ] = [ Ex
f(X,Ey)
] =M ([ Ex Ey ]) =M(x0), (5)
where x0,xT ∈ R2, their probability densities p0, pT ≥ 0, and M : R2 → R2. As an analogy to the mass transfer scenario (Monge, 1781), we consider the samples of independent noise Ex and Ey as the particles of materials in the mass transfer problem and regard the map M in Eqn. (5) as a special transformation of the independent noise samples. As shown in Fig. 1, one can consider the data points are transferred from the original positions (which are unmeasured) in the plane Ex–Ey at time 0 to the observed positions in the plane X–Y at time T . Such transformation considers the transfer as a causal process which moves the unmeasured independent noise x0 = [Ex, Ey]′ 1 and consequently leads to the observations xT = [X,Y ]′. From the perspective of FCM-based causal FIX (pBTi) discovery approaches, causal influences are represented by FCMs which represent the effect as a function of its direct cause and an unmeasured noise satisfying the FCM constraints:
1“′” denotes the transpose of vectors or matrices.
(i) The map constraint: the values of X are determined by the values of its corresponding noise, i.e., X = Ex, while the values of the effect depend on cause X and noise Ey;
(ii) The independence constraint: the noise terms are independent, i.e, Ex is independent of Ey .
How is causality defined with the FCM constraints connected with the optimality of maps in optimal transport? The optimal transport M∗ with the minimal Lp Wasserstein distance is not necessary to be the causal one given by Eqn. (5), because it has no information about the FCM constraints or the true data generation process. In other words, given two sample sets of x0 and xT , the couplings FIX (pBTi) given by optimal transport are not necessary to be the ones generated from the ground-truth FCM.

heading: 3.2 DYNAMICAL INTERPRETATION OF FCMS: OPTIMAL TRANSPORT UNDER THE FCM CONSTRAINTS
text: In this section we jointly consider causality and optimality of the maps in the Jacobian problem. It provides both a causal sense of the transformation and a dynamical view of FCMs. We first recap the dynamical formulation of theL2 Wasserstein distance, and then study such dynamical systems under the FCM constraints. Furthermore, we show their properties under the FCM and ANM constraints.
Dynamical interpretation of the L2 Wasserstein distance. Benamou and Brenier (2000) formulate the L2 Monge-Kantorovich problem as a convex space-time minimization problem in a continuum mechanics framework. Fixing a time interval [0, T ], they introduce the concepts of smooth time-dependent density ρ(t,xt) ≥ 0 and velocity fields v(t,xt). When they are clear from context, we denote them by ρ and v. Because we are considering the bivariate case, xt ∈ R2 and v ∈ R2. Then, they give the dynamical formulation of W 22 :
W 22 (p0, pT ) = inf ρ,v T ∫ R2 ∫ T 0 ρ(t,xt)|v(t,xt)|2dxtdt, (6)
s.t. { initial and final conditions: ρ(0, ·) = p0, ρ(T, ·) = pT the continuity equation: ∂tρ+∇ · (ρv) = 0 .
They prove that the solutions of (6) are pressureless potential flows, of which the fluid particles are FIX (pBTi)
not subject to any pressure or force and their trajectories are determined given their initial positions and velocities or given their initial and final positions. Suppose that M∗ is the solution given by W 22 . The corresponding flows follow time evolution equation,
xt = x0 + t
T v(t,xt), where v(t,xt) = v(0,x0) =M∗(x0)− x0 and t ∈ [0, T ]. (7)
The time evolution equation shows that xt is just a convex combination of x0 and M∗(x0).
As an analogy between the optimal transport problem and causal direction determination, the density ρ and the velocity v of moving particles can be considered as the probability density and the velocity of changing values of data points. Moreover, the dynamical interpretation of the L2 Wasserstein distance introduces a time variable and provides a natural time interpolation ρ(t,xt) of ρ0 and ρT together with the velocity field v(t,xt) under the least action principle. Similarly, we can also have the natural time interpolation p(t,xt) between p0 and pT as well as the velocity field v(t,xt) under the least action principle, which is the dynamical interpretation of FCMs.
Dynamical L2 Wasserstein distance under the FCMs constraints. First, we introduce FCM constraints in the context of the dynamical L2 Wasserstein distance. According to time evolution equation (7), we know that the velocity is fully determined by the initial and final values of xt. We first consider FCM constraint (i). For the initial and final values of the cause, the value of its observation X is equal to its noise value Ex in Eqn. (5). Consequently, denoting v = [vx, vy]′, where vx and vy represent the velocities along the Ex/X-axis and Ey/Y -axis respectively as shown in Fi.g 1, we know that ∀t ∈ [0, T ], xt ∈ R, vx(t, xt) = 0. As for FCM constraint (ii), it implies that the initial noise of cause and effect, corresponding to Ex and Ey at the time 0, are independent. Therefore, we have the FCM constraints for the dynamical L2 Wasserstein distance,
(I) The map constraints: v(t,xt) = [vx(t, xt), vy(t, yt)]′, ∀xt, yt ∈ R, t ∈ [0, T ], vx(t, xt) = 0;
(II) The independence constraint: Two random variables (cause and effect) at the initial time have the joint probability density function p0(x0) and they are independent.
Second, we characterize the properties of the dynamical optimal transport (L2 Wasserstein distance) under constraints (I) and (II). According to Eqn. (5), the form of M∗ is determined as
M∗(x0) =
[ Ex
f(Ex, Ey)
] ; (8)
otherwise, the FCM constraints will be violated. Moreover, the two-dimensionalW 22 under the FCM constraints can be computed with the one-dimensional W 22 as shown in Prop. 1 (the derivation is in App. D), Proposition 1. Under constraints (I) and (II) the square of the L2 Wasserstein distance between p0 and pT is
W 22 (p0, pT ) = EEx [ W 22 ( p(Ey), p(Y |Ex) )] . (9)
Furthermore, we consider the constraints of ANMs and characterize the corresponding dynamical systems for now, which can be directly extended to PNLs as mentioned in Sec. 4.2. Based on constraints (I) and (II), we further introduce the ANM constraint,
(III) the effect is the sum of noise and a nonlinear function of cause as defined in Eqn. (3). Theorem 1 ( Zero divergence of the velocity field ). Under constraints (I) and (II), the dynamical systems given by the L2 Wasserstein distance are pressureless flows. Further under ANM constraint (III), they become volume-preserving pressureless flows, of which the divergence of the velocity field, v(t,xt) = [vx(t, xt), vy(t, yt)] ′, satisfies
div v(t,xt) = ∂vx(t, xt) ∂xt + ∂vy(t, yt) ∂yt = 0,∀t ∈ [0, T ], xt, yt ∈ R,
where div is the divergence operator in vector calculus.
Thm. 1 determines the corresponding dynamical systems of FCMs and ANMs by analyzing their densities and velocity fields (the details of the proof are in App. D), and shows a property of the corresponding dynamical systems of ANMs. Importantly, the property indicates a potential criterion for causal direction determination, that in the causal direction the divergence of the velocity field is zero everywhere, while in the reverse direction it may not always hold. Next, we will verify the criterion rigorously, propose an algorithm based on it, and show the extension for the PNL cases.

heading: 4 CAUSAL DIRECTION DETERMINATION WITH OPTIMAL TRANSPORT
text: In this section we define a divergence measure as a criterion of determining causal direction between two variables for ANMs. Based on the criterion, we then provide an algorithm to identify the causal direction, named by the divergence measure with optimal transport (DIVOT).

heading: 4.1 DIVERGENCE MEASURE AS A CRITERION OF DETERMINING CAUSAL RELATION BETWEEN TWO VARIABLES
text: We first define the divergence measure and then show that it is a valid criterion for identifying the causal direction in the bivariate case under the identifiability conditions of ANMs, as shown in Prop. 2 (the proof is in App. D). Proposition 2 (Divergence measure as a causal discovery criterion). Define the divergence measure,
D(v) = ∫ R2 |div v|2p0(x0)dx0 = Ex0 [ |div v|2], (10)
where v = M∗(x0) − x0. Suppose that constraints (I), (II), (III), and the identifiability conditions of ANMs (Hoyer et al., 2008) are satisfied. The divergence measure of the corresponding dynamical system satisfies D(v) = 0 if and only if X is the direct cause of Y .
Nevertheless, there are some challenges to compute the divergence measure. For example, we need to solve a two-dimensional optimal transport problem and compute the derivative of a velocity field for all samples. Another challenge with computing M∗ is that we in general have no information about x0. Such issues are all solved with DIVOT.

heading: 4.2 PROPOSED METHOD: DIVOT
text: We first provide an overview of the algorithm, DIVOT, for determining the causal direction between two variables, and then introduce the four steps to compute the divergence measure.
Overview of DIVOT. Alg. 1 is based on the divergence-measure criterion of causal direction determination (Prop. 2) and determines causal direction between two variables. Given data of X and Y , we want to infer whether X causes Y (X → Y ) or Y causes X (Y → X). We compute the divergence measure of both directions. And the one with the smaller value of its divergence measure is the causal direction. Here, we present the algorithm for the experiments, and we also provide the
NEW (xS6T)
modification of the algorithm for including the independent case and the significance of the results in App. A.
Algorithm 1: DIVOT: Causal direction determination using the divergence measure with optimal transport
Input: data {(xi, yi)} and noise distribution p(E; θ0)
Output: X → Y or Y → X . 1 Def DIVOT({(xi, yi)}, p(E; θ0)): 2 data, θ = {(xi, yi)}, θ0 3 LX→Y = Div(data, p(E; θ)) 4 data, θ = {(yi, xi)}, θ0 5 LY→X = Div(data, p(E; θ))
6 if LY→X < LX→Y then 7 return Y → X 8 else 9 return X → Y
10 Def Div(data, p(E; θ)): 11 noise = Sampling(p(E; θ)) 12 datas, noises = OT(data, noise) 13 loss = VarDiv(datas, noises) 14 return minθ loss
Noise data generation: The first step of computing the divergence measure: To compute the divergence measure, we need to know the velocity field v as defined in time evolution equation (7). It requires the couplings of the data of x0 = [Ex, Ey]′ and xT . But in the bivariate causal discovery task, only the data of xT = [X,Y ]′ are given. Therefore, as shown in Line 11 of Alg. 1, we first deal with the issue due to the lack of the noise data of x0, denoted by {(eix, eiy)}. To obtain the noise data, we want to assume a multivariate probability distribution of x0 with the density p0(Ex, Ey), and then sample the data from it, represented by (eix, e i y) ∼ p0(Ex, Ey). Furthermore, due to the FCM constraints, we know that p0(Ex, Ey) = p0(X,Ey) = p(X)p(Ey). So we only need to assume the probability distribution of Ey and parameterize it with θ, denoted by p(Ey; θ). Suppose that the dataset of xT with N samples is given, denoted by {(xi, yi)}N . We first sample a data set of Ey with the sample size N , denoted by {eiy}N , e.g., in the experiments of this work, we use the simplified reparameterization trick,
eiy = f noise θ (e source y ) = θ × esourcey and esourcey ∼ N (0, 1)/U(0, 1), (11)
where fnoiseθ is a monotonic function, and e source y is sampled from a standard normal distribution or a uniform distribution. As the monotonic fnoiseθ can be very flexible, we can represent flexible noise distributions with monotonic neural networks as in (Huang et al., 2018). Next, we randomly match the data of X , denoted by {xi}N , with {eiy}N , which gives the {(xi, eiy)}N as the data of x0.
Optimal transport: Find couplings of the observation and the generated noise. Computing the divergence measure requires the couplings of the data of x0 and xT , because we need v =M∗(x0)− x0; in other words, given the data of x0 and xT , we need to solve a two-dimensional optimal transport problem, which gives M∗ as in Line 12 of Alg. 1. According to Prop. 1, solving the two-dimensional optimal transport problem under the FCM constraints is equivalent to solving onedimensional optimal transport problems. More specifically, the expectation in Eqn. (9) of Prop. 1 is computed by using N samples of Ex, denoted by {eix}N , and then the two-dimensional W 22 in Eqn. (9) is computed with N one-dimensional W 22 , i.e., W 2 2 (p0, pT ) = 1 N ∑ iW 2 2 (p(Ey), p(Y |Ex = eix)). Importantly, solving the one-dimensional optimal transport problem or computing the onedimensionalW 22 can be implemented with a sorting operation as in (Kolouri et al., 2019). Therefore, the couplings are found by which given each value of Ex we find the corresponding values of Ey and Y from {(eix, eiy)}N and {(xi, yi)}N , and then match such sorted values of Ey and Y .
Variance-based divergence measure. Suppose that the previous step has found the couplings, i.e., given any sample (eix, e i y) we have its corresponding (xi, yi). Then, we will compute the divergence of a velocity field as in Line 13 of Alg. 1. According to the definition, |div v|2 = (dvydy )2. Since all the eiy are matched with yi, we have v i y = yi − eiy . A straightforward way to approximate the derivative is using its nearest neighbour pair (eixb, e i yb) and (xbi, ybi) and then approximate it with ddy (vy)|y=yi = (yi−eiy)−(ybi−e i yb)
yi−ybi ; however, it suffers the following issues (especially in the few-sample scenario) : (a) the denominator is in general a small number and the distance to the nearest neighbour can be large in the few sample case, which makes the computation unstable and inaccurate; (b) the deviation on the X-axis makes the approximation a biased estimate especially when the gradient of g(X) + Ey at X = x is large. Therefore, we propose the variance-based divergence measure working better in practice. It is straightforward to see that under constraints (I) and (II), the value of the divergence measure of an ANM is zero if and only if the value of the variance-based divergence measure is zero, which is defined as
Dvar(v) = EX [V[Vy|X |X]], (12) where V[Vy|X |X] represents a conditional variance of the velocity field Vy (a random variable) at position X at the initial time. For example, V[Vy|x|X = x] represents the variance of all the velocities vy at the positions where X = x at the initial time. Let {(xi, yi)}N denote the causeeffect pair dataset with the sample size N , and {viy|x}Nx denote all the velocities viy at position X = x at the initial time, the sample size is Nx, and their mean value is vy|x. Then, V[V y|x|X = x] = ∑Nx i=1(v i y|x − vy|x)2/(Nx − 1), and
Dvar(v) ≈ 1
N ∑ x∈{xi}N ∥∥sort(−→yx)− sort(−→ey)− ave(−→yx −−→ey)∥∥22 Nx − 1 , (13)
where {xi}N represents the set of all the values ofX and its sample size isN ;−→yx is the vector of the Y samples where X = x; −→ey is the vector of the Ey samples where Ex = x; sort(·) sorts a vector; ave(·) computes the vector mean; and‖·‖22 is the square of a `2 norm.
Minimization w.r.t θ. Given the data {(xi, yi)}N and the density p(Ey; θ), we can compute the divergence measure with the variance-based method. Note that we only initialize θ with some random value, and p(Ey; θ) is not necessary to be the true distribution or even significantly different from the true one, which can lead to the wrong result of the divergence measure. Therefore, as shown in Line 14 of Alg. 1, we minimize the divergence measure w.r.t. θ. For the minimization, one can derive the gradient w.r.t θ in a simple parameterization case as (11), while in the complex case one can use auto-differentiation. According to Prop. 2, in the causal direction the divergence measure is zero if and only if p(Ey; θ∗) with the optimal parameter θ∗ is the true noise distribution a.e., implied by the identifiability of ANMs. In this work, we used autograd and RMSProp (gradient descent) of JAX for the minimization.
Extension to PNLs. The measure (13) can be directly extended to the PNL cases. Since h in PNL (4) is an invertible function, by considering h−1(Y ) as a new random variable, h−1(Y ) = g(X)+Ey is an ANM. Thus, for PNLs, under the identifiabililty conditions of PNLs (Zhang and Hyvärinen, 2009b), Prop. 2 still holds; and we only need replace −→yx in (13) with −→̃ yx = f PNL ω ( −→yx), where fPNLω is an invertible function, e.g., it can be the simplified version of (Zhang et al., 2015a), −→̃ yx = f PNL ω ( −→yx) = −→yx + ωa · tanh(ωb · −→yx + ωc), (14) where ωa, ωb, and ωc are positive scalars. Moreover, fPNLω can also be the monotonic neural networks, such as (Huang et al., 2018). NEW (GA77)Extension to the multivariate case. For simplicity and clarity of the paper, we focus on elaborating the connection between FCMs and dynamical systems and developing the theoretical basis and the method in the bivariate case, which are essential for the further development of FCM-based causal discovery methods. In the case of multiple variables, one can use a constraint-based method to find the causal skeleton (the undirected causal graph) and then use the extension of our method for the edge orientation similar as (Zhang and Hyvärinen, 2009b; Monti et al., 2020; Khemakhem et al., 2021). See App. B for details.

heading: 5 RELATED WORK
text: There are mainly two types of causal discovery methods for static causal direction determination between two variables. The first one introduces model assumptions of FCMs to achieve the identifiability of causal direction. Most of such methods are based on LiNGAM (Shimizu et al., 2006), ANMs (Hoyer et al., 2008; Mooij et al., 2009) and some of them are based on the more general models, PNLs useing MLP for representing nonlinear functions (PNL-MLP) (Zhang and Hyvärinen, 2009b), PNLs using warped Gaussian process and mixture of Gaussian noise (PNL-WGP-MoG) (Zhang et al., 2015a). Recently, Khemakhem et al. (2021) propose an autoregressive flow based model (CAREFL), of which the assumption is more general than ANMs and stricter than PNLs. They commonly apply (non)linear regression to learn the function g in (3) (or together with h in (4)), and then test the independence between the independent noise Ey in (3) (the residual) and the cause X . And the independence test is commonly Hilbert-Schmidt independence criterion (HSIC) (Gretton et al., 2005); however, as argued by Yamada and Sugiyama (2010), the kernel width limits its practical use and its common heuristic value limits the flexibility of the function approximation. Moreover, there are other criteria proposed in the first type of methods, such as the likelihood ratio (LLR) (Hyvärinen and Smith, 2013), maximum likelihood based criterion (MML), regression error based causal inference (RECI) (Blöbaum et al., 2018) and mutual information (Zhang and Hyvärinen, 2009a; Yamada and Sugiyama, 2010). Nevertheless, the common issue of the first type of methods is that they restrict the model classes and their performance depends on model assumptions. The second type of methods achieve idenifiability by proposing other principles instead of restricting model classes, such as exogeneity (Zhang et al., 2015b), randomness of exogeneous variables (Entropic) (Compton et al., 2021), and independent causal mechanisms, e.g., GPI (Stegle et al., 2010) and IGCI (Janzing et al., 2012).

heading: 6 EXPERIMENTS
text: We demonstrate and evaluate our method on the synthetic and real-world cause-effect pair data (Mooij et al., 2016). Moreover, we also provide the experiments and the discussion of our method
NEW (GA77, pBTi) in the presence of unknown confounding in App. C. The details of implementation and experiments are in App. E. Synthetic data. We evaluate DIVOT with Gaussian noise as in (11) on the datasets with different (non)linear functions and samples sizes. We generate synthetic data with the ANMs: 1) Y = X + Ey; 2) Y = 0.1(2.5X)3 − 0.1X +Ey; 3) Y = sin(4X) +Ey; 4) if x < 0, Y = 0.5X3 −X +Ey and if x ≥ 0, Y = 1 − 0.5X3 + X + Ey , where X ∼ U(−1, 1) and Ey ∼ U(0, 1) are uniform distribution. For each ANM, we generate datasets with the sample sizes 10, 25, 50, 100, 200, and 500. With each sample size, 100 different datasets are generated for computing the accuracy.
Fig. 2a shows that DIVOT consistently recovers the causal direction for all the cases. Since DIVOT has no smoothness constraint of functions, it can deal with the case 4) which is a discontinuous
function. We compared DIVOT with the results of CAREFL (Khemakhem et al., 2021), RECI (Blöbaum et al., 2018), and other benchmark methods in Appendix. As shown in Fig. 7 of App. E, our method performs better than the others. Moreover, we also show the robustness of DIVOT to the prior misspecification as the experiments in CAREFL (Khemakhem et al., 2021). The synthetic data are generated with uniform distribution noise Ey ∼ U(0, 1), where the p(Ey; θ) of DIVOT is either uniform distribution U(0, 1), beta distribution B(a = 0.5, b = 0.5), or standard normal distribution N (0, 1). As shown in Fig. 2b, DIVOT with the misspecified noise distributions has similar performance with the one using the correct class of distributions.
Tübingen cause-effect pair dataset. We apply DIVOT to the Tübingen cause-effect pair dataset (Mooij et al., 2016). This is a collection of real-world cause-effect pairs. We use the variancebased divergence measure (13) and parameterize p(Ey; θ) as in (11). As for the PNL extension, we implement (14) without the positivity constraint (in practice imposing no constraint on ω also performs well in this simple formulation). Moreover, we found that given a value of X or Y of the datasets, there are often few samples. Thus we consider a range of values ofX , which may introduce bias of the divergence measure; hence we also use a linear debiasing function for reducing the bias. Then, we minimize the variance-based divergence measure over parameters with autograd of JAX. More details about the debiasing function and optimization can be found in App. E. To compare with the results reported in other works, we use the maximum number of the sample size N = 500 as (Stegle et al., 2010; Zhang et al., 2015a) and run all the experiments with 3 random seeds as (Stegle et al., 2010), though DIVOT is efficient enough for datasets with larger sample sizes as shown in Appendix. In addition, we normalize data and select the ones within 2 standard deviation, because DIVOT is sensitive to the outliers as optimal transport. As shown in Tab. 1, our proposed DIVOT, especially the extension for PNLs, outperforms than the other methods. The reported results in Tab. 1 are taken from (Stegle et al., 2010) and (Zhang et al., 2015a). And our results are based on the same datasets as them, i.e., with 68 and 77 cause-effect pairs respectively. Moreover, the entropic causal inference (Compton et al., 2021) is reported with 64.21% accuracy; CAREFL (Khemakhem et al., 2021) is reported with 73% accuracy on 108 pairs; RECI (Blöbaum et al., 2018) is reported with 76% weighted accuracy on the 100 pairs. From Tab. 1, we can also see that the other ANM/PNL-based methods are more sensitive to the choice of noise distributions. As shown in Tab. 1, the ANM/PNL with Gaussian noise performs worse than the one with a more complex distribution (mixture of Gaussian) or the one combined with a more stable measure, such as HSIC independence test. In contrast, the PNL extension of DIVOT with the Gaussian noise has a state-of-the-art result.

heading: 7 CONCLUSION
text: In this paper, we provide a new dynamical-system perspective of FCMs in the context of identifying causal relationships in the bivariate case. We first demonstrate the connection between FCMs and optimal transport, and study the dynamical systems of the optimal transport under constraints of FCMs. We then show that FCMs correspond to pressureless potential flows and that ANMs as a special case corresponding to the pressureless potential flows with the velocity field divergence equal to zero. Based on such findings, we propose a divergence measure-based criterion for causal discovery and provide an efficient optimal transport-based algorithm, DIVOT, for identifying causal direction between two variables. The experimental results on both synthetic and real datasets show that compared with the state-of-the-art methods, DIVOT has state-of-the-art results and promising properties for flexibility, efficiency, and prior misspecification robustness. We hope that the connection between FCMs and optimal transport has the potential of helping understand general FCMs from the dynamical perspective and inspiring more generic causal discovery methods.

heading: B EXTENSION TO THE MULTIVARIATE CASE
text: The problem of the naive solution to the multivariate case. A direct extension to the multivariate case is as the extensions of other bivariate causal discovery methods, such as (Khemakhem et al., 2021), (Zhang and Hyvärinen, 2009b), and (Monti et al., 2020). One can first apply constraint/score-based methods to get a causal skeleton, an undirected causal graph, and then use the extension of our measure for finding all the causal directions. A problem of the extension in some of the other bivariate works assuming causal sufficiency is that they directly applied their bivariate methods to each edge of the causal skeleton without considering the DAG structure. This can lead to the wrong results, especially in the case where they disregard the common parent/confounder of an edge.
Extending our method directly to the multivariate case. We provided the extensions of Prop. 1 and the variance-based measure, which can be used for computing the measure value considering the DAG structure and then orienting the edges based on the causal skeleton in the multivariate case. Our extension to the multivariate case has the following properties: Given the causal skeleton,
1. it distinguishes Markov equivalent classes under identifiability conditions of ANMs (or PNLs);
2. it still benefits from the closed-form 1D optimal transport solution, and doesn’t have the computational issue as in the high-dimensional optimal transport methods because of the FCM constraints;
3. it provides a score, of which the causal structure has the minimum value compared with all DAGs of the causal skeleton; moreover, the measure value of a DAG is the summation of all measure values of the causal modules/conditionals/mechanisms.
First, we show the derivation of the extension of Prop. 1 for ANMs, given which the derivation in the PNL case is straightforward (one can consider what we did for Eqn. (14) in our paper). The proposition shows how to efficiently compute theL2 Wasserstein distance between high dimensional distributions under FCM constraints. Suppose that the general ANM is Xi = gi(PAi) + Ei, where i = 1, ...,m, Ei is the noise term of Xi, and PAi denotes the parent variables of Xi. The square of L2 Wasserstein distance is
W 22 (p0, pT ) = m∑ i=1 EPAi [W 22 (p(Ei), p(Xi|PAi)],
of which the derivation is shown in App. D.2.
Second, as a direct implication of Thm. 1, given the couplings of Ei and Xi, the corresponding dynamical system has zero divergence of its velocity field; in other words, the corresponding dynamical system which moves samples of x0 to samples of xT under ANM constraints has zero divergence on each dimension. And the variance-based measure is
Dvar(v) ≈ m∑ i=1 1 N ∑ k∈{samples of PAi }N ||sort(−−−−−→xi|PAi=k)− sort(−→ei )− ave(−−−−−→xi|PAi=k −−→ei )||22 Nk − 1 ,
where −−−−−→xi|PAi=k is the data vector of which the elements are the Xi values of the samples with PAi taking the value k; −→ei is the data vector of which elements are the generated noise samples; Nk represents the length of the vector −−−−−→xi|PAi=k or −→ei . Next, as the direct extension, one could enumerate all possible DAGs of the causal skeleton, and compute their measure values, of which the minimum value is corresponding to the causal graph. Because the causal skeleton is given, it must be the case that one of the two variables of an edge is the cause and the other one is the effect. So the enumerated graphs have two situations: 1) all the edges are correctly oriented; 2) the causal direction of at least one edge is wrong such that the measure value of at least one causal module is significantly larger than the correct one (note that considering a child as the direct cause leads to increasing the measure value, while omitting a cause is not necessary to increase the measure value of the causal module). Therefore, we can simply choose the graph with the minimum measure value as the causal one. As mentioned in the paper, it is also very important to develop practical algorithms for large-scale real-world problems, and there are some points for future works to further explore: 1. Testing the significance of the results in the multivariate case, i.e., whether the minimum one is significantly smaller than the others. One could apply bootstrapping to the dataset and have a pvalue for the measure value. One could also apply bootstrapping for each causal module, however, there may exist the problem of the multiple statistical test issue with family-wise errors in this way. 2. Developing an efficient search algorithm in the multivariate case without relying on constraint/score-based methods. This requires to analysing the measure value of a causal module in more situations (e.g., omitting a parent, involving an independent variable/non-child descendant/non-parent ancestor, or a case mixing the mentioned factors), and then considering the characteristics to develop an efficient search algorithm similar as the greedy search algorithm (Chickering, 2002). NEW (GA77, pBTi)C UNKNOWN CONFOUNDING
The unknown confounding has different influences on the results of our method in different situations: 1) independent case: when two variables are independent; 2) causal case: when there is a causal relationship between two variables. Depending on how the unknown confounding influences the pair of variables, it can make our method
1. reverse the direction of the result in the causal case;
2. disregard the causal relationship in the causal case;
3. introduce extraneous causal relationship in the independent case.
To understand the results, we could analyze some toy examples intuitively. One extreme independent case is that X = U and Y = Ey + U , where U is the unobserved confounder; Ey is the noise of Y ; X has no noise. Then our method will show that X is the cause of Y , even though there is no causal relationship. Similarly, suppose that U and the noise of X and Y have the same distribution with variance equal to 1, and that X = Ex + 1000U and Y = Ey + U . Then, the distribution of X can be dominated by U , which leads to the wrong result for the same reason as the extreme case.
In the causal case, it follows the same reason. For example,X = Ex+U and Y = X+Ey+1000U . Then the distribution of Y can be dominated by U . Therefore, the FCM is close to X = Ex + U and Y = 1000U , which leads to reversing the causal direction in the result. Moreover, when X = Ex+1000U and Y = X+Ey+1000U , then both distributions ofX and Y can be dominated
by U ; hence, it is close to X = 1000U and Y = 1000U , which leads to disregarding the causal relationship in the result.
Nevertheless, in practice, the impact of unknown confounding is more complex and all factors can be mixed together with the impact of finite samples and the function properties. Therefore, we provide the experiments based on the synthetic datasets. From the experimental results, we can also find that the strength of the confounding and the difference of the confounding strength on the two variables are two important factors resulting in the wrong results of our method.
Experimental results for the unknown confounder. We generate a dataset with 1000 samples for each FCM: FCM1) X = U and Y = U ; FCM2) X = Ex + wx × U and Y = Ey + wy × U ; FCM3) X = Ex + wx × U and Y = f(X) + Ey + wy × U , where Ex and Ey are the noise terms; U is the unknown confounder; wx and wy are the coefficients of U representing the confounding strength; we used f(X) = X and f(X) = sin(4X) in the experiments. As for the data generation, Ex, Ey , and U follow the uniform distribution, U(0, 1), and we vary the coefficients wx and wy .
Instead of testing whether a result is significantly close to zero, we use bootstrapping to resample 50 datasets for each generated dataset and then apply a two-sample test to test whether the measure values in the two different directions are significantly different. Because although the measure value, in theory, is zero in the causal direction when there is no unknown confounder, with finite samples the measure value will be larger than zero in practice. And to decide which measure value is close to zero or not, it requires a threshold which can vary in different applications, or deriving a statistical test which is nontrivial for the measure without assuming the type of noise distributions and can be the future work of our method. Thus, we test whether two measure values are significantly different, and if so, we then conclude that the smaller one is in the causal direction; if they are not significantly different, we then conclude that it is the independent case. As for the experimental results, when the p-value is close to zero, it means that the two measure values are significantly different. And in general, one could take the significant level at 0.05 to make a conclusion. We used 0.05 for the experiments.
The experimental result of FCM1 is that the p-value is equal to 1.0, which means that the measure values in the two directions are not significantly different.
As for the experimental results of FCM2, we found that in the independent case, when the coefficients wx and wy are the same, we can get the correct results; when the coefficients are different, the method will give the wrong results. In different applications/scenarios, there are different tolerance ranges for our methods such that when the difference of the confounding coefficients is within the range, even if the confounding coefficients are different, we can still get the correct results.
As for the experimental results of FCM3, we found that in the causal case, the confounding strength is a factor influencing the results more, which is different from the independent case. Moreover, we found that when the nonlinear function f(X) is non-monotonic, such as sin(4X), the number of correct results is larger than the one in the linear case. Because the non-monotonic function itself can introduce a type of asymmetry which can indicate the causal direction for our method, the non-monotonic functions can be easier than the monotonic functions.
Future work for dealing with the unknown confounders. Causal discovery in the presence of unknown confounding is still an open problem, but there are some promising results in recent years such as (Janzing and Schölkopf, 2018; Mastakouri et al., 2021). Especially, (Mastakouri et al.,
2021) shows the identifiability results on time-series data in the presence of memoryless unknown confounders where the confounder doesn’t have an autocorrelation effect. A potential research direction with our framework is that by recovering the corresponding dynamical process of a static causal discovery problem while considering the memoryless confounding, we can determine the causal direction between two variables in the presence of memoryless unknown confounders.

heading: D PROOFS AND DERIVATION

heading: D.1 DERIVATION OF PROP. 1
text: Under constraints (I) and (II),
W 22 (p0, pT ) = ∫ |M∗(x0)− x0|2p0(x0)dx0 definition of W 22
= ∫∫
(f(Ex, Ey)− Ey)2p(Ey)dEyp(Ex)dEx Eqn. (8), and the independence of Ex and Ey = EEx [W 22 (p(Ey), p(Y |Ex))]
NEW (GA77)

heading: D.2 PROOF OF THE EXTENSION OF PROP. 1
text: Suppose that in the multivariate ANM, x0 = [E1, ..., Em]′, and xT = [X1, ..., Xm]′, and Xi = gi(PAi) + Ei, where PAi denotes the parent of Xi.
W 22 (p0, pT ) = ∫ |M∗(x0)− x0|2p0(x0)dx0
= ∫ |M∗([E1, .., Em]′)− [E1, .., Em]′|2p(E1, ..., Em)dx0
= m∑ i=1 ∫ |M∗i ([E1, .., Em]′)− Ei|2p(E1, ..., Em)dx0,
(where M∗i represents the i-th element of M ∗)
= m∑ i=1 ∫ |M∗i ([E1, .., Em]′)− Ei|2p(Ei)dEi p(E1:m−i)dE1:m−i,
(where 1 : m− i represents 1, ..., i− 1, i+ 1, ...,m, and the independence of noise implies p(Ei)p(E1:m−i))
= m∑ i=1 ∫ |f(Ei, E1:m−i)− Ei|2p(Ei)dEi p(E1:m−i)dE1:m−i
= m∑ i=1 ∫ |f̃(Ei,M∗1:m−i(E1:m−i))− Ei|2p(Ei)dEi p(E1:m−i)dE1:m−i
= m∑ i=1 ∫ |f̃(Ei, X1:m−i))− Ei|2p(Ei)dEi p(E1:m−i)dE1:m−i,
(where applying the change of variable formula)
= m∑ i=1 ∫ |f̃(Ei, X1:m−i))− Ei|2p(Ei)dEi p(X1:m−i)dX1:m−i|det(J)|−1
= m∑ i=1 ∫ |f̃(Ei, X1:m−i))− Ei|2p(Ei)dEi p(X1:m−i)dX1:m−i,
(where |det(J)|−1 = 1 for ANMs)
= m∑ i=1 EPAi [ ∫ |f̃(Ei, PAi))− Ei|2p(Ei)dEi]
= m∑ i=1 EPAi [W 22 (p(Ei), p(Xi|PAi)]

heading: D.3 PROOF OF THM. 1
text: Part I. We first derive the time interpolation p(t,xt) of p0 and pT under dynamical FCM constraints. Suppose that the FCMs are Y = f(X,Ey). As for the dynamical formulation of the L2 Wasserstein distance under constraints (I) and (II), according to the Jacobian equation, p(t,xt) = p0(x0)/|det(J ∂xt
∂x0 )|, t ∈ (0, T ),where xt,x0 ∈ R2 and J ∂xt ∂x0 is the Jacobian matrix w.r.t
x0. Moreover, ∂ log p(t,xt) ∂t = − ∂ ∂t | det(J ∂xt ∂x0 )|
| det(J ∂xt ∂x0
)| , where J ∂xt ∂x0 = [ 1 + tT , 0 t T ∂f ∂X , t T ( ∂f ∂Ey − 1) + 1 ] .
According to Jacobi’s formula and (Tao, 2013),
d dt |det(A(t))|
= ±det(A(t))tr ( A(t)−1 d
dt A(t) ) = ±det(A(t))tr ( det(A(t))−1adj(A(t)) d
dt A(t) ) = ±
∑ ij ( adj(A(t))′ d dt A(t) ) .
where is element wise multiplication and adj(A) is adjugate matrix of A and the sign takes the same sign as det(J ∂xt
∂x0
). Furthermore, since we know that under the structural constraints (I) and
(II), M∗ is Eqn. (8). Thus, we replace A(t) with J ∂xt ∂x0 = I + tT (∇M∗ − I), and then we have
J ∂xt ∂x0 = [ 1 + tT , 0 t T ∂f ∂X , t T ( ∂f ∂Ey − 1) + 1 ] ;
adj(J ∂xt ∂x0
)′ =
[ t T ( ∂f ∂Ey − 1) + 1, tT ∂f ∂X
0, 1 + tT
] ;
d dt J ∂xt ∂x0 = 1 T (∇M∗ − I) =
[ 0, 0
1 T ∂f ∂X , 1 T ( ∂f ∂Ey − 1)
] .
Therefore,
d dt |det(J ∂xt ∂x0 )| = ± 1 T ( ∂f ∂Ey − 1)(1 + t T ),
which takes the same sign as det(J ∂xt ∂x0 ).
Part II. We then proof the property of the corresponding dynamical systems of ANMs in Thm. 1.
Because ∂ log p(t,xt)∂t = − ∂ ∂t | det(J ∂xt ∂x0 )|
| det(J ∂xt ∂x0
)| , under the conditions in Thm. 1, ∂f ∂Ey = 1 and it is obvious
that ddt |det(J ∂xt ∂x0 )| = 0 and |det(J ∂xt ∂x0 )| 6= 0, ∀t ∈ [0, T ]. Furthermore, according to the theorem of instantaneous change of variables (Chen et al., 2018), which is a variant of Fokker-Plank equation, we know that
∂ log p(x(t))
∂t = −tr(∇v) = −divv.
Therefore, div v = 0.

heading: D.4 PROOF OF PROP. 2
text: Necessary direction. We prove the necessary direction by showing that given X is the direct cause of Y in an ANM, D(v) = 0. In the causal direction, because of the time evolution equation (7), we know that it is sufficient to check the divergence of the velocity field at time 0, and that the divergence taking zero value everywhere at time 0 implies the divergence taking zero value everywhere for t ∈ [0, T ]. Because of Thm. 1, we know the div v = 0 at time 0. Therefore, D(v) = 0.
Sufficient direction. We prove the sufficient direction by showing the contradiction with the identifiability of ANMs. Given that D(v) = 0 under the constraints (I), (II), (III), and the identifiability conditions of ANMs, it tells us that there is an ANM in the form of Eqn. (3) which is consistent with the data distribution and has independent noise. Because of the identifiability of ANMs shown
by (Hoyer et al., 2008), there is no ANM with Y → X which is consistent with the data distribution and has independent noise at the same time. Therefore, it can only be the case the X is the cause of Y .
Moreover, there is a concern that it may happen that in the direction which is not the causal direction, the model is not an ANM and its divergence measure is equal to zero. This may be problematic if the divergence measure is applied to applications. However, we will show that in general, this will not happen by proving that under weak assumptions D(v) = 0 implies the model is an ANM. First, D(v) = 0 implies div v = 0 for all x0 with positive probability densities p(x0) > 0. In the following we say div v = 0 in short. Second, div v = 0 implies ∂p(t,xt)∂t = 0 a.e. according to Thm. 1. Consequently, d| det(J ∂xt ∂x0 )|
dt =0, which leads to ∂f ∂Ey = 1 where Y = f(X,Ey) a.e. according to the proof of Thm. 1. Therefore, under the assumptions: 1) p(Ey) is positive in a continuous range of Ey; 2) in the range, ∂f∂Ey = 1 a.e. implies that it holds everywhere, we have Y = Ey + C, where C is a quantity which doesn’t change with Ey , e.g., it can be a function of X or a constant. Therefore, D(v) = 0 implies Y = Ey + g(X) further under assumptions 1) and 2).

heading: E EXPERIMENTS AND DETAILS OF DIVOT
text: For the synthetic data experiments in Sec. 6, we visualized the generated data with the sample size 500 in Fig. 3. For determining causal direction, we used the variance-based divergence measure in
Eqn. (13). The experiment setup is the same for experiments with different ANMs. For each ANM, we run experiments with different sample sizes, 10, 25, 50, 100, 200, and 500. With each sample size, 100 different datasets are generated for testing the accuracy of DIVOT. For all the synthetic data experiments in Sec. 6, we use batches for computing Eqn. (13) without using debiasing functions. In the following, we introduce positions (App. E.1), batches (App. E.2), debiasing functions (App. E.3), and optimization details (App. E.4) of DIVOT; as well as show its robustness to prior misspecification (App. E.5), more comparison with other benchmark methods (App. E.6), and its efficiency (App. E.7).

heading: E.1 NUMBER OF POSITIONS FOR COMPUTING THE VARIANCE-BASED DIVERGENCE MEASURE
text: Suppose a dataset {(xi, yi)}N with sample size N . For computing the divergence measure, we call x ∈ {xi}N as a position and for each x ∈ {xi}N we need to specify the corresponding vector −→yx of which the number of elements is Nx and generated samples of −→ey . In practice, it is not necessary to use all the N positions for computing Eqn. (13). For the synthetic data experiments in Sec. 6, whenN > 50, we choose 50 positions out ofN ; otherwise, we choose all theN positions. To select a position (when N > 50), we first find the maximal and minimal values of {xi}N , compute the interval length, ( max({xi}N )−min({xi}N ) ) /50, and then choose a position every such length (or its nearest position available in the data). We can then compute Eqn. (13) for causal direction determination.

heading: E.2 BATCHES FOR THE FINITE (OR LIMITED) SAMPLE SCENARIO
text: In the finite sample scenario, it can be the case that there is no other data at a position x. Therefore, we use the neighbors of x as a batch of data for computing Eqn. (13). We consider all the data in a batch as having the same value x. We represent the batch size with the percentage of the total numberN of samples. For the synthetic data experiments in Sec. 6, we use the batch size 0.4 for the datasets with sample size 10; 0.2 for the datasets with sample sizes 25 and 50; 0.15 for the datasets with sample sizes 100 and 200; and 0.05 for the datasets with sample size 500.
Note that for the synthetic data experiments in Sec. 6, the batch size is not larger than 0.2 except the extremely small datasets with 10 samples. Because using larger batch sizes introduces bias by concatenating data at different positions in a batch. For the experiments in Fig. 4a, we generate a synthetic dataset with sample size 100 and use all the 100 positions for computing Eqn. (13). We can see that choosing a proper batch size can increase the accuracy of DIVOT by increasing the number of samples in a position; however, if further increasing the batch size such that it is larger than 0.2, the accuracy is decreased because concatenating data at different positions leads to a different corresponding noise distribution compared with the ground-true p(Ey). Thus, as shown in Fig. 4b, the estimation of the noise distribution becomes worse and worse with increasing the batch size. Such bias can be problematic especially in the few-sample case where we have to choose a larger batch size for computing the divergence measure. Therefore, we use a debiasing function such that we can use a larger sample size for the few-sample case without sacrificing the performance of DIVOT.

heading: E.3 DEBIASING FUNCTIONS FOR THE FEW-SAMPLE SCENARIO
text: Because the mentioned bias is introduced by omitting the position information, we introduce it back with a debiasing function gd(x) to reduce the bias by modifying Eqn. (13) as
Dvar(v) ≈ 1
N ∑ x∈{xi}N
∥∥∥∥sort(−−−−→ydebaisx )− sort(−→ey)− ave(−−−−→ydebaisx −−→ey)∥∥∥∥2 2
Nx − 1 , eiy ∼ p(Ey; θ),
−−−−→ ydebaisx =
−→yx − gd(x;w), (15) where gd(x;w) is a (non)linear function parameterized with ω, e.g., in our experiments the linear debiasing function is gd(x) = w × x. As shown in Fig. 4, using the linear debiasing function guarantees not only the accuracy while using a larger batch size than 0.2, but also the correctness of the density estimation.
For a more complicated scenario, e.g., the ANM uses a nonlinear function, one can apply a neural network to the debiasing function. We found that as a by-product, a sufficient flexible debiasing function can be used for estimating g(x) of ANMs in Eqn. (3). As shown in Fig. 5, we generate synthetic 500 data of an ANM with the nonlinear function Y = g(X)+Ey = 0.1(2.5X)3−0.1X+ Ey , where X ∼ U(−1, 1) and Ey ∼ U(−0.5, 0.5) satisfy uniform distribution. The gd(x) is close to the ground-truth g(x). Because this is not the main focus of the work, we would like to refer the readers to the Jupyter notebook in supplementary materials for the details of the implementation.
Nevertheless, the purpose of using debiasing functions is not to estimate the g(x) accurately but to reduce the bias introduced by using a large batch size. This means that when we use a restrictive class of gd(x), it can lead to a noticeable difference of gd(x) and g(x). But for the purpose of causal direction determination, the divergence measure does not require to specify the nonlinear functional form of ANMs. Moreover, DIVOT based on optimal transport does not require to estimate g(x) either. Therefore, DIVOT can still determine causal direction correctly in such case.

heading: E.4 OPTIMIZATION METHODS AND THE CONVEXITY OF THE OBJECTIVE FUNCTIONS OF
text: DIVOT
In this section, we introduce the optimization method of the synthetic experiments in Sec. 6 and the convexity of the objective function. Next, we introduce the optimization method of DIVOT in the real-world data experiments in Sec. 6 of which the divergence measure uses a linear debiasing function and standard Gaussian noise.
As for the synthetic data experiments, we minimize the variance-based divergence measure in Eqn. (13) as an objective function of which the parameter is θ. Because when θ > 0, the linear function fnoiseθ (e source y ) does not change the order of sorting results of −−−−→ esourcey , we have
Dvar(v)
≈ 1 N ∑ x∈{xi}N
∥∥∥sort(−→yx)− fnoiseθ (sort(−−−−→esourcey ))− ave(−→yx − fnoiseθ (−−−−→esourcey ))∥∥∥2 2
Nx − 1 ,
= 1
N ∑ x∈{xi}N
∥∥∥sort(−→yx)− θ × sort(−−−−→esourcey )− ave(−→yx − θ ×−−−−→esourcey )∥∥∥2 2
Nx − 1 ,
fnoiseθ ( −−−−→ esourcey ) = θ × −−−−→ esourcey , θ > 0, e source y ∼ U(0, 1). (16)
It is obvious that the divergence measure in Eqn. (16) is convex on θ > 0. The convexity is achieved by parameterizing the noise distribution with a linear function fnoise, and there are many exited toolboxes for solving the minimization problem.
For the synthetic data experiments, we used gradient descend for finding the optimal θ∗. The gradient jax.grad(loss) is computed with the autograd in JAX (Bradbury et al., 2018). We update θ by specifying a step size sz and θ := θ−jax.grad(loss)×sz. If after the update θ < 0 (which has never happened), we set the value of θ as a positive number close to zero. We used sz = 1 for all the synthetic data experiments. Moreover, because of the convexity, one can also use an one-step update method to directly get the best parameter θ∗ by finding the root of the gradient function of Eqn. (16). We simply give a range of θ and find θ∗ by a binary search method such that the gradient at θ∗ is equal to zero. The experiments of DIVOT for ANMs in the real-world data experiments also used the one-step update for updating θ and the range θ is specified as [0, 100]. In addition, we need to optimize over the parameter w of the linear debiasing function. We use gradient descent to find the best parameter, w∗. We update θ every 10 updates of w. For the PNL extension of DIVOT, we used the cyclic learning rate (Smith, 2017) with gradient descent. We update θ every 10 updates of w and ω in (14). The program of DIVOT is terminated when the divergence measure converges. If a more complicated scenario requires fnoiseθ to be a nonlinear function, one may need to use the gradient descend method instead of the one-step update. We find that it is sufficient to use the standard Gaussian distribution with the linear fnoiseθ for DIVOT to have promising results in the experiments on the Tübingen datasets, which indicates that DIVOT is robust to the choice of models. We then investigate the robustness of DIVOT to prior misspecification.

heading: E.5 ROBUSTNESS TO PRIOR MISSPECIFICATION
text: We use different hypothesized noise distributions in DIVOT to test the robustness to the misspecification of noise distribution. The synthetic data are generated as in Sec. 6 with uniform distrbution noise Ey ∼ U(0, 1). In DIVOT, we use one of the three hypothesized distributions: uniform distribution U(0, 1), beta distribution B(a = 0.5, b = 0.5), and standard normal distribution N (0, 1). As shown in Fig. 6, DIVOT with the misspecified noise distributions has similar performance with the one using the correct class of distributions, which shows the robustness of DIVOT.

heading: E.6 COMPARISON WITH RESULTS OF BENCHMARK METHODS
text: We compare DIVOT with other benchmark methods, such as ANM (Hoyer et al., 2008), CAREFL (Khemakhem et al., 2021), RECI (Blöbaum et al., 2018), and LLR (Hyvärinen and Smith, 2013). As in (Khemakhem et al., 2021), the synthetic data are generated with the Laplace distribution as the noise distribution. DIVOT uses batches without debiasing functions. And the data preprocessing is the same as in the real-world data experiments in Sec. 6. We run experiments on each dataset for 100 times. As shown in Fig. 7, for the sample size larger than 100 the accuracy of DIVOT is 100% which shows that DIVOT has promising results and performs better compared the other methods, especially in the linear case.
performs competitively and is able to robustly identify the underlying causal direction.
distribution, while the noise distribution is chosen to be either a Laplace, Student-t or Gaussian distribution.
Table 1: Percentage of correct causal variables identi-
fied over 108 pairs from the Cause Effect Pairs benchmark.
CAREFL Linear LR ANM RECI
73 % 66% 69 % 69%
sented in Table 1. We note that the proposed method performs better than alternative algorithms.
Arrow of time on EEG data Finally, we consider the performance of CAREFL in inferring the arrow of time from open-access electroencephalogram (EEG) time series (Dornhege et al., 2004). The data consists of 118 EEG channels for one subject. We only consider the first n time points, where n 2 {150, 500}, after which each of the channels is randomly reversed. More details on the preprocessing can be found in Appendix F.3. The goal is to correctly infer whether xt ! xt+1 or xt+1 ! xt for each channel. This is a useful test case for causal methods since the true direction is known to be from the past to the future. We report in Figure 3 the accuracy as a function of the percentage of channels considered, sorted from highest to lowest confidence (i.e. by how high the amplitude of the output of each algorithm is). For average to high confidence, CAREFL is comparable in performance to the baseline methods,
but performs better in the low confidence regime. We also note that the performance of CAREFL improves by increasing the sample-size, which is to be expected from a method based on deep learning.
5.2 Interventions
To demonstrate that CAREFL can answer interventional queries, we will consider both a synthetic controlled example, as well as real fMRI data.
0 20 40 60 80 100 DeFLsLRn rDte
0.0
0.2
0.4
0.6
0.8
1.0
Pr RS Rr tLR n FR rr eF t
n 150
0 20 40 60 80 100 DeFLsLRn rDte
0.0
0.2
0.4
0.6
0.8
1.0
Pr RS Rr tLR n FR rr eF t
n 500
AlgRrLthP CA5EFL CA5EFL-16 LLneDr L5 5ECI A10
Figure 3: Performance on finding the arrow of time of
EEG data, as a function of decision rate (percentage of channels — sorted by decreasing confidence — we have to classify).

heading: E.7 EFFICIENCY OF DIVOT
text: For common causal discovery tasks in the bivariate case, the sample size 10000 is a large and challengeable one. Causal discovery methods need to consider the efficiency of algorithms especially in the large sample size scenario. We apply DIVOT with batches and no debiasing functions to the synthetic data of which the sample size is 10000 generated with Y = 0.1((2.5X)3 − X) + Ey and Ey ∼ U(0, 1). As shown in Table 7, we test the running time of DIVOT to determine causal direction with different number of positions and different batch sizes. The experiments are based on MacBook Pro (15-inch, 2018) with 2.9 GHz 6-Core Intel Core i9. Our implementation is based on JAX (Bradbury et al., 2018) which uses Apache License and the running time is measured with the command %timeit in JAX.
F IDENTIFIABILITY CONDITIONS OF ANMS IN (HOYER ET AL., 2008)
Because the identifiability conditions of ANMs in (Hoyer et al., 2008) are important and necessary for our Thm. 1, we include them here:
Let the joint probability density of x and y be given by p(x, y) = pn(y − f(x))px(x), where pn, px are probability densities on R. If there is a backward model of the same form, i.e., p(x, y) = pñ(x − g(y))py(y), then, denoting ν := log pn and ξ := log px, the triple (f, px, pn) must satisfy the following differential equation for all x, y with ν′′(y − f(x))f ′(x) 6= 0:
ξ′′′ = ξ′′(−ν ′′′f ′ ν′′ + f ′′ f ′ )− 2ν′′f ′′f ′ + ν′f ′′′ + ν ′ν′′′f ′′f ′ ν′′ − ν ′(f ′′)2 f ′ ,
where we have skipped the arguments y − f(x), x, and x for ν, ξ, and f and their derivatives, respectively. Moreover, if for a fixed pair (f, ν) there exists y ∈ R such that ν′′(y−f(x))f ′(x) 6= 0 for all but a countable set of points x ∈ R, the set of all px for which p has a backward model is contained in a 3-dimensional affine space.

Review summaries:
The paper tackles the problem of causal discovery in the basic case where a pair of variables is considered. In particular, it is concerned with Functional Causal Models and how to establlish the causal direction for a pair of continuous variables. It motivates the contibution by poiting out thet the performance of availabel algorithms is sensitive to the model assumptions, which makes it difficult for practitioners to use them in real settings. Then, the main contibution is a novel dynamical-system view of Functional Causal Models which aims to identify the causal direction in the case of pairs of variables. The paper exploits connection between functional causal models and optimal transport. Then, the paper studies the problem of optimal transport by taking into account the constraint originating from the considered functional causal model.

The authors frame the bivariate causal discovery problem in terms of the analysis of a dynamical-system. They use results from the field of optimal transport to interpret additive noise models from this framework. They also develop a novel criterion and a causal discovery algorithm based thereupon, and compare their results with previous methods used in bivariate causal discovery.