Title: NEURAL CAUSAL DISCOVERY WITH LEARNABLE IN-PUT NOISE

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.

1. INTRODUCTION:

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.

2. RELATED WORK:

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 Θ(N 2 ) 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.

3. METHOD:



3.1. PROBLEM DEFINITION:

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 ∈ R M 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 X t-1 = {X (i) t-1 }, i = 1, 2, ...N . We also denote X ( ĵ) t-1 = X t-1 \X (j) t-1 , i.e. X t-1 ex- cluding 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 h i that are unknown to the learner:

         x (1) t := h 1 (X t-1 , u 1 ) x (2) t := h 2 (X t-1 , u 2 ) ... x (N ) t := h N (X t-1 , u N ) (1) for t = K + 1, K + 2, ... . Here u i ∈ R M , 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 h i . 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

( ĵ) t-1 , X (j) t-1 , x (i) t ) have probability density function P (X ( ĵ) 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 = i does not directly structurally cause x (i) t , if for all possible values of X ( ĵ) t-1 and u l , l ∈ 1, 2, ...N , the function X 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 = i as default.

(j) t-1 → h i (X t-1 , u i ) is constant in X (j) t-1 . Otherwise, we say X (j) t-1 directly structurally causes x (i) t .

Granger causality (Granger & Newbold (1986) ) We say X (j) t-1 , j = i does not Granger-cause x

(i) t , if P (x (i) t |X (j) t-1 , X ( ĵ) t-1 ) = P (x (i) t |X ( ĵ) t-1 ), i.e. the conditional distribution function of x (i) t given X (j) t-1 , X ( ĵ) t-1 is identical to the conditional distribution function of x (i) t given X ( ĵ) 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 di- rectly structurally causes x (i)

t for each j = 1, 2, ...N .

3.2. GRANGER CAUSALITY IMPLIES DIRECT STRUCTURAL CAUSALITY IN EQ. (1):

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 = 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.

3.3. OUR METHOD:

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 Schreiber (2000) ), may not give a good estimate. Alternatively, various works have resorted to testing whether the prediction of

(i) t ; X (j) t-1 |X ( ĵ) t-1 ), see

x (i) t given X (j) t-1 , X ( ĵ) t-1 is better than given only X ( ĵ)

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 , X t-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 ) = (η 0 1, η 0 1, ...η 0 1), where each element η 0 1 is a KM - dimensional vector, same dimension as X (j) t-1 . 4: (f θ * , η * ) ← Minimize (f θ ,η) RX,x (i) , [f θ , η] (Eq. 4

) with e.g. gradient descent 5:

W ji ← 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:

EQUATION

where

X(η) t-1 := X t-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 ∈ R KM , 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

EQUATION

as another measure for causality, where

(f θ * , η * ) = argmin (f θ ,η) R X,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:

EQUATION

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 ) η 2 j,l

, where l denotes the l th 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:

EQUATION

Under review as a conference paper at ICLR 2019 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 W ji , which we term Causal Inference with Learnable Noise (CILN).

4. EXPERIMENTS:

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 Rupper 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).

4.1. SYNTHETIC EXPERIMENT WITH LOG-NORMAL CAUSAL STRENGTHS:

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 h i in Eq. ( 1):

EQUATION

