Year: 2018

Abstract:
Learning causal relations from observational time series with nonlinear interactions and complex causal structures is a key component of human intelligence, and has a wide range of applications. Although neural nets have demonstrated their effectiveness in a variety of fields, their application in learning causal relations has been scarce. This is due to both a lack of theoretical results connecting risk minimization and causality (enabling function approximators like neural nets to apply), and a lack of scalability in prior causal measures to allow for expressive function approximators like neural nets to apply. In this work, we propose a novel causal measure and algorithm using risk minimization to infer causal relations from time series. We demonstrate the effectiveness and scalability of our algorithms to learn nonlinear causal models in synthetic datasets as comparing to other methods, and its effectiveness in inferring causal relations in a video game environment and real-world heart-rate vs. breath-rate and rat brain EEG datasets.

text: Learning causal relations from observational time series with nonlinear interactions and complex causal structures is a key component of human intelligence, and has a wide range of applications. Although neural nets have demonstrated their effectiveness in a variety of fields, their application in learning causal relations has been scarce. This is due to both a lack of theoretical results connecting risk minimization and causality (enabling function approximators like neural nets to apply), and a lack of scalability in prior causal measures to allow for expressive function approximators like neural nets to apply. In this work, we propose a novel causal measure and algorithm using risk minimization to infer causal relations from time series. We demonstrate the effectiveness and scalability of our algorithms to learn nonlinear causal models in synthetic datasets as comparing to other methods, and its effectiveness in inferring causal relations in a video game environment and real-world heart-rate vs. breath-rate and rat brain EEG datasets.

heading: 1 INTRODUCTION
text: From an early age, humans have a remarkable ability to infer causal relations from pure observations (White & Milne (1997); Scholl & Tremoulet (2000); Buchsbaum et al. (2012)). By observing that a left dot moving towards a right dot and the right dot moves correspondingly (the launching effect, see Michotte (1963)), a human can quickly infer that the left dot causes the right dot to move (Scholl & Tremoulet (2000)). In fact, much of our way of thinking is via cause and effect. Causal analysis also permits counterfactual reasoning, answering what would have happened if the cause had happened differently. The learning of causality also constitutes much of the scientific endeavor, for example finding the cause of a certain cancer (Bosch et al. (2002)), or discovering gene regulatory networks (Lozano et al. (2009)). In addition, causality plays a key role in neuroscience (Neves et al. (2008); Seth et al. (2015)), economics (e.g. Granger (1969); Stock & Watson (1989)) and finance (Hiemstra & Jones (1994); Granger et al. (2000)).
The study of causality has a long history, yet the application of neural nets to learning causal relations has been scarce. There have been various works that propose methods to infer causal structures with limited model space (e.g. Granger (1969)) which may not be able to model complex nonlinear causal relations, or propose measures to quantify causal strength (e.g. Schreiber (2000); Janzing et al. (2013)) that may not scale to high-dimensional data. For these, the use of universal function approximators like neural nets may be beneficial in inferring causality, which motivate us to propose causal measures that are not only theoretically founded, but also amenable to the learning of function approximators. On the other hand, in the deep learning community, the learning of causal models has not been prevalent, in part due to the lack of theoretical understanding between learning a prediction model and learning a causal model. This also motivates us to propose causal measures, obtained via learning a prediction model, that can deduce causality.
The contributions of this work are as follows:
• We propose a novel measure to quantify causality from observational data, and an effective algorithm, Causal Inference with Learnable Noise (CILN), to estimate it. It is based on minimizing a learnable noise risk of a prediction model, allowing function approximators such as neural nets to learn complex causal relationships.
• We demonstrate on nonlinear synthetic datasets that our method outperforms other causal measures by a large margin, and can scale to a larger number of time series. We also demonstrate that our models are effective on real-world datasets.

heading: 2 RELATED WORK
text: The study of causality has a long history, and has been approached from different perspectives. Pearl (e.g. Pearl (2002; 2009); Pearl et al. (2009)) defines causality in terms of intervention and structural dependence, under the structural equation models (SEM). Granger (Granger (1969); Granger & Newbold (1986)) defines causality via prediction: if the prediction of Y via a linear model can be improved by including the information of X, then X causes Y in the Granger sense. Since its proposal, Granger causality has been widely applied in economics (e.g. Joerding (1986)) and neuroscience (e.g. Deshpande et al. (2009); Seth et al. (2015)). To learn nonlinear causal relations, later works also extend Granger causality to kernel methods (Ancona et al. (2004); Marinazzo et al. (2008a;b); Sindhwani et al. (2012)). To clear up the relations between these two notions of causality, White et al. (2011) provide conditions under which Granger causality can deduce direct structural causality in a general settable system framework (White & Chalak (2009)), with direct structural causality as a natural extension of Pearl causality in settable systems. Our method also utilizes the high-level idea of inferring causal relations via prediction, and building on the work of White et al. (2011) we propose a novel measure that can likely uncover direct structural causality under certain conditions.
Numerous methods have been proposed to discover causal structures from data. One important class is constraint-based methods, for example PC (Spirtes et al. (2000)), rankPC (Harris & Drton (2013)), IC (Pearl (2002)), and FCI (Spirtes et al. (2000)), which require repetitive conditional independence tests. Hybrid methods, e.g. MMHC (Tsamardinos et al. (2006)), require repetitive estimation of conditional association score (e.g. using conditional mutual information). In comparison, under the scope of time series, our method can simultaneously discover multiple variables that directly cause the variable of interest. Score-based methods search for the structure that yields the optimal score w.r.t. the data, generally using greedy search methods, for example GES (Chickering (2002)), rankGES (Nandy et al. (2018)) and GIES (Hauser & Bühlmann (2012)). This in general requires Θ(N2) steps (N denoting the number of nodes in the graph), and the number of neighboring states may grow very large at each step. In comparison, our method only requires training N models for causal discovery with time series.
In addition, various measures have been proposed to quantify causality. Schreiber (2000) proposes transfer entropy as a measure for causality. It measures the mutual information between the current Y and the past of X, conditioned on the past of Y, to quantify the directional information transfer. As noted in Marinazzo et al. (2008a), Granger causality as defined in Granger & Newbold (1986) implies nonzero transfer entropy. Janzing et al. (2013) analyze several causal measures, and conclude that they are unsatisfactory measures of causal strength. They then propose causal influence, defined via the KL-divergence between the original joint distribution and the distribution with a set of causal arrows broken. We note that both the calculation of transfer entropy and causal influence requires density estimation of the full joint distributions of the input and output, which could easily become difficult in high dimensions. This motivates us to propose new causal measures based on risk minimization, which is much easier with high input dimensions. Besides, various other works have also proposed methods to infer causal structure under some specific conditions, for example, causal additive noise models (Hoyer et al. (2009)), information-geometric causal inference (Daniusis et al. (2012); Janzing et al. (2012)), dynamic causal modeling (Friston et al. (2003)), etc.
There has also been works that approach causality from a machine learning perspective, or have architectures that put causality to mind. Lopez-Paz et al. (2015) study causal inference as a supervised learning problem, where the pairwise causal directions are given as labels for training. Louizos et al. (2017) utilize a variational autoencoder (VAE) structure to estimate the unknown latent space summarizing the confounders and the causal effect. Kipf et al. (2018) propose a neural relational inference architecture, which simultaneously infers the interactions and learning the dynamics from observational data. Designing causal inference into the model architecture can also be beneficial in certain applications. For example, Kansky et al. (2017) propose a Schema Network that allows re-
gression planning from a goal through causal chains, and demonstrates zero-shot learning in a suite of variations of Atari Breakout games.
Our method also relates to sparse learning/feature selection methods, for example L1 regularization and group L1 regularization (e.g. Meier et al. (2008); Scardapane et al. (2017)). Although L1 or group L1 regularization encourage sparsity in the weights of neural nets, the regularization is model and input dependent, in contrast to our method’s invariance to model structure change and inputs rescaling. For example, a three-layer linear network and a collapsed one-layer linear network representing the same function can induce different L1 regularizations. For another example, if one of the input time series is scaled by a factor of 0.1 while other variables remain unchanged, the L1 regularization will be different due to a differently learned model. In comparison, our learnable noise risk when minimized is invariant to the above two kinds of changes, making it especially suitable to discover causality where the scale of different time series may span orders of magnitude, and the model structure may vary, as demonstrated in Section 4.1.

