Using machine learning to assess short term causal dependence and infer network links

The core goal of science is often described to be a generalization from observations to understanding,¹ commonly embodied in predictive theories. Related to this is the desire to use measured data to infer necessary properties and the structure of any description consistent with a given class of observations. On the other hand, it has recently emerged that machine learning (ML) is capable of effectively performing a wide range of interpretive and predictive tasks on data.² Thus, it is natural to ask whether machine learning might be useful for the common scientific goal of discovering structural properties of a system from data generated by this system. In this paper, we consider an important, widely applicable class of such tasks. Specifically, we consider the use of machine learning to address two goals.

• Goal (i): Determine whether or not a state variable of a time evolving system causally influences another state variable.

• Goal (ii): Determine the “strength” of such causal influences.

In the terminology of ML, Goal (i) is referred to as “classification ML,” and Goal (ii) is referred to as “regression ML.” These goals have previously been of great interest in many applications (e.g., economics,³ neuroscience,⁴ genomics,⁵ climate,⁶ etc.). Many past approaches have, for example, been based upon the concepts of prediction impact,^3,4 correlation,^7–9 information transfer,^10,11 and direct physical perturbations.^12,13 Other previous works have investigated the inference of network links from time series of node states assuming some prior knowledge of the form of the network system and using that knowledge in a fitting procedure to determine links.^9,14–17 In addition, some recent papers address network link inference from data via techniques based on delay coordinate embedding,¹⁵ random forest methods,¹⁸ network embedding algorithms,¹⁹ and feature ranking.²⁰ In this paper, we introduce a technique that makes the use of an ML training process in performing predictive and interpretive tasks and attempts to use it to extract information about causal dependencies. In particular, here, we use a particular type of machine learning (ML) called reservoir computing, an efficient method of time-series analysis, which has previously been successfully used for different tasks, e.g., prediction of chaotic dynamics^21–23 and speech recognition,^24,25 to mention a few. In our case, a “reservoir” dynamical system is trained such that it becomes synchronized to a training time series data set from the unknown system of interest. The trained reservoir system is then able to provide an estimation of the response to perturbations in different parts of the original system, thus yielding information about causal dependencies in the actual system. We will show that this ML-based technique offers a unique and potentially highly effective approach to determining causal dependencies. Furthermore, the presence of dynamical noise (either naturally present or intentionally injected) can very greatly improve the ability to infer causality,^14,15 while, in contrast, observational noise degrades inference.

We begin by considering the very general case of an evolving, deterministic, dynamical system whose state at time $t$ is represented by the $M$-dimensional vector $z(t)=[z1(t),z2(t),…,zM(t)]T$⁠, where $z(t)$ evolves via a system of $M$ differential equations, $dz(t)/dt=F(z(t))$⁠, and has reached a statistically steady dynamical state (perhaps chaotic). In this context, we frame the issue of causality as follows: Will a perturbation at time $t$ applied to a component $zi$ of the state vector $z(t)$ [i.e., $zi(t)→zi(t)+δzi(t)$] lead to a subsequent change at a slightly later time, $t+τ$⁠, of another scalar component $zj$ [i.e., $zj(t+τ)→zj(t+τ)+δzj(t+τ)$]; and how can we quantify the strength of this dependence? This formulation might suggest comparison of the evolutions of $z(t)$ that result from two identical systems, one with, and the other without, the application of the perturbation. However, we will be interested in the typical situation in which such a comparison is not possible, and one can only passively observe (measure) the state $z(t)$ of the (single) system of interest. Aside from knowing that the dynamics of interest evolves according to a system of the form $dz/dt=F(z)$⁠, we assume little or no additional knowledge of the system and that the available information is a limited-duration past time series of the state evolution $z(t)$⁠. Nevertheless, we still desire to deduce causal dependencies, where the meaning of causal is in terms of responses to perturbations as defined above. Since, as we will see, accomplishment of this task, in principle, is not always possible, our approach will be to first propose a heuristic solution and then numerically test its validity. The main message of this paper is that our proposed procedure can be extremely effective for a very large class of important problems. We will also delineate situations where our procedure is expected to fail. We emphasize that, as our method is conceptually based on consideration of responses to perturbations, in our opinion, it provides a more direct test of what is commonly of interest when determining causality than do tests based on prediction impact, correlation, or entropy metrics.

