Neural graphical modelling in continuous-time: consistency guarantees and algorithms
Authors: Alexis Bellot
Kim Branson
Mihaela van der Schaar
Venue: ICLR 2022 Poster (2022)
Keywords: Dynamical systems
graphical modelling
structure learning

Abstract: The discovery of structure from time series data is a key problem in fields of study working with complex systems. Most identifiability results and learning algorithms assume the underlying dynamics to be discrete in time. Comparatively few, in contrast, explicitly define dependencies in infinitesimal intervals of time, independently of the scale of observation and of the regularity of sampling. In this paper, we consider score-based structure learning for the study of dynamical systems. We prove that for vector fields parameterized in a large class of neural networks, least squares optimization with adaptive regularization schemes consistently recovers directed graphs of local independencies in systems of stochastic differential equations. Using this insight, we propose a score-based learning algorithm based on penalized Neural Ordinary Differential Equations (modelling the mean process) that we show to be applicable to the general setting of irregularly-sampled multivariate time series and to outperform the state of the art across a range of dynamical systems.

Published as a conference paper at ICLR 2022
NEURAL GRAPHICAL MODELLING IN CONTINUOUS- TIME: CONSISTENCY GUARANTEES AND ALGORIHHM
Alexis Bellot"
Kim Branson GlaxoSmithKlein, USA
Columbia University, USA
ab53050col 1e13005ecolubii.. Lumbia edu
k kim, m. www.brannon@gss.cco com
Mihaela van der Schaar University of Cambridge, UK
The Alan Turing Institute, UK University of California, Los Angeles, USA mv4720 cam. cw..cc.u....k..........
ABSTRACT
The discovery of structure from time series data is a key problem in fields of study working with complex systems. Most identifiability results and learning algorithms assume the underlying dynamics to be discrete in time. Comparatively few, contrast, explicitly define dependencies independently of the scale of observation and of the regularity of sampling. In this paper, we consider score-based structure learning for the study of dynamical systems. We prove that for vector fields parameterized in a large class of neural net- works, least squares optimization with adaptive regularization schemes consistently recovers directed graphs of local independencies in systems of stochastic differ- ential equations. Using this insight, we propose a score- score-based learning algorithm based on penalized Neural Ordinary Differential Equations (modelling the mean process) that we show to be applicable to the general setting of trreularly-sample multivariate time series and to outperform the state of the across a range of dynamical systems.
infinitesimal intervals of time,
INTRODUCTION
This paper deals with learning directed graphs from combination of temporal data and assumptions on the parameterization of the underlying structural dynamical system. Graphical models can offer parsimonious, interpretable representation of the dynamics of stochastic processes, and have proven to be especially useful in problems involving complex systems, non-linear associations and chaotic behaviour that are characteristic in a wide array of applications in biology (Trapnell et al., 2014; Qiu et al., 2017; Bracco et al., 2018; Raia, 2008; Qian et al., 2020), neuroscience (Friston et al., 2003; Friston, 2009) and climate science (Runge, 2018: 8; Runge et al., 2019). In these contexts, inferring graphical models from temporal data subject to practical limitations as how finely and regularly each variable can be measured over time a longstanding challenge.
Time series data is often assumed to be sequence of observations from underlying process evolving continuously in time. This underlying representation is fundamental to define the semantics of dependencies between sequences. Time defines an asymmetry between dependencies in dynamical systems, distinguishing between local, direct dependencies that occur over infinitesimal time intervals not mediated by other variables in the system and indirect dependencies that necessarily occur over longer time frames. In many applications, the underlying structural model is formalized as the state of random variables (e.g. x(t) E IR contemporaneously influencing the rate of change of the same or other variables (e.g. dx(())
dx(tt ((x(t))d + dw((),
x x(0) = Xp.
t E [0, T].
(1)
WWrk primarily conducted while the University of Cambridge and the Alan Turing Institute.

---

Published as a conference paper at ICLR 2022
Gycy
GNC
Gren
Time
(a) Data and adjacency matrix.
(b) Estimate with At = 0.05.
(c) Estimate with At = 0.25.
Figure 1: Visual comparison of the true and learned adjacency matrices G' of a 10-variable Lorenz system (see section 4.1). GwaM (the proposed continuous-time approach) and GNec (Neural Granger Causality (Tank et al., 2018)) are estimates of continuous and discrete-time algorithms respectively. Panel (a) shows data sample and the true adjacency matrix, panel (b) shows estimates with higher frequency of observation (At t, t4-1 0.05) and panel (c) with lower frequency of observation (At 0.25). The heat scale gives the strength of estimated functional interactions. An explicitly continuous-time model is more accurate and more robust to the sampling frequency than discrete-time alternatives for graphical modelling in dynamical systems.
where w(t) a d-dimensional standard Brownian motion and Xo is Gaussian random variable independent of w(t). The functional dependence structure of the vector field defines directed graph G and associated adjacency matrix G E {0, 1]axd, i.e., Gy if and only if x appears as an argument of f: [].. The problem of structure learning to search over the space graphs compatible with the data, but the pattern of observation in dynamical systems emphasize a number of differences with respect classical graphical modelling with static data or explicitly discrete-time stochastic process.
Observed data is sampled at a sequence of (often irregular) time points (f1,..., fn)) and is systemat- ically subsampled. Most work graphical modelling with time series data assume fundamentally discrete parameterization of the underlying structural model (e.g. based on vector autoregression models). Associations discrete-time in general do not correspond the structure of the underly- ing dynamical system and are highly dependent on the interval between observations. The same subsampled discrete model may disaggregate several continuous models, which are observation- ally equivalent at the subsampled frequency. see e.g. (Runge, 2018; Gong et al., 20155; Danks & Plis, 2013) and a worked example in Appendix A. The realm of problems that involve irregularly- sampled data are fundamentally out of scope in discrete-time in general. We complement this point Figure with an illustration of our performance results comparing a state of the art discrete-time graphical modelling method with our proposed continuous-time counterpart that is shown to be more accurate and more robust to the frequency and irregularity of sampling.
Diserete samples are not independent which can (and does) increase the sample complexity. An increasing sample frequency will produce an increasing number of distinct samples. However, samples become more dependent, and intuitively one expects that there limited information to be harnessed from given time interval [0, 7". Learning performance depends on the number of independent samples which a function both of the number of samples and the length of the observed interval T.
Non- parametrii graphical modelling in dynamical systems is relatively unexplored. Existing approaches rely on specific model assumptions (e.g. linearity, additivity) establish the consistency structure recovery even though flexible model families, such as neural networks, are increasingly used in related problems such as feature selection and graphical modelling with static data. In addition, consistent derivative approximations are typically required for consistency arguments which in practice involve choices on the smoothness of the interpolated curve and makes two-step strategies far from automatic and applicable in general dynamical systems.
Contributions. This paper establishes the consistency of score-based recovery of G when an analytic deep neural network model is imposed for f (such as feed-forwart networks with multiple hidden layers and convolutional neural networks) under general observation patterns including irregular sampling. In particular, consider penalized optimization problems of the form,
arg min fo
[|x(t)) x(1))
subject to (fo) S and dR(() fo (x(())) (2)
2

---

Published as a conference paper at ICLR 2022
where the observation process (x(()),... x(t,)) is given by irregular sequence of time points 0 S f1 <.... < tn S T. Pa,r ) is an adaptive group lasso constraint on the parameter space of fa. We analyze this problem with fixed dimension d and increasing sample size and horizon T the sample complexity of this problem depending both on the frequency of sampling n as well as on the time horizon T.
second contribution is to propose an instantiation of this method using differential equations with vector fields parameterized by neural networks (Chen al., 2018) model the mean process of (1) with the advantage of implicitly inferring variable derivatives instead of involving a separate approximation step (that common in the dynamical systems literature). This construction shows that, empirically, graphical models in continuous- tim can be inferred accurately in large range of settings despite rregularl-ssapplee multivariate time series data and non-linear underlying dependencies.
Code associated with this work may be found at http:: com/ axxxi weleexisbeee sbel lot and at https: /ttt:://gittuu com http:///nnderr.c..eaaa.aaaaaaaaa
2 RELATED WORK
A substantial amount of work devoted to graphical modelling has considered the analysis of penalized least squares and its variants, most prominently in the high-dimensional regression literature with i.i.d data, see e.g. (Friedman et al., 2008; Zou, 2006; Zhao & Yu, 2006). Closely related to our results are a number of extensions that have considered parameter identification neural networks, using for instance sparse one-to-on linear layers (Li et al., 2016), group lasso constraints of the input layer of parameters (Zhang et al., 2019) and input to output residual connections (Lemhadri et al., 2021). For large class of neural networks Dinh & Ho (2020) proved the consistency of adaptive regularization methods. The distinction with our formalism in (2) is that the observations are not corrupted by i.i.d. noise (since successive samples are correlated) and therefore standard concentration inequalities are not sufficient.
Learning graphical models with dependent noise terms is also a topic of significant literature in the context of Granger causality, proposed by Granger (1969) and also popularized by Sims (1980) within autoregressive models. Various authors have considered the consistency of penalized vector autoregression models and proposed tests of Granger causality using parameter estimates in these models, see e.g. (Nardi & Rinaldo, 2011; Kock & Callot, 2015; Adamek et al., 2020; Chernozhukov et al., 2019), and extended some of these approaches to models of neural networks, see e.g. (Tank et al., 2018; Khanna & Tan, 2019; Marcinkevičs & Vogt, 2021) (without however proving consistency of inference). Methods exist also using conditional independence tests such as those given by Runge et al. (2017) and transfer entropy principles originating in Schreiber (2000). The conceptual and statistical contrasts between discrete and continuous accounts of the underlying structural model substantial and are discussed in the Appendix A.
In the context of differential equations, penalized regression has been explored using two-stage collocation methods, first proposed by Varah (1982), by which derivatives are estimated on smoothed data and subsequently regressed on observed samples for inference. The consistency of parameter estimates has been established for linear models in parameters, as done for example in (Ramsay et al., 2007; Chen al., 2017; Wu et al., 2014; Brunton et al., 2016). From a modelling perspective, our approach in contrast is end-to-end, coupling the estimation of the underlying paths x and the vector field f. Graphical modelling has also been considered for linear stochastic differential equations by Bento et al. (2010). Similarly to the discrete-time literature, proposals exist for recovering non-linear vector fields via neural networks (see e.g. (Raissi et al., 2017: Bellot van der Schaar, 2021) I)) and Gaussian processes (see e.g. (Heinonen et al., 2018; Wenk et al., 2020)) but we not aware of any identifiability guarantees.
GRAPHICAL MODELLING IN CONTINUOUS
We consider the underlying structure of an evolving process to be described by multivariate dynamical system of d distinct stochastic processes x = (#),..., [0, T] + x with each instantiation in time (t) for j 1,...,d d and >> 0 defined in a bounded open set X' c Ik.
3

---

Published as a conference paper at ICLR 2022
Definition (Neural Dynamic Structural Model (NDSM)). We say that (Z1,.... F.) [0 X" follows a Neural Dynamic Structural Model if there exist functions f1..... E F such that fg x + R and,
dz,(t)= fj((())t du,,t),
x(ta)= xx. Xo.
t [0
(3)
with F defined the space of analytic defined in bounded, real-valued intervals and w,(t) (t) standard Brownian motion independently generated across processes j'.
neural networks with sets of parameters 0 E θ
We will write fep (Ji,.... Sa] for the true underlying vector field, parameterized by a set of parameter values lg. It will be useful define each layer of each network precisely. Let A E IRdxh denote the d x h weight matrix (we omit biases for clarity) in the input layer of fj.j = 1,..., dd
Let Ain E R for m = 2..... . M 1, denote the weight matrix of each hidden layer, and let Asrr E 12x11 be the h x 1 dimensional output layers of each sub-network such that,
f = o( o...- (0( ..hr)
j I.....d,
(4)
where o(-) is an analytic activation function (e.g. tanh, sigmoid, arctan, softplus, etc.) and X E Raxd is the sequence of n d-dimensional instantiations of x.
Assumption 1 (Observation process). The data in practice, is a partial sequence of observations of'x at n time points (f1....,tt, fn) sampled from a temporal point process with positive intensity such that,
(xI,...,x.) Xn.) ~ NN,, En),
(5)
with dependency structure encoded in En E R The closer in time two observations are the more closely correlated we can expect them be. We assume the data to normalized, i.e. diagonal elements of Σn to be equal to I. are the instantiations of the mean process that can be described by an system of ordinary differential equations dx(t)= = f(x(t))dtt x(0) = Xo, t E [0, T).
Time points at which observations are made are thus themselves assumed stochastic, driven by an independent temporal point process with intensity lima--o (rObserrvation [1, l + dij(He) > for any t >> 0 with respect to a filtration H, that denotes sigma algebras generated by any sequence of prior observations. Perfectly homogencous and systematic subsampling has measure zero under this probability model. This is important because it will enable, in principle, infer local conditional independencies (defined below) arbitrarily well with increasing sample size.
3.1 GRAPHICAL PRESENTATION
The stochastic process x by itself defines a local independence model that can be used to characterize asymmetric dependencies within stochastic processes (Eichler & Didelez, 2012; Eichler, 2013; Didelez, 2012).
Definition 2 (Local independence). A process is locally independent of given if. for each time point I, the past up until time of z gives us the same predictable information about as the past of x and y until time I, where Hily, 2) the filtration generated by y and up time l.
=))
This independence structure may be represented by a (cyclic) directed graph G = -V, E), where each process is associated with a distinct vertex in V and there is directed edge (I + xx) and only if there exist no conditioning subset of V such Iy E V is locally conditionally independent of z E V.
Lemma (Uniqueness of local independence graphs, Proposition 3.6 (Mogensen et al., 2020)). In the context of Neural Dynamic Structural models, two processes are locally dependent given any subset of other processes if and only if x& appears in the differential equation of 2; i.e. |10afs * O. Moreover, for any such that [0=f; L2 there exists an equivalent vector fielo f such that the euclidian norm of its column vectors [|[A[]= O.
Proof. All proofs are given in Appendix B.
For analytic function spaces F, the vector field is locally Lipschitz, i.c., the system is stable and the diffusion process has unique stationary measure that Gaussian. We assume unique solutions also T co.
4

---

Published as a conference paper at ICLR 2022
This Lemma specifies equivalence relation between functional dependence graphs given by the underlying dynamical system and local independence graphs'. is clear that enforcing [|[AA]] 0 will remove the local dependence of the j-th stochastic process on the k-th stochastic process but is not the case that all functions f with this particular local independence (10x,fyl L 0) necessarily have zero-valued k-th column in its parameters Ai. This proposition shows that such cases there exists equivalent vector field (i.e. that defines the exact same input-output map) that does have +||2 = 0. Local independence graphs may be recovered by searching for such a solution, in theory, if the stochastic process is fully observed.
In practice, with finite samples and complex functions, the map between model and data may not necessarily be identifiable, and a priori should not be expected for highly parameterized neural networks (e.g. simple rearrangement of the nodes in the same hidden layer leads to configuration that produces the same mapping as the generating network). We define next local consistency desirable and target property for estimators in practice.
Definition 3 (Local consistency). An estimator fa - -f Sa) is locally consistent for any ō >>0 there exists N& and T' such that for > Ns and T > Ts. we have |0x]z]| [g + 02, for all k.j E (1....,d) d such that xy is locally significant for 2j, and have 110xffllz, = 0 otherwise, with probability least 8.
3.2 SINITE-SAMPLE IDENTIHABILITY
Write R., (fo) =
|x(t,) ( (the dependence on fo implicit in x) and R(fa) for its population counterpart. The difficulty arises from the geometry of the loss function around the set of loss minimizers,
(0 E θ R(fs) R )),
where Ô is the parameter space. The @ (of all weight vectors that produce the same input-output map as the generating model) can be quite complex and the behavior of a generic estimator in this set not necessarily reflect local dependencies.
It is possible, however, to constrain the solution space to a subset of "well-behaved" optima for which G is uniquely identifiable even if the full set of parameters 0 not. It is sufficient for identifiablity of the local independence model to recover the group structure of input layer parameters [A]] exactly. In the context of analytic vector fields f. we define the adjacency matrix G {0, 1 1]dxd associated with G such that Gaj * O if and only if [j0=/ * 0. If this pattern can be recovered exactly, then the inferred corresponds to the functional structure. Optimization with this group structure desiderata has been often considered before as penalized optimization problem,
arg min fo
Rx, (fa). subject to dx(()) (t) and p(.)) S n
where we have suppressed the dependence of p on the sample size and time horizon for readability. Two popular constraints are the group lasso (GL) and the adaptive group lasso (AGL), see e.g. (Zou, 2006; Zhao & Yu, 2006), defined as,
PCL(fo) == AGL Σ |[[A]]-ale, k
PAGL(fo) = AAGL
[[[A)]+*12.
A 11
respectively, where A is the GL estimate to problem (7), λOE AAGL determine the regularization strength and may vary with n and T, 7 > 0 and || ll2 is the Euclidian norm. As with other adaptive lasso estimators, AGL uses its base estimator to provide rough data-dependent estimate to shrink groups of parameters with different regularization strengths. As n and grow, the weights for non-significant features get inflated while the weights for significant ones remain bounded, allowing AGL to exactly identify significant parameters.
2 Note, however, that this is not true in general for marginalization of the local independence graph (that we do not consider), i.c. the context of unobserved processes (Mogensen et al., 2020), in which case several graphs may encode the same set of local independence relations and only equivalence classes may be identifiable. Os denotes the partial derivative with respect the k-th argument of fz, Il !|e, is the functional Ta norn.
5

---

Published as a conference paper at ICLR 2022
The following generalization bound will be useful to define convergence rates for penalized solutions.
Lemma (Generalization bound). Assume En to be invertible and let a = (ai...., .n)) such that
a1 > an > 0 are its eigenvalues. For any d > O, there exists C3 >> 0 such that,
||al22
(fo - -R((fo)
S C.
V
log
n1
(8)
||a1|12
with probability at least 1- 8.
Remark on interpretation. Note that ||a||2 S ||a||1 = n (the data is assumed to be scaled to have variance for all observations) and that larger values of ||a||2 occur with a greater difference in magnitude in the entries of a. The difference in magnitude in the principal components of X is defined by the dependence between samples. A strong dependence between samples (as would be expected with frequently observed time series) leads to proportionally larger magnitude of (the first component of a)) as more variance the data is explained by a single direction variation and thus decreases the effective sample size making ||||| closer to n. The bound formalizes the trade-off between the number of samples and their dependence. By increasing the observation frequency one can produce arbitrarily large number of distinct samples but samples become more dependent and therefore less useful for concentration of the empirical error around its population value (unless one simultaneously increases the time horizon T').
We now show convergence and local consistency of the adaptive group lasso penalized estimator. The following lemmas use a similar proof technique Dinh & Ho (2020) with the difference that the convergence speeds differ due to sample dependency.
Lemma 3 (Convergence of Adaptive Group Lasso). Let On E θ be the parameter solution of (7) with adaptive group lasso constraint. For any ô O, assuming that AAGL. there exists > 0, 0g > 0, No > 0 and T3 > O such that,
)) 1
(MaCe +
[|al12) n
min 6E
10 011 s C's
log
||a||2
with probability at least
Lemma 4 (Local consistency of Adaptive Group Lasso) Let > 0, > 0, > 0. AAGL 21( 1 )>//w=)), and MaG. ((A2t)), then the adaptive group lasso (solution to problem (7)) is locally consistent.
Remark on interpretation. There exists a well defined time complexity i.e., a minimum time interval such that, observing the system at appropriate frequency enables us to reconstruct the network with high probability. The sample complexity is inversely proportional to the time spacing between samples. Lemma 3 implies that with increasing sample size and time horizon the graph G detined such that [G]i if and only if [|]]]+]] = the local independence graph with high probability. The group lasso will generally not be locally consistent because it forces all parameters be equally penalized. Some evidence for this claim in the context of feature selection was provided by Zou (2006).
Lemma 4 gives an asymptotic guarantee on structure learning. With finite samples, for the estimated structure to have good accuracy for G. we have to require that the local dependencies between processes (that define the non-zero entries of G) are "sufficiently large" for a given effective sample size. We make a minimum restricted strength assumption whose form reads as,
A min > Co (MaiL +
(1cell2 n
n
V log
(10)
[|a||2
where we have defined |min := min [[]]]][]]] j,k = I,...,d, 1|0rS + 0} to be the minimum column norm of first layer parameters among all locally dependent stochastic processes. This condition on the design of the problem is not testable but versions of are essentially necessary in the context of structure learning via parameter estimation (Van de Geer et al., 2011; Zhao & Yu, 2006), and allows us to describe next a guarantee on the finite sample consistency of structure learning with the Adaptive Group Lasso.

