Optimal Transport for Causal Discovery
Authors: Ruibo Tu
Kun Zhang
Hedvig Kjellstrom
Cheng Zhang
Venue: ICLR 2022 Spotlight (2022)
Keywords: causal discovery
optimal transport
functional causal model

Abstract: To determine causal relationships between two variables, approaches based on Functional Causal Models (FCMs) have been proposed by properly restricting model classes; however, the performance is sensitive to the model assumptions, which makes it difficult 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. 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 for determining 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-nonlinear models. Our method demonstrated state-of-the-art results on both synthetic and causal discovery benchmark datasets.

Published as a conference paper at ICLR 2022
OPTIMAL TRANSPORT FOR CAUSAL DISCOVERY
Ruibo Tu
Kun Zhang
KTH Royal Institute of Technology ruiboêkth.se
Carnegie Mellon University
Mohamed bin Zayed University of Artificial Intelligence kunz @cmu. edu
Hedvig Kjellström
Cheng Zhang
KTH Royal Institute of Technology Silo AI
Microsoft Research
Cheng. Zhangemicrosoffi croso com
hedvigekth.se
ABSTRACT
To determine causal relationships between two variables, approaches based on Functional Causal Models (FCMs) have been proposed by properly restricting model classes; however, the performance is sensitive to the model assumptions, which makes difficult to use. In this paper, we provide novel dynamical-syss 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. provides new dimension for describing static causal discovery tasks while enjoying more freedom for modeling the quantitative causal influences. In particular, we we show that Additive Noise Models (ANMs) cor- respond poleme-preserving pressureless flows. Consequently, based their velocity field divergence, we introduce a criterion for determining causal direc- tion. With this criterion, we propose novel optimal transport-based algorithm for ANMs which is robust to the choice of models and extend it to post-nonlinear models. Our method demonstrated ofate-of-the-ard results on both synthetic and causal discovery benchmark datasets.
INTRODUCTION
Determining causal relationships between two variables is a fundamental and challenging causal discovery task (Janzing al. 2012). Conventional constraint-based and ccore-based causal discov- ery methods identify causal structures only up to Markov equivalent classes (Spirtes et al.) (000). in which some causal relationships are undetermined. To address this challenge, properly constrained functional çausal 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 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- nonlinonl models (PNLs) (Zhang and Hyvärinen, 20096). While some of the models, such as PNLs, are highly flexible, the constraints are still re- strictive and difficult to interpret and relax. Inevitably, the performance of these methods 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 time dimension for the static causal discovery task, interpret FCMs with dynamical causal processes under the least action principle (Amool d (203)). The new interpretation connects FCMs with 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.

---