heading: 3 METHOD

heading: 3.1 PROBLEM DEFINITION
text: Although causal data can be inferred from observations at individual time points, we are primarily concerned with time series. Time series have extra structure that is particularly useful for causal inference: causes must precede their effects. Therefore, we consider N time series x(1), x(2), ...x(N), where each time series x(i) = (x(i)1 , x (i) 2 , ...x (i) t , ...) and each x (i) t ∈ RM is an M -dimensional vector. Denote X(i)t−1 = (x (i) t−K , x (i) t−K+1, ...x (i) t−1) with maximum time horizon of K, and Xt−1 = {X(i)t−1}, i = 1, 2, ...N . We also denote X (ĵ) t−1 = Xt−1\X (j) t−1, i.e. Xt−1 excluding X(j)t−1, to notationally differentiate with the variable of interest X (j) t−1. We assume that the time series is generated by a canonical settable system (White et al. (2011)). We adopt the settable system (White & Chalak (2009)) paradigm due to the following reasons: (1) as a natural extension to SEMs, it facilitates optimization, equilibrium, and learning; (2) it can formally link Granger causality with Pearl causality; (3) it is general enough to encompass a large number of practical scenarios, including time series. Intuitively, each variable in a settable system can either be determined by direct setting (which formalizes interventions), or by a response function on the settings of other variables. A canonical settable system is a settable system where each variable’s setting equals its response, allowing the system to evolve naturally without intervention, formalizing time series. In this paper, we also assume causal sufficiency (Peters et al. (2017)), i.e., each time series x(i) can only be structurally caused by the time series from x(1), x(2), ...x(N), and generated by stationary response functions hi that are unknown to the learner: x (1) t := h1(Xt−1, u1) x (2) t := h2(Xt−1, u2) ...
x (N) t := hN (Xt−1, uN )
(1)
for t = K + 1,K + 2, ... . Here ui ∈ RM , i = 1, 2, ...N are noise variables that are mutually independent, are independent of any X(i)t−1, x (i) t , i ∈ {1, 2, ...N}, and are effective arguments of the response functions hi. Also in this paper, we only consider “causality in mean", i.e. the causal relations, if exists, influence the mean value of other variables. For any i, j ∈ {1, 2, ...N}, we assume that the variables (X(ĵ)t−1, X (j) t−1, x (i) t ) have probability density function P (X (ĵ) t−1, X (j) t−1, x (i) t ). We will leave time series with hidden variables (therefore confounding can occur) for future work. We note that even without considering hidden variables, Eq. 1 is very general and already encompasses a wide range of scenarios.
In order to state the goal of the learner who wants to learn causality from the observations x(1), x(2), ...x(N), we need to rigorously define causality. Here, we restate the definitions of direct structural causality (White et al. (2011)) and Granger causality (Granger & Newbold (1986)) using our notations of the system Eq. 1, the former being a natural extension to Pearl causality (Pearl (2009)) in the settable systems.
Direct structural causality (White et al. (2011)) We say X(j)t−1, j 6= i does not directly structurally cause x(i)t , if for all possible values of X (ĵ) t−1 and ul, l ∈ 1, 2, ...N , the function X (j) t−1 → hi(Xt−1, ui) is constant in X (j) t−1. Otherwise, we say X (j) t−1 directly structurally causes x (i) t .
Granger (1969) defines causality in terms of conditional expectations. Later, Granger & Newbold (1986) define it using conditional distributions. We use the latter definition, since it is more general, and also facilitates connection with direct structural causality, and the algorithms proposed in this paper. In this work, we only consider causality between different time series, i.e. requiring j 6= i as default.
Granger causality (Granger & Newbold (1986)) We say X(j)t−1, j 6= i does not Granger-cause x (i) t , if P (x(i)t |X (j) t−1,X (ĵ) t−1) = P (x (i) t |X (ĵ) t−1), i.e. the conditional distribution function of x (i) t given X (j) t−1,X (ĵ) t−1 is identical to the conditional distribution function of x (i) t given X (ĵ) t−1. Otherwise, we say X(j)t−1 Granger-causes x (i) t .
The goal of the learner is, given only x(1), x(2), ...x(N), determine for any x(i)t , whether X (j) t−1 directly structurally causes x(i)t for each j = 1, 2, ...N .

heading: 3.2 GRANGER CAUSALITY IMPLIES DIRECT STRUCTURAL CAUSALITY IN EQ. (1)
text: For our system Eq. (1), applying the results by White et al. (2011), we have that Granger causality is a sufficient condition for direct structural causality. See Appendix A for the detailed theorem and proof.
For system Eq. 1, for any i, j ∈ {1, 2, ...N}, i 6= j, if X(j)t−1 Granger-causes x (i) t , then X (j) t−1 directly structurally causes x(i)t .
The reason that here Granger causality can deduce direct structural causality is in part due to the fact that for system Eq. (1), conditional exogeneity (White et al. (2011)) is automatically satisfied by the assumptions of the system.
Note that the reverse of the statement is not true, i.e. Granger non-causality does not necessarily imply direct structural non-causality (White & Lu (2010) give several examples). They also note that these instances are exceptional, in that mild perturbations to their structures destroy the Granger non-causality.