---

Published as a conference paper at ICLR 2022
(Finite sample local consistency of Adaptive Group Lasso). Under the conditions of Lemma with the additional minimum restricted strength assumption in (10) on the problem design for particular values of n and a. the Adaptive Group Lasso recovers the structure G exactly with high probability.
Lemma 5
3.3 ALGORITHM: NEURAL GRAPHICAL MODELLING
Neural ODEs proposed by Chen et al. (2018) are a family of continuous-time models which can be used to define the mean process of stochastic differential equation explicitly to be the solution to an ODE initial-value problem in which fg is free parameter specified by a neural network. For each estimate of fø, the forward trajectory can be computed using any numerical ODE solver:
8(()). ODESSoo x(to), (1,...,tn))
(11)
and thus implicitly enforcing the constraint dx(tt) = fa (x(t)) while optimizing for R. (fo) and thus solving for (2). Gradients with respect 0 may be computed with adjoint sensitivities and gradient descent algorithm can be used to backpropagate through the ODE solver and the continuous state dynamics update the parameters of fo, shown by Chen et al. (2018).
Remark on optimization. The adaptive group lasso constraint not differentiable and non-separable which precludes applying coordinate optimization algorithms. However, wide variety of techniques from optimization theory have been developed to tackle this case. A general way of doing so is through proximal optimization, see e.g. (Parikh & Boyd, 2014) that leads exact zeros in the columns of the input matrices without having to use cut-off value for selection. The proximal step for the group lasso penalty given by a group wof-thhresoolding operation on the input weights and can be interleaved with conventional gradient update steps. Please find all details in the Appendix.
We call this algorithm for structure learning and graphical modelling the Neural Graphical Model (NGM).
4 EXPERIMENTS
This section makes performance comparisons on controlled experiments designed to analyzed important challenges for graphical modelling with time series data: the irregularity of observation times, the sparsity of observation times, the non-linearity of dynamics, and the differing scale of processes in system.
We benchmark NGM against a variety of algorithmss, namely: Three representative vector autoregres- sion models: Neural Granger causality (Tank et al., 2018) in two instantiations, one based on feed forward neural networks (NGC-MLP) and one based on recurrent neural networks (NGC-LSTM),, and the Structural Vector Autoregression Model (SVAM, (Hyvärinen al., 2010), an extension of the LiNGAM algorithm to time series). A representative independence- based approach to structure learning with time series data: PCMCI (Runge et al., 2017), extending the PC algorithm. A represen- lative two-stage collocation method we call Dynamic Causal Modelling (DCM) in which derivatives first estimated on interpolations of the data and penalized neural network is learned to infer G (extending the linear models of (Ramsay et al., 2007; Wu et al., 2014; Brunton et al., 2016)).
Metric. We seek recover the adjacency matrix of local dependencies G between the state of all variables and their variation. All experiments are repeated 100 times and we report mean and standard deviations of the false discovery rate (FDR) and true positive rate (TPR) in recovery performance of G. Thresholds for determining the presence and absence of edges in G were chosen for maximum Fi score. For applications in biology, false positives and false negatives can have very different failure interpretations; we choose to report both TPR and FDR explicitly to emphasize the trade-offss of each method when used in practice. Comparisons based on the area under the ROC curve (evaluating the whole range of possible thresholds), experiments comparing different regularization schemes, hyperparameter configurations, dimensionality of processes and run times, as well as details regarding neural network architectures and implementation software may be found in Appendix C.
7

