Title: VARIABLE COMPLEXITY IN THE UNIVARIATE AND MULTIVARIATE STRUCTURAL CAUSAL MODEL

Abstract: We provide a critical view of the univariate Structural Causal Model. We show that by comparing the individual complexities of univariate cause and effect, in the Structural Causal Model, one can identify the cause and the effect, without considering their interaction at all. The entropy of each variable is ineffective in measuring the complexity, and We propose We suggest to capture it the complexity by the reconstruction error of an autoencoder that operates on the list of sorted samples the percentiles of the distribution. Comparing the reconstruction errors of the two autoencoders, one for each variable, is shown to perform well on the accepted benchmarks of the field. Hence, the decision as to which of the two is the cause and which is the effect may not be based on causality but on complexity. In the multivariate case, where one can ensure that the complexities of the cause and effect are balanced, we propose a new adversarial training method that mimics the disentangled structure of the causal model. We extend the results of Zhang & Hyvrinen (2010) to the multidimensional case, showing that such modeling is only likely in the direction of causality. Furthermore, the learned model is shown theoretically to perform the separation to the causal component and to the residual (noise) component. Our multidimensional method obtains a significantly higher accuracy than the literature methods.

1. INTRODUCTION:

A long standing debate in the causality literature, is whether causality can be inferred without intervention (Pearl, 2009; Spirtes et al., 2000) . The Structural Causal Model (SCM) (Spirtes et al., 2000) is a simple causative model for which many results demonstrate the possibility of such inference (Stegle et al., 2010; Bloebaum et al., 2018; Goudet et al., 2018; Lopez-Paz et al., 2017; 2015) . In this model, the effect (Y ) is a function of the cause (X) and some independent random noise (E).

In this work, we take a critical perspective of the univariate SCM. We demonstrate empirically that for the 1D univariate case, which is the dominant case in the existing literature, the SCM model leads to an effect that has a lower complexity than the cause. Therefore, one can identify the cause and the effect by measuring their individual complexities, with no need to make the inference based on both variables simultaneously. Thus, the decision as to which of the two is the cause and which is the effect may not be based on causality but on complexity.

Since we are dealing with unordered 1D univariate random variables, the complexity measure has to be based on the probability distribution function. As we show empirically, comparing the entropies of the distribution two random variables is ineffective for inferring the causal direction. We, therefore, consider the quantiles, i.e, fixed sized vectors that are obtained as sub-sequences of the sorted sampled values of the variable.

We consider suitable complexity scores for these vectors. In our analysis, we show that the reconstruction error of an autoencoder of a multivariate random variable is a valid complexity measure. In addition, we link the reconstruction error based complexity, in the case of variational autoencoders, to the differential entropy of the input random variable. Hence, by computing the reconstruction errors of trained autoencoders on these vectors we estimate the entropies of the quantile vectors of X and Y .

The challenges of measuring causality independently of complexity in the 1D case lead us to consider the multidimensional case, where the complexity can be controlled by, e.g., manipulating the dimension of the noise signal in the SCM model. Note that unlike Goudet et al. (2018) , we consider pairs of multivariate vectors and not many univariate variables in a graph structure. We demonstrate that for the multidimensional case, any method that is based on comparing the complexity of the individual random variables X and Y fails to infer causality of random variables. Furthermore, we extend the result of Zhang & Hyvrinen (2010) to the multidimensional case and prove that an SCM is unlikely to hold in both directions X → Y and Y → X for reasonable conditions.

Based on our observations, we propose a new causality inference method for multidimensional cause and effect. The algorithm learns three networks in a way that mimics the parts of the SCM model. The noise part is unknown and is replaced by a function that is constrained to be independent of the cause, as captured by an adversarial loss. However, we show analytically and empirically that even without the explicit constraint, in several cases, such an independence emerges.

Our empirical results support our analysis and demonstrate that in the 1D univariate case, assigning cause and effect based on complexity is competitive with the state of the art methods. In the multidimensional case, we show that the proposed method outperforms existing and new extensions of the literature methods to the multidimensional case.

1.1. PROBLEM SETUP:

We investigate the problem of causal inference from observational data. A non-linear structural causal model (SCM for short) is a generative process of the following form:

EQUATION

The functions g : R d f +de → R n and f : R n → R d f are fixed and unknown. In general, g and f are non-linear. Here, X is the input random variable and E is the environment random variable that is independent of X. We say that X causes Y if they satisfy a generative process such as Eq. 1.

We present methods for inferring whether X causes Y (denoted by X → Y ) or Y causes X, or neither. The algorithm is provided with i.i.d samples {(x i , y i )} m i=1 ∼ P m X,Y (the distribution of m i.i.d samples from the joint distribution P X,Y ) from the generative process of Eq. 1. In general, by (cf. Prop 4.8, (Peters et al., 2017) ), for any joint distribution P X,Y of two random variables X and Y , there is an SCM, Y = g(f (X), E), where E is a noise variable, such that, X |= E and h is some (measurable) function. Therefore, in general, deciding whether X causes Y or vice versa is ill-posed when only provided with samples from the joint distribution. However, Zhang & Hyvrinen (2010) showed for the one dimensional case (i.e., X, Y ∈ R) that under reasonable conditions, a representation Y = g(f (X)+E) holds only in one direction. In Sec. 3.2 we extend this theorem and show that a representation Y = g(f (X), E) holds only in one direction when g and f are assumed to be neural networks and X, Y are multidimensional (we call such SCMs neural SCMs).

Throughout the paper, we denote by P U [u] := P[U ≤ u] the cumulative distribution function of a uni/multi-variate real valued random variable U and P is a standard Lebesgue measure. Additionally, we denote by p U (u) = d du P U [u] the probability density function of U (if exists, i.e., P U [u] is absolutely continuous). We denote by E u∼U [f (u)] or by E p U [f (u)] the expected value of f (u) for u that is distributed by P U [u], depending on the context. The identity matrix of dimension n × n is denoted by I n . We recall a few notations from information theory (Shannon, 1948; Cover & Thomas, 2006) . The entropy of a random variable U , is defined as h(U

) := -E p U [log(p U (u))],

where U is the domain of U . The mutual information between two random variables U and V is denoted by

I(U ; V ) := E p (U,V ) log p U,V (u,v) p U (u)p V (v)

. Sometimes, we will denote I p (U ; V ), to specify the PDF, p, of the variables (U, V ).

1.2. RELATED WORK:

In causal inference, the algorithm is provided with a dataset of matched samples (x, y) of two random variables X and Y and decides whether X causes Y or vice versa. The early wisdom in this area asserted that this asymmetry of the data generating process (i.e., that Y is computed from X and not vice versa) is not apparent from looking at P X,Y alone. That is, in general, provided with samples from the joint distribution P X,Y of two variables X, Y does tell us whether it has been induced by an SCM from X to Y or from Y to X.