Furthermore, although the setting motivating our procedure is for deterministic systems, $dz/dt=F(z)$⁠, we will also investigate performance of our procedure in the presence of both dynamical noise [i.e., noise added to the state evolution equation, $dz/dt=F(z)$] and observational noise [i.e., noise added to observations of $z(t)$ used as training data for the machine learning]. Both types of noise are, in practice, invariably present. An important result from our study is that the presence of dynamical noise can very greatly enhance the accuracy and applicability of our method (a similar point has been made in Refs. 14 and 15), while observational noise degrades the ability to infer causal dependence.

To more precisely define causal dependence, we consider the effect of a perturbation on one variable on the other variables as follows. Taking the $jth$ component of $dz/dt=F(z)$⁠, we have

for $j=1,2,…,M$⁠. Perturbing $zi(t)$ by $δzi(t)$⁠, we obtain, for small $τ$⁠, that the component of the orbit perturbation of $zj$ at time $(t+τ)$ due to $δzi$ is

We define the Short Term Causal Dependence (STCD) metric, $fji$⁠, of $zj$ on $zi$ by

where $⟨(…)⟩$ denotes a long time average of the quantity $(…)$ over an orbit and the function $G$ is to be chosen in a situation-dependent manner. For example, later in this paper, we consider examples addressing Goal (i) [where we want to distinguish whether or not $∂Fj(z)/∂zi$ is always zero] for which we use $G(q)=|q|$⁠, while, when we consider an example addressing Goal (ii) and are concerned with the time-averaged signed value of the interaction strength, we then use $G(q)=q$⁠. In either case, we view $fji$ as quantifying the causal dependence of $zj$ on $zi$⁠, and the key goal of this paper will be to obtain and test a machine learning procedure for estimating $fji$ from observations of the state evolution $z(t)$⁠. For future reference, we will henceforth denote our machine learning estimate of $fji$ by $f^ji$⁠. In the case of our Goal (i) experiments, where $G(q)=|q|$⁠, we note that $fji$ defined by (1) is an average of a non-negative quantity and thus, $fji≥0$⁠, as will be our estimate, $f^ji≥0$⁠. Furthermore, for this case, we will define STCD of $zi$ on $zj$ by the condition, $fji>0$⁠, and, when using our machine learning estimate $f^ji$⁠, we shall judge STCD to likely apply when $f^ji>ϵ$ where we call $ϵ>0$ the discrimination threshold. In the ideal case $(f^ji=fji)$⁠, the discrimination threshold $ϵ$ can be set to zero, but, in practice, due to the error in our estimate, we consider $ϵ$ to be a suitably chosen positive number. We note that, in the ideal case, $ϵ=0$ can be regarded as a test for whether or not $Fj(z)$ is independent of $zi$⁠.

As a demonstration of a situation for which the determination of STCD from observations of the motion of $z(t)$ on its attractor is not possible, we note the case where the attractor is a fixed point (a zero-dimensional attractor). Here, the measured available information is the $M$ numbers that are the coordinates of the fixed point, and this information is clearly insufficient for determining STCD. As another problematic example, we note that in certain cases, one is interested in a dynamical system that is a connected network of identical dynamical subsystems and that such a network system can exhibit exact synchronization of its component subsystems²⁶ (including cases where the subsystem orbits are chaotic). In the case where such a synchronized state is stable, observations of the individual subsystems are indistinguishable, and it is then impossible, in principle, for one to infer causal relationships between state variables belonging to different subsystems. More generally, in addition to the above fixed point and synchronization examples, we note that the dimension of the tangent space at a given point $z∗$ on the attractor is, at most, the smallest embedding dimension of the part of the attractor in a small neighborhood of $z∗$⁠. Thus, the full $M×M$ Jacobian of $F(z)$ at $z∗$ cannot be precisely determined from data on the attractor when the local attractor embedding dimension at $z∗$ is less than $M$⁠, which is commonly the case. Thus, these examples motivate the conjecture that to efficiently and accurately infer STCD, the orbital complexity of the dynamics should be large enough so as to encode the information that we seek. Note that these considerations of cases where inference of STCD is problematic do not apply to situations with dynamical noise, e.g., $dz/dt=F(z)+(noise)$⁠, as the addition of noise may be roughly thought of as introducing an infinite amount of orbital complexity. Alternatively, the addition of noise increases the embedding dimension of the data to that of the full state space, i.e., $M$⁠.