heading: 3.3 OUR METHOD
text: Based on section 3.2, if we have an algorithm that can deduce Granger causality in system Eq. (1), then we can immediately deduce direct structural causality. In general, a direct test of Granger causality is difficult. In particular, when the number of time series is large or the dimension M of each time series is large, density or mutual information-based methods like transfer entropy (I(x(i)t ;X (j) t−1|X (ĵ) t−1), see Schreiber (2000)), may not give a good estimate. Alternatively, various works have resorted to testing whether the prediction of x(i)t given X (j) t−1,X (ĵ) t−1 is better than given only X(ĵ)t−1, under some limited functional space. For example, in his original work, Granger (1969) investigates causality with linear function predictors. Later works have extended it to kernel methods, e.g., Ancona et al. (2004); Marinazzo et al. (2008a;b); Sindhwani et al. (2012), which essentially estimate linear Granger causality on the feature space of the kernel.
To enable causal learning with potentially highly nonlinear response functions, it may be desirable to use universal function approximators (Hornik (1991)) such as neural nets. As the main contribution of this paper, we will provide a novel causal measure and corresponding algorithm to estimate it. The measure is based on optimizing an objective containing a function approximator, thus allowing neural nets to apply.
Our algorithm is inspired by asking a counterfactual question during the learning of a prediction model. Specifically, it asks:
Algorithm 1 Causal Inference with Learnable Noise (CILN)
Require x(i)t ,Xt−1, for i ∈ {1, 2, ...N}, t ∈ T = {K + 1,K + 2, ...}. Require η0: a small value for initialization of η. Require λ: coefficient for the mutual information term. 1: for i in {1, 2, ...N} do: 2: Initialize function approximator fθ. 3: Initialize η = (η1, η2, ...ηN ) = (η01, η01, ...η01), where each element η01 is a KM− dimensional vector, same dimension as X(j)t−1. 4: (fθ∗ ,η∗)← Minimize(fθ,η)R̂X,x(i), [fθ,η] (Eq. 4) with e.g. gradient descent 5: Wji ← I(X̃ (j)(η∗j ) t−1 ;X (j) t−1), for j = 1, 2, ...N . 6: end for 7: return W
How much noise can I add to X(j)t−1, without making the best prediction of x (i) t worse?
To give a quantitative answer to this question, we define a learnable noise risk:
RX,x(i) [fθ,η] = EXt−1,x(i)t ,
[( x
(i) t − fθ(X̃ (η) t−1)
)2] + λ · N∑ j=1 I(X̃ (j)(ηj) t−1 ;X (j) t−1) (2)
where X̃(η)t−1 := Xt−1 + η (or element-wise, X̃ (j)(ηj) t−1 := X (j) t−1 + ηj · j , j = 1, 2, ...N ) is the noise-corrupted inputs with learnable noise amplitudes ηj ∈ RKM , and j ∼ N(0, I). λ > 0 is a positive hyperparameter for the mutual information I(·, ·). Intuitively, the minimization of the second term I(X̃(j)(ηj)t−1 ;X (j) t−1) requires the noise amplitude ηj to go up. The minimization of the first term requires the noise amplitude ηj to go down, and the larger causal strength from X (j) t−1 to x (i) t , the larger this force. The minimization of the two terms strikes a balance, at which point the I(X̃
(j)(ηj) t−1 ;X (j) t−1) measures how many bits of information does time series j need to provide to the
learner, without making the prediction worsened. Thus, we propose to use
Wji = I(X̃ (j)(η∗j ) t−1 ;X (j) t−1) (3)
as another measure for causality, where (fθ∗ ,η∗) = argmin(fθ,η)RX,x(i) [fθ,η] 1. In Appendix B, we analyze the qualitative and quantitative properties of the learnable noise risk, and give intuitions why it is likely to select the variables that directly structurally causes x(i)t .
Empirically, we minimize the following empirical risk:
R̂X,x(i), [fθ,η] = 1 |T| ∑ t∈T ( x (i) t − fθ(X̃ (η) t−1) )2 + λ N∑ j=1 I(X̃ (j)(ηj) t−1 ;X (j) t−1) (4)
In general, it may be difficult to estimate the mutual information I(X̃(j)(ηj)t−1 ;X (j) t−1) with large dimension of X(j)t−1 such that the expression is also differentiable w.r.t. ηj . Utilizing the property of Gaussian channels, in Appendix C we prove that I(X̃(j)(ηj)t−1 ;X (j) t−1) ≤
1 2 ∑KM l=1 log ( 1 + Var(X(j)t−1,l) η2j,l ) , where l denotes the lth element of a vector, and Var(X(j)t−1,l) is the
variance of X(j)t−1,l across t. Therefore, in practice to improve efficiency, we can optimize an upper bound of the learnable noise risk:
R̂upper X,x(i), [fθ,η] = 1 |T| ∑ t∈T ( x (i) t − fθ(X̃ (η) t−1) )2 + λ 2 N∑ j=1 KM∑ l=1 log ( 1 + Var(X(j)t−1,l) η2j,l ) (5)
1Note also that throughout this paper, when talking about causal matrices, the (j, i)th element always denotes the causal strength from j to i.
When the dimension of X(j)t−1 is large, differentiable estimate of the mutual information, such as MINE (Belghazi et al. (2018)), can be applied. We provide Algorithm 1 to empirically estimate Wji, which we term Causal Inference with Learnable Noise (CILN).

heading: 4 EXPERIMENTS
text: To demonstrate that our proposed method works, we test it on both synthetic and real datasets. We first use synthetic datasets, where we know the underlying causal structure and compare with previous causal measures. We then test whether our algorithm can infer causal structure from watching an agent playing video games. Finally, we apply our algorithm to real-world heart-rate vs. breath-rate and rat EEG datasets to test its effectiveness. We use the R̂upper
X,x(i), [fθ,η] (Eq. 5) for optimization
for all experiments. The metrics we use are the standard metrics of area under the precision-recall curve (AUC-PR) (Davis & Goadrich (2006)), and area under the ROC curve (AUC-ROC).

heading: 4.1 SYNTHETIC EXPERIMENT WITH LOG-NORMAL CAUSAL STRENGTHS
text: In this experiment, we evaluate our method together with other methods with a nonlinear synthetic dataset generated to have a known causal structure (hidden to the methods being compared). We study how they perform with varying number N of time series, with N up to 30. To generate the data, we let each x(i)t have dimension M = 1, and also set the maximum time horizon K = 3, so each X(j)t−1 is a K ×M = 3 × 1 matrix. We use the following realization of the response function hi in Eq. (1):
x (i) t = hi(Xt−1, ut) = H1  N∑ j=1 [ Aji H2(Bj X(j)t−1) ]+ ut, i = 1, 2, ...N (6) where ut ∼ N(0, I) ∈ RM , denoting element-wise multiplication, and H1 and H2 are two nonlinear functions to make the response functions nonlinear. In this experiment, we use H1(x) = softplus(x) = log(1 + ex), and H2(x) = tanh(x), we also find similar performance with other choices of nonlinear functions. Bj is a K ×M random matrix, whose element is sampled from U [−1, 1]. Aji is a K ×M matrix, with 0.5 probability of being a zero matrix and 0.5 probability of being a nonzero random matrix, characterizing the underlying causal strength from j to i. Crucially, to reflect that the causal strength may span different orders of magnitude, if Aji is sampled to be a nonzero matrix, then the amplitude of each of its element is sampled from a log-normal distribution with µ = 1, σ = 0, their sign sampling from U{−1, 1}. Denote Aindi as the 0-1 indicator matrix of causality (Aindi,ji = 1 if |Aji| > 0; 0 otherwise). The goal of each algorithm being evaluated is to produce an N ×N causal matrix Ã, where each entry Ãji characterizes the causal strength from j to i. Then the flattened Ã is evaluated against the flattened Aindi (excluding diagonal elements of the matrices), producing the metrics of AUC-PR and AUC-ROC.
In general, for a large N , the number of possible causal graphs grows double exponentially: there are 2N 2
possible matrix of Aindi. To give an estimate, for N = 3, 4, 5, 8, 10, 20, 30, there are 512, 6.6 × 104, 3.3 × 107, 1.8 × 1019, 1.2 × 1030, 2.6 × 10120, 8.5 × 10270 number of possible graphs, respectively. Therefore, estimating the causal graph is in general a non-trivial task when N is large. We compare our algorithm, with previous methods including transfer entropy (Schreiber (2000)), causal influence (Janzing et al. (2013)), linear Granger causality (Granger (1969); Ding et al. (2006)), and a baseline of mutual information Ãji = I(X (j) t−1;x (i) t ) (which gives Ãji = Ãij). The implementation details for each method and each experiment are provided in Appendix D and Appendix E, respectively. Table 1 and 2 shows the average AUC-PR and AUC-ROC with each N , each with four random initializations of the true underlying causal matrix A and dataset.
We see that not only does our method outperform other methods by a large margin across all Ns, it also shows good performance when N is as large as 30, demonstrating our method’s capability and scalability to infer complex causal structures from interacting time series. For the Causal Influence method, although it has very good mathematical properties, it may be impractical in practice, as is also shown in the table. This is due to that it is defined as the KL-divergence between (Xt−1, x (i) t−1)
and (X(j̃)t−1, x (i) t−1), each of which is an (NK + 1)M−dimensional vector, which can quickly go to high dimensions, where density estimation required to calculate KL-divergence is in general data-hungry and difficult. In comparison, our method that estimates causal strength via minimizing prediction errors is comparatively easier in high dimensions (only have to predict a M dimensional vector conditioned on the inputs), which contributes to a better performance when N is large.
Since in practice, we do not know the underlying causal structure a priori, it presents a greater challenge to select the model capacity for fθ, as compared with supervised learning method where we can do cross-validation. To see how the capacity of the function approximator fθ influences our method, we vary the number of layers and the number of neurons in each layer at N = 10. Table 3 summarizes the result. We see that our method’s performance here is hardly influenced by the model capacity, with only a slight degradation at very low capacity. This shows that our method is quite tolerant and stable with model capacity variations.