---

Published as a conference paper at ICLR 2022
DCM LSTM
10
06
MGC
04
PCMEI
POMO
SHAM
SWAM
NGH
+
MGA
0.2
04
60
30
30
001
0.5
01
015
0.5
001
0.55
0.1
0.15
0.5
% observations removed
observations removed
Interval observation
Interval observation
Performance with irregular observation times.
(b) Performance with varying interval of observation.
True positive (higher better) and false discovery (lower better) performance comparisons on Lorenz's model. Thresholds are chosen for maximum Fi score. We omitted plotting NGC-MLP which gave very similar results to NGC-LSTM. NGM the proposed approach.
Figure 2:
4.1 LORENZ'S S CHAOTIC MODEL
We begin by considering irregularly sampled data and sparsely sampled data to investigate the benefit of modelling dynamics continuously in time.
We use Lorenz's model (Lorenz, 1996) as an example of the kind of chaotic systems observed in biology e.g.. electrodynamics of cardiac tissue (Goldberger West, 1987), gene regulatory networks (Heltherg et al., 2019), etc. The continuous dynamics in a d-dimensional Lorenz model are,
d
dt 2;( (1)+1(t) z.-2(1)) 11) (t) F odw;
i= 1=....., d
where xt) == (t), (t) == (t), xd+1 (t) := (t), F is forcing constant which determines the level of non- non-lineariyy and chaos the series and (t) is standard independent Brownian motion across = 1..... d. The initial state of each variable is sampled from a standard
Gaussian distribution, dl set to 10, to 10 and a to 0.5. We illustrate sample trajectories of this system in Figure 1.
Irregularly sampled data is generated by removing randomly percentage of the regularly sampled data (with 0.1 time interval between observations and 1000 observations). For consistency, discrete-time methods here use cubic spline interpolations evaluated at regular time intervals.
Frequently and sparsely sampled data is generated by varying the time interval of observation (the data always being regularly sampled).
Results. Performance results given in Figure 2. It is important here two look at both TPR and FDR panels for each experiment. For instance, a model returning always a fully connected graph G' will have TPR= 1 but FDR= O. Thus looking at both measures together, NGM significantly improves performance over competing approaches. And moreover, NGM's 's performance more robust (worsens less) with increasingly irregular sampling and with increasingly sparse data. The behaviour of discrete-time methods is highly heterogeneous. NGC-LSTM, SVAM and PCMCI are highly dependent the interval of observation, both too frequent and too sparse measurement times leading to poor FDR. Similarly, as we introduce more irregular sampling TPR decreases and FDR increases which we hypothesize is due eгror introduced in the interpolation. This pattern is consistent with our intuition (illustrated in Figure 1) that direct and indirect effects become indistinguishable in discrete time and thus have high FDR with irregular or sparse data. In the case of DCM, it has good performance with frequently observed time series and regular data but rapidly deteriorates otherwise which we hypothesize is due to worsening approximations of derivatives in those cases.
4.2 RÖSSLER'S ' HIYPERCHAOTIC MODEL
Next, we consider data with non-linear dynamics, said to exhibit hyperchaotic behaviour, to demon- strate the flexibility of learned functional relationships.
Chaotic systems, such as Lorenz's 's model, are characterized by one direction exponential spreading. If the number of directions of spreading is greater than one the hehavior of the system is hyperchaotic
8

---

Published as a conference paper at ICLR 2022
Rössler (d = 10)
Rössler (d = 50)
Glycolysis
TPR
FDR
TPR
FDR
TPR f
FDR 4
NGC-MLP
(.05)
(.06)
.31 (.04)
.67 (.06)
.60 (.04)
.51 (.05)
NGC-LSTM
.49 (.04)
.53 (.04)
.38 (.04)
.64 (.08)
.69 (.04)
.40 (.04)
SVAM
.17 (.03)
.84 (.04)
.03 (.08)
.95 (.09)
.61 (.02)
.53 (.07)
PCMCI
.10 (.03)
.92 (.02)
.09 (.06)
.89 (.06)
.56 (.05)
.43 (.06)
DCM
(.01)
.10 (.04)
.97 (.01)
.31 (.07)
.67 (.04)
.49 (.05)
NGM (ours)
.6 (.01)
.02 (.01)
.95 (.01)
.04 (.01)
.84 (.04)
.44 (.09)
Table 1: Performance comparisons Rössler's model and the yeast Glycolysis model.
(see e.g. (Barrio et al., 2015)), and much more complicated to predict. In practice this has been observed in chemical reactions (Eiswirth et al., 1992) and EEG models of the brain (Dafilis et al., 2013). Rossler (1979) the first hyperchaotic system of differential equations and here consider generalization of this model to arbitrary dimensions and non-linear vector fields as in (Meyer et al., 1997). The d' dimeesnsional generalized Rössler model is given by,
d
di 21 (I) a.i (1) =1) + odwi (4),
d
dt (t) = sin(r;-1 ( odw; i= 2, 2....,
d dt =a(t) = bta(t) (I)- + adwa(t).
We use typical parameters for hyperchaotic behaviour: 0, = 0.1, 4, 2, as in (Meyer et al., 1997). This system is observed over a sequence of 1000 time points with 0.1 time unit interval after randomly initializing each variable to sample from standard Gaussian distribution and d is set to 10 and 50 to evaluate also performance higher dimensions. a is set 0.1.
Results. Performance results are given in Table and the contrast between non-linear and linear methods is stark. NGM continues to strongly outperform other methods with almost perfect recovery of the graph G' in both low and dighediimesional regimes. The strongest baseline is DCM which also models the underlying graph in continuous-time.
4.3 MODEL OF OSCILLATIONS IN YEAST GLYCOLYSIS
We conclude with experiment from the pharmacology literature which makes extensive of dynamical models to determine the interaction patterns of drugs in the body. In these systems, is the difference in scale of different biochemicals that makes structure recovery difficult.
The glycolytic oscillator model is standard benchmark for this kind of systems. It simulates the cycles the metabolic pathway that breaks down glucose in cells. We simulate the system presented in equation (19) by Daniels & Nemenman (2015) defined by 7 biochemical components and fully described in Appendix C.5.
For biology in particular, under our as- sumptions, feature of NGM is that it plicitly discovers and consistent graph- ical structure such that the model can be
Predicted processes
True processes
25 20
used to simulate the expected effect of in- terventions by modifying the weights of trained network which may have a large impact on the design of (laboratory) exper- iments. We show, for illustration, NGM's simulation of the mean behaviour of the glycolytic oscillator in Figure 3, which suc- cessfully recovers the true mean dynamies.
20 10
10
100
50
75
100
Time
Time
Figure 3: True and NGM-estimatee Glycolytic oscillations.
Results. Table shows that performance on this data is more heterogencous. Although NGM outperforms competitive with other methods, FDR is high. Better understanding how model

---

Published as a conference paper at ICLR 2022
structure in systems with different scales (see this explicitly in Figure with numerous variables with values equal to a small fraction of the largest observations) thus remains an important challenge.
5 DISCUSSION
An emerging trend in biology, climate science and healthcare is the measurement of increasing amounts of datatypes (individual gene transcript levels, protein abundances, molecular concentrations of pollution in the air, biomarkers in the hospital, etc.) neous observation patterns. A graphical model that stable over a large number of variables, over a large number of sampling frequencies and observation patterns has attractive properties scientists in all these domains. In this paper, we have discussed graphical modelling from a continuous-tii perspective and have shown the consistency of penalized least squares problems in general models of differential equations with analytic vector fields. As an instantiation of this method, we propose novel graphical modelling algorithm, Neural Graphical Model, that models the latent vector field explicitly with penalized extensions to Neural ODEs and is applicable to general tregularly-samplle multivariate time series. We conclude with some additional remarks.
increasing resolution but with heteroge-
Marginalized local independence graphs. In general, we cannot expect that local independence graphs uniquely identify the underlying functional dependency structure if unobserved processes influence the dynamics of the system. In such cases one can define an equivalence class of local independence graphs which have been characterized recently by Mogensen et al. (2020). We have not considered this setting. Anecdotally however, form of violation of this assumption was included in the Lorenz system (perturbed by constant factor F which not modelled) for which NGM achieves almost perfect structure recovery. Some robustness to violations of the assumption no unobserved processes may thus be expected although we not perform formal investigation of this phenomenon.
Switching dynamical systems. In the context of Neural ODE, f was defined using a continuous neural network and x(I) is always continuous l. Trajectories modeled by an ODE can thus have limited representation capabilities when discontinuities in the state occur. Examples dynamical systems which are hybrid discrete and continuous systems that choose which dynamics follow based on discrete switch, and are popular neuroscience and finance (see e.g. Friston (2009)). For these systems, a practical extension would be to follow Chen et al. (2020) and modify the gradient update to include discrete changes in the event state which case, infer the causal graph, partial derivatives 0kfj would be defined piece-wise.
switching
ACKNOWLEDGEMENT
This work was supported by the Alan Turing Institute under the EPSRC grant EP/N510129/1, the ONR and the NSF grants number 1462245 and number 1533983.
REFERENCES
Robert Adamek, Stephan Smeekes, and Ines Wilms. Lasso inference for high-dimensional time series. arXiv preprint arXiv:2007.1095: 2020.
Roberto Barrio, M Angeles Martínez, Sergio Serrano, and Daniel Wilezak. When chaos meets hyperchaos: 4d rössler model. Physics Letters A, 379033):2300-2305 2015.
Alexis Bellot and Mihaela Van Der Schaar. Flexible modelling of longitudinal medical data: A bayesian nonparametric approach. ACM Transactions on Computing for Healthcare, 1(1):1-15. 2020.
Alexis Bellot and Mihaela van der Schaar. Policy analysis using synthetic controls continuous-time. arXiv preprini urXiv:2102.0 2021.
José Bento, Morteza Ibrahimi, and Andrea Montanari. Learning networks of stochastic differential equations, arXiv preprint 1011.44 2010.
10

---

Published as a conference paper at ICLR 2022
Annalisa Bracco, Fabrizio Falasca, Athanasios Nenes, Ilias Fountalis, and Constantine Dovrolis. Advancing climate science with knowledge-discoveyy through data mining. upj Climate and Amosphericc Science, 1(1):166, 2018.
Steven L Brunton, Joshua L Proctor, and Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the national academy of sciences, 113155)333223333, 2016.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. In Advances in neural information processing systems, pp. 6571-6583 2018.
Ricky TQ Chen, Brandon Amos, and Maximilian Nickel. Learning neural event functions for ordinary differential equations. arXiv preprint 2XX1:20011.03002 2020.
Shizhe Chen, Ali Shojaie, and Daniela M Witten. Network reconstruction from high-dimensional ordinary differential equations. Journal of the American Statistical Association, 112(520):1697- 1707, 2017.
Yonghong Chen, Govindan Rangarajan, Jianfeng Feng. and Mingzhou Ding. Analyzing multiple nonlinear time series with extended granger causality. Physics Letters A, 324(1):26-33 2004.
Victor Cheroozuukov, Wolfgang Härdle, Chen Huang, and Weining Wang. Lasso-driven inference in time and space. Available SSRN 3/88362, 2019.
Mathew Dafilis, Federico Frascoli, Peter J Cadusch, and David TJ Liley. Four dimensional chaos and intermittency in mesoscopic model of the electroenceehallgram. Chaos: An Interdisciplinary Journal of Nonlinear Science, 23(2):02311 2013.
Bryan C Daniels and Ilya Nemenman. Efficient inference of parsimonious phenomenological models of cellular dynamics using s-systems and alternating regression. PloS one, 10(3)=e01 19821, 2015.
David Danks and Sergey Plis. Learning causal structure from undersampled time series. 2013.
Denver Dash. Restructuring dynamic causal systems in equilibrium. In AISTATS. Citeseer, 2005.
Itai Dattner, Chris AJ Klaassen, et al. Optimal rate of direct estimators in systems of ordinary differential equations lincar in functions of the parameters. Electronic Journal of Statisties, 9(2): 1939-1973, 2015.
Vanessa Didelez. Asymmetric separation for local independence graphs.
arXiv preprint
arXiv:1206.6841, 2012.
Vu Dinh and Lam Si Tung Ho. Consistent feature selection for analytic deep neural networks. arXiv preprint 1=Xi:22000000097, 2020.
Phil Dowe. Wesley salmon's process theory of causality and the conserved quantity theory. Philosophy of science, 59(2):195-216, 1992.
Michael Eichler. Causal inference with multiple time series: principles and problems. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 371(1997): 40110613, 2013.
Michael Eichler and Vanessa Didelez. Causal reasoning in graphical time series models. arXiv preprint arXiv:/206.55
2012.
M Eiswirth, Th-M Kruel, G Ertl, and FW Schneider. Hyperchaos in chemical reaction. Chemical physics letters, 193(4):305-310, 1992.
Jerome Friedman, Trevor Hastie, and Robert Tibshirani. Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 9(3):432 441, 2008.
Karl Friston. Causal modelling and brain connectivity in functional magnetic resonance imaging. PLo8 biol,
2009.
11