We base our considerations on a type of machine learning called reservoir computing, originally put forward in Refs. 27 and 28 (for a review, see Ref. 29). We assume that we can sample the time-series data $z(t)$ from our system at regular time intervals of length $τ$ so that we have a discrete set of observations ${z(0),z(τ),z(2τ),…}$⁠. To begin, we first describe a reservoir-computer-based machine learning procedure in which the reservoir computer is trained to give an output $z^(t+τ)$ in response to an $M$-dimensional input $z(t)$ as illustrated in Fig. 1.

For our numerical tests, we consider a specific reservoir computer implementation (Fig. 1) in which the reservoir consists of a network of $R≫M$ nodes whose scalar states, $r1(nτ),r2(nτ),…,rR(nτ)$⁠, are the components of the $R$-dimensional vector $r(nτ)$⁠.

The nodes interact dynamically with each other through an $R×R$ network adjacency matrix $A$⁠, and their evolution is also influenced by coupling of the $M$-dimensional input $z(nτ)$ to the individual nodes of the reservoir network by the $M×R$ input coupling matrix $Win$ according to the neural-network type of evolution equation (e.g., Refs. 21–23 and 29–31),

where $tanh⁡(v)$ for a vector $v=(v1,v2,v3,…)T$ is defined as $(tanh⁡v1,tanh⁡v2,tanh⁡v3,…)T$⁠. For proper operation of the reservoir computer, it is important that Eq. (2) satisfies the “echo state property”^21,27,29 (in nonlinear dynamics, this condition is also known as “generalized synchronization”^32–34): given two different initial reservoir states, $r1∗$ and $r2∗$⁠, for the same input time series of $z$⁠, the difference between the two corresponding reservoir states converges to zero as they evolve in time [that is, $|r1(t)−r2(t)|→0$ as $t→∞$⁠, implying that, after a transient initial period, $r(t)$ essentially depends only on the past history of $z$⁠, $z(t′)$ for $t′≤t$⁠, and not on the initial condition for $r$].

Using measured input training data over a training interval of length $Tτ$⁠, which begins after the initial transient period mentioned above, we use Eq. (2) to generate $r(τ),r(2τ),…,r(Tτ)$⁠. We also record and store these determined values $r(nτ)$ along with the corresponding inputs, $z(nτ)$⁠, that created them. The matrices $A$ and $Win$ are regarded as fixed and are typically chosen randomly. In contrast, the $R×M$ output coupling matrix $Wout$⁠, shown in Fig. 1, is regarded as an adjustable linear mapping from the reservoir states $r$ to an $M$-dimensional output vector $z^$⁠,

“Training” of the machine learning reservoir computer then consists of choosing the $RM$ adjustable matrix elements (“weights”) of $Wout$ so as to make $z^(nτ)$ a very good approximation to $z(nτ)$ over the time duration $(τ,2τ,…,Tτ)$ of the training data. This is done by minimization with respect to $Wout$ of the quantity, ${∑n=1T∥z(nτ)−Woutr(nτ)∥2}+β∥Wout∥2$⁠. Here, $β∥Wout∥2$⁠, with $β$ small, is a “ridge” regularization term³⁵ added to prevent overfitting, and $(r(nτ),z(nτ))$ are the previously recorded and stored training data. In general, $R≫M$ is required in order to obtain a good fit of $z^$ to $z(t)$⁠. For illustrative purposes, we now consider the ideal case where $z^=z$ (i.e., the training perfectly achieves its goal).