where u t ∼ N (0, I) ∈ R M , denoting element-wise multiplication, and H 1 and H 2 are two nonlinear functions to make the response functions nonlinear. In this experiment, we use H 1 (x) = softplus(x) = log(1 + e x ), and H 2 (x) = tanh(x), we also find similar performance with other choices of nonlinear functions. B j is a K × M random matrix, whose element is sampled from U [-1, 1]. A ji 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 A ji 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 A indi as the 0-1 indicator matrix of causality (A indi,ji = 1 if |A ji | > 0; 0 otherwise). The goal of each algorithm being evaluated is to produce an N × N causal matrix Ã, where each entry Ãji characterizes the causal strength from j to i. Then the flattened Ã is evaluated against the flattened A indi (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 2 N 2 possible matrix of A indi . To give an estimate, for N = 3, 4, 5, 8, 10, 20, 30, there are 512, 6.6 × 10 4 , 3.3 × 10 7 , 1.8 × 10 19 , 1.2 × 10 30 , 2.6 × 10 120 , 8.5 × 10 270 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 Ãji = I(X

(j) t-1 ; x (i) t ) (which gives Ãji = Ã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 N s, 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 (X t-1 , x 

t-1 ), each of which is an (N K + 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. 

AUC-PR AUC-ROC:



4.2. EXPERIMENTS WITH VIDEO GAMES:

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 W ji 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 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.

4.3. HEART-RATE VS. BREATH-RATE AND RAT BRAIN EEG DATASETS:

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 W ji 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 W ji 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. 

5. DISCUSSION AND CONCLUSION:

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 W ji 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.

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 |X t-1 )x (i)

t . Let j = 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 .

B.1 PROVING A LEMMA:

Here we prove Lemmas 1.1. Lemma 1.1. Suppose that Assumption 1 holds, we have

argmin f θ R X,x (i) [f θ ] = dx (i) t P (x (i) t |X t-1 )x (i) t (7) and min f θ R X,x (i) [f θ ] = E Xt-1,x (i) t x (i) t -dx (i) t P (x (i) t |X t-1 )x (i) t 2 (8)

In other words, for the MSE risk, its minimum is attained when f θ (X t-1 ) is the expectation of x

(i) t conditioned on X t-1 .

Proof. The proof of the lemma is adapted from Papoulis (1985) . The risk

R X,x (i) [f θ ] = E Xt-1,x (i) t x (i) t -f θ (X t-1 ) 2 = dX t-1 dx (i) t • P (X t-1 , x (i) t ) x (i) t -f θ (X t-1 ) 2 = dX t-1 P (X t-1 ) dx (i) t P (x (i) t |X t-1 ) x (i) t -f θ (X t-1 ) 2 Note that here (x (i) t -f θ (X t-1 )) 2 ≡ x (i) t -f θ (X t-1 ), x (i) t -f θ (X t-1 ) is an inner product in R M .

For any X t-1 , treating f θ (X t-1 ) ∈ R M as a vector, let's calculate its value such that the integral

F (f θ (X t-1 )) := dx (i) t P (x (i) t |X t-1 ) x (i) t -f θ (X t-1 ) 2 attains its minimum. Let 0 = ∂ ∂f θ (X t-1 ) F (f θ (X t-1 )) = ∂ ∂f θ (X t-1 ) dx (i) t P (x (i) t |X t-1 ) x (i) t -f θ (X t-1 ) 2 = -2 dx (i) t P (x (i) t |X t-1 ) x (i) t -f θ (X t-1 )

we have dx

(i) t P (x (i) t |X t-1 )x (i) t = dx (i) t P (x (i) t |X t-1 )f θ (X t-1 ) = f θ (X t-1 ) dx (i) t P (x (i) t |X t-1 ) = f θ (X t-1 )

Therefore, for any X t-1 , f θ (X t-1 ) = dx

(i) t P (x (i) t |X t-1 )x (i)

t is the only stationary point for F (f θ (X t-1 )).

Taking the second derivative, we have

∂ 2 (∂f θ (X t-1 )) 2 F (f θ (X t-1 )) = 2 dx (i) t P (x (i) t |X t-1 )I = 2I

where I is an M × M identity matrix, which is always positive definite.

Therefore, for any X t-1 , f θ (X t-1 ) = dx

(i) t P (x (i) t |X t-1 )x (i) t is the only global minimum of F (f θ (X t-1 )) w.r.t. f θ (X t-1 ). Since R X,x (i) [f θ ] = dX t-1 P (X t-1 )F (f θ (X t-1

))

The minimum of the risk R X,

x (i) [f θ ] is attained iff F (f θ (X t-1 )) attains minimum at every X t-1 , i.e., f θ (X t-1 ) = dx (i) t P (x (i) t |X t-1 )x (i) t

is true for any X t-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 |X t-1 )x (i) t . Therefore, argmin f θ R X,x (i) [f θ ] = dx (i) t P (x (i) t |X t-1 )x (i) t and min f θ R X,x (i) [f θ ] = E Xt-1,x (i) t x (i) t -dx (i) t P (x (i) t |X t-1 )x (i) t 2 B.2 MINIMUM MSE WITH DIFFERENT VARIABLES

Lemma 1.2. Suppose that Assumption 1 holds, and

X (U ) t-1 ,X (V ) t-1 , X (W )

t-1 ⊂ X t-1 are mutually exclusive sets of variables satisfying

EQUATION

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-1 dX (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-1 dX (V ) t-1 dX (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-1 dX (V ) t-1 dX (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-1 dX (V ) t-1 dX (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-1 dX (V ) t-1 P (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 , lead- ing 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 u i are effective arguments of the response functions h i , and that each h i is "causality in mean". Therefore, dx

EQUATION

, hence the inequality.

Lemma 1.3. Suppose that Assumption 1 holds, and X (D)

t-1 ⊆ X t-1 are the set of variables that directly structurally causes x

(i) t . Then ∀X (S) t-1 ⊆ X t-1 with X (S) t-1 = X (D) t-1 , we have min f θ E X (D) t-1 x (i) t -f θ (X (D) t-1 ) 2 < min f θ E X (S) t-1 x (i) t -f θ (X (S) t-1 ) 2 Specifically, we have min f θ E X (D) t-1 x (i) t -f θ (X (D) t-1 ) 2 < min f θ E X ( D) t-1 x (i) t -f θ (X ( D) t-1 ) 2 where X ( D) t-1 = X t-1 \X (D) t-1 . Proof. For any X (S) t-1 , let X (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, and

X (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 ⊆ X t-1 with X (S) t-1 = X (D) t-1 , the corresponding X (U ) t-1 , X (V ) t-1 , X (W ) t-1 , x (i) t sat- isfy the condition for Lemma 1.2. Since X (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 ⊥ ⊥ 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 Theo- rem 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 = X t-1 \X (D) t-1 = X (D) t-1 and letting X (S) t-1 = X ( D) t-1 .

B.3 QUALITATIVE AND QUANTITATIVE BEHAVIORS OF THE LEARNABLE NOISE RISK:

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 ∈ X t-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

EQUATION

where

MMSE (i) (ρ) = min η,f θ E Xt-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 ⊆ X t-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 = X t-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) . 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 R X,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)

EQUATION

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 R X,x (3) [f θ , η x , η y ]:

EQUATION

[f θ , η x , η y ] w.r.t. a and b, we get

a * = η 2 y Σ x η 2 x η 2 y + η 2 x Σ x + η 2 y Σ x + η 2 x Ω x + Ω x Σ x b * = η 2 x (Σ x + Ω x ) + Σ x Ω x η 2 x η 2 y + η 2 x Σ x + η 2 y Σ x + η 2 x Ω x + Ω x Σ x

Substituting into R X,x (3) [f θ , η x , η y ], we have

R X,x (3) [η x , η y ] = min f θ R X,x (3) [f θ , η x , η y ] = η 2 y (Σ x Ω x + η 2 x (Σ x + Ω x )) η 2 x η 2 y + η 2 x Σ x + η 2 y Σ x + η 2 x Ω x + Ω x Σ x + λ 2 log 1 + Σ x η 2 x + λ 2 log 1 + Σ x + Ω x η 2 y

Here we have neglected the constant Ω y . To obtain R X,x (3) [ρ], let ρ 1 = tanh 1 2 log 1 + Σx η 2

x , ρ 2 = tanh 1 2 log 1 + Σx+Ωx η 2

x , we have η 2 x = 1-ρ1 2ρ1 Σ x , η 2 y = 1-ρ2 2ρ2 (Σ x + Ω x ). Substituting, we have 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 R X,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 R X,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 R X,x (i) [ρ] behave qualitatively similar. When λ → +∞, the global minimum for R X,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.

C UPPER BOUND FOR THE LEARNABLE NOISE RISK:

In this section, 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 ) η 2 j,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 ∈ R KM , and j ∼ N (0, I). We have

EQUATION

where l is the l th 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 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 . 

:

Note also that throughout this paper, when talking about causal matrices, the (j, i) th element always denotes the causal strength from j to i.

A video showing the game playing can be seen at https://goo.gl/XGzppc

Figures and Tables:

Figure1: Causal strength W ji inferred by our method with watching a trained CNN playing Breakout. The (j, i) element denotes the inferred causal strength from j to i.

Figure 3: (a) Causal strength W ji inferred by our method with the heart-rate vs. breath-rate dataset, averaged over 50 initializations of f θ . The shaded areas in this figure and Fig. 4 are the 95% confidence interval. (b) Upper: the filtered causality index vs. varying width of Gaussian kernel σ (Marinazzo et al. (2008a)); lower: transfer entropy vs. r, the length scale(Schreiber (2000)); (c) The causality index for breath→heart (lower) and heart→breath (upper) inAncona et al. (2004).

Figure 4: Causal strength inferred by our method with the EEG datasets, for different maximum time horizon K, averaged over 50 initializations of f θ . (Left) Causal strength W ji for the EEG signal for a normal rat. (Right) Causal strength W ji of EEG signal from the same rat, after brain lesion.

Fig.5shows the landscape of MMSE (3) (ρ) and R X,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 R X,x (3) [ρ] has global minimum along ρ 1 = 0 largely due to this property. Therefore, for this particular example, whenR X,x (3) [ρ] is minimized, ρ 1 = 0, i.e. I(x (1) t-2 ,x(1)(η * 1 ) t-2 ) = 0.

Figure 5: (a) MMSE (3) (ρ) and (b) R X,x (3) [ρ] in section B.3, for Σ x = 1, Ω x = 2, λ = 1.

Figure 6: (a) R X,x (3) [ρ] for (a) λ = 0.01, (b) λ = 0.5, (c) λ = 2 and (d) λ = 10 in section B.3, for Σ x = 1, Ω x = 2.

Average AUC-PR vs. N , with random sampling of A indi . Bold font marks the top method for each N .

Table content: <table><tr><td>(i) t-1 )</td></tr></table>

Average AUC-ROC vs. N , with random sampling of A indi . Bold font marks the top method for each N .

Table content: <table><tr><td>and (X</td><td>( j) t-1 , x</td></tr></table>

Average AUC-PR and AUC-ROC for different network structures for N = 10 with our method. Here for example,(8, 8, 8)  means that the f θ has 3 hidden layers, each with 8 neurons.

Table content: <table><tr><td>Neurons in hidden layers</td><td/><td/></tr><tr><td>(8)</td><td>0.8969</td><td>0.9086</td></tr><tr><td>(8, 8)</td><td>0.9309</td><td>0.9388</td></tr><tr><td>(16, 16)</td><td>0.9404</td><td>0.9455</td></tr><tr><td>(8, 8, 8)</td><td>0.9339</td><td>0.9410</td></tr><tr><td>(16, 16, 16)</td><td>0.9284</td><td>0.9288</td></tr><tr><td>(8, 8, 8, 8)</td><td>0.9312</td><td>0.9350</td></tr><tr><td>(16, 16, 16, 16)</td><td>0.9199</td><td>0.9202</td></tr></table>

Peer review (iclr2019):
 The paper proposes an approach to learn nonlinear causal relationship from time series data that is based on empirical risk minimization regularized by mutual information. The mutual information at the minimizer of the objective function is used as causal measure. The paper is well written and the proposed method well motivate and intuitive. However I am concerned by the assumption that the lagged variables X_{t-1}^{(j)} follow a diagonal gaussian distribution. This appears to be very restrictive, since typically the values of time series j at time t-1 are typically depending say of those that time t-2, t-3 etc. Another key concern concerns scalability. The authors mention gene regulatory networks , neuroscience etc as key applications. Yet the experiments considered in the paper are limited to very few time series. For instance the simulation experiments use N=30, which is much smaller than the number of time series usually involved say in gene regulatory network data. The real data experiments use N= 6 or N=2. This is way to small. The real data experiments (sections 4.2 and 4.3) are not very convincing, not only because of the very small size of N, but also because there is no comparison with the other approaches. How do these compare? Does the proposed approach offer insights on these datasets which are not captured by the comparison methods?