---

Published as a conference paper at ICLR 2022
Karl J Friston, Lee Harrison, and Will Penny. Dynamic causal modelling. Neuroimage, 19(4): 1273-1302, 2003.
Stuart Glennan. Rethinking mechanistic explanation. Philosophy of science, (9(S3):S342-S33) 2002.
Stuart Glennan. Mechanisms, causes, and the layered model of the world. Philosophy and Phe- nomenological Research, 81(2):362-38) 2010
2010.
Ary Goldberger and Bruce West. Applications of nonlinear dynamics to clinical cardiology. Annals of the New York Academy of Sciences, 504:195-213, 1987.
Mingming Gong, Kun Zhang, Bernhard Schoelkopf, Dacheng Tao, and Philipp Geiger. Discovering temporal causal relations from subsampled data. In International Conference on Machine Learning, pp. 1898-1906. PMLR, 2015.
Clive WJ Granger. Investigating causal relations by econometric models and cross-spectral methods. Econometrica: journal of the Econometric Society, pp. 124-438, 1969.
Markus Heinonen, Cagatay Yildiz, Henrik Mannerström, Jukka Intosalmi, and Harri Lähdesmäki. Learning unknown ode models with gaussian processes. In International Conference on Machine Learning, pp. 1959-1968. PMLR, 2018.
Mathias Heltberg, Sandeep Krishna, and Mogens Jensen. On chaotic dynamics in transcription factors and the associated effects in differential gene regulation. Nature communicationns 10(1): 1-10, 2019.
Vân Anh Huynh-Thne Thu and Guido Sanguinetti. Combining rreed and dynamical systems for the inference of gene regulatory networks. Bioinformatics, 31 110):1611-166, 2015.
Aapo Hyvärinen, Kun Zhang, Shohei Shimizu, and Patrik O Hoyer. Estimation a structural vector autoregression model using non-gaussianity. Journal of Machine Learning Research, 1 11(5), 2010.
Saurabh Khanna and Vincent YF Tan. Economy statistical recurrent units for inferring nonlinear granger causality. arXiv preprint 1rXiv:191..09879, 2019.
Anders Bredahl Kock and Laurent Callot. Oracle inequalities for high dimensional vector autoregres- sions. Journal of Econometrics, 186(2):325-344, 2015.
Rebecca M Kuiper and Oisín Ryan. Drawing conclusions from cross-lagged relationships: Re- considering the role of the time-interval. Structural Equation Modeling: Multidisciplinary Journal, 25(5):809-82) 809-823,, 2018
Ismael Lemhadri, Feng Ruan, and Rob Tibshirani. Lassonet: Neural networks with feature sparsity. In International Conference on Artificial Intelligence and Statistics, pp. 10-18. PMLR, 2021.
David Lewis. Counterfactual dependence and time's arrow. Noñs, pp. 455 476, 1979.
David K Lewis. Void and object. Department of Philosophy, University of Melbourne, 1998.
Yifeng Li, Chih- Yu Chen, and Wyeth W Wasserman. Deep feature selection: theory and application to identify enhancers and promoters. Journal of Computational Biology, 23(5):322-33) 2016.
Edward N Lorenz. Predictability: A problem partly solved. In Proc. Seminar on predictability, volume 1, 1996.
Peter Machamer. Activities and causation: The metaphysics and epistemology of mechanisms. International studies in the philosophy of science, 18(1):27-39, 2004.
Ričards Marcinkevičs and Julia E Vogt. Interpretable models for granger causality using self- explaining neural networks. arXiv preprint arXiv:2101.07600, 2021.
Th Meyer, MJ Bünner, A Kittel, and J Parisi. Hyperchaos the generalized rössler system. Physical Review E, 56(5):0069, 1997.
12

---