For the purpose of estimating STCD, we now wish to eliminate the quantity $r$ from the basic reservoir computer system [Eqs. (2) and (3)] to obtain an evolution equation solely for the state variable $z$⁠. To do this, we would like to solve (3) for $r$ in terms of $z$⁠. However, since $R$⁠, the dimension of $r$⁠, is much larger than $M$⁠, the dimension of $z$⁠, there are typically an infinite number of solutions of (3) for $r$⁠. To proceed, we hypothesize that it may be useful to eliminate $r$ by choosing it to be the solution of (3) with the smallest $L2$ norm. This condition defines the so-called Moore-Penrose inverse³⁶ of $Wout$⁠, which we denote $W^out−1$⁠; i.e., the minimum $L2$ norm solution for $r$ is written as $r=W^out−1z$⁠. We emphasize that $W^out−1z$ is not necessarily expected to give the correct $r$ obtained by solving systems (2) and (3). However, from numerical results to follow, our choice will be supported by the fact that it often yields very useful estimates of $fji$⁠.

Now, applying $Wout$ to both sides of Eq. (2) and employing $r=W^out−1z$ to eliminate $r(nτ)$ from the argument of the $tanh$ function in Eq. (2), we obtain a surrogate time $−τ$ map for the evolution of $z$⁠, $z((n+1)τ)=H[z(nτ)]$⁠, where $H(z)=Wouttanh⁡[(AW^out−1+Win)z]$ . Here, we note that we do not claim that this map in itself can be used for time-series prediction in place of Eqs. (2) and (3), which were commonly used in previous works (e.g., Refs. 21–23,30, and 31). Rather, we use it as a symbolic representation of the result obtained after eliminating the reservoir state vector $r$ from Eqs. (2) and (3). In particular, the prediction recipe using Eqs. (2) and (3) is always unique and well defined, in contrast to the above map, where $Wout−1$ is clearly nonunique. Therefore, we use this map only for causality estimation purposes, as described below. Differentiating $H(z)$ with respect to $zi$⁠, we have

where $δij$ is the Kronecker delta, and we propose to use Eqs. (1) and (4) to determine STCD.

In our numerical experiments, the number of training time steps is $T=6×104$ for Figs. 2 and 3 and $T=2×104$ for Fig. 4. In each case, the actual training data are obtained after discarding a transient part of $2×104$ time steps, and the reservoir system sampling time is $τ=0.02$⁠. The elements of the input matrix $Win$ are randomly chosen in the interval $[−0.1,0.1]$⁠. The reservoir is a sparse random network of $R=5000$ nodes for Figs. 2 and 3 and of $R=1000$ nodes for Fig. 4. In each case, the average number of incoming links per node is $3$⁠. Each nonzero element of the reservoir adjacency matric $A$ is randomly chosen from the interval $[−a,a]$⁠, and $a>0$ is then adjusted so that the maximum magnitude eigenvalue of $A$ is $0.9$⁠. The regularization parameter is $β=10−4$⁠. These parameters are adapted from Ref. 23. The average indicated in Eq. (1) is over $1000$ time steps. The chosen time step $τ$ is sufficiently small compared to the time scale over which $z(t)$ evolves that the discrete time series $z(nτ)$ is a good representation of the continuous variation of $z(t)$⁠.

Although we use a specific reservoir computing implementation, we expect that, with suitable modifications, our approach can be adapted to “deep” types of machine learning,² as well as to other implementations of reservoir computing^24,25,37,38 [notably, implementations involving photonics,²⁴ electronics,³⁷ and field programable gate arrays (FPGAs)²⁵].