heading: 4.2 EXPERIMENTS WITH VIDEO GAMES
text: To see how our method can infer causal relations in real videos games, and potentially helping reinforcement learning (RL) or imitation learning (IL), we apply our method to the causal inference between the trajectories of different objects from a trained CNN RL-agent playing Atari Breakout
.
games2 (Bellemare et al. (2013), implementation details see Appendix F). Fig 1 shows the inferred Wji matrix, with the (j, i)th element denoting the inferred causal strength from j to i. We see that there is a prominent causal direction from the ball’s y position to the reward, which correctly summarizes that the ball’s y position has a large influence on the reward. Additionally, the other discovered causal directions with causal strength greater than 1 include: brick→ reward, ball-y→ brick, reward → paddle, ball-x → action, ball-y → action. The former two correctly summarize causal chains from ball-y → brick → reward, the latter two show that the ball’s x and y positions also have big causal influences on the trained agent’s action: in order that the ball does not fall to the bottom, the agent has to position itself at the right position depending on the x and y positions of the ball.
In comparison, mutual information (Fig. 2 (a)) gives a symmetric matrix that does not differentiate the two possible causal directions. Moreover, it discovers a high mutual information between brick and paddle, which does not related to the underlying causal mechanism, and also missed the causal arrows ball-y→brick→reward. For transfer entropy, although it discovers a number of prominent causal arrows (e.g. ball-x→action, ball-y→action, ball-y→reward, brick→reward), it also gives relative high scores for some incorrect causal arrows: brick→ action, ball-y→ball-x, and missed an important causal arrow ball-y→brick. For linear Granger, it gives most arrows 0 or negative score, failing to discover any useful causal arrows. For causal influence, it also fails to discover useful causal arrows.
The inferred causal matrix and learned models may be useful for downstream learning in RL/IL. The learned models fθ can serve as a succinct subset of the environment model p(s′, r|s, a) that predicts future state s′ and reward r based on current state s and action a. Moreover, the agent can utilize the causal matrix to learn policy much more quickly, by reasoning with a causal chain from the reward backward, and from the action forward and backward, and focusing learning policies that can maximally influence the entities in the causal chains. Take the current Breakout game for example. After discovering the causal matrix (Fig. 1) by watching a teacher agent playing, the agent can reason from the reward backward on the causal chain ball-y→brick→reward, and understand that in order to maximize its reward, it has to influence the ball’s y position and the number of
2A video showing the game playing can be seen at https://goo.gl/XGzppc
bricks. Reasoning from the action backward on causal chains ball-x→action and ball-y→action, it understands that the teacher agent can obtain high rewards by making actions based on the ball’s x and y position (there are no causal chains from the action forward in Fig. 1 so the agent as to learn on its own). The agent can then utilize the above reasoning to develop policies π(a|s) where the state s can focus on the ball’s x and y position instead of the entire environment state, and the action a can focus on those that can maximally influence not only the reward, but also the ball-y and brick that is on the reward’s backward causal chain. This significantly reduces the search space for the policies, and the agent may probably require significantly less number of episodes to learn a good policy. It will be an exciting future research direction to incorporate this causal reasoning into RL and IL to improve sample efficiency.

heading: 4.3 HEART-RATE VS. BREATH-RATE AND RAT BRAIN EEG DATASETS
text: Now we test our algorithm with real-world datasets. As a common dataset studied in previous causal works, we use the time-series of the breathing rate and instantaneous heart rate of a sleeping patient suffering from sleep apnea (samples 2350-3550 of data set B from Santa Fe Institute time series contest held in 1991, available in PhysioNet). We apply our method to infer the causal strength between the breathing rate and heart rate, with different maximum time horizon K. The result is shown in Fig. 3. The causal strength from heart to breath is significantly higher than the reverse direction, consistent with the results from previous causal inference methods (Schreiber (2000); Ancona et al. (2004); Marinazzo et al. (2008a)) as also shown in Fig. 3. Notably, the causal strength remains at roughly the same level for different Ks, in contrast to the decaying causality index w.r.t. increasing history length in (Ancona et al. (2004), Fig. 3 (c)) showing a merit of our method in estimating causal strength across different time-horizons, aided by the flexibility of fθ in extracting the right information to predict the future. The implementation details in this section is provided in Appendix G.
As a second real-world example, we apply our algorithm in estimating the causal strength of the EEG signals between the right and left cortical intracranial electrodes (Quiroga), also studied in Ancona et al. (2004); Quiroga et al. (2002); Marinazzo et al. (2008a). Figure 4 (left) shows the inferred causal strength Wji for the EEG signals of a normal rat. We see that there is only a slight asymmetry, with the right channel having a slightly stronger influence on the left channel than the reverse direction. Figure 4 (right) shows Wji for the EEG signals with unilateral lesion in the rostral pole of the reticular thalamic nucleus. We see that there is stronger causal influence from the left to the right channels. Compared with the result of previous works Ancona et al. (2004); Marinazzo et al. (2008a) as also shown in Fig. 7, we see that all methods correctly infer the causal relations before and after brain lesion. In addition, our method shows non-decaying causal strength with increasing history length, in contrast to the decaying causality index in Ancona et al. (2004), again demonstrating our method’s insensitivity against history length, due to its flexibility in extracting the right amount of information in order to predict the future. The above two applications demonstrate our method’s capability in inferring the causal relations from noisy, real-world data.

heading: 5 DISCUSSION AND CONCLUSION
text: In this paper, we have addressed causal inference by proposing a novel causal measure, defined via a novel learnable noise risk, that allows function approximators like neural nets to learn complex causal relations via the Causal Inference with Learnable Noise (CILN) algorithm. We provide intuitions that the algorithm is likely to discover variables that directly structurally cause the variable of interest. We demonstrated in synthetic nonlinear datasets that our method outperforms previous methods by a large margin in inferring causal relations with complex structures, can scale to a large number of time series, and is hardly influenced by the capacity of the function approximator. Our method also correctly infers the causal arrows by watching a trained CNN playing Breakout, and give causal directions consistent with prior works in real-world heart-rate vs. breath-rate and rat EEG datasets.
Our CILN algorithm provides many opportunities for downstream tasks. As discussed in Section 4.2, we can incorporate CILN into reinforcement learning or imitation learning, and by reasoning backward and forward with the discovered causal chains, the agent may be able to learn useful policies with much less number of episodes. It may also help interpret the learned neural nets, by quantifying the causal strength from the input to different hidden neurons to the output neurons. Since Wji is based on mutual information which is scale-free, it may provide additional useful information on the internal mechanisms of neural net, in addition to the weights and gradients. CILN can also help decipher complex real-world systems, for example in neuroscience, economics and finance. As for the error effect of the algorithm, it depends on how the inferred causal matrix is utilized qualitatively or quantitatively downstream. For application in RL/IL, a falsely discovered causal arrow will increase the policy search space and thus reduce sample efficiency, but may not influence final performance since the RL/IL algorithm will eventually learn a (sub)optimal policy within the larger search space. Missed causal arrows (false negatives) may lead to the agent’s negligence of certain good policies, therefore it is important for the RL/IL algorithm to always allow a certain amount of exploration in addition to focusing policies on the causal chain. For helping interpretability of neural nets and deciphering complex real-world systems, the false positives/negatives of causal relations will influence the qualitative understanding and potential decision making downstream, and it is important to set a higher/lower threshold for the causal score depending on whether false positives or negatives are deemed less desirable.
Above all, we believe our work can not only pave the way for future exciting advancements in enabling machine learning models to understand causality, a key component of human intelligence, but also endow researchers with a useful tool for deciphering the causal relations in complex systems.