Published as a conference paper at ICLR 2022
Søren Wengel Mogensen, Niels Richard Hansen, et al. Markov equivalence of marginalized local independence graphs. The Annals of Statistics, 48(1):539-555 559, 2020.
Joris M Mooij, Dominik Janzing, and Bernhard Schölkopf. From ordinary differential equations to structural causal models: the deterministic case. arXiv preprint arXiv: 2013.
Yuval Nardi and Alessandro Rinaldo. Autoregressive process modeling via the lasso procedure.
Journal of Multivariate Analysis, 102(3):528-54, 2011.
Roxana Pamfil, Nisara Sriwattanaworachai. Shaan Desai, Philip Pilgerstorfer, Konstantinos Geor- gatzis, Paul Beaumont, and Bryon Aragam. Dynotears: Structure learning from time-series data. In International Conference on Artificial Intelligence and Statistics, pp. 1595-1605. PMLR, 2020.
Neal Parikh and Stephen Boyd. Proximal algorithms. Foundations and Trends in optimization, 1(3): 127-239, 2014.
Judea Pearl. Causality. Cambridge university press, 2009.
Huw Price. Agency and probabilistic causality. The British Journal for the Philosophy of Science, 42 (2):157-176, 1991.
Zhaozhi Qian, Ahmed Alaa, Alexis Bellot, Mihaela Schaar, and Jem Rashbass. Learning dynamic and personalized comorbidity networks from event data using deep diffusion processes. In International Conference on Artificial Intelligence and Statistics, pp. 3295-3305. PMLR, 2020.
Xiaojie Qiu, Qi Mao, Ying Tang, Li Wang, Raghav Chawla, Hannah A Pliner, and Cole Trapnell. Reversed graph embedding resolves complex single-cell trajectories. Nature methods, 14(10):979, 2017.
Federica Raia. Causality in complex dynamic systems: A challenge in earth systems science education. Journal of Geoscience Education, 56(1):81-94 2008.
Maziar Raissi. Deep hidden physics models: Deep learning of nonlinear partial differential equations. The Journal of Machine Learning Research, 19(1):93--955, 2018.
Maziar Raissi, Paris Perdikaris, and George Em Kariiadaiss. Physics informed deep learning (part i): Data-driven solutions of nonlinear partial differential equations. arXiv preprint arXiv: 171.:1671.. 2017.
Jim O Ramsay, Giles Hooker, David Camphell, and Jiguo Cao. Parameter estimation for differential equations: generalized smoothing approach. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 69(5):741-796. 2007.
OE Rossler. An equation for hyperchaos. Physics Letters A, 71(2-3):155-157, 1979.
Paul K Rubenstein, Stephan Bongers, Bernhard Schölkopf, and Joris M Mooij. From deterministic odes to dynamic structural causal models. arXiv preprint arXiv: 1608. 080288 2016.
Donald B Rubin. Causal inference using potential outcomes: Design, modeling, decisions. Journal of the American Statistical Association, 100(469):322 331, 2005.
J Runge, Sejdinovic, and Flaxman. Detecting causal associations in large nonlinear time series datasets. arXiv preprint arXiv: (eX2iv/702.070000 2017.
Jakob Runge. Causal network reconstruction from time series: From theoretical assumptions to practical estimation. Chaos: An Interdisciplinary Journal of Nonlinear Science, 28(7)075310, 2018.
Jakob Runge, Sebastian Bathiany, Erik Bollt, Gustau Camps- Calps-Valls, Dim Coumou, Ethan Deyle, Clark Glymour, Marlene Kretschmer, Miguel D Mahecha, Jordi Muñoz-Mari, et al. Inferring causation from time series in earth system sciences. Nature communications, 10(1):1-13, 2019.
Wesley Salmon. Scientific explanation and the causal structure of the world. Princeton University Press, 1984.
13

---

Published as a conference paper at ICLR 2022
Bernhard Schölkopf. Causality for machine learning. arXiv preprint arXiv:1911.10500 2019.
Thomas Schreiber. Measuring information transfer. Physical review letters, 85(2):461, 2000.
Peter Schulam and Suchi Saria. Reliable decision support using counterfactual models. In Advances in Neural Information Processing Systems, 1697-1708, 2017.
Shai Shalev-Shwartz and Shai Ben-David. Understanding machine learning: From theory to algorithms. Cambridge university press, 2014.
Christopher A Sims. Macroeconomics and reality. Econometrica: journal of the Econometric Society, pp. 1-48, 1980.
Hossein Soleimani, Adarsh Subbaswamy, and Suchi Saria. Treatment-response models for coun- terfactual reasoning with continuousttmme, continuous-valued interventions. arXiv preprint 02X00 2017.
Peter Spirtes, Clark N Glymour, Richard Scheines, and David Heckerman. Causation, prediction, and search. MIT press, 2000.
Alex Tank, Ian Covert, Nicholas Foti, Ali Shojaic, and Emily Fox. Neural granger causality for nonlinear time series. arXiv preprint arXiv: 1802.05842 2018.
Cole Trapnell, Davide Cacchiarelli, Jonna Grimsby, Prapti Pokharel, Shuqiang Li, Michael Morse, Niall J Lennon, Kenneth Livak, Tarjei S Mikkelsen, and John L Rinn. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nature biotechnology. 32(4):381, 2014.
Sara Van de Geer, Peter Bühlmann, and Shuheng Zhou. The adaptive and the thresholded lasso for potentially misspecified models (and a lower bound for the lasso). Electronic Journal of Statistics, 5:688-7455 -298,,749100 2011.
James M Varah. A spline least squares method for numerical parameter estimation in differential equations. SIAM Journal Scientific and Statistical Computing, 3(1):28-46, 1982.
Philippe Wenk, Gabriele Abbati, Michael Osborne, Bernhard Schölkopf, Andreas Krause, and Stefan Bauer. Odin: Ode-informed regression for parameter and state inference in time-continuous dynamical systems. Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pp. 6364-6371, 2020.
Jon Williamson. Mechanistic theories of causality part i. Philosophy Compass, 6(6) 2011.
James Woodward. A theory of singular causal explanation. Erkenntnis, 21(3):231-266 1984.
Hulin Wu, Tao Lu, Hongqi Xue, and Hua Liang. Sparse additive ordinary differential equations for dynamic gene regulatory network modeling. Journal of the American Statistical Association, 109 (506):700-716, 2014.
Zuogon Yue, Johan Thunberg, and Jorge Goncalves. Inverse problems for matrix exponential in system identification: system aliasing. arXiv preprint arXiv:1605.06973, 2016.
Anru R Zhang and Yuchen Zhou. On the non-asymptotic and sharp lower tail bounds of random variables. Stat, 91))=331, 2020.
Huaqing Zhang, Jian Wang, Zhanquan Sun, Jacek M Zurada, and Nikhil R Pal. Feature selection for neural networks using group lasso regularization. IEEE Transactions on Knowledge and Data Engineering. 32(4):659 -63, 2019.
Peng Zhao and Bin Yu. On model selection consistency of lasso. Journal of Machine learning research, 7(Nov):5541-2566 2006.
Hui Zou. The adaptive lasso and its oracle properties. Journal of the American statistical association, 141877) 2006.
14

---

Published as a conference paper at ICLR 2022
Appendix
This appendix is outlined as follows:
Section A discusses additional related work.
Section A. contrasts discrete-time and continuous-time structural models considering in- consistencies in structure learning by modelling the underlying continuous-tima process in discrete-time.
Section A.2 discusses other work modelling differential equations.
Section A.3 comments the philosophical debate around the nature of causality in dynamical systems and its ramifications in machine learning.
Section proves lemmas discussed in the main body of this paper.
Section C includes additional experiments and all experimental details.
Section C.1 gives performance results function of feature dimensionality, neural network parameterization and gives run time comparisons.
Section C.2 gives performance results in terms of the area under the ROC curve.
Section C.3 presents purely synthetic data generating mechanism and tests different regular- ization schemes to justify the adaptive group lasso empirically.
Section C,4 includes algorithm implementation details.
Section C.5 gives the data generating system for the Glycolysis experiment.
15

---

Published as a conference paper at ICLR 2022
A RELATED WORK
A.1 GRAPHICAL MODELLING IN DISCRETE-TIME
Graphical modelling with time series has been driven by applications in causality. A first practical definition of causality for inference with time series data was given by Granger (Granger, 1969). A time series 2; is said to Granger-cause Ij if omitting the past of a time series I; in time series model including sj 's own and other covariates' past increases the prediction error of the next time step of z.. Implementations of this principle are variants of vector auto-regressive (VAR) models and its extensions (see e.g. (Sims, 1980; Chen et al., 2004; Tank al., 2018; Pamfil al., 2020)), for example in the linear case assuming,
x(1 + AC) Bix(() +... Bax(1 - kat) + e())
~1
(12)
for some distribution P. Each one of the matrices BI...., B2 then describe lagged causal relation- ships with different lags. Subsampling occurs when the frequency of observation is lower than At and renders VAR models generally unidentifiable although specific exceptions exist and have been explored in (Gong et al., 2015; Danks Plis, 2013). Alternatives to VAR models include the PC algorithm with conditional independence testing methods accounting for auto-corredd correlallst between successive observations as done by Runge (2018); Runge et al. (2017; 2019) and transfer entropy principles (Schreiber, 2000).
We illustrate next with an example why, from a learning perspective, graphical modelling in discrete- time cannot be consistently applied for the purpose of graphical modelling in dynamical systems without strong assumptions on the observation process or parameterization of the underlying structural model.
A.1.1 EXAMPLE
Assume the true dynamics of two processes x(t) = (z1(t), to be given by,
d x(t) Ax(t), dt
- 1 2 0.5
A =
(13)
The underlying dynamics are given by A: causing an increase in the rate of change of æ1 (with value 2) while 12 having a large negative causal effect on the rate of change of αj (with value -4).
The corresponding discrete-time model (with one time lag for simplicity) may be defined as,
"Causal" strength
X; x2
X1 X22
x(t + At) Ba: x(t).
(14)
where At is the time interval of observation. The two models are connected by a simple relation, given by Bat exp(A- A/), that can be used to uniquely compute the off-diagonal entries [Bar]12 and [Bai]21 indicating the strength of "causal" ef- fects x2 + 1 and æ + 22 respectively, under the Granger-causality paradigm (Granger, 1969).
20 33
25
45
Time interval
Figure Discrete-time model inferred causal in- teraction with underlying continuous- process.
As first observation, irrespective of the time interval At note that single discrete- eiscrete-time model Ba may describe multiple underlying mechanisms A as matrix logarithms log BA are not identified Bar has complex eigenvalues (see e.g. (Yue et al., 2016) for a formal result). Second, for a single model defined by A the discrete-time causal interpretation (i.e. Bat) may change dramatically as a function of the time interval At which the process is observed, as shown in Figure 4. Describing causality within dynamical systems in discrete-tima is inherently an ill-posed problem. In fact, only for dynamical systems that are bivariate, stable, and non-oscillating can expect consistent conclusions at all measurement intervals.
Proposition 1 (Causal inconsistency in discrete time (Kuiper & Ryan, 2018)). The sign of off- diagonal entries of B and A agree for all time intervals At > 0 if defines a bivariate, stable, and noc-lscilgatiig system of differential equations.
16

---

Published as a conference paper at ICLR 2022
Proof. (Given here for completeness can also be found in (Kuiper & Ryan, 2018)) Using the relationship Bai = exp(A- Ar} we can write,
Bar V-1 exp(DA-ac)V
(15)
where A. represents the diagonal matrix containing the scalar exponentials of the eigenvalues of multiplied by the scalar At and V the matrix of eigenvectors. This can be used to relate two estimated discrete-timg matrix coefficients with each other,
Ba
B(Agt/A11) A,t
(16)
Consider bivariate process with eigenvalues of an estimated Ba, denoted by Å and X2, and let Aet = n At. Then, misleading causal conclusions are obtained if the sign of the entries of Ba, differ from Batt. In this case they may be computed explicitly. holds that the sign of off-diagonal entries are equal if and only if (Λ1 X2) (determining the sign of entries of Ba,)) and (Xī - Xς) (determining the sign of entries of Ba,t) are equal. This is the case if both A and X2 lie between and i..., stable, bi-variate, non-oscillating system, but will not hold in general otherwise.
Counter-examples of discrepancies can be found for many other types of differential equations in (Kuiper & Ryan, 2018).
A.2 RELATED WORK ON MODELLING DIFFERENTIAL EQUATIONS
Two-stage collocation methods were first proposed by Varah (1982). The authors proposed to fit smoothing estimate &(: h, () to the observations 1/:....: Un with smoothing parameter h, and using it and its derivative with respect to / in order to estimate the vector field f in system of differential equations,
min 1 n
[dx(t;;)) h) (x( h))
(17)
where ((;;h) = arg min Σ lly: x(1:; h)ll is the smooth interpolation function and most x(=)))E
approaches assume f to be linear although variations of this principle have been developed using non-linear functions of observations (e.g. see (Dattner et al., 2015)) and using bases of functions without explicitly determining their form (e.g. see (Chen et al., 2017)). There are two important differences between the proposed approach and two-stage collocation methods.
The fact that one must choose the interpolation and smoothing function leads to a very different analysis of the properties of estimators. Estimates rely heavily on the smoothing estimates obtained consistency has only been shown for certain values of the smoothing parameter h that are hard to choose in practice (Chen et al., 2017).
The objectives of two-stage collocation methods is the statistical estimation of the parameters 0 of fo rather than of the parameters of the graphical structure induced by fg. As a contrast, in the particular case of analytic neural networks the exact set of parameters in the data generating mechanism is not uniquely identified multiple alternatives define the same input-output map of the vector field fg. even in the infinite-sample regime.
Broader literature. Modelling time series is a wide and varied research topic. One can train feed-forward recurrent neural networks approximate differential equation (Chen al., 2018; Raissi, 2018) and model interventions (Bellot & van der Schaar, 2021). Gaussian Processes have been adapted to fit differential equations (Raissi et al., 2017) and have been proposed to model continuous-time interventions (Schulam & Saria, 2017; Soleimani et al., 2017). And tree-based methods are also popular in biology to model irregularly sampled data (Huynh- Thu & Sanguinetti, 2015; Bellot & Schaar, 2020). The objective in these papers however is to extrapolate the latent state of the process i.e., the forward problem: inferring the latent state of the trajectory from time series. This description has not yet involved structure learning which, in contrast, is concerned with the inverse problem i.c., using discrete measurements of trajectories to infer the underlying structural model.
17