In publications such as (Pearl, 2009; Spirtes et al., 2000) it is argued that in order to decide whether X causes Y or vice versa, one needs to observe the influence of interventions on the environment parameter. To avoid employing interventions, most publications assume prior knowledge on the generating process and/or independence between the cause and the mechanism.

Various methods for causal inference under the SCM model have been suggested. Many of these methods are based on independence testing, where the algorithm models the data as Y = g(f (X), E) (and vice versa) and decides upon the side that provides a better fitting in terms of mapping accuracy and independence between f (X) and E = r(X, Y ). The LiNGAM (Shimizu et al., 2006) algorithm assumes that the SCM takes the form Y = βX + E, where X |= E and β ∈ R and E is non-Gaussian. The algorithm learns β, such that, X and Y -βX are independent by applying independent component analysis (ICA for short). The Direct-LiNGAM (Shimizu et al., 2011) extends this method and replaces the mutual information minimization with a non-parametric kernel based loss (Bach & Jordan, 2003) . However, the computation of this loss is of order Θ(m 2 ) in the the worst case (m is the number of samples).

The ANM approach (Hoyer et al., 2009 ) extends LiNGAM's modeling and assumes that Y = f (X) + E, where X |= E. In their modeling they employ a Gaussian process for the learned mechanism between the two random variable. The function f is trained to map between X and Y (and vice versa) and the method then tests whether, X and f (X) -Y are independent. The independence test is based on kernels (Gretton et al., 2005) .

A different extension of LiNGAM is the PNL algorithm by Zhang & Hyvrinen (2010) . This algorithm learns a mapping between X and Y (and vice versa) of the form Y = g(f (X) + E), where f (X) and E are restricted to be independent. To do so, PNL trains two neural networks f and g to minimize the mutual information between f (X) and E = g -1 (Y ) -f (X). The main disadvantage of this methods is the reliance on the minimization of the mutual information. It is often hard to measure and optimize the mutual information directly, especially in higher dimensions. In many cases, it requires having an explicit modeling of the density functions, because of the computation of expected log-probability within the formulation of the entropy measure.

In our method, we take a similar approach to the above methods. However, our independence constraint is non-parametric and is applied on the observations rather on an explicit modeling of the density functions. In addition, we do not assume restrictive structural assumptions and treat the generic case where the effect is of the form Y = g(f (X), E). In addition, our method is computationally efficient.

Several other methods take different approaches. For example, the Information Geometric Causal Inference algorithm (IGCI) (Daniusis et al., 2012) approach determines the causal relationship in a deterministic setting Y = f (X) under an independence assumption between the cause X and the mechanism f , Cov(log f (x), p X ) = 0.

The Conditional Distribution Similarity Statistic (CDS) (Fonollosa, 2016) measures the standard deviation of the values of Y (resp. X) after binning in the X (resp. Y ) direction. The lower the standard deivation, the more likely the pair to be X → Y . The CURE algorithm (Sgouritsa et al., 2015) compares between X → Y and Y → X directions in the following manner: is given as follows: if we can estimate p X|Y based on samples from p Y more accurately than p Y |X based on samples from p X , then X → Y is inferred.

The BivariateFit method learns a Gaussian process regressor in both directions and decides upon the side that had the lowest error. The RECI method (Bloebaum et al., 2018 ) trains a regression model in both directions, and returns the side that produced a lower MSE loss. The utilized regression models were: a logistic function, polynomial functions, support vector regression and neural networks. The CGNN algorithm (Goudet et al., 2018) uses the Maximum Mean Discrepancy (MMD) distance between the distribution produced by modeling Y as an effect of X, (X, g(X, E)) (and vice versa), and the ground truth distribution. The algorithm compares the two distances and returns the direction that led to a smaller distance. The Gaussian Process Inference model (GPI) (Stegle et al., 2010) builds two generative models, one for X → Y and one for Y → X. The distribution of the candidate cause variable is modelled as a Gaussian mixture model, and the mechanism f is a Gaussian process. The causal direction is determined from the generative model that best fits the data.

Finally, it is worth mentioning that several other methods, such as (Heinze-Deml et al., 2017; Zhang et al., 2011) assume a different type of SCM, where the algorithm is provided with separate datasets that correspond to different environments, i.e., sampled i.i.d from P X,Y |E , where the value of E is fixed for all samples in the dataset. In these publications, a different independence condition is assumed: Y is independent of E given X. This assumption typically fails in our setting. In our paper, we focus on a vanilla SCM, where the algorithm is only provided with observational samples of X and Y = g(f (X), E) (i.i.d samples). The samples are not divided into subsets that are invariant w.r.t E.

2. THE UNIVARIATE CASE:

In this section, we intend to show that the univariate case is a bit simplistic. For this purpose, we start by describing a method for identifying the cause and the effect, which considers each of the two variables independently without considering the mapping between them. The method is based on the statistics of the quantiles of each of the two distributions. In Sec. 2.2, we show that the reconstruction error of an autoencoder is directly linked to the differential entropy of the distribution it models.

2.1. A METHOD FOR INFERRING DIRECTIONALITY:

We hypothesize that due to the information processing inequality, the entropy of the cause h(X) in the SCM model is higher than the entropy h(Y ) of the effect. However, estimating the entropy of each random variable from its samples does not present a consistent difference between the entropies h(X) and h(Y ).

Our method, therefore,The proposed method computes an alternative a complexity score for X and, independently, for Y . It then compares the scores and decides that the cause is the random variable that has the larger score among them.

Our scoring function has a few stages. As a first step, it converts the random variable at hand (say, X) into a multivariate random variable. This is done by sorting the samples of the random variable, and then cutting the obtained list into fixed sized vectors of length k. We discard the largest measurements in the case where the number of samples is not a multiple of k. We denote the random variable obtained this way by U .

At the second stage, the method trains an autoencoder A : R k → R k on these quantile vectors. Formally, A is trained to minimize the following objective:

EQUATION

where (a, b) is some loss function. In our implementation, we employ the L 2 -loss function, defined as (a, b) = a -b 2 2 . Finally, the method uses the value of L recon (A), which we refer to as the AEQ score, as a proxy for the complexity of X (smaller loss means lower complexity). It decides that X or Y is the cause, based on which side provides a higher AEQ.

2.2. RECONSTRUCTION ERROR AS AN ESTIMATOR OF THE ENTROPY:

Next, we would like to provide better understanding of the AEQ method. Since we use the reconstruction error as a measure of complexity, we next link it to entropy, which is an acceptable complexity measure. This will explain that the AEQ method estimates and compares the entropies of the percentiles of X and Y . This does not imply that the AEQ method compares between the entropies of X and Y .

The analysis This is done in the context of a learned variational autoencoder (VAE) (Kingma & Welling, 2014) , with a fixed latent space dimension.