heading: Appendix

heading: A THEOREM 1 AND PROOF
text: Here, we formally propose a theorem corresponding to the statement in section 3.2, and provide proof afterwords.
Theorem 1. For system Eq. 1, for any i, j ∈ {1, 2, ...N}, i 6= j, if X(j)t−1 Granger-causes x (i) t , then X (j) t−1 directly structurally causes x (i) t .
Proof. We base the proof on the Theorem 5.6 in White et al. (2011). Firstly, by definition, the system Eq. (1) belongs to the canonical settable system (Def. 3.3 in White et al. (2011)), on which their Theorem 5.6 is based. To prove that in our system Granger causality can deduce direct structural causality, we only have to prove that the assumption A.1 and assumption A.2 in White et al. (2011) are satisfied by our system. If we identify our x(i)t with their Y1,t, our Xt−1 with their Yt−1, our x(j)t with their Y2,t, our ui,t (our ui at time t) with their U1,t, our uj,t with their U2,t, their Zt = ∅, Wt = ∅, then our system Eq. (1) satisfies their Assumption A.1. Additionally, by definition, our ui ∈ RM , i = 1, 2, ...N are random variables that are mutually independent, and also independent of any X(i)t−1, x (i) t , i ∈ {1, 2, ...N}. Therefore, our system satisfies their strict exogeneity (Yt−1,Zt) ⊥ U1,t (in our representation (Xt−1,∅) ⊥ ui,t), which is a sufficient condition for Assumption A.2. Therefore, both their Assumption A.1 and Assumption A.2 are satisfied by our system Eq. (1). Applying their Theorem 5.6, we prove Theorem 1.

heading: B ANALYZING THE PROPERTY OF THE LEARNABLE NOISE RISK
text: Firstly we state the assumption that will be used throughout this section: Assumption 1. Assume that fθ ∈ F is a continuous function and has enough capacity so that it can approximate any ∫ dx
(i) t P (x (i) t |Xt−1)x (i) t . Let j 6= i and assume that P (X (j) t−1) has support with
intrinsic dimension of KM .
Also we emphasize that in this paper, the expected risks (with symbol R) are w.r.t. the distributions, and the empirical risks (with symbol R̂) are w.r.t. a dataset drawn from the distribution, with finite number of examples. The theorems in this paper are all proved w.r.t. distributions (assuming infinite number of examples). Sample complexity results will be left for future work.
The structure of this section is as follows. First in subsection B.1 and B.2, we prove three lemmas that will be helpful for the following analysis. Then in subsection B.3, we analyze the property of the learnable noise risk both qualitatively and quantitatively, and argue why the learnable noise risk is likely to select the variables that directly structurally causes x(i)t .

heading: B.1 PROVING A LEMMA
text: Here we prove Lemmas 1.1. Lemma 1.1. Suppose that Assumption 1 holds, we have
argminfθRX,x(i) [fθ] = ∫ dx (i) t P (x (i) t |Xt−1)x (i) t (7)
and
minfθRX,x(i) [fθ] = EXt−1,x(i)t
[( x
(i) t −
∫ dx
(i) t P (x (i) t |Xt−1)x (i) t
)2] (8)
In other words, for the MSE risk, its minimum is attained when fθ(Xt−1) is the expectation of x (i) t conditioned on Xt−1.
Proof. The proof of the lemma is adapted from Papoulis (1985). The risk
RX,x(i) [fθ] = EXt−1,x(i)t
[( x
(i) t − fθ(Xt−1) )2] = ∫ dXt−1dx (i) t · P (Xt−1, x (i) t ) ( x (i) t − fθ(Xt−1)
)2 = ∫ dXt−1P (Xt−1) ∫ dx (i) t P (x (i) t |Xt−1) ( x (i) t − fθ(Xt−1)
)2 Note that here (x(i)t − fθ(Xt−1))2 ≡ 〈 x (i) t − fθ(Xt−1), x (i) t − fθ(Xt−1) 〉 is an inner product in RM .
For any Xt−1, treating fθ(Xt−1) ∈ RM as a vector, let’s calculate its value such that the integral F (fθ(Xt−1)) := ∫ dx (i) t P (x (i) t |Xt−1) ( x (i) t − fθ(Xt−1) )2 attains its minimum.
Let 0 =
∂
∂fθ(Xt−1) F (fθ(Xt−1))
= ∂
∂fθ(Xt−1)
∫ dx
(i) t P (x (i) t |Xt−1)
( x
(i) t − fθ(Xt−1) )2 = −2 ∫ dx
(i) t P (x (i) t |Xt−1)
( x
(i) t − fθ(Xt−1) ) we have ∫
dx (i) t P (x (i) t |Xt−1)x (i) t =
∫ dx
(i) t P (x (i) t |Xt−1)fθ(Xt−1)
= fθ(Xt−1)
∫ dx
(i) t P (x (i) t |Xt−1)
= fθ(Xt−1) Therefore, for any Xt−1, fθ(Xt−1) = ∫ dx (i) t P (x (i) t |Xt−1)x (i) t is the only stationary point for F (fθ(Xt−1)).
Taking the second derivative, we have ∂2
(∂fθ(Xt−1))2 F (fθ(Xt−1)) = 2
∫ dx
(i) t P (x (i) t |Xt−1)I = 2I
where I is an M ×M identity matrix, which is always positive definite. Therefore, for any Xt−1, fθ(Xt−1) = ∫ dx (i) t P (x (i) t |Xt−1)x (i) t is the only global minimum of

heading: F (fθ(Xt−1)) w.r.t. fθ(Xt−1).
text: Since RX,x(i) [fθ] = ∫ dXt−1P (Xt−1)F (fθ(Xt−1))
The minimum of the risk RX,x(i) [fθ] is attained iff F (fθ(Xt−1)) attains minimum at every Xt−1, i.e.,
fθ(Xt−1) =
∫ dx
(i) t P (x (i) t |Xt−1)x (i) t
is true for any Xt−1. Given Assumption 1, we know that fθ ∈ F has enough capacity such that it can approximate any ∫ dx
(i) t P (x (i) t |Xt−1)x (i) t . Therefore, argminfθRX,x(i) [fθ] = ∫ dx (i) t P (x (i) t |Xt−1)x (i) t
and
minfθRX,x(i) [fθ] = EXt−1,x(i)t
[( x
(i) t −
∫ dx
(i) t P (x (i) t |Xt−1)x (i) t
)2]