---

Published as a conference paper at ICLR 2022
A.3 PHILOSOPHICAL DEBATE ABOUT THE NATURE OF CAUSALITY
A mechanistic interpretation of causality views causal claims as claims about the existence of a mechanism or process that mediates events one time with events later time and has been formalized mathematically as system of structural differential equations. Proponents of this view hold that the metaphysical connection between mechanisms and causality very close: two events are causally connected if and only if they are connected by an underlying physical mechanism. Glennan (Glennan, 2010; 2002) and Machamer (Machamer, 2004) define mechanism as complex system of interacting parts, but process accounts have also been proposed (Salmon, 1984; Dowe, 1992) that focus on the fact that causal processes manifest conserved quantity. the latter, an interaction between two processes is causal if there is exchange of a conserved quantity between them. Mechanistic theories causality, from a philosophical standpoint, are normally contrasted with probabiliittic (Pearl, 2009; Spirtes et al., 2000), counterfactual (Rubin, 2005; Lewis, 1979; 1998) and manipulationist (Pearl, 2009; Price, 1991; Woodward, 1984) theories of causality. These define a connection be causal if and only if change in makes difference to the other. This distinction a useful conceptual contrast even though it has been noted that these accounts overlap in some measure: for instance, mechanisms can be given a counterfactual analysis and thus would be form of difference-making theory (Williamson, 2011).
Differential equations and their causal semantics (e.g. formalizing the meaning of interventions or counterfactuals) are often derived from those of (static or time- (imeinndpenddent Structural Causal Models (SCMs, Definition 7.1.1 (Pearl, 2009)) that, in dynamical systems, have been interpreted as defining an intermediate layer of expressiveness between the underlying model of differential equations and statistical models based on associations (Schölkopf, 2019). The connection between these conceptual layers is of ongoing interest. For instance, SCMs have been shown to describe changes in equilibrium states (if they arise) (Dash, 2005; Mooij et al., 2013), and have been shown to describe changes in asymptotic dynamics (Rubenstein et al., 2016) under carefully defined interventions, and has been the focus of most work at the intersection of dynamical systems and causality.
In any case, the underlying nature and semantics of causality does not influence the theory and algorithms presented this paper as long as we assume that the causal system of interest may be represented as Neural Dynamic Structural model. In such model, causation across time in dynamical systems due to a derivative (e.g., velocity dx) causing a change in its integral (e.g.. position x). All other causation is contemporaneous, occurring between two variables on a time-scale that is smaller than the time-step of observation and defined by the sparsity pattern of f.
18

---

Published as a conference paper at ICLR 2022
B PROOFS
We will use RR, Rn,(fo) (fo) and (0) interchangeably. With the notation introduced in the main body of this paper, will use the following Lemma that extends the standard Taylor's inequality around local optimum to unbounded sets 0"
Lemma (Lemma 3.2. (Dinh & Ho, 2020)) There exist c22 v 0 and such that [R(0) RR(00) 2 (2d(0, @r)" for all E O.
d(0, θ*) is the minimum Fuclidian norm between and any element of O
B.
PROOF OF LEMMA
The uniqueness of local independence graphs in systems of fully observed stochastic processes was given Proposition 3.6 by Mogensen al. (2020). We restate the Lemma for convenience.
Lemma (Uniqueness of local independence graphs, Proposition 3.6 (Mogensen et al., 2020)). In the context of Neural Dynamic Structural models, two processes are locally dependent given any subset of other processes if and only if I. appears in the differential equation of æ.. i.e. [10x fill,lll,, + O. Moreover, for any such that [|0x.[]] L2 = there exists an equivalent vector field such that the euclidian norm of its column vectors [[[Aq]=*112 = O.
Proof. Consider the class of NDSMs with F such that fy is independent of process Xg, that is [50&, 0, and the class of NDSMs Fo with fj parameterized such that [[[4]]=*[12 = 0. We will how that J Fo-
It is clear that the class of NDSMs with Fo is contained in the class of NDSMs with F as process X interacts with the vector field /3 of X; only if the corresponding entries the first layer of the analytic neural network are non-zero: Fo C.F.
For the converse consider NSDM whose vector field f, of x; is independent of the k-th process x; and consider a set of processes X and such that x except for the k-th entry of & which is set to the zero function ãR [0,T] T' - 0. Because of independence, and because the two processes differ only in their k-th entry: fy (x))) (x))..)
Now define À₁ be the matrix such that [Āi]æ 0 and [Ai] [Ai]ak for all k'' * k. Then it holds that À₃ x1)= = A1)(1). And therefore also (x(((1;; A1) (x(((, A1) (x(1;;; By definition of Fo. (: A1) E F and thus FC Fo- FCF.
B.2 PROOF LEMMA 2
We restate the Lemma for convenience.
Lemma (Generalization bound). Assume En, to he invertible and let a = (aa,..., a.) such that 01 > > On > O be its eigenvalues. For any 8 O, [|ala 3, there exists C such that,
||0|12
(Ro R(fo)] S
Clog (|al|2)
(18)
with probability at least
Proof. Without loss of generality, we consider each differential equation separately. nRR,, (fa) =:
X is sum of dependent squared normal random variables. With this notation X (Xi...., X,) E R where X; = (Y,j 4 (1+) E R is a univariate random variable and Yı- E R is the j-th random variable at time 4 given of the observation model and underlying dynamical system. X (X1,..., Xx.) has joint distribution defined by a mean μ and covariance matrix Σ,, as defined the main body of this paper. We may write Z = Σx 1/2 (X- - (1) and,
x (Z + - 1/2 E)
(19)
Let Σ,, = AV be the eigendecomposition of Σ, where V is an orthogonal basis of eigenvectors and A is diagonal matrix of positive eigenvalues (o1,..., an). U VZ then is also
19

---

Published as a conference paper at ICLR 2022
multivariate normal, with expectation and identity covariance matrix (since VTV VVT In). For VEE 1/2 rewriting the decomposition above in terms of U and
x (U + u)² A(U u) - Σ a: au u:)
(20)
The distribution of (fo (fo) = X" X is thus weighted non-central x² random variable.
Write W VV RRx= By applying the concentration results in Theorem 6 and 7 in (Zhang & Zhou, 2020) we have,
=2n
P RR (0) - [ 8) S exp
||a|13 + >/= (E -
2
(21)
2
-02 n18
S exp
(22)
||a||2
for 8 such that,
|al|
2 2Σ
プラ=
Wo)l((((():0) - (((())))))))))(())))))
O<S<
(23)
n2|a 1
n
C1 C2 are two scalars not depending on a. The remaining of the proof follows the argument of (Dinh & Ho, 2020). We define events,
A(O, s) = (0) R(0) > s},
(24)
B(0, 8) = (30' E θ such that 110 - 0 S
4M, and RR (0) R((0) 8}
5
(25)
- (I) -
R., (., S Ma 12,,,, v0, 0' E e),
(26)
where the last event is defined with respect to the Lipschitz constant Ms of R., (assumed to be Lipschitz with probability at least - ó, that is IP(C) 2 1 - 8). Let m dim(O), there exist (y(m) 2 and finite set H θ such that,
ec U v(O,e),
|7| Ca/em,
(27)
where we choose = s/(MMs). V(O,e) denotes the open ball centered at with radius and [H|] denotes the cardinality of H. In other words, e-covers θ and the incquality involving the cardinality of H. follows because θ is bounded subset of Euclidian space, see e.g. section 27.1 in (Shalev-Shwartz & Ben-David, 2014). By a union bound over all elements H.
P ( E :: (R) R(0) > 8) S Ca(4M)) exp -C2
n
2
22
(28)
||00112
Since B(0,s) s) Cc C A(0, s) and 76 C we have,
2
P(30 E θ R (0) R(0) > s) < Ca(4MM) exp -C2
n
2
+
8.
(29)
[|a||2
/C log(k)
Now let k | and let s
for notational simplicity. To complete the proof we need to
choose C such that,
m
4M5k
C
exp C- C. log(k)) 8.
(30)
(loo((k)
20

---

Published as a conference paper at ICLR 2022
This inequality holds if,
Ca(4Ma) km C2 1m-Cx-C S
(31)
since (C)) 2 1, which can be obtained if m - C2 0 and C2(4M)) 3m-C2-0 v
since & by assumption, so that,
C 2 C2 log(3) (log(C's) + m (4M)) + log (//)))
(32)
B.3 PROOF LEMMA
To traverse this result, will start by considering the convergence of the group lasso.
Lemma 7 (Convergence of Group Lasso). For any > 0, assuming that AGL + there exists
> 0, Cg > Ns 0 and Ts > O such that,
min EEB*
||0, - 011 S Cs
cL
1/vv-1
||0112 n
V
log
||ox|12
n2
(33)
with probability at least 1- 5.
Proof. Recall the group lasso and adaptive group lasso penalty terms,
PGr (0) AGE. Σ
||A])a|l
and
PAGI. (0) = AAGE
k,j=1 [[[A]] [[A]]
k,==1
By definition, we have,
R., + p0a) S Rr + (0));
(34)
where 0,, arg min
(0) + (0) the parameter solution to the group lasso.
age
It holds then that,
min l0e*
011 s R((a) R(6o) C2 00,
(35)
R((m)) R, (dw)) (R((0) R.,(oo)l + (Q,) R. (0o)
(36)
||a||2
2
Clog
m ||all2
+ (PG) - (0n))
(37)
n
||ll2
<
C log
+ AGI. K 1100 - 0m.12,
(38)
||a||a
where the first inequality is due Lemma 6 (for some 02, 1/ 0), the second inequality is due the triangle inequality, the third inequality is due to equation (34) and Lemma 1, and the fourth inequality comes from the Lipschitzness of PGL. We have used K The last step is given by Youngs's 's inequality, e.g. as stated in section 5.1 (Dinh & Ho, 2020), to conclude that,
denote the Lipschitz constant of PoL-
min O0
1|0, - 011 Cs
[Iaell2
V |log
[|ol12
+ >
(39)
for some constant C3.
We will now state and prove Lemma 3.
Lemma (Convergence of Adaptive Group Lasso). For any 8 > 0, assuming that XAGL. 0 there exists v > 0, Cg > 0. N₂ > and T3 > 0 such that,
10, 01 C's Pscc + (1a|122
min 8EB.
log ||a/12/
(40)
21

---

Published as a conference paper at ICLR 2022
with probability at least 1 - 5.
Proof. By the convergence of the group lasso [[[A]] is bounded away from zero for any process that causally significant for process j. k, E {L...., d}. Let the set of causally significant pairs (k, j) be denoted S. Then can define,
M(O) Σ (RJ)ES
[|l < oo.
(41)
[[<9]
since ||[A]]=11 > O.
Let On = arg min (0) + PAGL(O) be the solution parameters of the adaptive group lasso problem. deG
By a similar derivation to that used (35),
|||||2
R(Ōx)) R(00)) 2C
V log
+ AAGE. (M(0.)-M(00))
(42)
||c||12
|||||
S 2C
log
+ MOGL)M(On).
(43)
||cx||2
And since M(0) is finite positive scalar, again by a similar derivation to that used in (35),
1/v
(11call2 n
min SE8
||0, - 0112 S Cs
V log
+ dAca
(44)
||a|12
for some constant Cs > 0.
B.4 PROOF LEMMA
We restate the Lemma for convenience.
Lemma 4 (Local consistency of Adaptive Group Lasso). Let 7 > > 0, AaGL.
2(( || la )->/u=+), and AAGL (((7T then the adaptive group lasso locally consistent.
Proof. By the convergence of the Group Lasso, for any pair (j. k) of non-significant processes,
1 + ||a||2
[[[A]]=]l2 S Cs
log ||al|2
(45)
with probability at least 8-8. 8.
holds therefore that,
lim 1
2 xc.
(46)
Tala Do
Now assume for contradiction that there exists pair (k,j) of non-locally significant processes (that is, Fg is not locally significant for æ.y) such that [[[A5][*11 # (the notation above the matrix A denotes estimation with the adaptive group lasso) and define ((Ōm) to be equal to O, except that non-significant parameters are set to zero. By the definition of 0,, as minimizing the empirical risk regularized by the adaptive group lasso constraint,
R2 (O,) + AAGL
||A|] Ro,O)(
(47)
[[[Aa]
Since the regularization term on the right-hand side is zero by the definition of o(0). Then by Assumption 3,
AaGe ||LAA]=.11
|1A//1l12 I (n,()(m)) RR,(
(48)
< Msllo((m)) - 0
(49)
(50)
22

---

Published as a conference paper at ICLR 2022
But since have assumed |[[A]]+112 it follows from above that AAGL [| S Ms which is a contradiction of Lemma that proved the convergence of the group lasso and in particular that 3G for non-locally significant pairs of processes (k.
lim T-
B.5 PROOF OF LEMMA 5
We restate the Lemma for convenience.
Lemma 5 (Finite sample local consistency of Adaptive Group Lasso). Under the conditions of Lemma with the additional minimum restricted strength assumption in (10) the problem design for particular values of n and a, the Adaptive Group Lasso recovers the structure G exactly with high probability.
Proof. First, note that for the set of loss minimizers O defined in eq. (6) and by using the fact that neural networks are analytic, does hold that for any two locally dependent processes 2y and Ty the first layer parameters of any model with minimum loss are bounded away from zero, i.e.
[[[A]] x|12 for some 0.
To see this, assume for a contradiction that no such e exists, and therefore that there exists [A]]+ O such that ||[A]]-+ = since neural networks are analytic and each one of the parameters is defined in bounded intervals. This would imply that there exists neural network with the same input-output relationship as the true model fj that does not depend on its k-th input, which a contradiction because ||0=ff |l, O.
Next, given the minimum strength condition on the column norms of first layer parameters related to locally dependent processes,
|A Iiin > Os (AaGL +
||a| n
n
log
(51)
[|a+ ll
where recall that [1|miin min [[[A{] j.k I,...,d,
||0 * 0}, have that by Lemma 3 that estimated parameters [|Al] are bounded away from zero for specific values of and n since,
1/v
|| [Al1-ol12 Co (AxCL + (1lall2) n
n 11a1112
min 80
log
(52)
with high probability.
23

---

Published as a conference paper at ICLR 2022
C EXPERIMENTAL DETAILS
C.I RESULTS AS FUNCTION OF FEATURE DIMENSIONALITY, NEURAL NETWORK PARAMETERIZATION AND RUN TIME COMPARISONS
Performance with increasing number of variables is monitored to some extent with the Rössler experiment. We extent this analysis in this section to include performance comparisons as function of more variables. We also include run-time comparisons and performance comparisons function of different model parameterizations to understand the practical use of the proposed approach. We limit our comparisons here to NGC-LSTM, DCM and NGM.
Rössler = 10)
Rössler (d = 50)
Rössler (d 100)
Rössler = 100)
TPR f
FDR 4
TPR
FDR 4
TPR
FDR L
TPR
FDR
NGC-LSTM
.49 (.04)
.53 (.04)
.38 (.04)
.64 (.08)
DCM
.87 (.01)
.10 (.04)
.97 (.01)
.31 (.07)
.94 (.04)
.35 (.05)
.90 (.04)
.40 (.05)
NGM (ours)
.96 (.01)
.02 (.01)
.95 (.01)
.04 ((01)
.95 (.04)
.05 (.02)
.89 (.04)
.05 (.02)
Table 2: Performance comparisons on Rössler's model.
Rössler (d = 10)
Rössler (d 50)
Rössler (d 1000
Rössler (d 200)
NGC-LSTM
2523
5321
DCM
12
74
320
991
NGM (ours)
291
480
923
2289
Table Learning time seconds.
Rössler (d = 10)
Rössler (d = 50)
TPR t
FDR J
TPR t
FDR +
(1,10) Default
.96 (.01)
.02 (.01)
(.01)
.04 (.01)
(1,50)
.95 (.01)
.02 (.01)
.95 (.01)
.04 (.01)
(2,10)
.96 (.01)
.02 (.01)
.96 (.01)
.03 (.01)
(2,50)
.92 (.01)
(.01)
(.01)
(.01)
(5,10)
.95 (.01)
.03 (.01)
.95 (.01)
.05 (.01)
(5,50)
.94 (:01)
.04 (.01)
.96 (.01)
.06 (.01)
Table 4: Performance comparisons with different neural network architectures. Notation: (number of hidden layers, number hidden units each hidden layer).
Note that the parameter of interest is a norm over the columns of parameter matrices in the first NN layer which we found be sufficiently coarse not be sensitive small variations in architecture choices. Our default parameterization worked well across all datasets we analyzed. In Table 4, we show additional experiments that cover possible configurations that span reasonable choices that practitioner may make. We find that performance results are largely invariant to differences in the number layers and layer size.
C.2 RESULTS USING THE AREA UNDER THE ROC CURVE
In this section report all experiments in the main body of this paper using the area under the ROC curve (AUC) as performance metric. Figure 5 contains comparisons for experiments using the Lorenz model and Table contains comparisons for the Rössler and Glycolytic experiments.
24

---

Published as a conference paper at ICLR 2022
10
DCM
0.9
NGC-LSTM
PCMCI
AUC
AUC
0.8
NGM
0.7
0.6
0.5
0.6
0.01
0.05
0.1
0.15
0.25
Interval observation
% observations removed
Figure Experiments on Lorenz's model.
Rössler (d = 10)
Rössler (d = 50)
Glycolytic
NGC-MLP
.73 (.02)
.70 (.02)
.57 (.04)
NGC-LSTM
.75 (.03)
.74 (.02)
.65 (.04)
SVAM
.56 (.03)
.50 (.05)
.60 (.03)
PCMCI
.51 (.04)
.50 (.05)
.62 (.03)
DCM
.95 (.01)
.90 (.02)
.70 (.02)
NGM (ours)
.99 (.01)
.98 (.01)
.78 (.03)
Table 5: AUC on Rössler and Glycolytic data. Numbers in parenthesis standard deviations.
C.3 PURELY SYNTHETIC EXPERIMENT
This experiment investigates the behaviour of different regularization schemes with purely synthetic data generating mechanism. We generate 1000 observations of 50 variables X E 121000x50 using regular evaluations of the process defined by,
d dt =0 (,(x(() + dw,((), x(0) Xo:
(53)
for = I,..., 50, where each vector field component 9, 50 E is parameterized by neural network with three hidden layers of 10 nodes, such that depends only on 5 (locally significant) of its 50 arguments. The initial state Xp as well as all weights and biases are independently drawn from standard Gaussian random variables before setting the locally not significant first layer columns to zero. The time interval between observations is fixed 0.1 units.
The results are given in Table 6.
We write NGMØ1. for NGM with group lasso regularization, NGM AGI. for NGM with adaptive group lasso regularization and NGM,, for NGM with conventional lasso regularization. True positive rates are comparable across regularization schemes but false discovery rates are much lower for adaptive regularization which suggest that standard lasso and group lasso algorithms may not be aggressive enough to enforee sparsity strictly. As consequence, we may need cut-off values to interpret NGM&L locally which are difficult to specify in practice, while NGMAGL sets most non-local processes exactly without further post-processing
TPR
FDR +
NGM
.63 (.04)
.30 (.08)
NGMgL
.71 (.05)
.25 (.08)
NGMACL
.70 (.04)
.17 (.05)
Table 6: Regularization choices.
25

---

Published as a conference paper at ICLR 2022
C.4 ALGORITHM IMPLEMENTATION
C.4.1 NEURAL GRAPHICAL MODELLING (NGM)
Proximal gradient descent. The proximal step for the group lasso penalty is given by group soft-thresolding operation on the input weights.
In each iteration, proximal gradiient steps make two update computations to the relevant parameters, denoted here {.. = .....,d), d). They are.
1. 0 t - vV)
2. 0 + arggin 0112 adG., Ek.z |||4)]=]]
where the second part is the proximal operator with respect to the constraint, which for each 0 E [A]]a equivalent soft-threshold group-wise update,
[A]]. |[[4]]-011
[ & <-
max [0,]][]]] - a
(54)
Regularizing constants are chosen from the set (0.001, 0.01, 0.05, 0.1, 1, with y = using average test errors from random train-test splits of the corresponding dataset.
Threshold selection. With this optimization procedure we do not need to set a threshold for convert-
ing the weights to the presence / absence of edges in the graph. A non-zero estimate of [[[AA]=x]] is considered as presence of an edge in the underlying graph and zero estimate of [|[A]] is considered as absence of an edge.
Architecture. The integrand fo was taken to be feed-forward neural network as described with single hidden layer of size 10 and elu activation functions after each layer except after the output layer. In each case used the Adam optimiser implemented by PyTorch. Starting learning rates varied between experiments (with values between 0.001 and 0.01) before being reduced by half metrics failed to improve for certain number of epochs. It was enough in all experiments to consider the final parameter configuration (instead of the one with best validation performance) as only norms of first layer parameters of interest which we found not to be sensitive to the exact epoch choice. The same architecture was used for all experiments.
C.4.2 NEURAL GRANGER CAUSALITY
We implement Neural Granger Causality (Tank et al., 2018) with the code provided by the authors at https: //t:: com/ wewoovveerrNeuaaaaeee
We considered two architectures: an MLP to fit the lagged time series explicitly to its next value in time, and a LSTM model the hidden state capturing the relevant history information. We follow the author's 's implementation and use single hidden layer with 5 nodes. 5 lagged variables, relu activation function and hierarchical penalty optimized with their GISTA training procedure.
Threshold selection. NGC similarly uses an adaptive procedure to optimize parameters and presence / absencee edges in the underlying graph can simply be read off the non-zero parameters.
C.4.3 DYNAMIC CAUSAL MODELLING
Dynamic causal modelling attempts to recover the vector field f explicitly by fitting a multivariate linear model to map the set of variables X(t) Rd with estimated derivatives dX((t) E Rd (in our case computed separately with smoothing spline approximations).
Architecture. In our implementation, hidden layer of size 10 and elu activation functions after each layer except after the output layer. In each case used the Adam optimiser implemented by Porch..
parameterize f with neural networks. We use single
The variant of Dynamic Causal Modelling we used, with neural networks to approximate the mapping between variable states and their derívatives, is own implementation. The architecture is similar
26

---

Published as a conference paper at ICLR 2022
to that used in NGM with the exception that derivatives are approximated priori with the derivatives of natural cubic splines taken to interpolate the observed data. The smoothing hyperparameter in the spline fit was chosen for visual inspection preserve the trajectory of the curves in each data separately.
The optimization problem thus consisted in fitting the vector field explicitly such fit the approximated derivatives dx (t) evaluated at a given time t as well as possible. Derivatives are computed from cubic spline interpolation of the time series. We manually interpolation hyperparameters to obtain a visually faithful approximation of the observed trajectory for each dataset.
The optimization objective is given as,
T2
dx (1) (x(((),)
arg min EF
+ PPGA
(55)
n
dt
1=1
X here denotes the interpolation over time of the time series. We use similar regularization arguments for consistency with NGM and also proximal gradient descent for optimization.
Threshold selection. For DCM, in our implementation, we use the same proximal gradient descent method with the adaptive group lasso constraint and thus we do not require a threshold determine the presence / absence of edges the underlying graph. Presence absence of edges is defined by non-zero parameter values.
C.4.4 PCMCI
PCMCI (Runge et al., 2017; Runge, 2018; Runge et al., 2019) is a discrete-ting two-step approach that uses version of the PC-algorithm with the momentary conditional independence test account for autocorrelation the time series.
PCMCI implemented with the python package provided by the authors github. io/tigramite//
https:
Threshold selection. We chose to adjust for multiple testing with Benjamini-oocheer''s procedure and considered associations significant, determining presence / absence of edges, with p-values below 0.00001 (chosen here because gave good trade-off between TPR and FDR i.e., the values with maximum Fi score).
C.4.5 SVAM
SVAM (Hyvärinen et al., 2010) was implemented with the implementation provided by the authors at ht https::/ /tttts::/giii //gtthub...//////// com/cdt151 cd/1//ing ingam with the BIC model selection criterion and sin- gle lagged variable. Including more lagged variables did not alter our results much but adds an additional choice as to how to define causality we would have multiple estimated matrices of inter-relationships.
Threshold selection. We chose the threshold for converting the weights to the presence / absence of edges in the graph based on F1 scores on validation data and was consistent (around 0.1) across datasets.
C.5 DETAILS ON GLYCOLYSIS EXPERIMENT
Pharmacology is a branch medicine that makes extensive use dynamical models to determine the interaction patterns of drugs in the body. Like many other systems, the nonlinearity of dynamics in biology makes it hard to infer drug interactions from experimental data. Simple linear models are computationally efficient, but cannot incorporate these important nonlincarities. The glycolytic oscillator model a standard benchmark for these kinds of systems. It simulates the cycles of the metabolic pathway that breaks down glucose in cells. We simulate the system presented in Daniels and Nemenman (Daniels & Nemenman, 2015) in their equation 19, mimicking glycolytic oscillations in yeast cells.
27

---

Published as a conference paper at ICLR 2022
10
15
-10
1000
(a) Yeast glycolysis model.
(b) Rössler model with 10 variables.
Figure Sample trajectories.
The dynamics, defined by biochemical components denoted 21,22:---- .......,,0000 27, following system of equations (Daniels & Nemenman, 2015),
given by the
d
100 (t) =t(()
(t) = 2.5 1 (zo(t)/0.52)4
0.01dw1 (t)
dt
d
100 x (t) zt(
(t)
6x2(() (1 zt)))) 12x=(t)xs 0.01dw
dt
1 (z))5)))5))
d
dt (t) (t) (1 ==((t)) 16 2z()) (4 =t(()) 0.01dwyyt)
d
di (t) = 16 (t) (4 ze(()) 100 z1(() =t(() - 13 (za xt)( + 0.01dwa(t)
d
dt (() = (1) (1 -=s(t)) 100 (1) (1( (1) (1) 0.01
dt (t)()) - 100 21 (t) (t(
d
1 (ze (1))0.52) + 32 =z(() (4 =o(()) - 1.28 =o(+) 0.01dw;(t)
d
dt T7(() = 1.3 (t 2 (t)) 1.8 x7 (t) + 0.01dwy(t)
A sample of the trajectories is given in Figure 6, where after the first few time units, the system settles down onto simple limit-cycle behavior.
As in previous examples, the system is observed over a sequence of T time points with a 0.1 time unit interval after randomly initializing each variable in the ranges provided in Table I of (Daniels & Nemenman, 2015). The data is stacked into two matrices for X E RTx7 (and dX E TTT for methods using approximated derivatives) where each row of X a snapshot of the state in time.
28


Reviews:
[{"review_id": "vabFnzUPB4b2", "reviewer": "Reviewer_TPnG", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Strengths:\n\n-- This is a very well written paper that introduces a new framework for graphical modelling. Reading the paper was a real pleasure. Because I was asked to review the paper as an emergency reviewer, I couldn't afford time to check the proofs. \n\n-- The framework is very pleasant and novel and paves the way for other researchers to build on. \n\nWeaknesses:\n\n-- Unfortunately, my enthusiasm for this paper is somewhat tempered by the theoretical results, which do have potential, but also limitations. In traditional graphical modeling, one is interested in the recoverability of the underlying graph. Here, unless I missed something, in Lemma 3, the authors are proving some results concerning the parameters of the model. Why would we be interested in the parameters? It seems like learning the structure would be able to provide practitioners with more useful qualitative information about the relationships between the processes $x_j$. \n\n-- It's also not clear to me why one would be interested in the generalization bound as stated in Lemma 2. Generalization bounds are useful in supervised learning settings, but here one is interested (as I mentioned), the quality of the graph that is recovered from the observations. \n\n-- Even if we accept that the risk and difference between the estimated parameter and the true parameter are important quantities, it is unclear what the results imply. Can the authors comment about tightness of the upper bounds on the generalization error (8) or the consistency gap in (9)? I realize that a single paper cannot contain too many results, but I fail to appreciate how the bounds in (8) and (9) depend on some hardness parameter of the learning task. I am familiar with graphical models, and usually the results depend on the minimum strength $\\beta$ in the exponential family parameterization. Here, this dependence seems to be completely absent. ", "summary_of_the_paper": "This paper introduces a brand new graphical modeling framework from the perspective of neural ODEs. Traditionally structure learning involves using sampled data to learn the structure of graphs. This paper, however, looks at the graph structure learning problem from a different viewpoint, using continuous-time dynamics inspired from neural ODEs. Theoretical guarantees on parameter estimation are provided. Some experiments on benchmark time-series datasets are also conducted. ", "summary_of_the_review": "See main review above. The paper does have a lot of potential, but perhaps it is not fully mature at this point. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "uDTqGKXH4Wq", "reviewer": "Reviewer_xAV7", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper is masterfully written: the presentation is rigorous and yet reads effortlessly. I enjoyed reading the paper throughout. \n\nThe presented method is straightforward extension of neural ODEs with adaptive group lasso with proximal methods. The main contribution is the theoretical consistency and convergence analyses of such a system. The paper discusses learning bounds in terms of data fit, and in terms of sparsity pattern. The main insight is to show that one only needs to sparsify the first layer weight of the neural ODE. \n\nThe proposed method has outstanding performance in the experiments that cover observation irregularity and sparsity, and a biological and chaotic system. The method is compared against lots of competing methods. \n\nTechnical comments:\n* Lemma 1: If xk appears in differential of xj, then surely one expects to see non-zero jacobian between them, while here it is implied to be zero instead. \n", "summary_of_the_paper": "The paper proposes to learn Jacobian-sparse neural network ODEs from irregular trajectories of a dynamical system. The main contribution is the sparsity of the ODE Jacobian, which results in learning of differential covariate causalities. Learning the differential structure is an important real-world problem. The proposed method is elegant, simple and effective, although incremental.\n", "summary_of_the_review": "This is a fantastically written paper with an effective model for an important problem with solid theory and outstanding empirical performance. \n", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "c2hlzYXs-2", "reviewer": "Reviewer_4Lyp", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The proposed method is a straightforward application of the neural ODE and group Lasso method. The novelty is low. In addition, it is unclear how the proposed method can be used to handle heterogeneous dynamic systems, which are often the case of practical systems and might make the proposed method less robust.  ", "summary_of_the_paper": "This paper proposes a penalized neural ODE method for learning graphical models for multivariate dynamic systems. ", "summary_of_the_review": "This paper provides a reasonable method for graphical modeling of dynamic systems, although the technical novelty is limited.   ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Meta-review: This main focus of this paper is graph modeling. Specifically, this paper considers a setting in which data is generated under continuous time dynamics based on neural ODE. Theoretical results regarding parameter estimation are provided. The results are also supported by experiments.

The reviewers appreciate a thorough response to their questions and think that this paper would be of interest to ICLR and ML community. Please address reviewers comments in your final version.
Decision: Accept (Poster)