We are given a real-valued vector-random-variable U (not necessarily quantiles) with a PDF p U (u). We recall the setting of VAEs and the analysis of Zhao et al. (2017) . The VAE model employs a d dimensional latent variable z ∼ N (0, I d ), and considers the family of the distribution p θ (u|z) parameterized by θ ∈ Θ. The objective of this framework is to select θ that maximizes the likelihood of the samples of U within p θ (u|z) • p(z),

E u∼U [log p θ (u)] = E u∼U R d p θ (u|z) • p(z) dz (3) where p(z) = exp(-z 2 2 /2) • (2π) -k/2 is the PDF of z and p θ (u) = p θ (u|z) • p(z).

In typical situations, we select p θ to be the PDF of a normal distribution  and σ 2 1 is a variance hyperparameter of the framework. We denote by P = {p θ (u|z)|z ∈ R d , θ ∈ Θ} the class of density functions over X .

N (f θ (z), σ 2 1 • I n ), where f θ (z) ∈ F is a learned neural network

To optimize this objective, it is required to evaluate p -1 θ (u), which involves the computation of a high dimensional integral. Even though this integral can be approximated by sampling, this is a costy process that has to be done for every sample u. To solve this problem, we have a second distribution that is responsible to model the latent codes z that are most probable given u. This PDF is denoted by q φ (z|u) and corresponds to a normal distribution

N (g φ (u), σ 2 2 • I d ).

Here, g φ (u) is a neural network with parameters φ ∈ Φ.

Suppose we are given a fixed inference distribution q φ (z|u), which maps (probabilistically) inputs u to features z. We denote by q φ (u, z) := p U (u) • q φ (z|u) the joint distribution of u and z, the marginal q φ (z) := u q φ (u, z)du and the posterior q φ (u|z) = q φ (z|u) • p U (u)/q φ (z).

The optimization criterion is redefined as follows, where for each z we use a member of a class of distributions, P, to fit a different q φ (u|z) rather than the entire p U .

EQUATION

The following lemma links the reconstruction loss and the entropy (all proofs are in the appendix) Lemma 1. In the setting of Sec. 2.2. Let θ * be the global optimum of Eq. 4 for a sufficiently large F. Assume that P is sufficiently large, such that:

∀z ∈ Z, φ ∈ Φ : q φ (u|z) ∈ P (5)

Then,

EQUATION