heading: B.2 MINIMUM MSE WITH DIFFERENT VARIABLES
text: Lemma 1.2. Suppose that Assumption 1 holds, and X(U)t−1,X (V ) t−1, X (W ) t−1 ⊂ Xt−1 are mutually exclusive sets of variables satisfying
X (W ) t−1 ⊥ x (i) t |X (U) t−1, X (V ) t−1, X (V ) t−1 6⊥⊥ x (i) t |X (U) t−1, X (W ) t−1
Then
minfθEX(U)t−1,X(V )t−1,x(i)t
[( x
(i) t − fθ(X (U) t−1, X (V ) t−1) )2] < minfθEX(U)t−1,X(V )t−1,x(i)t [( x (i) t − fθ(X (U) t−1, X (W ) t−1 ) )2] Proof. Since Assumption 1 holds, according to Lemma 1.1, Lemma 1.2 is equivalent to
E X
(U) t−1,X (V ) t−1,x (i) t
[( x
(i) t −
∫ dx
(i) t P (x (i) t |X (U) t−1, X (V ) t−1)x (i) t
)2]
< E X
(U) t−1,X (W ) t−1 ,x (i) t
[( x
(i) t −
∫ dx
(i) t P (x (i) t |X (U) t−1, X (W ) t−1 )x (i) t )2] We have
E X
(U) t−1,X (W ) t−1 ,x (i) t
[( x
(i) t −
∫ dx
(i) t P (x (i) t |X (U) t−1, X (W ) t−1 )x (i) t
)2]
= ∫ dX
(U) t−1dX (W ) t−1 dx (i) t P (X (U) t−1, X (W ) t−1 , x (i) t )
( x
(i) t −
∫ dx
(i) t P (x (i) t |X (U) t−1, X (W ) t−1 )x (i) t )2 = ∫ dX
(U) t−1dX (V ) t−1dX (W ) t−1 dx (i) t P (X (U) t−1, X (V ) t−1, X (W ) t−1 , x (i) t )
( x
(i) t −
∫ dx
(i) t P (x (i) t |X (U) t−1, X (W ) t−1 )x (i) t )2 = ∫ dX
(U) t−1dX (V ) t−1dX (W ) t−1 P (X (U) t−1, X (V ) t−1)P (X (W ) t−1 |X (U) t−1, X (V ) t−1)·∫
dx (i) t P (x (i) t |X (U) t−1, X (V ) t−1)
( x
(i) t −
∫ dx
(i) t P (x (i) t |X (U) t−1, X (W ) t−1 )x (i) t )2 > ∫ dX
(U) t−1dX (V ) t−1dX (W ) t−1 P (X (U) t−1, X (V ) t−1)P (X (W ) t−1 |X (U) t−1, X (V ) t−1)·∫
dx (i) t P (x (i) t |X (U) t−1, X (V ) t−1)
( x
(i) t −
∫ dx
(i) t P (x (i) t |X (U) t−1, X (V ) t−1)x (i) t )2 = ∫ dX
(U) t−1dX (V ) t−1P (X (U) t−1, X (V ) t−1)
∫ dx
(i) t P (x (i) t |X (U) t−1, X (V ) t−1)
( x
(i) t −
∫ dx
(i) t P (x (i) t |X (U) t−1, X (V ) t−1)x (i) t )2 = E
X (U) t−1,X (V ) t−1,x (i) t
[( x
(i) t −
∫ dx
(i) t P (x (i) t |X (U) t−1, X (V ) t−1)x (i) t
)2]
The third equality (the one before the inequality) is due to that X(W )t−1 ⊥ x (i) t |X (U) t−1, X (V ) t−1, leading to P (X(U)t−1, X (V ) t−1, X (W ) t−1 , x (i) t ) = P (X (U) t−1, X (V ) t−1)P (X (W ) t−1 |X (U) t−1, X (V ) t−1)P (x (i) t |X (U) t−1, X (V ) t−1). The inequality step first uses the setting in Eq. (1) that the noise variables ui are effective arguments of the response functions hi, and that each hi is “causality in mean". Therefore, ∫ dx
(i) t P (x (i) t |X (U) t−1, X (V ) t−1)x (i) t 6=
∫ dx
(i) t P (x (i) t |X (U) t−1, X (W ) t−1 )x (i) t . Us-
ing Lemma 1.1, we have fθ(X (U) t−1, X (V ) t−1) =
∫ dx
(i) t P (x (i) t |X (U) t−1, X (V ) t−1)x (i) t minimizes∫
dx (i) t P (x (i) t |X (U) t−1, X (V ) t−1)
( x
(i) t − fθ(X (U) t−1, X (V ) t−1)
)2 , hence the inequality.
Lemma 1.3. Suppose that Assumption 1 holds, and X(D)t−1 ⊆ Xt−1 are the set of variables that directly structurally causes x(i)t . Then ∀X (S) t−1 ⊆ Xt−1 with X (S) t−1 6= X (D) t−1 , we have
minfθEX(D)t−1
[( x
(i) t − fθ(X (D) t−1) )2] < minfθEX(S)t−1 [( x (i) t − fθ(X (S) t−1) )2] Specifically, we have
minfθEX(D)t−1
[( x
(i) t − fθ(X (D) t−1) )2] < minfθEX(D̂)t−1 [( x (i) t − fθ(X (D̂) t−1) )2]
where X(D̂)t−1 = Xt−1\X (D) t−1 .
Proof. For anyX(S)t−1, letX (U) t−1 = X (D) t−1∩X (S) t−1,X (V ) t−1 = X (D) t−1\X (S) t−1,X (W ) t−1 = X (S) t−1\X (D) t−1 . Then X (U) t−1,X (V ) t−1,X (W ) t−1 are mutually exclusive, andX (D) t−1 = X (U) t−1∪X (V ) t−1,X (S) t−1 = X (U) t−1∪X (W ) t−1 . Now we prove that ∀X(S)t−1 ⊆ Xt−1 with X (S) t−1 6= X (D) t−1 , the corresponding X (U) t−1, X (V ) t−1, X (W ) t−1 , x (i) t satisfy the condition for Lemma 1.2. SinceX(D)t−1 are the set of variables that directly structurally causes x
(i) t , there does not exist a X (S) t−1 such that the corresponding X (V ) t−1 ⊥ x (i) t |X (U) t−1, X (W ) t−1 (otherwise
it violates the definition of direct structural causality). Thus X(V )t−1 6⊥⊥ x (i) t |X (U) t−1, X (W ) t−1 . To prove X (W ) t−1 ⊥ x (i) t |X (U) t−1, X (V ) t−1, note that X (W ) t−1 does not directly structurally cause x (i) t , then by Theorem 1, X(W )t−1 does not Granger-cause x (i) t , i.e. P (x (i) t |X (U) t−1, X (V ) t−1) = P (x (i) t |X (U) t−1, X (V ) t−1, X (W ) t−1 ), which is equivalent to X(W )t−1 ⊥ x (i) t |X (U) t−1, X (V ) t−1. The special case of X (D̂) t−1 follows directly that X (D̂) t−1 = Xt−1\X (D) t−1 6= X (D) t−1 and letting X (S) t−1 = X (D̂) t−1 .