In order to evaluate the effectiveness of our proposed method, we introduce mathematical model test systems that we use as proxies for the unknown system of interest for whose state variables we wish to determine STCD. We next use the test systems to generate simulated training data from which we determine STCD by our ML technique. We then assess the performance of the technique by the correctness of its results determined from the known properties of the test systems.

We first consider examples addressing our Goal (i) [$G(q)=|q|$ in Eq. (1)], and for our simulation test systems, we consider the case of a network of $N$ nodes and $L$ links, where each node is a classical Lorenz system³⁹ with heterogeneity from node to node, additive dynamical noise, and internode coupling,

The state space dimension of this system is $M=3N$⁠. The coupling of the $N$ nodes is taken to be only from the $y$ variable of one node to the $x$ variable of another node with coupling constant $c$⁠, and $akl(x,y)$ is either 1 or 0 depending on whether or not there is a link from $l$ to $k$⁠. The adjacency matrix $akl(x,y)$ of our Lorenz network (not to be confused with the adjacency matrix $A$ of the reservoir) is constructed by placing directed links between $L$ distinct randomly chosen node pairs. For each node $k$⁠, $hk$ is randomly chosen in the interval $[−h,+h]$⁠, and we call $h$ the heterogeneity parameter. Independent white noise terms of equal variance $σDyn2$ are added to the left-hand sides of the equations for $dx/dt$⁠, $dy/dt$⁠, and $dz/dt$⁠, where, for example, $⟨nkx(t)nk′x(t′)⟩=2δkk′δ(t−t′)$⁠. For $σ=c=h=0$⁠, each node obeys the classical chaotic Lorenz equation with the parameter values originally studied by Lorenz.³⁹ Furthermore, denoting the right-hand side of Eq. (5) by $Fxk$⁠, we have $∂Fxk/∂yl=10c$ or 0, depending on whether there is, or is not, a link from $yl$ to $xk$⁠.

Since in this case, the derivative $∂Fxk/∂yl$ is time independent, $⟨|∂Fxk/∂yl|⟩$ is also either $10c$ or 0, and adopting the notation $fkl(x,y)=⟨|∂Fxk/∂yl|⟩$⁠, we denote its machine learning estimate by our previously described procedure by $f^kl(x,y)$⁠. For a reasonably large network, the number $N2−N$ of ordered node pairs $(k,l)$ of distinct nodes is large, and we consequently have many values of $f^kl(x,y)$⁠. Bayesian techniques (see Ref. 40 and references therein) can be applied to such data to obtain an estimate $L^$ for the total number of links $L$⁠, and one can then set the value of $ϵ$ so that there are $L^$ values of $f^kl(x,y)$ that are greater than $ϵ$⁠. Less formally, we find that making a histogram of the values of $f^kl(x,y)$ often reveals a peak at zero and another peak at a higher positive value with a large gap or discernible minimum in between. One can then estimate $ϵ$ by a value in the gap or by the location of the minimum between the peaks, respectively. For simplicity, in our illustrative numerical simulations to follow, we assume that $L$ is known [approximately equivalent to the case that $L$ is unknown, but that a very good estimate (⁠$L^$⁠) has been obtained].

In this paper, we have formulated and tested a new, highly effective, machine-learning-based approach for inferring causal dependencies of state variables of an unknown system from time-series observations of these state variables. A key finding is that the effectiveness of our approach is greatly enhanced in the presence of sufficient dynamical noise, provided that the deleterious effect of observational noise is not too great. The competition between the opposing effects of these two types of noise will likely be the essential key factor determining the success or failure of causality inference in many of the most important situations of interest (e.g., in neuroscience and genomics). Much work remains to be done to more fully address the utility of our method. In particular, further numerical tests on diverse systems, and, especially, experimental studies in real world applications, will ultimately determine the circumstances under which the method developed here will be useful.

This work was supported by the U.S. National Science Foundation (NSF) (Grant No. DMS 1813027). The authors acknowledge useful discussion with Sarthak Chandra, Amitabha Sen, Woodrow Shew, Nuno Martins, Adrian Papamarcou, Erik Bollt, and, especially, Brian Hunt.