Published as a conference paper at ICLR 2022
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 applied in large range of upplications, not only because it is natural way to describe moving particles (Ambrosio et al., 2012) but also because of its recent improvement in 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 al. 2018). Among different optimal transport definitions, the L Wasserstein distance got extensive applications in statistics (Rachev and Rüschendori [998), functional analysis (Barthe 1998), et al. (McCann) 1997, Otto 1997). The dynamical formulation of the L Wasserstein distance is introduced by [Benamou and Brenier (2000) for relaxing the computational costs. We find that 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-systen view, which makes 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 without learning a regression model. 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 (Arool d 2013) for the static bivariate causal discovery task, and characterize properties of the corresponding dynamical systems (See. 3.1 and Sec. 3.2
2. A criterion for determining causal relationships between two variables. We study the correspond- ing 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 for deter- mining causal relationships and show that under the identifiability conditions of ANMs it is valid criterion for ANMs. which can be extended to PNLs directly (Sec. 3.2.
An optimal transport-basee basse approach (DIVOT) for distinguishing cause from effect between two variables. DIVOT inherits the advantages of oimeennnnioonal optimal transport. It can be computed efficiently and does not require independence tests, learning a regression model, or deriving likeli- hood 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. and Sec. 6
PRELIMINARIES
Optimal transport: the underdetermined Jacobian problem. We mainly follow the notations and the definitions of (Benamou and Brenier (000).. Suppose that two (probability) density func- tions, Px(x) and px(x) where x E IR, are non-negative and bounded with total mass one. The transfer of (x()) to (x() is realized with smooth one-to-one map M Ed + Rd. The Jacobian problem is to find M that satisfies the Jacobian equation, (x(x)) = Pr (M det(V (xn))) where X M ). V is the gradient in vector calculus, and det(-) denotes determinant. This is an underdetermined problem as many maps can be the solutions. A natural way is to choose the optimal one, e.g., the one with the lowest cost, A common cost function is the T,2 Wasserstein distance.
LP Wasserstein distance and its one-dimensional closed-form solution. The LP Wasser- stein distance between Pu and pr, denoted by We(ps. Pr), is defined by W Wp(Pu,Pr) = inf' Mr I [xx)))) p (x))xxx) where p 2 (Kantorovich) 1948). In this work, mainly use the square of the _ Wasserstein distance, denoted by W... Moreover, the one-dimensional (1D) LP Wasserstein distance has a closed-form solution, e.g., the 1D optimal solution of W3 is M' PT o Po. where Po and P1 are the cumulative distribution functions for Pu and pr, and "o" represents the function composition. In practice, the ID optimal solution can be computed with the average square distance between the sorted samples from po and pr (Kolouri et al. 2019).
Functional causal models. FCMs represent the effect Y as a function f()) of the direct cause X and independent noise Ey. where function ff describes the causal influence of X on Y, and E, is the exogenous variable/noise. Without any additional assumption on the functional classes, the causal direction is not identifiable (Hyvärinen and Pajunen| 1999 Zhang et al. 2015a). Roughly speaking.
2

---

Published as a conference paper at ICLR 2022
(a) Pressureless potential flow corresponding to FCMs
(b) Volume-preserving pressureless potential flow correspouding to ANMs
Figure 1: Panel (a) illustrates the trajectories of the pressureless potential flow corresponding to gen- eral FCMs, while panel (b) shows the trajectories of the volume-preserving pressureless potential flow corresponding to ANMs. In panel (b), the trajectories have velocities the Ez/X axis and move along straight lines that are parallel to each other when having the same values of X/Eg.
because given variable pair (X,Y). one always construct Y f(X, Ey) and another different FCM. J(Y, Ex), such that both of them have independent "noise" (Hyvärinen and Pajunen 1999 Zhang et al.,, (015a). Several works further introduce proper assumptions on model classes, which guarantees that the independence of cause and noise only holds the causal direction, e.g..
ANM: Y = 9(X) + Ey;
(1)
PNL: Y = h(g(X)+ + Ey),
(2)
where g and h are nonlinear functions and is invertible.
DYNAMICAL INTERPRETATION OF FUNCTIONAL CAUSAL MODELS
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 deter- mination problem. then study the optimal transport under the constraints of FCMs, show the corresponding dynamical systems of FCMs, and characterize the properties of such systems.
3.1 A REFORMULATION OF FUNCTIONAL CAUSAL MODELS
As introduced in Sec. 2 FCMs are used to approximate the true data generation process. Given the FCM, Y [(X, Ey), rewrite it in the vector form,
x
Ee [(X,E,)) M
1 M(xo), E
Ex
X
(3)
where Xo, XT IR2, their probability densities Po. Pr O, and M - As analogy to the mass transfer scenario (Monge 1781), we consider the samples of independent noise E, and E, the particles of materials and regard the map M in Eqn. 3 a special transformation of the independent noise samples. As shown in Fig. [ one can consider the data points are transferred from the original positions (which are unmeasured) in the plane Ez-Ey at time the observed positions in the plane X-Y at time T. Such transformation considers the transfer as dynamical process which moves the unmeasured independent noise X = EEg E,, and consequently leads to [X, yy. From the perspective of FCM-hased causal discovery approaches,
the observations causal influences are represented by FCMs which represent the effect function of its direct cause and an unmeasured noise satisfying the FCM constraints:
The map constraint: the values of X are determined by the values its corresponding noise, i.e., x Ez, while the values of the effect depend cause X and noise Ey: (it) The independence constraint: the noise terms are independent, i.e, Ez is independent of E-y-
denotes the transpose of vectors matrices.
3

---

Published as a conference paper at ICLR 2022
Note that the optimal transport M with the minimal L" Wasserstein distance is not necessary to be the one in Eqn. 3 because it has no information about the FCM constraints or the true data generation process. In other words, given two sample sets of Xp and XT, the couplings given by optimal transport are not necessary to be the ones generated from the ground-truth FCM.
3.2 DYNAMICAL INTERPRETATION OF FCMs: OPTIMAL TRANSPORT UNDER THE FCM CONSTRAINTS
In this section, we jointly consider the causality and the 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 the 7,2 Wasserstein distance, study such dynamical systems under the FCM constraints, and then show their properties under the FCM and ANM constraints.
Dynamical L2 Wasserstein distance. [enaamou and Brenier (000) formulate the L Monge- Kantorovich problem as a convex space-time minimization problem in continuum mechanics framework. Fixing time interval [0, I]. they introduce the concepts of the smooth time-dependent density ptt, xx) 2 and the velocity field v(t, xt). When they are clear from context, we denote them by / and Because give the dynamical formulation of W3:
considering the bivariate case, R and RR? Then, they
inf T ( [ P.v. J22
W2
initial and final conditions: p(0,-) Po, ((T,-) pr the continuity equation: 0,p+ + V (pv) =
s.t.
where V. denotes the divergence in vector calculus. They show that minimizing the objective func- tion in the optimization problem 4 is equivalent to finding the dynamical system with the least action (Amoll d 2013) and prove that the solutions of ④ are pressureless potential flows, of which the fluid particles are not subject to any pressure force and the trajectories are determined given their initial positions and velocities or given their initial and final positions. Suppose that M is the solution given by W2. The corresponding flows follow the time evolution equation,
X Xo + 7 ((,,,,,, where vt,,,,) v(0, M Xo and & [0,2].
(5)
The time evolution equation shows that x, is a convex combination of xo and M (xx) and that the velocity fields do not depend on time.
As an analogy between the optimal transport problem and causal direction determination, the density p and the velocity 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 L Wasserstein distance introduces time variable and provides a natural time interpolation p(t, x²) of po and pr together with the velocity field v(t. xt). Similarly, we can also have the natural time interpolation p(t, x,) between Po and pr well as the velocity field (t,, under the least action principle Amol'd 2013), which is the dynamical interpretation of FCMs.
Dynamical 2 Wasserstein distance under the FCMs constraints. First, we introduce FCM constraints in the context of the dynamical L Wasserstein distance. According to the time evolution equation (5). we know that the velocity is fully determined by the initial and final values of x,- We first consider FCM constraint ( For the initial and final values of the cause, the value of its observation X equal to its noise value E. in Eqn. (3). Consequently, denoting v [May v,,l' where Un and Vy represent the velocities along the E -axis and -a,iis respectively as shown Fi.g [ we know that At E [0, T]. Zt E I, Uz (1, xr.) = O. As for FCM constraint (i). it implies that the initial noise of cause and effect, corresponding to E, and E, at the time 0, are independent. Therefore, we have the FCM constraints for the dynamical I,2 Wasserstein distance,
(I The map constraints: vlt, xx) = (t, x4), "y (1, v.).. Ver, Wh E RR, E [0. T), (1, xx)) = 0:
(II) The independence constraint: two random variables (cause and effect) at the initial time have the joint probability density function (xo) and they are independent.
4

---

Published as a conference paper at ICLR 2022
Second, characterize the properties of the dynamical L Wasserstein distance under constraints and (II. According to Eqn. B). the form of M is determined as
Ex Eny
(xo))
[s
(6)
otherwise, the FCM constraints will be violated. Moreover, the W under the FCM constraints can be computed with the one me-dimensional W3 as shown in Prop. [
(the derivation is in App- D
Proposition 1. Under constraints ①1 and MI the square of the L Wasserstein distance between po and pr is
W (Po.Pr) = EE, W3
Furthermore, we consider the constraints of ANMs and characterize the corresponding dynamical systems for now, which can be directly extended to PNLs mentioned in Sec. 4.2 Based on constraints I and M.. we further introduce the ANM constraint,
(III) the effect is the sum of noise and nonlinear function of cause as defined in Eqn.
Theorem Zero divergence of the velocity field ). Under constraints ① and (L), the dynamical systems given by the L Wasserstein distance are pressureless flows. Further under ANM constraint (III), they become celume-preserving
p of which the divergence of the velocity field, (t, xx) [va(t, 21), 2y yol,' satisfies
div vt,, (t,
(t, vt) O, Vt E [0, T), Fz. Ur I,
Ox,
By.
where div is the divergence operator in vector calculus.
Thm. 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 an essential property of the corresponding dynamical systems of ANMs. The property indicates a potential criterion for causal direction determination, which the divergence of the velocity field is zero everywhere in the causal direction, while it may not always hold in the reverse direction. Next, will verify the criterion rigorously, propose an algorithm based on it, and show extension for the PNL cases.
CAUSAL DIRECTION DETERMINATION WITH OPTIMAL TRANSPORT
In this section, we define a divergence measure as a criterion for determining causal direction be- tween two variables for ANMs. Based on the criterion, then provide algorithm identify the causal direction, named by the divergence measure with optimal transport (DIVOT).
4.1 DIVERGENCE MEASURE AS A CAUSAL DISCOVERY CRITERION
We first define the divergence measure and then show that it is valid criterion for identifying the causal direction in the bivariate case under the identifiability conditions of ANMs, as shown in Prop.
2 (the proof App. D.
Proposition 2 (Divergence measure causal discovery criterion). Define the divergence measure,
D(v) = / Idiv v/ Po Exal [div v v
where v = (Xo) Xg. Suppose that constraints
(III), and the identifiability conditions
of ANMs (Hoyer al. 2008) are satisfied. The divergence measure of the corresponding dynamical system satisfies D(v) 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 two-dimensional optimal transport problem and compute the derivative of a velocity field for all samples. Another challenge of computing M is that we in general have no information about Xg- Such issues all solved with DIVOT.
5

---

Published as a conference paper at ICLR 2022
4.2 PROPOSED METHOD: DIVOT
We first provide an overview of the algorithm, DIVOT, for determining the causal direction between two variables and then introduce the four steps compute the divergence measure.
Overview of DIVOT. Alg. I is based on the Livergence-measury criterion for causal direction determination 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 in both directions. In practice, noise can have different variance in different applications, and the larger variance can lead to the larger measure value with finite samples. So we normalize the variance-basee measure value with the estimated noise variance. And then the one with the smaller normalized measure value is the causal direction. In App. A we provide the modified algorithms for including the independent case and the significance of the results with bootstrapping.
Algorithm 1: DIVOT: divergence measure with optimal transport for causal direction determi-
nation.
Input: data {(ri, ve) and noise distribution p(E; (b)
6 if Ly- < 4x-y then return Y + X
Output: X + Y or Y + X.
8 else
Def (1 (1(x;; y:) }, p(E; 0o)):
return X Y
data, A= {(13 y:)) 0g
Def Div (data, pl (E;
Lx_ . Div(data, 0));
data, 0= = (vr,, 0o
11
noise Sampling(p( E; 0))
data", noise" = OT(aaa, noise)
LY x Div(data, p( (E;
13 loss ming VarDiv(data", noise*) return loss/Var(p(
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 the time evolution equation 5 It requires the couplings of the data of x [Eg, Eyl' and xy [X, YJ'. But in the bivariate causal discovery task, only the data of XT are given. Therefore, as shown in Line 11 of Alg. [ we first deal with the issue due to the lack of the noise data of xg, denoted by ((e,,,;)) To obtain the noise data, we may assume multivariate probability distribution of x. with the density Po (Eg, Ey) and then sample data from it, represented by (e egy) ~ Po (Faz, Fyy) Fortunately, due to the FCM constraints, we know that Po E,,, Fy) = PX, E.y) = p(X)p((,,) (.y)) So we only need to assume the probability distribution of E, and parameterize it with 0, denoted by p(Ey; 0). Suppose that the dataset of Xr with N samples is given, denoted by M.). We first sample a data set of E, with the sample size N, denoted by fei,a,, e.g., in the experiments of this work, we use the simplified reparameterization trick,
Igoise (ey (egource)= = 0 ey gouuce and ey gourcee ~ NoO0 1)/2((0,1), 1),
(9)
where fgoise is a monotonic function, and eoouree is sampled from a standard normal distribution or a uniform distribution. As the monotonic fhaise can be very flexible, we can represent flexible noise distributions with monotonic neural networks in (Huang et al., 2018). Next, we randomly match the data of X, denoted by (z;) N. with {ee) N. which gives the { e n as the data of xn-
Com- Optimal transport finds the couplings of the observation and the generated noise.
puting the divergence measure requires the couplings of the data of Xo and XT because of v = M² Xo; in other words, given the data of x and x7, we need to solve a two-dimensional optimal transport problem, which gives M and the couplings as in Line 12 of Alg. [ According to Prop. [ solving the wwo-dimensional optimal transport problem under the FCM constraints equiv- alent to solving one- -icddimnsio optimal transport problems. More specifically, the expectation in Eqn. (7) can computed with the Monte Carlo estimator by using N samples of F.,, denoted by (ez N. and then the two-dimensional W? in Eqn. (7) is computed with N one-dimensional W?. ie., W3 (Po. & * E, (p( (p( Ey), p(Y E. = )) Importantly, solving the one-dimensional optimal transport problem or computing the miedenenniooal W3 can be implemented with a sorting operation in (Kolouri et al., 2019). Therefore, the couplings are found by which given each value of E,, we find the corresponding values of E, and Y from {(e), et e and ((x, yi) hy and then match such sorted values of F.y and Y

---

Published as a conference paper at ICLR 2022
Variance-based divergence measure. Suppose that the previous step has found the couplings, i.e,, given any sample (ea, f.y), we have its corresponding (zt, yr). Then, we will compute the divergence of a velocity field as in Line 13 of Alg. According to the definition, |div v22 all the ey are matched with Vi. we have oy yi dy- A straightforward way to approximate the derivative is using its nearest neighbour pair (cao: Cyb) and (abb, yb:) and then approximate it
dy )² Since
with d (wy))o=u, [w--e)- however, it suffers the followiing issues (especially in the dy 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) + F.g at X = x is large. Therefore, we propose the variance- baaedd divergence measure working better in practice. It is straightforward to see that under constraints ! 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
y.-
(10) where V(V, x [X] represents conditional variance of the velocity field Vy (a random variable) at position X at the initial time. For example, V/Vy|=/X z| represents the variance of all the velocities By at the positions where X = at the initial time. { Bs) denotes the cause-effect pair dataset with the sample size N. for hN, denotes all the velocities vi at position X at the
D Dvar (v) = Ex [V[V [X]],
initial time, where the sample size is Nz, and their mean value is Vyle- Then, y[Vy|æ/X x]
No (oy) / - 1). and
(ort(U..) -sootssossss soot(E,) - ave - -ct)
Dvar (v) 20 N
1
Σ
(11)
N₂
=E(z:)N
where { represents the set of all the values of X and its sample size is N; We is the vector of the Y samples where X = x;; Ey is the vector of the Ey samples where Eg = x; sort(.) sorts a vector;
ave(-) computes the vector mean; and 11- is the square of a & norm.
Minimization w.r.t 0. Given the data {(ra, yi) y:) n and the density p Eyy 0), we can compute the divergence measure with the variance- method. Note that we only initialize 0 with some ran- dom value, and p( 0) 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 13 of Alg. [. we minimize the divergence measure w.r.t. O. For the minimization, one can derive the gradient w.r.t 9 in simple parameterization case as ( while in the complex case one can use mito-differentiationn. According to Prop. _ the divergence measure in the causal direction is zero if and only if p( Ey; 0*) with the optimal parameter 0 is the true noise distribution a.e., implied by the identifiability of ANMs. In this work, we used aut r ad and RMSP rop (gradient descent) of JAX for the minimization.
Extension to PNLs. The measure II can be directly extended to the PNL cases. Since in PNL 2 is an invertible function, by considering / (YY) as new random variable, h- (Y) xXX)+ +E is an ANM. Thus, for PNLs, under the identifiabilityy conditions of PNLs (Zhang and Hyvärinen 2009b), Prop. still holds; and we only need replace V in (TT) with J. an invertible function, e.g., it can be the simplified version of (Zhang et al. 2015a),
=
fPLL (Tz). where fENL is
J= = (TN)((T)) = Th + wa tanh(wp Vz + we),
(12)
where Was Wh, and win are positive scalars. Moreover, fE8L can also be the monotonic neural net- works, such as (Huang et al. 2018).
Extension to the multivariate case. For simplicity and clarity of the paper, we focus on elabo- rating 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-basee 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 Hyvärinen] 2009b; Monti et al., 20200 Khemakhem et al. 2021). See App. B for details.
7

---

Published as a conference paper at ICLR 2022
0.1((2.5x) +
Y [aore E
1.4
1.0
1.0
0.4
0.4
y- 0.112.5xx' -0.1x+
0.7
0.5
y sintdkld +,
Uniform
Uniform
y-x+ 6
0.6
6.6
Normal
Normal
100
202
551
200
200
400
Sample
Sample size
Sample
(a) Accuracy different ANM datasets.
(b) Robustness to prior misspecification.
Figure 2: Performance on synthetic data generated from different ANMs with different sample sizes.
fpiece is discontinuous function: fpicce (z) 0.5.r I, S O: fpicco()) 1 0.5.r I, O. In panel (a), each result represents the percentage of correct results on datasets with sample size generated from an ANM. Panel (b) shows the results of applying DIVOT using different types of distributions, p( Ey: 8), to the datasets with uniform distribution noise E, ~ (...
5 RELATED WORK
There are mainly two types of causal discovery methods for static causal direction determination between two variables. The first one introduces model assumptions to achieve the identifiability of causal direction. Most of such methods are based on LiNGAM (Shimizu al. 2006) and ANMs (Hoyer et al.] 2008 Mooij et al. 2009). Some of them are based on the more general models, e.g., PNLs using MLP for representing nonlinear functions (PNL-MLP) (Zhang and Hyvärinen, 2009b) and PNLs using warped Gaussian process and mixture of Gaussian noise (PNL-WGP-MoG) Zhang et al., 2015a). Recently, Khhemakhem et al. (021) 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 in (T) (or together with h in 2). and then test the independence between the independent noise E, in (I) (the residual) and the Schhmdd independence criterion (HSIC) (Gretton 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 approxima- tion. Moreover, there are other criteria proposed in the first type of methods, such as the likelihood criterion (MML). regression error-based causal inference (RECI) (Blöbaum et al., 20188,, and mutual information (Zhang and Hyvärinen, 2009a; Yamada and Sugiyamal (010). Nevertheless, the common issue of the first type of methods is that they restrict the model classes and that their performance sensitive to model assumptions. The second type of methods achieves identifiability by proposing other principles instead of restricting model classes, such as exogeneity (Zhang et al. (2015)), the randomness of exogeneous variables (Entropic) (Compton et al. 2021 and independent causal mechanisms, e.g., GPI (Stegle et al. 2010) and IGCI (Janzing el al.
cause X. And the independence test is commonly
gatio (LLR) (Hyvärinen and Smith 2013), maximum likelihood-baaee
2012
6 EXPERIMENTS
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 in the presence of unknown confounding in App. C The details about experiments are in App. F
Synthetic data. We evaluate DIVOT with Gaussian noise as in ( 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 < 0, = 3 X + E, and if 2 O,Y = - 0.5X3 + X E, where ~ 1(-1,1) and E, ~ 10 uniform distribution. For each ANM, we generate datasets with the sample sizes 10, 25, 50. 100, 200, and 500. For each sample size, 100 different datasets are generated. Fig. 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 a discontinuous function. We compared DIVOT with the results of CAREFL (Keemkhhem al. 2021). RECI (Blöbaum al. 2018). and
8

---

Published as a conference paper at ICLR 2022
Table 1: Percentage (%) of recovering the true causal direction on the Tübingen cause-effect pair datasets Mooij et al. 2016) compared with the reported results of (Stegle et al. 2010) in the upper table and
Zhang et al.
2013a) in the lower table.
Ours (PNL)
Ours (ANM)
LINGAM
ANM-Gauss
ANM-MML
ANM-HSIC
PNL
GFI-HSIC
GPIMML
IGCT
78.923.9
7124
6223
1523
GNET
0823
6724
6244
7242
7621
Ours (PNL)
Ours (ANM)
ANM
WGP-Gaust
WGP-MoG
PNL-MLP
GPI
IGCI
76.6+4.5
67#3
63
67
73
70
72
73
other benchmark methods in Appendix. As shown in Fig. 7f App. F our method performs better than the others. Moreover, we also show the robustness of DIVOT to the prior misspecification the experiments in CAREFI. (Khemakhem et al. 202T). The synthetic data are generated with uniform distribution noise E, (0 1), where the pl Eyi 0) of DIVOT is either uniform distribution (1 1), beta distribution B(a = 0.5, b = 0.5), or standard normal distribution (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.
Tübingen cause-effect pair dataset. We apply DIVOT to the Tübingen cause-effect pair dataset [Mooij et al. 2016). This is a collection of real-world cause-effect pairs. We use the variance- based divergence measure [1]] and parameterize p( F,y; O) as in (9). As for the PNI. extension, we implement (12) without the positivity constraint (in practice imposing no constraint on w also per- forms 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 range of values of X, which may intro- duce bias of the divergence measure; hence we also use linear debiasing function for reducing the bias. Then, we minimize the variance- of JAX. More details about the debiasing function and optimization can be found in App. F To compare with the results reported in other works, we use the maximum number of the sample size N = 500 as (Stegle et al. [2100] 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. deviation, because DIVOT is sensitive the outliers as optimal transport. As shown in Tab. T our proposed DIVOT, especially the extension for PNLs, outperforms than the other methods. The reported results in Tab. 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 [Keemahhem et al.) [221 is reported with 73% accuracy on 108 pairs; RECI (Blöbaum
divergence measure over parameters with autogra ouoorad
addition, we normalize data and select the ones within 2 standard
jtt al. 2018) is reported with 76% weighted accuracy on the 100 pairs. From Tab. D we can also see that the other A ANM/PNL- -based methods are more sensitive to the choice of noise distributions. As shown Tab. [ the ANM/PNL with Gaussian noise performs worse than the one with a more complex distribution (mixture of Gaussian) or the one combined with more stable measure, such as HSIC independence test. In contrast, the PNL extension of DIVOT with the Gaussian noise has a state-of-the ofrtt--artttttte result.
7 CONCLUSION
In this paper, we provide new dynamical-system perspective of FCMs in the context of identify- ing 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 divergence measure-basee 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-an methods, DIVOT has ofate-of-the-an art results and promising properties for flexibility, efficiency, and prior misspecification robustness. We hope that the con- nection between FCMs and optimal transport has the potential of helping understand general FCMs from the dynamical perspective and inspiring more generic causal discovery methods.
9

---

Published as a conference paper at ICLR 2022
ACKNOWLEDGEMENTS
KZ would like to acknowledge the support by the National Institutes Health (NIH) under Contract RO1HL by the NSF-Convergence Accelerator Track-D award #2134901, and by the United States Air Force under Contract No. support of the Swedish e-Science Research Centre.
RT would like to acknowledge the funding
REFERENCES
L. Ambrosio, M. Colombo, G. De Philippis, and A. Figalli. Existence of eulerian solutions to the semigeostrophic equations in physical space: the 2-dimensional periodic case. Communications in Partial Differential Equations, 37(12):229-22277, 2012.
M. Arjovsky, S. Chintala, and L. Bottou. Wasserstein generative adversarial networks. In Proceed- ings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pages 214-223. PMLR, 06-11 Aug 2017.
V. I. Arnol'd. Mathematical methods of classical mechanics, volume 60. Springer Science & Busi- ness Media, 2013.
F. Barthe. Optimal young's inequality and its converse: a simple proof. Analysis GAFA, 8(2)22442422, 1998.
Geometric & Functional
J.-D. Benamou and Y. Brenier. A computational fluid mechanics solution the monge-kantorovich mass transfer problem. Numerische Mathematik, 84(3):375-393, 2000.
Blöbaum, D. Janzing, T. Washio, S. Shimizu, and B. Schölkopf. Cause-effect inference by com- paring regression errors. In International Conference on Artificial Intelligence and Statistics, pages 900-909. PMLR, 2018.
J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Miine and Q. Zhang. JAX: composable transformations of Python+NumPy programs, 2018. URL http::
com/ com//cogle/jjat
R. T. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud. Neural ordinary differential equations. arXiv preprini (tf::::80000777 2018
D. M. Chickering. Optimal structure identification with greedy search. Journal of machine learning research, 50755550 2002.
S. Compton, M. Kocaoglu, K. Greenewald, and D. Katz. Entropic causal inference: Identifiability and finite sample results. arXiv preprint arXiv:2101.03500 2021.
M. Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. In Advances in Neural Information Processing Systems, volume 26. Curran Associates, Inc., 2013.
A. Genevay, G. Peyré, and M. Cuturi. Learning generative models with sinkhorn divergences. In International Conference on Artificial Intelligence and Statistics, pages 1608-1617. PMLR, 2018.
A. Gretton, O. Bousquet, A. Smola, and B. Schölkopf. Measuring statistical dependence with hilbert-schmiid norms. In International conference on algorithmic learning theory, pages 63- 77. Springer, 2005.
Hoyer, D. Janzing, M. Mooij, J. Peters, and B. Schölkopf. Nonlinear causal discovery with additive noise models. Advances in neural information processing systems, 21 21:689-696 2008
C.W. Huang, D. Krueger, A. Lacoste, and A. Courville. Neural autoregressive flows. Interna- tional Conference on Machine Learning, pages 2078-2087. PMLR, 2018.
A. Hyvärinen and P. Pajunen. Nonlinear independent component analysis: Existence and uniqueness results. Neural networks, 12(3):429 4399, 1999.
A. Hyvärinen and S. M. Smith. Pairwise likelihood ratios for estimation of non- non-gaussian structural equation models. Journal of Machine Learning Research, 14(Jan): 11-152, 2013.
10

---

Published as a conference paper at ICLR 2022
D. Janzing and B. Schölkopf. Detecting non-causal artifacts in multivariate linear regression models. In Interaational Conference on Machine Learning, pages 1245-22555 PMLR, 2018.
D. Janzing, J. Mooij, K. Zhang, J. Lemeire, J. Zscheischler, P. Daniušis, B. Steudel, and B. Schölkopf. Infornation-geommrric approach to inferring causal directions. Artificial Intel- ligence, 182:1-31, 2012.
L. V. Kantorovich. On a problem of monge. Uspekhi Mat. Nauk., 1948.
1. Khemakhem, R. Monti, R. Leech, and A. Hyvarinen. Causal autoregressive flows. In International Conference on Artificial Intelligence and Statistics, pages 5520-3528, PMLR, 2021.
Kolouri, G. K. Rohde, and H. Hoffmann. Sliced wasserstein distance for learning gaussian mix- ture models. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 3427-3436, 2018.
S. Kolouri, K. Nadjahi, U. Simsekli, R. Badeau, and G. Rohde. Generalized sliced wasserstein distances, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
A. A. Mastakouri, B. Schölkopf, and D. Janzing. Necessary and sufficient conditions for causal fea- ture selection in time series with latent common causes. In International Conference Machine Learning, pages 7502-7511. PMLR, 2021.
R. McCann. A convexity principle for interacting gases. Advances in mathematics, 128(1):153- 179, 1997.
G. Monge. Mémoire sur théorie des déblais et des remblais. Histoire de l'Académie Royale des Sciences de Paris, 1781.
R. Monti, K. Zhang, and Hyvärinen. Causal discovery with general non-linear relationships using non-linear ica. In Uncertainty in Artificial Intelligence, pages 186-195. PMLR, 2020.
Mooij, D. Janzing, Peters, and B. Schölkopf. Regression by dependence minimization and its application to causal inference in additive noise models. In Proceedings of the 26th annual international conference on machine learning, pages 745-752, 2009.
J. M. Mooij, J. Peters, D. Janzing, J. Zscheischler, and B. Schölkopf. Distinguishing cause from effect using observational data: methods and benchmarks. The Journal of Machine Learning Research, 17(1):1103-125 1204, 2016.
Otto. Viscous fingering: an optimal bound on the growth rate the mixing zone. SIAM journal on applied mathematics, 57(4):982-990, 1997.
S. T. Rachev and L. Rüschendorf. Mass Transportation Problems: Volume I: Theory, volume 1. Springer Science & Business Media, 1998.
S. Shimizu, P. O. Hoyer, A. Hyvärinen, A. Kerminen, and M. Jordan. A linear non- non- gaussian acyclic model for causal discovery. Journal of Machine Learning Research, 7(10), 2006.
L. N. Smith. Cyclical learning rates for training neural networks. In 2017 IEEE winter conference on applications of computer vision (WACV), pages 464-472. IEEE, 2017.
P. Spirtes, C. Glymour, and R. Scheines. Causation, prediction, and search, 2nd edn., vol. 1, 2001.
O. Stegle, D. Janzing, Zhang, M. Mooij, and B. Schölkopf. Probabilistic latent variable models for distinguishing between cause and effect. Advances in neural information processing systems, 23:1687-1695, 2010.
T. Tаo. 2013.
Matrix URL
identities
derivatives of determinant identities,
www.aro.oooddersr...oo http://wwess..eseeeeeeeeeeeeee/e wwtp://1////////111/333333 mntrix-idente i tieesiissseee ote.rt/eeeeteeeeeeeeeeeeoooooene htt.:mmmiitt---dedeitiiiee

---

Published as a conference paper at ICLR 2022
M. Yamada and M. Sugiyama. Dependence minimizing regression with model selection for non- linear causal inference under non- non-gaus gausgasssian noise. Proceedings of the AAAI Conference on Artificial Intelligence, volume 24, 2010.
K. Zhang and A. Hyvärinen. Causality discovery with additive disturbances: An information- theoretical perspective. In Joint European Conference on Machine Learning and Knowledge Discovery Databases, pages 570-585. Springer, 2009a.
K. Zhang and Hyvärinen. On the identifiability of the post-nonlinear causal model. In Pro- ceedings of the Twenty-Fifth Conference Uncertainty in Artificial Intelligence, pages 647-655, 2009b.
K. Zhang, Z. Wang, J. Zhang, and B. Schölkopf. On estimation of functional causal models: gen- eral results and application the post-nonlinear causal model. ACM Transactions Intelligent Systems and Technology (TIST), 7(2):1-22, 2015a.
K. Zhang, J. Zhang, and B. Schölkopf. Distinguishing cause from effect based on exogeneity. arXiv preprint arviv: 504 0, 2015b.
/504:
12

---

and the Published as a conference paper at ICLR 2022
APPENDIX
In App. A provide the modification of Alg. [ for including the independent significance of the results. In App. we provide the outlook of the extension to the multivariate case. In App. we provide the analysis, the discussion, and the experimental results of our method in the presence of unknown confounding. In App. D we provide the proofs of Prop. [ the extension of Prop. [I] Prop. 2 and Thm. In App. E we include the identifiability conditions of ANMs in (Hoyer et al. 2008 In App. we introduce the details of the DIVOT implementation and the
experiments:
App- introduce the concept, position, for computing the divergence measure
App. show a potential problem of computing the variance- based divergence measure case and introduce batches of data to deal with the problem; App. F3 show that using batches can introduce bias and lead wrong causal relationships
in the finttelimmieddsampee
(especially in the few-sample case) and then introduce debiasing functions reduce the bias, which used for the experiments on real-world data Sec. 6
App. F44 introduce the optimization methods used in the experiments synthetic and
real-world data in Sec. 6 and show the convexity of the objective function, the divergence
measure D. under the construction in Sec. [6]
App. F.5
show the robustness of DIVOT to prior misspecification;
App. 166
compare DIVOT with other benchmark methods;
App.1
F.7
show the efficiency of DIVOT with its running time.
A MODIFICATION OF ALG. FOR INCLUDING THE INDEPENDENT CASE AND THE SIGNIFICANCE OF THE RESULTS
The proposed measure is able to deal with the independent case without relying on other methods or tests after adapting the output conditions of the algorithm. The values of proposed measure in the independent case are zero in the two directions, while the measure values of the causal case are a zero value in direction and non-zero value in another direction. The modified output conditions are shown in Tab. 2
To determine that the value of our measure is zero or not, one could follow a
Table 2: Modified output conditions for the independent case.
Output of DIVOT
Div(X =
Div(X < Y) =
X, Y independent
True
True
X-+
True
False
X+-Y
False
True
similar way as (Blöbaum et af. 2018) choosing a threshold for real-world applications practically. However, different applications may have different thresholds. For example, in practice, when the noise distributions have different variance in the finite/few sample scenario, the larger variance can lead to the larger measure value. Although one can handle the problem well by normalizing the variance-based measure value with the estimated noise variance, deriving a statistical test for the finite/few sample case is still the ideal way, which is another nontrivial task without assuming the type of noise distributions and will be the future work of our method.
Bootstrapping for the significance of the results. Instead of testing whether a measure value is significantly zero or non ner-zero, we suggest using a bootstrapping method and testing whether the two measure values are significantly different (one could also test whether the difference between the two measure values are significantly zero or not):
Step 1. use bootstrapping (resampling with replacement) to get (e.g.. 50) bootstrapping datasets; Step 2. compute two measure values in the two directions for each bootstrapping dataset;
Step 3. apply two-sample test (T-test) to check whether the mean of one measure value in one direction is significantly different from another one in the other direction.
Step 4. if two measure values are significantly different, we conclude that the smaller one is in the causal direction; if they not significantly different, conclude that it is the independent case.
13

---

Published as a conference paper at ICLR 2022
Moreover, did experiments for the where the causal relationship is weak that the FCM with a causal relationship is similar to the independent case (e.g. when the coefficient of the direct cause is extremely small in the linear case, it is close to the independent case where the measure values are close). We generate four datasets with 1000 samples for the four FCMs in the synthetic data experiments in Sec. 6 with a weighting factor w: Y
f(X) F.y. and then use the bootstrapping method and T-test, Tab. 3 shows the p-values, of which small values indicate the sig- nificantly different measure values. If the significant level is 0.05, our method can tell the difference between two measures in the causal case when > 0.02.
Table Experimental results of testing the significance of the results.
w
0.01
0.02
0.03
0.04
0.05
M1
0.155
0.038
0,0006
2.429e-08
3.091e-10 11
M2
0.251
0.070
0.0003
0.0001
1.935e-07
M3
0.390
0.119
0.0196
0.0003
9.578e-08
M4
0.408
0.098
0.0037
2.57
2.346e-11
W
0.01
0.02
0.03
0.04
0.05
MI
indep
X cause
X cause
X cause
X cause
M2
indep
indep
X cause
X cause
X cause
M3
indep indep
X cause
X cause
X cause
M4
indep indep
X cause
X cause
X cause
B EXTENSION TO THE MULTIVARIATE CASE
The problem of the naive solution to the multivariate case. A direct extension to the multi- variate case as the extensions of other bivariate causal discovery methods, such as (Khemakhem et al. 202T). (Zhang and Hyvärinen| 2009b), and (Monti et al. 2020). One can first apply constraint- /score-based methods to get 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 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- 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 the 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 modues/cnnditoollseeeaiissss
First, show the derivation of the extension of Prop. for ANMs, given which the derivation in the PNL case is straightforward (one can consider what we did for Eqn. (12) in our paper). The proposition shows how efficiently compute the 2 Wasserstein distance between high dimensional distributions under FCM constraints. Suppose that the general ANM is X; gi( (P.A.) + E,, where i i= I, m, E, is the noise term of X,, and PA; denotes the parent variables of X,, The square of L Wasserstein distance is
W (P0. PT) = Σ EEA
14

---

Published as a conference paper at ICLR 2022
of which the derivation is shown in App. D.2
Second, as a direct implication of Thm. | given the couplings of E, and X,, the corresponding dynamical system has zero divergence of its velocity field; in other words, the corresponding dy- namical system which moves the samples of X to the samples of Xy under ANM constraints has zero divergence on each dimension. And the variance-based measure is
||sort [.IPAA(=))
sort( et) ave avo(T:)PAA 2
Dvar (v)
Σ
N RE (samples of PA; )
where T. is the data vector, of which the elements are the X; values of the samples with P.A, taking the value k; &j is the data vector, of which the elements are the generated noise samples; Nx represents the length of the vector TiJPA;=k or o
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 where one of the two variables of an edge 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 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 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 p- value 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 this way, 2. Developing an efficient search algorithm in the multivariate case without relying on constraint- Ascore-based methods. This requires to analysing the measure value of causal module in more sit- uations (e.g., omitting parent, involving an independent variable/non- descendant/non- -preent ancestor, or a case mixing the mentioned factors) and considering the characteristics to develop an efficient search algorithm similar as the greedy search algorithm (Chickering 2002).
C UNKNOWN CONFOUNDING
The unknown confounding has different influences on the results of our method in different situ- ations: 1) independent case: when two variables are independent; 2) causal case: when there is causal relationship between two variables. Depending on how the unknown confounding influences the pair of variables, 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 E, + U, where U is the unobserved confounder; E, 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 Eg + 1000U and Y E, 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 Eg +U and Y Then the distribution of Y can be dominated by U. Therefore, the FCM is close to X = E, + U and Y = 1000U, which leads to reversing the causal direction in the result. Moreover, when x Feg 10000 and Y = X F.y 1000U, then both distributions of X and Y can be dominated
15

---

Published as a conference paper at ICLR 2022
by U; hence, it is close X 10000 and Y 1000U, which leads to disregarding the causal relationship in the result.
Nevertheless, in practice, the impact of unknown confounding is 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, also find that the strength 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. for each FCM:
We generate a dataset with 1000 samples
FCM1) X = U and Y U;
FCM2) X = E2 + Wp x U and Y Ey + wy x U:
FCM3) X E + WG X U and Y f(X) Ey w U.
where Eg and Ey are the noise terms; U is the unknown confounder; wa and wy are the coefficients of U representing the confounding strength; used f(X) X and f(X) sin(4X) ) in the experiments. As for the data generation, Eg, Eg. and follow the uniform distribution, U((0, 1), and we vary the coefficients Wz 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 two-sample test 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 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 different applications, or deriving 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, 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 FCMI is that the p-value is equal to 1.0, which means that the measure values in the two directions not significantly different.
As for the experimental results of FCM2, we found that in the independent case, when the coeffi- cients wa and Wy are the same, 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.
get the correct results; when the coefficients are different, the
Table 4: The experimental results of FCM2 (the independent case) in the presence of unknown confounding.
W2 Awy 0.1
0.1
1.0
10
Wx My
0.1
1.0 causal
10 causal causal indep
0.354
1.57e-39
9.20e-05
0.1
indep
1.0
0.159
3.20e-49
1.0
indep
10
0.451
10
As for the experimental results of FCM3, we found that in the causal case, the confounding strength factor influencing the results more, which 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 type of asymmetry which can indicate the causal direction for our method, the non-monotonic functions can 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 Schölkopf 2018 Mastakouri al. 2021). Especially, (Mastakouri et al.
16

---

Published as a conference paper at ICLR 2022
Table 5: The experimental results of FCM3 (the causal case) with f(X) unknown confounding.
X in the presence of
WX wx WW
0.1
1.0
10
100
0.1
8.56e-56
2.98e-36
0.157
1.42e-06
1.0
1.22e-49
0.00027
9.60e-49
3.26e-51
10
0.002
7.02e-05
0.0018
0.007
100
0.385
0.398
0.577
0.70
w wy
WX
0.1
1.0
10
100
0.1
X cause
X cause
indep
X cause
1.0
X cause
X cause
Y cause
Y cause
10
Y cause
Y cause
Y cause
X cause
100
indep
indep
indep
indep
Table 6: The experimental results of FCM3 (the causal case) with S(X) = sin(4X) in the presence of unknown confounding.
Wa yyy
0.1
1.0
10
100
0.1
1.42--06
1.66e-68
5.62c-16
2.91e-09
1.0
5.68e-76
4.70e-81
5.62e-16
6.94e-51
10
2.17e-55
2.70e-47
1.36e-51
0.50
100
2.93e-06
2.70c-47
9.81e-24
0.08
W
w
0.1
1.0
10
100
0.1
X cause
X cause
X cause
X cause
1.0
X cause
X cause
X cause
Y cause
10
Y cause
Y cause
X cause
indep
100
Y cause
Y cause
X
cause
indep
2021) shows the identifiability results on time-series data in the presence of memoryless unknown confounders where the confounder doesn't 't have autocorrelation effect. A potential research direction with our framework that by recovering the corresponding dynamical process of a static causal discovery problem while considering the memoryless confounding, can determine the causal direction between two variables in the presence memoryless unknown confounders.
D PROOFS AND DERIVATION
D. DERIVATION OF PROP.
Under constraints and
W2 (Po, Pr)
f M (xo) xol (x0))dxx
definition of W
[J Ey) - (pE,) )dpp (EE)dEE
Eqn. (6. and the independence of E,, and Ey
EE,
D.2 PROOF OF THE EXTENSION OF PROP.
Suppose that in the multivariate ANM, Xy = [E,,... .... Emi's and Xy = [X,,.., Xmn',, and X; ge Ax) E,, where P P, denotes the parent of X..
17

---

Published as a conference paper at ICLR 2022
(po.o) f IM (xx) -xolpu((xo
/ ((E1(E,,.. Emll" - Eml.., ((1, .... Em dbxo
Σ / |M Em_ Ed pd) ..... Em )dxx.
(where Mt represents the i-th element of M*)
/ [M E,.l)) (EE E1.m ()
(where 1 m. - represents 1,..., - 1, + 1,..., m. and the independence of noise implies pl Ee (()m))(Emm-.))
[5(Fs, Fism-i) - p( (d p( F1:m (d
/ [F(E, - E.l pEI'dE,) pl E1 -d)
Σ [f(E,, (E,'dE,, p (Eumm
(where applying the change of variable formula)
[[(E, Xi) E. p( ildeee()
[f(Fs, X1mm Egl2 ((F)) deE pl 1)
(where [det(J)
for ANMs)
EPA.1 / |F(E E. E.P(EE)) m
Σ EPA.
D.3 PROOF OF THM. 1
Part We first derive the time interpolation p(t, x,) of Po and pr under dynamical FCM con- straints. Suppose that the FCMs are Y f(X, Fey)) As for the dynamical formulation of the 2 Wasserstein distance under constraints ( and )), E T), where X1, Xu E IR and J Gxg the Jacobian matrix w.r.t
according to the Jacobian equation,
p(t, x1) (x(x))/1 (J
0%0
dt det([J yxr )l
1+
Xo- Moreover, Bog p(t.)))
axi det(J
380
where Jasc
T Tax:
1 05
1) +
0%0
18

---

Published as a conference paper at ICLR 2022
According to Jacobi's formula and (Tao 2013).
d
dt (e)()(t))
dett(A(t))tr
[A(t) dt A(t))
d
+ det((((t)tr
(de((A(t)- adj( di A(1)
d + ΣΣ adj(A(0)) @ de A(t))
where is element wise multiplication and adj(A) adjugate matrix of A and the sign takes the same sign as det(J dxa Furthermore, since we know that under the structural constraints ① and
M is Eqn. (6) Thus, replace A(L) with J Dxg 8x0
= =I {VM I), and then have
1 1+
J Jas
T (
1) +
T 8x
1 {
0f
1) +1,
$ af
adj (J oxs y
BE
T. ax 1+ T
0,
d J Jas dt ong
0. 1 a
0 81 DE
(VM* - =
T
TEX
T
1)
Therefore,
d det(J oxs
1 0f
t 1)(1 + +
dt
U30
T
DE,
which takes the same sign as det
oxx
Part II
We then proof the property of the corresponding dynamical systems of ANMs in Thm.
det(J ax,
Because 0 Blog p(t.x.)
under the conditions in Thm. 1 0f and is obvious
Ot
de
aE,
that d det(Jamm Jm ) and
[J
+ O, Vt E [0, 7]. Furthermore, according to the theorem
of instantancous change of variables Chen et al. 2018), which is variant of Fokker-Plant equation, we know that
log p(x(t))
tr( tr(Vv) = -div v.
Øt
Therefore, div v = 0.
D.4 PROOF OF PROP.
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 ().. we know that it is sufficient check the divergence of the velocity field at time O, and that the divergence taking zero value everywhere at time 0 implies the divergence taking zero value everywhere for i [). Because of Thm. we know the div 0 at time Therefore, D(v)= 0.
Sufficient direction. We prove the sufficient direction by showing the contradiction with the iden- tifiability of ANMs. Given that D(v) 0 under the constraints 4. M. (III). and the identifiability conditions of ANMs, it tells us that there is an ANM in the form Equ. (I) which is consistent with the data distribution and has independent noise. Because of the identifiability of ANMs shown
19

---

Published as a conference paper at ICLR 2022
by (Hoyer al. 2008), there is 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 a concern that may happen that in the direction which is not the causal direction, the model is not 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 general, this will happen by proving that under weak assumptions D(v)
implies the model is ANM. First, D(v) = 0 implies div v = 0 for all Xo with positive probability densities p(xo) >> 0. In the following say div v in short. Second, div 0 implies @p(t.xx) 0 a.e. according Thm.
di det(J sx )l 750
0f where Y = f(x, Fy) a.e. according to the proof of Thm. Therefore, under the assumptions: 1) EEy) is positive in continuous range of Ey; 2) in the range, aE, is a quantity which doesn't change with Fy, e.g., it can be function of X or constant. Therefore, D(v) = implies Y F.g g(X) further under assumptions 1) and 2).
Consequently,
=0, which leads to dE
af = 1 a.e. implies that it holds everywhere, we have = Ey + C, where C
E IDENTIFIABILITY CONDITIONS OF ANMS IN (HOYER ET AL., 2008
Because the identifiability conditions of ANMs in (Hoyer et al. for our Thm. 1 we include them here:
2008 important and necessary
Pa ((=)))=(xx), where Pao pz are probability densities on R. If there is a backward model of the same form, i.e., p(æ, y) log Pn. and == log Pa. Po. the triple (J.Pa,Po) must satisfy the following differential equation for all x, with v" f(z))f # 0:
Let the joint probability density of and y be given by p(x,
pa (2 - s(v) (y,)), then, denoting
5" +"" 20" vn v/ v/v' f" 1f(1")2 = -- 10 g
1/'1/" 1/
f
where have skipped the arguments y/ f(x), x, and for v, &. and I and their derivatives, respectively. Moreover, for fixed pair (ʃ, v) there exists y E R such that v" - (y-[(x)))['(e + O for all but a countable set of points E I, the set of all Pa for which p has a backward model is contained in a 3-dimensional affine space.
F EXPERIMENTS AND DETAILS OF DIVOT
For the synthetic data experiments in Sec. we visualized the generated data with the sample size 500 in Fig. For determining causal direction, we used the suriance-base divergence measure in
30 42 14
de dh
Figure 3: Synthetic data generated from ANMs with different (non)linear functions. From left to right, they are (1) Y X + Eyg (2) Y 0.1(2.5X) - 0.1X + Egi (3) sin(4X) + Egi (4) Y == 0.5X3 - + E, if < 0, and Y = 0.5x3 + X + Ey, otherwise, where X ~ U(-1,1) and E, ~ U(, 1) are uniform distributions.
Eqn. 11 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. For each sample size, 100 different datasets are generated. For all the synthetic data experiments in Sec. 6 we use batches for computing Eqn. (11) without using debiasing functions. In the following, we introduce positions (App. F.I. batches (App. F.2), debiasing functions (App. F.3), and optimization details (App. [44) of DIVOT; as well as show its robustness to prior misspecification (App. F.5). more comparison with other benchmark methods (App. [66,, and its efficiency (App. F.7
20

---

Published as a conference paper at ICLR 2022
0.7
Deeliiidbbaasnng Dorlv)
2.2
0.6
Estimated parameter:
Accuracy
0.4-
0.3
0.2
debbbaiigg
D
Batch size
Batch size
(a) The influence of batch sizes determining causal direction the improvement with linear debiasing function.
(b) The influence of batch sizes estimating p(E,; 0) and improvement with linear debiasing function. The ground-truth 0
Figure 4: The influence of the batch size and the debiasing function causal direction determina- tion panel (a) and density estimation in panel (b). The synthetic data (of which the sample size is 100) are generated from the linear ANM, Y X Ey, where
U) and E, (4(O, 1) are uniform distributions. In panels (a) and bb). the batch size is represented as the ratio of the sample size, e.g., 0.1 represents that the batch size 0.1 100; the red triangular represents the experiments with only using the variance-base based divergence measure; the blue star represents using the the srriance-based divergence measure and linear debiasing function. The error bar in panel (b) represents standard deviation.
F.1 NUMBER OF POSITIONS FOR COMPUTING THE VARIANCE-BASEN DIVERGENCE MEASURE
Suppose a dataset ( y:) with sample size N. For computing the divergence measure, we call {ar) N position and for each { we need to specify the corresponding vector Te of which the number of elements is N₂ and generated samples of ey. In practice, it is not necessary use all the N positions for computing Eqn. (IT). For the synthetic data experiments in Sec. 6 when N > 50, we choose 50 positions out of N: otherwise, we choose all the N positions. To select a position (when N > 50), we first find the maximal and minimal values of {a;) N. compute the interval length, (max( (x) minn( /50, and then choose a position every such length (or its nearest position available in the data). We can then compute Eqn. (11 for causal direction determination.
F.2 BATCHES FOR THE FINITE (OR LIMITED) SAMPLE SCENARIO
In the finite sample scenario, it can be the case that there is no other data at position r. Therefore, the neighbors of a batch as having the same value r. We represent the batch size with the percentage of the total number N of samples. For the synthetic data experiments in Sec. 6 we 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.
batch of data for computing Eqn. (1).. We consider all the data in
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 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. (TT). 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, 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 different corresponding noise distribution compared with the ground souund-ruue true EEyy) 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 faw-sample case where we have to choose a larger batch size for computing the divergence measure. Therefore, use debiasing function such
21

---

Published as a conference paper at ICLR 2022
1.0655056
9c(x)
1.5
1.5
g(x) in ANM
1.0
1.0
0.5
0.0
0.5
-0.5
-1.0
-1.5-
-1.00
-0.75
-0.50
0.00
0.25
0.50
0.75
1.00
-2.0
-1.00
-0.50
0.00
0.75
X
X
(a) The debiasing function and ( in the ANM.
(b) The samples of (a1, vi (
Figure 5: As by-product, debiasing functions can be used for estimating the nonlinear function of an ANM in Eqn (T). We optimize the parameters
neural network representing nonlinear debiasing function and the parameter of uniform distribution representing the hypothesized noise distribution. As for panel (a). the ground-true true noise distribution is corresponding to = and the estimated parameter 0 is 1.066. As for panel (b), given dataset {(c;; y:)), y ydabies = - ga(m:)
that we can use a larger sample size for the few-sample case without sacrificing the performance of DIVOT.
F.3 DEBIASING FUNCTIONS FOR THE FEW- SCENARIO
Because the mentioned bias is introduced by omitting the position information, we introduce it back with debiasing function ga( to reduce the bias by modifying Eqn. 11 as
sort (ort(y sort( (g)) - ave debala
Ey) cg ~ p(Er;; 0),
Dvas(V)
Σ
N₂
2E(=;tN
debais = V/ ga(z; w),
(13)
where (a(;; w) a (non)linear function parameterized with w, e.g., in our experiments the linear debiasing function is ga(x) guarantees not only the accuracy while using larger batch size than 0.2, but also the correctness of the density estimation.
x As shown in Fig. [ using the linear debiasing function
For more complicated scenario, e.g., the ANM uses a nonlinear function, one can apply neural network the debiasing function. We found that as a by-product, sufficient flexible debiasing function can be used for estimating g(z) ) of ANMs in Eqn. ( As shown in Fig. 5 we generate synthetic 500 data of an ANM with the nonlinear function Y x(X)++++ E, 0.1 (2.5x) a 0.1X+ Ey. where X ~ (1-1,1) 1) and Ey (..5,, 0.5) satisfy uniform distribution. The (1( is close the ground-truth (x)). Because this not the main focus of the work, we would like 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 large batch size. This means that when we use restrictive class of ga(z), it can lead to noticeable difference of ga(z) and g(x). But for the purpose of causal direction determination, the divergence measure does not require specify the nonlinear functional form of ANMs. Moreover, DIVOT based on optimal transport does not require to estimate gxx) either. Therefore, DIVOT still determine causal direction correctly in such case.
22

---

Published as a conference paper at ICLR 2022
F.4 OPTIMIZATION METHODS AND THE CONVEXITY OF THE OBJECTIVE FUNCTIONS OF 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 Sec. c 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. I1 as an objective function, of which the parameter is A. Because when 0 > O, the linear function fg fooLse (eg (eoorce)) does not change the order of sorting results of ey egqurcé, have
(v)r(v)
(Ta) - fg (sort( @oourod)) - ave fnotsn se esouros) (w
y
y
N
N₂
(ort(TE) - x sort( aourcb
ave Th - 0 x E
aourch
N
N
fo foo1se
(Cy
Sourrd 06 0
x
0,
rou c
1(0
(14)
It is obvious that the divergence measure in Eqn. (14) convex on 1 > 0. The convexity is achieved by parameterizing the noise distribution with 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 8". The gradi- ent jax. . glass) is computed with the autograd in JAX (Bradbury et al., 2018). We update 0 by specifying step size sz and 0 = 0 0 jax grad sz. If after the update 0 < 0 (which has never happened), we set the value of as positive number close to zero. We used s2 = 1 for all the synthetic data experiments. Moreover, because of the convexity. one can also use one-step update method to directly get the best parameter 6 by finding the root of the gradient function of Eqn. (4). We simply give a range of f and find 6 by a binary search method such that the gradient at 8 is equal to zero. The experiments of DIVOT for ANMs in the real-world data experiments also [0, 100). In addition, we need
used the one-step update for updating 0 and the range 0 is specified 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 w in (12). The program of DIVOT terminated when the divergence measure converges. If more complicated scenario requires froise to be nonlinear function, one may need to use the gra- dient descend method instead of the one-step update. We find that it is sufficient to use the standard Gaussian distribution with the linear fø fnoise for DIVOT to have promising results in the experiments on the Tübingen datasets, which indicates that DIVOT is robust to the choice of models. We then investigate the robustness of DIVOT to prior misspecification.
F.5 ROBUSTNESS TO PRIOR MISSPECIHICATION
We use different hypothesized noise distributions in DIVOT to test the robustness to the misspecifi- cation of noise distribution. The synthetic data are generated in Sec. with uniform distrbution noise E, ~ (0, 1). In DIVOT, we use one of the three hypothesized distributions: uniform distri- bution (0, 1). beta distribution B(a 0.5, b = 0.5), and standard normal distribution (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.
F.6 COMPARISON WITH RESULTS OF BENCHMARK METHODS
We compare DIVOT with other benchmark methods, such as ANM (Hoyer et al., 2008), CAREFL (Khemakhem et al. 202I). RECI (Blöbaum et al., 2018), and LLR (Hyvärinen and Smith 2013).
As in (Khemakhem al. 2021). the synthetic data are generated with the Laplace distribution as
23

---

Published as a conference paper at ICLR 2022
Y=x+E, Ey
0.1((2.5X)) +
[piece(X) + Ey
1.0
1.0
1.0
0.9
0.9
0.88
Accuracy
0.8
0.7
0.7
Uniform
Uniferm
Uniform
Beta
Beta
Beta
0.4
Normal
Normal
Normal
0.5
0.5
100
200
300
400
500
200
400
200
400
Sample size
Sample size
Sample size
Figure 6: Robustness to prior misspecification: The performance of DIVOT with different hypothe- sized distributions of Eyy 0) on the datasets with uniform distribution noise Ey ~ L(, 1).
Y= + Ey
Y= Ey
1.0
1.0
0.8
0.8
Accuracy
Accuracy
0.6
0.6
100
200
300
400
500
0.0
100
200
300
400
500
Sample size
Sample size
Linear SEM
Additive noise SEM
X₂ = + 22
X2 X1 + x + Z2
Proportion correct
0.8
Proportion correct
0.6
0.4
Algorithm
CAREFL
0.2
CAREFL-NS
Linear
0.0
100
200
300
400
500
100
200
300
400
500
Sample
Sample
Figure 7: Comparison of the results of DIVOT (the first row) with the reported results in (Khe- makhem et al. 2021) (the second row). The figures in the second row are taken from the screen- shot of (Khemakhem et al.] (202I). The noise (Z, and E,y) distribution of the synthetic data is a Laplace distribution. CAREFL represents the causal autoregressive model (Khemakhem et al., 2021;; CAREFL-NS represents the causal autoregressive model without scaling: RECI represents the method, regression eπror causal inference (Blöbaum et al., (018); ANM (Hoyer et al. 20088) uses Gaussian process while ANM-NN uses a neural network: linear LR the linear likelihood ratio method (Hyvärinen and Smith 2013).
the noise distribution. DIVOT uses batches without debiasing functions. And the data preprocessing 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.
24

---

Published as a conference paper at ICLR 2022
F.7 EFFICIENCY OF DIVOT
For common causal discovery tasks in the bivariate case, the sample size 10000 is a large and chal- lengeable Causal discovery methods need 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 10000 generated with and Ey ~ (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 19. Our implementation is based on JAX (Bradbury et al., 2018 which uses Apache License and the running time is measured with command Stimeit JAX.
0.1((2.5XX) -X) Ey
Table Running time of DIVOT on synthetic data, of which the sample size is 10000. The DIVOT with batches and no debiasing function is tested with different number of positions and batch sizes.
batch size
50 positions
100 positions
0.001
240 ms + 3.48 ms
431 ms 4.33 ms
0,01
7.71 + 475 ms
14.3 + 262 ms
0.1
Imin 7s 1.81
2min 19s + 2.11 8
25


Reviews:
[{"review_id": "uSMtXTwci4i", "reviewer": "Reviewer_pbTi", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "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": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Good and very interesting paper on a well-known but challenging problem. In particular the derivation of the new measure through the connection with dynamical systems is both intriguing and inspiring. \n\nPerhaps main drawback is that the information density in the paper is quite high and technical, and assumes familiarity of the reader with literature from other fields that may not be well-known to the wider ICLR audience, making it hard to understand certain parts of the paper in isolation. Even with a background in physics I had think back hard to remember details of vector field divergences and incompressible flow mechanics. That said, the ideas in the paper are well motivated and supported by rigorous theoretical derivations, making it sufficiently worthwhile for an interested reader to spend the extra effort.\nI briefly entertained the hope the approach might also apply in the presence of confounding, but as with similar approaches it does still require causal sufficiency and additive noise components. Performance shows improvement on existing methods, though still in the somewhat disappointing range of around 75% accuracy for a collection of 50-50 decisions. But that just shows the fundamental difficulty of the problem.\n\nIn short: good paper on challenging subject, bit technical but clear accept.\n\nFew remarks/questions:\np1, bottom: ‘.. of which the modern form was initiated by ..’, ‘.. by augmenting a time-dimension …’, ‘ … in the context of the Jacobian problem ..’ => these references are unlikely to be known by the majority of readers, and it would help to give just a minor hint of what these might entail \np2,top, point 2. ‘pressureless flow’ => do you mean the more standard notion of ‘incompressible flow’? or else explain the difference\np2,mid: ‘.. from p_0(x) to p_T(x) …’ => technically the random variables do not need to obtain their value after a constant time T, e.g. different coin throws end up heads or tails, but in the actual underlying dynamical system different throws can take different times to come to rest. Does this impact the interpretation in terms of the map M or is the actual value of T essentially irrelevant? If so then explain or perhaps remove it form the description.\np3,bot: ‘The sense of causality is enforced .. FCMs’ => this sentence is overly vague / ambiguous: please clarify\np4, final sentence above 3.2 .’Therefore optimality .. problem.’ => the implication of this sentence is unclear\np5, Prop2: can the authors comment/speculate on whether in the presence of unobserved confounders, in almost all cases we would expect to find D(v) =/= 0 in both the causal and anti-causal direction ?\nidem: ‘ .. and the identifiability conditions of ANMs (Hoyer,2008) …’ => these conditions are crucial and should be explicitly stated in the current paper\np7, Ext. to PNL: it is not entirely clear how this invertible function should be chosen based on a given data set\np7, ‘Minimisation’ , ‘Note .. is not necessary to be the true distribution or even significantly different from the true one’ => this sentence is unclear\n", "summary_of_the_paper": "The paper considers the problem of establishing causal direction in a bivariate system under the assumptions of additive noise (possibly post-nonlinear) and no unobserved confounding. It establishes a connection between functional causal models (FCM) and optimal transport problems in dynamical systems to derive a new divergence metric that can be efficiently computed as a conditional variance estimate, and shows this quantity is zero iff in the causal direction..Efficacy of the metric is evaluated on synthetic data as well as the real-world Tuebingen data stand found to outperform existing measures.", "summary_of_the_review": "Good and inspiring paper on challenging subject, bit technical but clear accept.\n", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "XRLTSIbh4BA", "reviewer": "Reviewer_GA77", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "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": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "I believe that the authors' attempt of reformulating bivariate causal discovery problem is a novel and inspiring endeavor. Given that the problem of causal discovery is a conceptual problem as much as it is a statistical problem, the introduction of such new perspectives are certainly helpful for the discussion in this field. The authors' presentation is concise and clear, and the paper is easy to follow (modulo the additional comments I provide below).\n\nOne potential avenue for improvement is the motivation of the present reframing. What shortcoming of the information geometric, information theoretic, probabilistic, or more recently, topological analyses of the problem of causal discovery motivated the authors to present an optimal transport-based perspective? What additional advantages do the authors' formulation (or its potential sequels) provide to our understanding and methodology regarding causality? I find their work interesting regardless, but the preceding point is crucial for any work that formulates a novel mathematical framework causal discovery.\n\nThe authors limit their presentation to bivariate causal discovery problem. I do not think this is necessarily a crippling limitation, given the novelty of their theoretical perspective. However, the authors still need to provide an outlook for other cases as the authors do not seem to position their work as a purely conceptual development. If practical implications are to be considered, at least initial discussions of the multivariate case, as well as other standard challenges such as latent confounding, selection bias etc. should be provided.\n\nAnother potential development of the current work is to improve its presentation. If accepted, I encourage the authors to use the extra space conferred by the camera-ready version to provide a better introduction to the field of optimal transport, its prevalent use cases, as well as various concepts that are used in the field. The authors could also consider moving the related work section to the end of the introduction to further strengthen their initial setup.", "summary_of_the_paper": "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.", "summary_of_the_review": "The current paper proposes a novel and well-presented formulation of the causal discovery problem, and thus is of interest to the conference audience. Improved comparison with existing frameworks and expanded discussion regarding practical use cases are recommended to improve the paper.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "G8Gsa4D-9WC", "reviewer": "Reviewer_xS6T", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "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": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "As strenghts of the proposed method I see that is is nicely grounded on theorethical results, the paper is well structured, organized and written, and I personally enjoyed reading it.\n\nWeaknesses, as weaknesses are cocerned I see the following main weak points:\n1) it is not clear what happens, or at least I was not able to understand what happens in the case when the considered pair of variables are not related at all. Indeed, te proposed algorithm does not take into account this possibility. Lines 6-9 of Algorithm 1, which furthermore ignores that all we have is a sample and that which is which should be at least addressesd, what about the quantity you use are the same or nearly the same?\n2) the decision to remove outliers is in itself something that practically could also be understood but from the theoretical point of view opens some doubts about the effectiveness of the proposed apporach. Indeed, standardization and subsequent removal of observations acconding to the 2 standard deviation limits brings bias in my opinion\n3) numerical results seem to be in line with state of the art algorithms, however it seems that the proposed approach suffers from greater uncertainty about accuracy when compared to other literature approaches.\n\nI also found some minor issues with typos.", "summary_of_the_paper": "The paper tackles the problem of causal discovery in the basic case where a pair of variables is considered.\nIn particular, it is concerned with Functional Causal Models and how to establlish the causal direction for a pair of continuous variables.\nIt motivates the contibution by poiting out thet the performance of availabel algorithms is sensitive to the model assumptions, which\nmakes it difficult for practitioners to use them in real settings.\nThen, 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.", "summary_of_the_review": "The paper tackles an interesting and well studied problem under a general setting, i.e. assuming uncounfondedness.\nIt is well written and proposes a nice method to discover the direction of causal relationship when considering pair of variables in isolation.\nI have expected some more discussion about what to do in the case where the terms used in Algorithm 1 were not too different or extremely close. Furthermore, the results reported in Table 1 witness that the proposed methods achieves accuracy values similar to state of the art algorithms, while bringing more uncertainty about the true value of accuracy than other methods. Furthermore, for ANM the proposed method seems to be inferior to other state of the art algorithms.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper is a solid contribution to researchers in this field, as it provides a new idea for the basic problem of determining the direction of causality between two variables, using the functional causal model as a dynamical system and optimal transport.
Decision: Accept (Spotlight)