heading: B.3 QUALITATIVE AND QUANTITATIVE BEHAVIORS OF THE LEARNABLE NOISE RISK
text: In this section, we analyze the qualitative and quantitative behaviors of the learnable noise risk (Eq. 2), with varying noise levels ηj . For each variable X (j) t−1 ∈ Xt−1, j = 1, 2, ...N , define
ρj = tanh ( I(X (j) t−1; X̃ (j)(ηj) t−1 ) ) ∈ [0, 1] as a “rescaled" mutual information between X(j)t−1 and
X̃ (j)(ηj) t−1 . When ηj = 0 so that X̃ (j)(ηj) t−1 = X (j) t−1, ρj = 1. When all elements of ηj → ∞, ρj = 0. Denoting ρ = (ρ1, ρ2, ...ρN ), we can then rewrite the learnable noise risk (Eq. 2) as
RX,x(i) [fθ,ρ] = MMSE (i)(ρ) + λ · N∑ j=1 arctanh(ρj) (9)
where MMSE(i)(ρ) = minη,fθ EXt−1,x(i)t ,
[( x
(i) t − fθ(X̃ (η) t−1) )2] subject to ρj =
tanh ( I(X
(j) t−1; X̃ (j)(ηj) t−1 ) , j = 1, 2, ...N . Let X(D)t−1 ⊆ Xt−1 be the set of variables that
directly structurally causes x(i)t , and denote the corresponding set of ρj as ρ (D). Denote X (D̂) t−1 = Xt−1\X (D) t−1 and the corresponding set of ρj as ρ
(D̂). For any i = 1, 2, ...N , we have the following properties:
1. MMSE(i)(ρ) attains maximum at ρ = 0.
2. MMSE(i)(ρ) is monotonically decreasing w.r.t. each ρj . 3. MMSE(i)(ρ) ∣∣ ρ(D)=1,ρ(D̂)=0 < MMSE(i)(ρ) ∣∣ ρ(D)=0,ρ(D̂)=1 (using Lemma 1.3).
4. MMSE(i)(ρ) attains minimum at ρ(D) = 1. MMSE(i)(ρ) ∣∣ ρ(D)=1 is constant w.r.t. ρ(D̂).
To get a better intuition of the landscape of RX,x(i) [fθ,ρ], let’s investigate a simple example. Let the response function be:  x (1) t := h1(u1) = √ Σx · u1 x (2) t := h2(x (1) t−1, u2) = x (1) t−1 + √ Ωx · u2
x (3) t := h3(x (2) t−1, u3) = x (2) t−1 + √ Ωy · u3
(10)
where u1, u2, u3 are independent unit Gaussian variables, and Xt−1 = (X (1) t−1, X (2) t−1, X (3) t−1) =(
(x (1) t−2, x (1) t−1), (x (2) t−2, x (2) t−1), (x (3) t−2, x (3) t−1) ) . For RX,x(3) [fθ,ρ] = MMSE
(3)(ρ) + λ ·∑3 j=1 arctanh(ρj), since only x (1) t−2 and x (2) t−1 are d-connected to x (3) t , at the minimization of RX,x(3) [fθ,ρ], only x (1) t−2 and x (2) t−1 may have a finite η ∗ j,l (the other η ∗ j,l are all infinite). Therefore, setting the ηj,l not corresponding to x (1) t−2 and x (2) t−1 as infinity, and let x̃ (1) t−2 = x (1) t−2 + ηx · x, x̃ (2) t−1 = x (2) t−1 + ηy · y , x and x being independent unit Gaussian variables. Let fθ(x (1) t−2, x (2) t−1) = a · x(1)t−2 + b · x (2) t−1, then we can get an analytic expression for RX,x(3) [fθ, ηx, ηy]:
RX,x(3) [fθ, ηx, ηy]
= a2Σx + (b− 1)2(Σx + Ωx) + a2η2x + b2η2y + 2a(b− 1)Σx + Ωy + λ 2 log ( 1 + Σx η2x ) + λ 2 log ( 1 + Σx + Ωx η2y ) Minimizing RX,x(3) [fθ, ηx, ηy] w.r.t. a and b, we get
a∗ = η2yΣx
η2xη 2 y + η 2 xΣx + η 2 yΣx + η 2 xΩx + ΩxΣx
b∗ = η2x(Σx + Ωx) + ΣxΩx
η2xη 2 y + η 2 xΣx + η 2 yΣx + η 2 xΩx + ΩxΣx
Substituting into RX,x(3) [fθ, ηx, ηy], we have
RX,x(3) [ηx, ηy]
= min fθ RX,x(3) [fθ, ηx, ηy]
= η2y(ΣxΩx + η 2 x(Σx + Ωx))
η2xη 2 y + η 2 xΣx + η 2 yΣx + η 2 xΩx + ΩxΣx
+ λ 2 log ( 1 + Σx η2x ) + λ 2 log ( 1 + Σx + Ωx η2y )
Here we have neglected the constant Ωy . To obtain RX,x(3) [ρ], let ρ1 = tanh ( 1 2 log ( 1 + Σxη2x )) ,
ρ2 = tanh ( 1 2 log ( 1 + Σx+Ωxη2x )) , we have η2x = 1−ρ1 2ρ1 Σx, η2y = 1−ρ2 2ρ2
(Σx + Ωx). Substituting, we have
RX,x(3) [ρ] = MMSE (3)(ρ) + λ · 2∑ j=1 arctanh(ρj)
= (ρ2 − 1)(Σx + Ωx)((ρ1 − 1)Σx − (ρ1 + 1)Ωx)
(1 + ρ1 + ρ2 − 3ρ1ρ2)Σx + (1 + ρ1)(1 + ρ2)Ωx + λ · arctanh(ρ1) + λ · arctanh(ρ2)
Fig. 5 shows the landscape of MMSE(3)(ρ) andRX,x(3) [ρ], for Σx = 1,Ωx = 2, λ = 1. We see that MMSE(3)(ρ) satisfies the above mentioned four properties. Particularly, MMSE(3)(ρ) ∣∣ ρ1=1,ρ2=0 >
MMSE(3)(ρ) ∣∣ ρ1=0,ρ2=1
. After adding λ · arctanh(ρ1) + λ · arctanh(ρ2), the RX,x(3) [ρ] has global minimum along ρ1 = 0 largely due to this property. Therefore, for this particular example, when RX,x(3) [ρ] is minimized, ρ1 = 0, i.e. I(x (1) t−2, x̃ (1)(η∗1 ) t−2 ) = 0.
By varying the value of λ, we can tune the relative influence of the two terms MMSE(3)(ρ) and∑2 j=1 arctanh(ρj). The landscape corresponding to λ = 0.01, 0.5, 2, 10 are plotted in Fig. 6. We see that when λ 1, the MMSE term dominates, and it is possible that the global minimum of RX,x(3) [ρ] is not at ρ1 = 0. This is similar to the effect of a L1 regularization, where if the coefficient λ for the L1 is vanishingly small, the L1 regularization will barely influence the loss landscape. When λ is not vanishingly small, as in Fig. 6 (b), we see that the global minimum of RX,x(3) [ρ] lies on ρ1 = 0. When λ → +∞, the ∑2 j=1 arctanh(ρj) term dominates and the global minimum is at ρ1 = 0, ρ2 = 0.
In general, we expect RX,x(i) [ρ] behave qualitatively similar. When λ→ +∞, the global minimum for RX,x(i) [ρ] is at ρ∗ = 0. As we ramp down λ, the dimension that has largest influence on MMSE will first host the global minimum with nonzero ρ∗j , which is most likely the variable that directly structurally causes x(i)i . When λ is further ramping down, we expect that the variables that host the global minimum with nonzero ρj will more likely be those that directly structurally causes x (i) i , due to the landscape influenced by the four properties of MMSE. This can justify the learnable noise risk as a good objective for causal discovery/variable selection. The experiments in the paper will empirically test the performance of the learnable noise risk.

heading: C UPPER BOUND FOR THE LEARNABLE NOISE RISK
text: In this section, we prove that I(X̃(j)(ηj)t−1 ;X (j) t−1) ≤ 12 ∑KM l=1 log ( 1 + Var(X(j)t−1,l) η2j,l ) . We formally state the theorem as follows:
Theorem 2. Let X̃(j)(ηj)t−1 := X (j) t−1 + ηj · j , j = 1, 2, ...N be the noise-corrupted inputs with learnable noise amplitudes ηj ∈ RKM , and j ∼ N(0, I). We have
I(X̃ (j)(ηj) t−1 ;X (j) t−1) ≤
1
2 KM∑ l=1 log
( 1 + Var(X(j)t−1,l)
η2j,l
) (11)
where l is the lth element of a vector, std(X(j)t−1,l) is the standard deviation of X (j) t−1,l across t. The equality is reached whenX(j)t−1 obeys a multivariate Gaussian distribution with diagonal covariance matrix Σ satisfying Σl,l = Var(X (j) t−1,l) + η 2 j,l.
Proof. We have I(X̃
(j)(ηj) t−1 ;X (j) t−1) = H(X̃ (j)(ηj) t−1 )−H(ηj · j)
= H(X̃ (j)(ηj) t−1 )−
( KM
2 log(2πe) + KM∑ l=1 1 2 log(η2j,l)
)
Here H(·) is differential entropy. For X̃(j)(ηj)t−1 , its variance at the lth dimension is
Var(X̃(j)(ηj)t−1,l ) = Var(X (j) t−1,l + ηj · j)
= Var(X(j)t−1,l) + Var(ηj,l · j,l)
= Var(X(j)t−1,l) + η 2 j,l
The second equality is due to that X(j)t−1 is independent of j . Using the principle of maximum entropy, the distribution that maximizes H(X̃(j)(ηj)t−1 ) subject to the constraint of Var(X̃ (j)(ηj) t−1,l ) = Var(X(j)t−1,l) + η 2 j,l, l = 1, 2, ...KM is a Gaussian distribution whose diagonal covariance matrix Σ satisfies Σl,l = Var(X (j) t−1,l)+η 2 j,l. Its entropy is H(X̃ (j)(ηj) t−1 ) = KM 2 log(2πe)+ ∑KM l=1 1 2 log(η 2 j,l+ Var(X(j)t−1,l)). Therefore,
I(X̃ (j)(ηj) t−1 ;X (j) t−1)
≤
( KM
2 log(2πe) + KM∑ l=1 1 2 log(η2j,l + Var(X (j) t−1,l))
) − ( KM
2 log(2πe) + KM∑ l=1 1 2 log(η2j,l)
)
= 1
2 KM∑ l=1 log
( 1 + Var(X(j)t−1,l)
η2j,l
)
The equality is reached when X(j)t−1 obeys a multivariate Gaussian distribution with diagonal covariance matrix Σ satisfying Σl,l = Var(X (j) t−1,l) + η 2 j,l.
D IMPLEMENTATION DETAILS FOR THE METHODS
Here we state the implementation details for our method, as well as other methods being compared. Throughout this paper, unless otherwise specified, we use the standard technique in Kraskov et al. (2004) to estimate the KL-divergence and mutual information, which is used in our implementations of Mutual information, Transfer Entropy and Causal Influence.

heading: D.1 OUR METHOD
text: Without stating otherwise, our method (Algorithm 1) as a default uses a three layer neural net, with two hidden layers having 8 neurons and SELU (Klambauer et al. (2017)) activation, and the last layer having linear activation. Adam (Kingma & Ba (2014)) optimizer with learning rate = 10−4 is used as default throughout this paper. We set η0 = 0.01 and λ = 0.01. We use 10000 epochs, with earlystopping such that if the best validation loss does not go down for 40 monitoring points (we monitor the validation loss every 20 epochs), do early-stop. It also has a 400 epoch warm-up period where the mutual information term is turned off, to allow fθ to find a good initial model as a start. We use the approximation I(X̃ (j)(η∗j )
t−1 ;X (j) t−1) ' 12 ∑KM l=1 log(1 + 1/χ 2 j,l) in the risk and also in estimating
Wji, as discussed in the main text in 3.3. Here χj,l = ηj,l
std(X(j)t−1,l) is the relative noise scale w.r.t. the
standard deviation of each element, l denoting the lth element of the KM -dimensional vector. In this work, we fix χj,l to be the same for each j, and let χj be a single parameter instead of a vector. This simplifies the risk calculation, and also to a first order invariant to the reparameterizaiton of each time series X(j)t−1.

heading: D.2 TRANSFER ENTROPY
text: We use the definition of transfer entropy as defined in Schreiber (2000). In that work the transfer entropy is defined for two time series. To deal with multiple time series, we let X(ĵ)t−1 also include other time series, similar to the extension of transfer entropy as in Lizier et al. (2008).

heading: D.3 CAUSAL INFLUENCE
text: For causal influence, we use the same network architecture as in our method, to learn a prediction model. Then the KL divergence is estimated via the technique in Kraskov et al. (2004).

heading: D.4 LINEAR GRANGER
text: We use the definition of linear Granger causality (Eq. (7) and (8) in Ding et al. (2006)) to calculate linear Granger causality. Specifically, we estimate the variance of the residue of a linear predictor of x(i)t−1 with and without X (j) t−1 (also conditioned on X (ĵ) t−1), using Levinson-Whittle(-Wiggins) and Robinson algorithm (Morf et al. (1978)). Then the linear Granger causality equals the log of the ratio of the two variances.
E IMPLEMENTATION DETAILS FOR SYNTHETIC EXPERIMENTS
For all experiments in this section, each metrics is obtained by performing the experiments (including generation of the dataset and the training) four times with seed = 0, 30, 60, 90 and averaging the resulting metrics. For the ground-truth causal tensor A, each element Aji is a K ×M matrix, with 0.5 probability of being an all-zero matrix, and 0.5 probability of being a nonzero matrix. If Aji is a nonzero matrix, its each element is sampled from a log-normal distribution with µ = 0 and σ = 1. For B, each Bj is also a K × N matrix, with each element sampling from U [−1, 1]. We use H1(x) = softplus(x) = log(1 + ex), and H2(x) = tanh(x) in equation (6). As a default, 500 time series each with length of 22 are generated from Eq. (6), each of which is wrapped into 19 (Xt−1, x (i) t ) pairs (since K = 3), so there are in total 500× 19 = 9500 examples for each dataset. The train-test-split is 4:1 for all experiments in this paper.

heading: F DETAILS FOR THE VIDEO GAME DATASET
text: Here, we implement a custom Atari Breakout game in the OpenAI Gym (Brockman et al. (2016)) environment, mimicking the original game, where we can access the state of the ball, paddle and bricks, etc. This representation is also used in the OO-MDP (Diuk et al. (2008)) paradigm for a more efficient representation of the environment state. We use the DQN algorithm, the same CNN architecture as in Mnih et al. (2015) to train an RL agent. Then we let it play the game for 20000 steps, obtaining a dataset with time-length of 20000 steps (if the agent dies, we restart the game) and 6 time series: action, paddle’s x position, ball’s x position, ball’s y position, number of bricks and reward. We then feed the time series (each time series normalized to mean of 0 and variance of 1) to our method, the same procedure as performed in the synthetic experiment, to let it produce an inferred causal matrix Wji, which is shown in Fig. 1. All the datasets used in this paper and code will be open-sourced upon publication of the paper.
G IMPLEMENTATION DETAILS FOR THE REAL-WORLD DATASET
For the two real-world datasets, we obtain the data with the same procedure as in Ancona et al. (2004). Then the data (each time series normalized to mean of 0 and variance of 1) are fed into our algorithm to infer the causal strength Wji with the default settings as described in Appendix D and E. For each K = 1, 2, ...20, the experiments are run for 50 times with seed from 0 to 49, and Fig. 3 and Fig. 4 are obtained by averaging over the inferred W matrix.
H CAUSALITY RESULTS FOR THE RAT-EEG DATASET BY PREVIOUS WORKS

Review summaries:
In the manuscript entitled "Neural Causal Discovery with Learnable Input Noise" the authors describe a method for automated causal inference under the scenario of a stream of temporally structured random variables (with no missingness and a look-back window of given size).  The proposed approach combines a novel measure of the importance of fidelty in each variable to predictive accuracy of the future system state ("learnable noise risk") with a flexible functional approximation (neural network).  Although the setting (informative temporal data) is relatively restricted with respect to the general problem of causal inference, this is not unreasonable given the proposed direction of application to automated reasoning in machine learning.  The simulation and real data experiments are interesting and seem well applied.