Therefore, from this lemma, we observe that if g φ is learned, such that, the expected norm of the encodings of u, E u∼U g φ (u) 2 2 is bounded, then, the differential entropy h(U ) is proportional reconstruction error of the learned autoencoder. In particular, for two random variables U 1 and U 2 , their entropies can be compared, by considering the reconstruction error of their autoencoders with the same hyperparameters (latent space dimension d, σ 2 1 , σ 2 2 , etc').

3. THE MULTIVARIATE CASE:

For the univariate case, one can consider the complexity of the X and Y variables of the SCM model and infer directionality. We propose the AEQ complexity for this case, since more conventional complexities are ill-defined for unordered 1D data or, in the case of entropy, found to be ineffective.

In the multivariate case, one can consider the complexity of the random variables in various ways. We focus on the family of complexity measures C(X) that satisfy the assumption that when X and Y are independent

C(X, Y ) ≥ max(C(X), C(Y )).

(7) Sample complexity measures that satisfy this condition are the Shannon Entropy and the Kolmogorov Complexity. Lem. 3 in Appendix. B.2 shows that the complexity that is based on autoencoder modeling is also in this family.

The next result shows that for any complexity measure C in this family, one cannot infer directionality in the multivariate SCM based on C. Lemma 2. Let C be a complexity measure of multivariate random variables (i.e, non-negative and satisfies Eq. 7). Then, there are triplets of random variables (X, E, Y ) and ( X, E, Y ), such that, Y = g(X, E), Y = g ( X, E), C(X) < C(Y ) and C( X) > C(Y ). Therefore, C cannot serve as a score for causal inference.

3.1. A METHOD FOR CAUSAL INFERENCE:

The causality inference algorithm trains neural networks G, F, R and D. The success of fitting these networks serves as the score for the causality test.

The function F models the function f , G models g and R(Y ) aims to model the environment parameter E. In general, our method aims at solving the following objective:

EQUATION

where A is either X or Y and B is the other option, and

a i = x i , b i = y i or a i = y i , b i = x i accordingly.

To decide whether X → Y or vice versa, we train a different triplet G, F, R for each direction and see if we can minimize the translation error subject to independence. We decide upon a specified direction if the loss can be minimized subject to independence. In general, searching within the space of functions that satisfy A |= R(B) is an intractable problem. However, we can replace it with a loss term that is minimized when A |= R(B).

Independence loss We would like R(B) to capture the information encoded in E. Therefore, restrict R(B) and A to be independent in each other. We propose an adversarial loss for this purpose, which is a modified version of a loss proposed by Brakel & Bengio (2017) and later analyzed by (Press et al., 2019) .

This loss measures the discrepancy between the joint distribution P A,R(B) and the product of the marginal distributions P A × P R(B) . Let d F (d R ) be the dimension of F 's output (R). To measure the discrepancy, we make use of a discriminator D : R n+d R → [0, 1] from a class C that maximizes the following term:

L D (D; R) := 1 m m i=1 (D(a i , R(b i )), 1) + 1 m m i=1 (D(â i , R( bi )), 0) (9)

where D is a discriminator network, and l(p, q) = -(q log(p) + (1 -q) log(1 -p)) is the binary cross entropy loss for p ∈ [0, 1] and q ∈ {0, 1}. In addition, {(â i , bi )} m i=1 are i.i.d samples from P A × P B . To create such samples, we can simply generate samples from âi ∼ P A and bi ∼ P B independently and then arbitrarily match them into couples (â i , bi ).

To restrict that R(B) and A are independent, we train R to fool the discriminator D and to convince it that the two sets of samples are from the same distribution,

EQUATION

Full objective The full objective of our method is then translated into the following program: min

G,F,R L err (G, F, R) + λ 1 • L indep (R; D) min D L D (D; R) (11)

Where λ 1 is some positive constant. The discriminator D minimizes the loss L D (D; R) concurrently with the other networks. We denote by c real the percentage of samples (a i , b i ) that the discriminator classifies as 1 and by c fake the percentage of samples (â i , bi ) that are classified as 0.

We note that when c real + c fake ≈ 1, the discriminator is unable to discriminate between the two distributions, i.e., it is wrong in classifying half of the samples.

Our method decides if X causes Y or vice versa by comparing the score L err (G, F, R). A lower error means a better fit. The full description of the architecture employed for the encoders, generator and discriminator is given in Appendix. A.

3.2. ANALYSIS:

In this section, we analyze the proposed method. In Thm. 1, we show that if X and Y admit a SCM in one direction, then it admits a SCM in the opposite direction only if the involved functions satisfy a specific partial differential equation.

Theorem 1 (Identifiability of neural SCMs). Let P X,Y admit a neural SCM from X to Y as in Eq. 1, such that p X , and the activation functions of f and g are three-times differentiable. Then it admits a neural SCM from Y to X only if p X , f , g satisfy Eq. 34 in the Appendix.

This result generalizes the one-dimensional case presented in (Zhang & Hyvrinen, 2010) , where a one-dimensional version of this differential equation is shown to hold in the analog case.

In the following theorem, we show that minimizing the proposed losses is sufficient to recover the different components, i.e., F (X) ∝ f (X) and R(Y ) ∝ E, where A ∝ B means that A = f (B) for some invertible function f . Theorem 2 (Uniqueness of Representation). Let P X,Y admit a nonlinear model from X to Y as in Eq. 1, i.e., Y = g(f (X), E) for some random variable E |= X. Assume that f and g are invertible. Let G, F and R be functions, such that,

L err := E (x,y)∼(X,Y ) [ G(F (x), R(y)) -y 2 2 ] = 0 and G and F are invertible functions and X |= R(Y ). Then, F (X) ∝ f (X) and R(Y ) ∝ E.

Here,

L err := E (x,y)∼(X,Y ) [ G(F (x), R(y)) -y 2 2 ]

is the mapping error proposed in Eq. 8. In addition, the assumption X |= R(Y ) is sufficed by the independence loss. In Lem. 6 in Appendix. B.3 we extend this lemma to the case where the mapping loss is not necessarily zero and get rid of the assumption that G is invertible. In this case, instead of showing that R(Y ) ∝ E, we provide an upper bound on the reconstruction of E out of R(Y ) (and vice versa) that improves as L err decreases.

We note that when the entropy of R(Y ) is smaller than the entropy of E, then, the independence X |= R(Y ) follows, even without the discriminator. Typically, when R has a limited capacity (i.e., small number of layers or output dimension), then, the entropy of R(Y ) is limited as well. Theorem 3 (Emergence of independent representations). Let P X,Y admits a nonlinear model from X to Y as in Eq. 1, i.e., Y = g(f (X), E) for some random variable E |= X. Assume that X and Y are discrete random variables. Let G, F and R be functions, such that, L err := E (x,y)∼(X,Y ) [ G(F (x), R(y)) -y 2 2 ] = 0 and G is an invertible function. Assume that h(R(Y )) ≤ h(E) and that g is invertible. Then, we have:

F (X) |= R(Y ), F (X) ∝ f (X) and R(Y ) ∝ E.

To conclude this section, by Thm. 1, under reasonable assumptions, if X and Y admit a multivariate SCM in direction X → Y , then, there is no such representation in the other direction. By Thm. 2, by training our method in both directions, one is able to capture the causal model in the correct direction. This is something that is impossible to do in the other direction by Thm. 1. By Thm. 3, we show that in cases where h(R(Y )) ≤ h(E), the independence between F (X) and R(Y ) emerges implicitly. This is the typical case, when R's capacity is limited.

4. EXPERIMENTS:

This section is divided into two parts. The first, is devoted to showing that causal inference in the one-dimensional case highly depends on the complexities of the distributions of X and Y . In the second part of this section, we show that our multivariate causal inference method outperforms existing baselines. Most of the baseline implementations were taken from the Causality Discovery Toolbox of Kalainathan & Goudet (2019) . The experiments with PNL (Zhang & Hyvrinen, 2010) , LiNGAM (Shimizu et al., 2006) and GPI (Stegle et al., 2010) are based on their original matlab code.

4.1. ONE-DIMENSIONAL DATA:

We compared the autoencoder method on several well-known one dimensional cause-effect pairs datasets. Each dataset consists of a list of pairs of real valued random variables (X, Y ) with their direction 1 or 0 depending on X → Y or Y → X (resp.). For each pair, we have a dataset of samples {(x i , y i )} m i=1 .

Five cause-effect inference datasets, covering a wide range of associations, are used. CE-Net (Goudet et al., 2018) contains 300 artificial cause-effect pairs generated using random distributions as causes, and neural networks as causal mechanisms. CE-Gauss contains 300 artificial causeeffect pairs as generated by Mooij et al. (2016) , using random mixtures of Gaussians as causes, and Gaussian process priors as causal mechanisms. CE-Multi (Goudet et al., 2018) contains 300 artificial cause-effect pairs built with random linear and polynomial causal mechanisms. In this dataset, simulated additive or multiplicative noise is applied before or after the causal mechanism.

The real-world datasets include the diabetes dataset by Frank & Asuncion (2010) , where causality is from Insulin → Glucose. Glucose curves and Insulin dose were analysed for 69 patients, each serves as a separate dataset. To match the literature protocols, the pairs are taken in an orderless manner, ignoring the time series aspect of the problem. Finally, the Tübingen cause-effect pairs dataset by Mooij et al. ( 2016) is employed. This dataset is a collection of 100 heterogeneous, handcollected, real-world cause-effect samples.

The autoencoder employed in our method is a fully-connected 5-layered neural network with 3 layers for the encoder and 2 layers for the decoder. The hyperparameters of this algorithm are the sizes of each layer, the activation function and the input dimension, i.e., length of sorted cuts (denoted by k in Sec. 2). Throughout the experiments, we noticed that the hyperparameter that has the highest influence is the input dimension. For all datasets results are stable in the range of 200 ≤ k ≤ 300, and we therefore use k = 250 throughout the experiments. For all dataset we employed the ReLU activation function, except the Tübingen dataset, where the sigmoid activation function produced better results (results are also reasonable with ReLU, but not state of the art).

In addition to our method, we also employ the entropy of each individual variable as a complexity method. This is done by binning the values of the variables into 50 bins. Other numbers of bins produce similar results.

Tab. 1 presents the mean AUC for each literature benchmark. As can be seen, the AEQ complexity measure produces reasonable results in comparison to the state of the art methods, indicating that the 1D SCM model can be overcome by comparing per-variable scores. On the popular Tübingen dataset, the one-sided AEQ computation outperforms all literature methods.

Tab. 2 presents accuracy rates for various methods on the Tübingen dataset, where such results are often reported in the literature. As can be seen, our simple method outperforms almost all other methods, including methods that employ supervised learning of the cause-effect relation.

4.2. MULTIVARIATE DATA:

We compared our method on several synthetic datasets. Each dataset consists of a list of pairs of real multivariate random variables (X, Y ) with their direction 1 or 0, depending on X → Y or Y → X (resp.). For each pair, we have a dataset of samples {(x i , y i )} m i=1 . We employ five datasets, covering multiple associations. Each dataset contains 300 artificial causeeffect pairs. The cause random variable is of the form X = h(z), where h is some function and

z ∼ N (0, σ 2 1 • I n ). The effect is of the form Y = g(u(X, E)), where E ∼ N (0, σ 2 2 • I n ) is independent of X,

u is a fixed function that combined the cause X and the noise term E and g is the causal mechanism. For each dataset, the functions h and g are taken from the same family of causal mechanisms H. Each pair of random variables is specified by randomly selected functions h and g.

MCE-:

Poly is generated element-wise polynomials composed on linear transformations as mechanisms and u(X, E) = X + E. MCE-Net pairs are generated using neural networks as causal mechanisms and u is the concatenation operator. The mechanism in MCE-SigMix consists of linear transformation followed by element wise application of q a,b,c (x) := ab(x + c)/1 + |b • (x + c)|, where a, b, c are random real valued numbers, which are sampled for each pair and x = x + e, where e is the environment random variable. In this case, u(X, E) = X + E.

In order to create each one of the synthetic datasets we took the standard synthetic data generators of Kalainathan & Goudet (2019) and extended them to produce multivariate causal pairs. We noticed that a-priori, the produced datasets are imbalanced in a way that the reconstruction error of a standard autoencoder on each random variable can be employed as a score that predicts the cause variable with a high accuracy. Therefore, in order to create balanced datasets, we varied the amount of noise (Hoyer et al., 2009) -59.5% CURE (Sgouritsa et al., 2015) -60.0% a GPI (Stegle et al., 2010) -62.6% PNL (Zhang & Hyvrinen, 2010) -66.2% CGNN (Goudet et al., 2018) -74.4% RECI (Bloebaum et al., 2018) -77.5% SLOPE (Marx & Vreeken, 2017 a The accuracy of CURE is reported on version 0.8 of the dataset in (Sgouritsa et al., 2015) as 75%. In Bloebaum et al. ( 2018 dimensions and their intensity until the autoencoder reconstruction error of both X and Y became similar. Note that for these multivariate variables, we do not use quantiles and use the variables themsevles. As the AutoEncoder reconstruction results in Tab. 3 show, in the MCE-SigMix dataset balancing was only partly successful.

In Tab. 3 we also compare our method to BivariateFit and ANM (Hoyer et al., 2009) algorithms on these datasets and present favorable performance. The other literature methods are designed for the univariate case and therefore, we were unable to extend them. The only baseline methods that are extendable to multi dimensions are PNL (Zhang & Hyvrinen, 2010) , CGNN (Goudet et al., 2018) and RECI (Bloebaum et al., 2018) . However, their runtime becomes too long. Specifically, However, RECI's runtime is of order O(n 3 ), where n is the input dimension. In addition to the synthetic dataset, we also employ the MOUS-MEG real world dataset, provided to us by the authors of (Anonymous, 2020) . This dataset is part of Mother Of Unification Studies (MOUS) dataset (Schoffelen et al., 2019) . This dataset contains magneto-encephalography (MEG) recordings of 102 healthy dutch-speaking subjects performing a reading task (9 of them were excluded due to corrupted data). Each subject was asked to read 120 sentences in Dutch, both in the right order and randomly mixed order, which adds up to a total of over 1000 words. Each word was presented on the computer screen for 351ms on average and was separated from the next word by 3-4 seconds. Each time step consists of 301 MEG readings of the magnetometers, attached to different parts of the head. For more information see (Schoffelen et al., 2019) . For each pair (X, Y ), X is the interval [-1.5s, -0.5s] relative to the word onset concatenated with the word embedding (using the spaCy python module with Dutch language model), this presents the subject in his "rest" state (i.e. the cause). Y is the interval [0, 1.0s] relative to the word onset, which presents the subject in his "active" state (i.e. the effect). The results, also given in Tab. 3 show a clear advantage over the literature methods.

To validate the soundness of the dataset, we ran made a few experiments on variations of the dataset and report the results in Tab. 4. As can be seen, a dataset where the cause consists of the word embedding and the effect consists of the subject's "active" state is highly imbalanced. This is reasonable since the word embedding and the MEG readings are encoded differently and are of different dimensions. In addition, when the cause is selected to be the "rest" state and the effect is the "active" state, the various algorithms are unable to infer which side is the cause and which one is the effect since the word is missing.

Emergence of independence:

To check the importance of our adversarial loss in identifying the direction of causality and capturing the implicit independent representation f (X) and E, we applied our method without training R against the discriminator. Therefore, in this case, the discriminator only serves as a test whether X and R(Y ) are independent or not. Recall that our analysis shows that the discriminator is not necessary under reasonable conditions. As can be seen in Tab. 5, the adversarial loss improves the results when there is no implicit emergence of independence, however, in cases where there is emergence of independence, the results are similar. This is verified in Tab. 5, which shows similar performance with and without the discriminator.

As mentioned in Sec. 3.1, the distance between c real + c fake to 1 indicates the amount of dependence between X and R(Y ). We denote by Ind C the mean values of |c real + c fake -1| over all pairs of random variables, epochs and samples when training our method in the causal direction. The same mean score when training in the anti-causal direction is denoted Ind E. As is evident from Tab. 5, the independence is similar between the two directions, emphasizing the importance of the reconstruction error in the score.

We noticed that the values of Ind C and Ind E are smaller for the full method. However, in MCE-Poly and MCE-SigMix they are still very small, and therefore, there is implicit emergence of independence between X and R(Y ), even without explicitly training R(Y ) to be independent of X. This can explain the fact that the drop in the AUC is significant but not drastic. 

5. SUMMARY:

We identify an inbalance in the complexities of cause and effect in the unidimensional SCM model and suggest a method to exploit it. Since the method does not consider the interactions between the variables, its success in predicting cause and effect indicates an inherent bias in the unidimensional datasets. Turning our attention to the multivariate case, where the complexity can be actively balanced, we propose a new method in which the learned networks models the underlying SCM model itself. Since the noise term E is unknown, we replace it by a function of Y that is enforced to be independent of X. We also show that under reasonable conditions, the independence emerges even without explicitly enforcing it.

A ARCHITECTURE:

The function G, F , R and D in the multivariate method in Sec. 3 are fully connected neural networks and their architectures are as follows: F is a 2-layered network with dimensions 100 → 60 → 50, R is a 3-layered network with dimensions 100 → 50 → 50 → 20, G is a 2-layers neural network with dimensions 50 + 20 → 80 → 100 (the input has 50 dimensions for F (X) and 20 for R(Y )).

The discriminator is a 3-layers network with dimensions 100 + 20 → 60 → 50 → 2 (the input is X and R(Y )). The activation function in all networks is the sigmoid function except the discriminator that applies the leaky ReLU activation. For all networks, the activation is not applied at the output layer.

B ANALYSIS B.1 TERMINOLOGY AND NOTATIONS:

We recall some relevant notations and terminology. The identity matrix of dimension n is denoted by I n . For a vector x = (x 1 , . . . , x n ) ∈ R n we denote x 2 := n i=1 x 2 i the Euclidean norm of x. For a differentiable function f : R m → R n and x ∈ R m , we denote by

EQUATION

the Jacobian matrix of f in x. For a twice differentiable function f : R m → R, we denote by

EQUATION

the Hessian matrix of f in x. Additionally, for a twice differentiable function f : R m → R n , f (x) = (f 1 (x), . . . , f n (x)), we denote the Hessian of f by H(f (x)) := (H(f 1 (x)), . . . , H(f n (x))).

For a scalar function f : R m → R instead of using the Jacobian notation, the gradient notation will be employed, ∇(f (x)) := J(f (x)).

B.2 AUTOENCODER AS A COMPLEXITY MEASURE:

Let F = {H d } ∞ d=1 be a family of classes of autoencoders A : R d → R d . Assume that the family F is closed to fixations, i.e., for any autoencoder A ∈ H d1+d2 and any fixed vector y * ∈ R d2 (x * ∈ R d1 ), we have: A(x, y * ) 1:d1 ∈ H d1 (A(x * , y) d1+1:d2 ∈ H d2 ). Here, v i:j = (v i , . . . , v j ). We consider that this is the typical situation when considering neural networks with biases.

Let X be a random variable. Let X be a multivariate random variable dimension d. We define the autoencoding complexity of X as follows:

EQUATION

Lemma 3. Let {H d } ∞ d=1 be a family of classes of autoencoders that is closed to fixations. The function C(X; F) is a complexity measure according to Sec. 3.

Proof. First, since (a, b) ≥ 0 for all a, b ∈ R k , this function is non-negative. Next, we would like to show that C α (X, Y ) ≥ max(C(X), C(Y )). Let A * be the minimizer of E x∼X [ (A * (x), x)] within H d1+d2 . We consider that there is a vector y * , such that,

E (x,y)∼(X,Y ) [ (A(x, y), (x, y))] ≥E y∼Y E x∼X [ (A(x, y), (x, y))] ≥E x∼X [ (A(x, y * ), (x, y * ))] ≥ E x∼X [ (A(x, y * ) 1:d1 , x)] (15)

We note that A(x, y * ) 1:d1 ∈ H d1 . Therefore,

EQUATION

By similar considerations, C(X, Y ; F).

B.3 PROOFS FOR THE RESULTS:

In this section we provide the proofs of the main results in the paper.

Lemma 1. In the setting of Sec. 2.2. Let θ * be the global optimum of Eq. 4 for a sufficiently large F. Assume that P is sufficiently large, such that:

∀z ∈ Z, φ ∈ Φ : q φ (u|z) ∈ P (5)

Then,

EQUATION

Proof. This is a corollary of the analysis of Zhao et al. (2017) . By their analysis, we have:

EQUATION

Since we selected p θ * (u|z) to be the PDF of a normal distribution N (f θ * (z), σ 2 1 • I n ), we have:

EQUATION

Therefore, since z distributed by q φ (z|u) can be represented as z|u

EQUATION

In addition, we have:

I q φ (z; u) = E u∼U [D KL (q φ (z|u) p(z))] = E u∼U [D KL (N (g φ (u), σ 2 2 • I d ) N (0, I d ))] = E u∼U 1 2 d • σ 2 2 + g φ (u) 2 2 + d log(σ 2 2 ) = d • (σ 2 2 -1 -log(σ 2 2 )) 2 + 1 2 E u g φ (u) 2 2 (20)

Therefore, we obtained the desired equation.

Lemma 2. Let C be a complexity measure of multivariate random variables (i.e, non-negative and satisfies Eq. 7). Then, there are triplets of random variables (X, E, Y ) and ( X, E, Y ), such that, Y = g(X, E), Y = g ( X, E), C(X) < C(Y ) and C( X) > C(Y ). Therefore, C cannot serve as a score for causal inference.

Proof. Let X be a random variable and E |= X, such that, Y = g(X, E). Assume that C(X) < C(Y ). Then, let X be a random variable independent of X, such that, C(X ) > C(Y ). Then, we have: C(X, X ) > C(Y ) and we have: Y = g (X, X , E), for g (a, b, c) = g(a, c).

The following lemma is an extension of Thm. 1 in (Zhang & Hyvrinen, 2010) to real valued random variables of dimension > 1.

Lemma 4. Assume that (X, Y ) can be described by both:

Y = g 1 (f 1 (X) + E 1 ), s.t: X |= E 1 and g 1 is invertible (21) and X = g 2 (f 2 (Y ) + E 2 ), s.t: Y |= E 2 and g 2 is invertible (22)

Assume that g 1 and g 2 are invertible and let:

T 1 := g -1 1 (Y ) and h 1 := f 2 • g 1 T 2 := g -1 2 (X) and h 2 := f 1 • g 2 (23)

Assume that the involved densities p T2 , p E1 and nonlinear functions f 1 , g 1 and f 2 , g 2 are third order differentiable. We then have the following equations for all (X, Y ) satisfying:

EQUATION

where η 1 (t 2 ) := log p T2 (t 2 ) and η 2 (e 1 ) := log p E1 (e 1 ).

Proof. The proof is an extension of the proof of Thm. 1 in (Zhang & Hyvrinen, 2010) . We define:

EQUATION

Since g 2 is invertible, the independence between X and E 1 is equivalent to the independence between T 2 and E 1 . Similarly, the independence between Y and E 2 is equivalent to the independence between T 1 and E 2 . Consider the transformation F : (E 2 , T 1 ) → (E 1 , T 2 ):

EQUATION

The Jacobian matrix of this transformation is given by:

EQUATION

Since I commutes with any matrix, by Thm. 3 in (Silvester, 1999) , we have:

det(J(F (E 2 , T 1 ))) = det (-J(h 2 (T 2 )) • J(h 1 (T 1 )) -I • (I -J(h 2 (T 2 )) • J(h 1 (T 1 )))) = 1

(28) Therefore, we have: p T2 (t 2 ) • p E1 (e 1 ) = p T1,E2 (t 1 , e 2 )/| det J| = p T1,E2 (t 1 , e 2 ). Hence, log(p T1,E2 (t 1 , e 2 )) = η 1 (t 2 ) + η 2 (e 1 ) and we have:

EQUATION

Therefore,

EQUATION

The independence between T 1 and E 2 implies that for every possible (t 1 , e 2 ), we have:

∂ 2 log p T 1 ,E 2 (t1,e2) ∂e2∂t1 = 0.

Lemma 5 (Reduction to post-linear models). Let f

(x) = σ 1 (W d . . . σ 1 (W 1 x)) and g(u, v) = σ 2 (U k . . . σ 2 (U 1 (u, v))) be two neural networks. Then, if Y = g(f (X), E) for some E |= X, we can represent Y = ĝ( f (X) + N ) for some N |= X. Proof. Let f (x) = σ 1 (W d . . . σ 1 (W 1 x)) and g(u, v) = σ 2 (U k . . . σ 2 (U 1 (u, v)

)) be two neural networks. Here, (u, v) is the concatenation of the vectors u and v. We consider that U

k (f (X), E) = U 1 1 f (X) + U 2 1 E.

We define a noise variable N := U 2 1 E and have: X |= N . In addition, let f (x) := U 1 1 f (x) and ĝ(z) := σ 2 (U k . . . σ 2 (U 2 σ 2 (z))). We consider that: Y = ĝ( f (X) + N ) as desired.

Theorem 1 (Identifiability of neural SCMs). Let P X,Y admit a neural SCM from X to Y as in Eq. 1, such that p X , and the activation functions of f and g are three-times differentiable. Then it admits a neural SCM from Y to X only if p X , f , g satisfy Eq. 34 in the Appendix.

Proof.

Let f i (z) = σ 1 (W i,d . . . σ 1 (W i,1 z)) and g i (u, v) = σ 2 (U i,k . . . σ 2 (U i,1 (u, v))) (where i = 1,

2) be pairs of neural networks, such that, σ 1 and σ 2 are three-times differentiable. Assume that:

EQUATION

for some E 1 |= X and E 2 |= Y . By Lem. 5, we can represent

Y = ĝ1 ( f1 (X) + N 1 ),

where

EQUATION

and also,

X = ĝ2 ( f2 (Y ) + N 2 ),

where

EQUATION

From the proof of Lem. 5, it is evident that the constructed ĝ1 , f1 and ĝ1 , f2 are three-times differentiable whenever σ 1 and σ 2 are. Therefore, by Thm. 1, the following differential equation holds:

EQUATION

where T 1 := ĝ-1 1 (Y ) and h 1 := f2 • ĝ1 T 2 := ĝ-1 2 (X) and h 2 := f1 • ĝ2

and η 1 (t 2 ) := log p T2 (t 2 ) and η 2 (n 1 ) := log p N1 (n 1 ).

Theorem 2 (Uniqueness of Representation). Let P X,Y admit a nonlinear model from X to Y as in Eq. 1, i.e., Y = g(f (X), E) for some random variable E |= X. Assume that f and g are invertible. Let G, F and R be functions, such that, L err := E (x,y)∼(X,Y ) [ G(F (x), R(y)) -y 2 2 ] = 0 and G and F are invertible functions and X |= R(Y ). Then, F (X) ∝ f (X) and R(Y ) ∝ E.

Proof. Since F and f are invertible, one can represent: F (X) = F (f -1 (f (X))) and f (X) = f (F -1 (F (X))). Similarly, since G and g are invertible, we also have: (F (X), R(Y )) ∝ (f (X), E). Since (F (X), R(Y )) ∝ (f (X), E) and F (X) ∝ f (X), we have: R(Y ) = Q(F (X), E). However, R(Y ) |= F (X) and therefore, we can represent R(Y ) = P (E) and vice versa.

Theorem 3 (Emergence of independent representations). Let P X,Y admits a nonlinear model from X to Y as in Eq. 1, i.e., Y = g(f (X), E) for some random variable E |= X. Assume that X and Y are discrete random variables. Let G, F and R be functions, such that, L err := E (x,y)∼(X,Y ) [ G(F (x), R(y)) -y 2 2 ] = 0 and G is an invertible function. Assume that h(R(Y )) ≤ h(E) and that g is invertible. Then, we have: F (X) |= R(Y ), F (X) ∝ f (X) and R(Y ) ∝ E.

Proof. Since g and f are invertible, we have: h(Y ) = h(f (X)) + h(E) ≤ h(X) + h(E). In addition, we have: h(R(Y )) ≤ h(E). Therefore, since G is invertible, h(Y ) = h(F (X), R(Y )) ≤ h(F (X)) + h(R(Y )) ≤ h(X) + h(E) = h(Y ). Hence, F (X) and R(Y ) are independent.

In the following three lemmas we extend Thms. 2 and 3. Lemma 6. Let P X,Y admit a nonlinear model from X to Y as in Eq. 1, i.e., Y = g(f (X), E), for some random variable E |= X. Assume that f and g are invertible. Assume that X and Y are discrete and Y is over a set Y, such that, ∀y 1 = y 2 ∈ Y, we have, y 1 -y 2 2 > ∆ > 0. Let G, F and R be functions, such that, L err := E (x,y)∼(X,Y ) y -G(F (x), R(y)) 2 2 < ∆/2 and F is invertible. Then, F (X) ∝ f (X) and

I(R(Y ); E) ≥ 1 -Lerr ∆ • h(E) -Lerr ∆ .

Proof. First, since both F and f are invertible, F (X) ∝ f (X). By the proof of Lem. 10 in (Press et al., 2019) , there is a function r, such that: Therefore, by Lem. 6 in (Press et al., 2019) , we have:

P (x,

EQUATION

By the analysis in the proof of Lem. 10 in (Press et al., 2019) , we derive that:

EQUATION

Finally, since F (X) is a function of X which is independent of E, we obtain the desired inequality:

EQUATION

We mention that in Lem. 6, the function R(Y ) can hold all of the information present in Y . Therefore, in order to suffice that R(Y ) holds only the information present in E, one can restrict that h(R(Y )) ≤ h(E) + as will be shown in the next lemma. We note that under the conditions of Lem. 6, h(R(Y )) ≤ h(E) + 1 is equivalent to I(F (X); R(Y )) ≤ 2 for 1 , 2 that are functions of each other.

Lemma 7. In the setting of Lem. 6. Assume that h(R(Y )) ≤ h(E) + , for some constant > 0.

Then, I(F (X); R(Y )) ≤ Lerr ∆ • h(Y ) + Lerr ∆ + . In addition, there are functions r 1 , r 2 , such that,

EQUATION

and

EQUATION

Proof. By Lem. 10 in Press et al. (2019) , we have:

EQUATION

In addition, by the data processing inequality and the assumption h(R(Y )) ≤ h(E) + , we have:

EQUATION

Finally, by Feder & Merhav (1994) , there is a function r 1 , such that,

EQUATION

In addition, there is a function r 2 , such that,

EQUATION

as desired.

:



Figures and Tables:

) they reran this algorithm and achieved an accuracy rate around 60%. b The accuracy scores reported in (Lopez-Paz et al., 2015) are for version 0.8 of the dataset, in (Lopez-Paz et al., 2017) they reran RCC (Lopez-Paz et al., 2015) on version 1.0 of the dataset.

y)∼(X,Y ) [r(F (x), R(y)) = y] ≥ 1 -L err ∆ ≥ 0.5(36)Since g is invertible, there is a function u(Y ) = E. In particular,P (x,y)∼(X,Y ) [u(r(F (x), R(y))) = u(y)] ≥ 1 -L err ∆ ≥ 0.5(37)

Mean AUC rates of various baselines on different one dimensional cause-effect pairs datasets. Our simple AEQ algorithm achieves competitive results on most datasets.

Table content: <table><tr><td>Dataset</td></tr></table>

Accuracy rates of various baselines on the CE-Tüb dataset. Our simple algorithm AEQ achieves almost SOTA accuracy.

Table content: <table><tr><td>Method</td><td>Supervised</td><td>Acc</td></tr><tr><td>LiNGAM (Shimizu et al., 2006)</td><td>-</td><td>44.3%</td></tr><tr><td>BivariateFit</td><td>-</td><td>44.9%</td></tr><tr><td>Entropy as a complexity measure</td><td>-</td><td>52.5%</td></tr><tr><td>IGCI (Daniusis et al., 2012)</td><td>-</td><td>62.6%</td></tr><tr><td>CDS (Fonollosa, 2016)</td><td>-</td><td>65.5%</td></tr><tr><td>ANM</td><td/><td/></tr></table>

Mean AUC rates of various baselines on different multidimensional dimensional causeeffect pairs datasets. The datasets are designed and balanced, such that an autoencoder method would fail. Our method achieves SOTA results.

Table content: <table><tr><td>Method</td><td colspan="4">MCE-Poly MCE-Net MCE-SigMix MOUS-MEG</td></tr><tr><td>BivariateFit</td><td>54.7</td><td>48.4</td><td>48.2</td><td>44.2</td></tr><tr><td>ANM (Hoyer et al., 2009)</td><td>52.2</td><td>51.1</td><td>46.4</td><td>52.4</td></tr><tr><td>PNL (Zhang &amp; Hyvrinen, 2010)</td><td>76.4</td><td>54.7</td><td>16.8</td><td>56.3</td></tr><tr><td>CGNN (Goudet et al., 2018)</td><td>47.8</td><td>67.8</td><td>58.8</td><td>40.9</td></tr><tr><td>AE reconstruction</td><td>57.2</td><td>42.4</td><td>22.3</td><td>41.2</td></tr><tr><td>Our method</td><td>95.3</td><td>70.0</td><td>98.5</td><td>87.6</td></tr></table>

Results of various methods on different variations of the MOUS-MEG dataset.

Table content: <table><tr><td>Dataset</td></tr></table>

Emergence of independence. Ind C (Ind E) is the mean of |c real + c fake -1| over all pairs of random variables, epochs and samples, when training the method from X to Y (vice versa). w/o backprop means without backpropagating gradients from D to R.

Table content: <table><tr><td/><td/><td>Full method</td><td/><td colspan="3">w/o backprop</td></tr><tr><td>Dataset</td><td colspan="6">AUC Ind C Ind E AUC Ind C Ind E</td></tr><tr><td>MCE-Poly</td><td>95.3</td><td>0.06</td><td>0.05</td><td>95.1</td><td>0.10</td><td>0.10</td></tr><tr><td>MCE-Net</td><td>70.0</td><td>0.16</td><td>0.20</td><td>65.1</td><td>0.55</td><td>0.55</td></tr><tr><td colspan="2">MCE-SigMix 98.5</td><td>0.05</td><td>0.06</td><td>98.8</td><td>0.16</td><td>0.20</td></tr><tr><td colspan="2">MOUS-MEG 87.6</td><td>0.61</td><td>0.61</td><td>80.7</td><td>0.74</td><td>0.75</td></tr></table>

Peer review (iclr2020):
Edit after author rebuttal and author additions: I have updated my score from a weak reject (3) to a weak accept (6). Justification: 1. The authors have pointed out that I misunderstood one of their contributions: pointing out that they are demonstrating that the univariate case is an insufficient setting to test causal discovery methods because it can be done without even looking at the conditional distributions (just the marginals). This contribution seems important to orient future work. They have made this more clear in their recent upload and would likely make it even more clear in a camera-ready version. 2. The authors have made their contribution to the multivariate setting more substantial by adding evaluation on the MOUS-MEG real-world dataset and have better positioned their work relative to others by adding comparisions to multivariate extensions of PNL and CGNN. ==================================================================================================== Original Review: Summary: The authors focus on the problem of inferring whether the causal structure X —> Y or Y —> X. They first consider the case where both X and Y are scalar (univariate) random variables and then consider the case where X and Y are vector-valued (multivariate) random variables. In the scalar case, motivated by the idea that the effect could be less entropic than the cause (due to data processing inequality), they introduce a method based on comparing reconstruction losses of X and Y and show competitive results in Tables 1 and 2. They establish that this method is not sufficient for the multivariate case in Lemma 2 and move to a new method for the multivariate case. They prove identifiability for this new method in for the multivariate case in Section 4.2 and claim state-of-the-art (SOTA) results in Table 3.? Main contributions: - Presents a causal discovery technique for the univariate cases that only examines the marginal distributions of X and Y and seems fairly competitive (Tables 1 and 2) - Extends the post-nonlinear identifiability analysis of Zhang & Hyv¨arinen (2009) from scalars to vectors and proved that their method will actually identify the correct causal direction - Demonstrates competitive experimental results for both their univariate method - Claims SOTA results for their multivariate method Decision: I lean toward rejecting this paper because 1) I have several questions about the univariate case (see below) that would need to be resolved before I lean toward accept, 2) although I am not too familiar with the literature, I believe that this paper may be missing key related work that also uses independence testing for causal discovery (see, e.g., Heinze-Deml et al. (2017)’s Invariant Causal Prediction for Nonlinear Models), and 3) I am not yet convinced that the comparison done in Table 3 is fair and exhaustive.?? Sufficient reason to accept: If the theorems in Section 4.2 are found to checkout, and the SOTA results in Table 3 are found to be fair, exhaustive comparisons to the previous SOTA, their contribution to the multivariate case would seem to be sufficient for acceptance. I believe more discussion between the authors and reviewers is necessary here. ?Questions about univariate case: ?1. The motivation for the first method (entropy decreasing along a Markov chain due the data processing inequality) seems to only be valid when Y := f(X), but not necessarily when Y := f(X) + E. For example, let f be the identity function and E be independent to X. How did you resolve this argument against the intuition?? 2. Also, I thought the data processing inequality relates mutual information between variables, not necessarily their entropies. Can you make this connection more clear?? Context for questions 3 and 4: In Section 3.1, you write, “estimating the entropy of each random variable from its samples does not present a consistent difference between the entropies h(X) and h(Y). Our method, therefore, computes an alternative complexity score for X and, independently, for Y.” You then go on to link the entropy to the reconstruction error (your method) in Lemma 1 and show competitive results in Tables 1 and 2. 3. Why do you want to link the reconstruction error to entropy if you found a purely entropy-based method did not work? 4. Why did the purely entropy-based method not work while your method worked if the two are linked? Questions about multivariate case: 5. Are you certain that BivariateFit and ANM are the only models that you should be comparing against for this multivariate setting? 6. What is CGNN’s runtime? Would you be able to compare against CGNN in time for a potential camera-ready version of this paper?