Year: 2022

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.

heading: 1 INTRODUCTION
text: This paper deals with learning directed graphs from a combination of temporal data and assumptions on the parameterization of the underlying structural dynamical system. Graphical models can offer a 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; Runge et al., 2019). In these contexts, inferring graphical models from temporal data subject to practical limitations as to how finely and regularly each variable can be measured over time is a longstanding challenge.
Time series data is often assumed to be a sequence of observations from an 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) ∈ Rd) contemporaneously influencing the rate of change of the same or other variables (e.g. dx(t)),
dx(t) = f(x(t))dt+ dw(t), x(0) = x0, t ∈ [0, T ], (1) ∗Work primarily conducted while at the University of Cambridge and at the Alan Turing Institute.
where w(t) a d-dimensional standard Brownian motion and x0 is a Gaussian random variable independent of w(t). The functional dependence structure of the vector field f defines a directed graph G and associated adjacency matrix G ∈ {0, 1}d×d, i.e., Gij = 1 if and only if xj appears as an argument of fi = [f ]i. The problem of structure learning is to search over the space of graphs compatible with the data, but the pattern of observation in dynamical systems emphasize a number of differences with respect to classical graphical modelling with static data or explicitly discrete-time stochastic process.
• Observed data is sampled at a sequence of (often irregular) time points (t1, . . . , tn) and is systematically subsampled. Most work on graphical modelling with time series data assume a fundamentally discrete parameterization of the underlying structural model (e.g. based on vector autoregression models). Associations in discrete-time in general do not correspond to the structure of the underlying dynamical system and are highly dependent on the interval between observations. The same subsampled discrete model may disaggregate to several continuous models, which are observationally equivalent at the subsampled frequency, see e.g. (Runge, 2018; Gong et al., 2015; Danks & Plis, 2013) and a worked example in Appendix A. The realm of problems that involve irregularlysampled data are fundamentally out of scope in discrete-time in general. We complement this point in Figure 1 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.
• Discrete 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 is limited information to be harnessed from a given time interval [0, T ]. Learning performance depends on the number of independent samples which is a function both of the number of samples n and the length of the observed interval T .
• Non-parametric graphical modelling in dynamical systems is relatively unexplored. Existing approaches rely on specific model assumptions (e.g. linearity, additivity) to establish the consistency of 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-forward networks with multiple hidden layers and convolutional neural networks) under general observation patterns including irregular sampling. In particular, we consider penalized optimization problems of the form,
arg min fθ
1
n n∑ i=1 ||x(ti)− x̂(ti))||22, subject to ρn,T (fθ) ≤ η and dx̂(t) = fθ(x̂(t))dt, (2)
where the observation process (x(t1), . . . ,x(tn)) is given by an irregular sequence of time points 0 ≤ t1 < · · · < tn ≤ T . ρn,T (fθ) is an adaptive group lasso constraint on the parameter space of fθ. We analyze this problem with fixed dimension d and increasing sample size n 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 .
A second contribution is to propose an instantiation of this method using differential equations with vector fields parameterized by neural networks (Chen et al., 2018) to model the mean process of (1) with the advantage of implicitly inferring variable derivatives instead of involving a separate approximation step (that is common in the dynamical systems literature). This construction shows that, empirically, graphical models in continuous-time can be inferred accurately in a large range of settings despite irregularly-sampled multivariate time series data and non-linear underlying dependencies.
Code associated with this work may be found at https://github.com/alexisbellot and at https://github.com/vanderschaarlab/mlforhealthlabpub.

heading: 2 RELATED WORK
text: 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 in neural networks, using for instance a sparse one-to-one 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 a 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; Marcinkevič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 are 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 et 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)) and Gaussian processes (see e.g. (Heinonen et al., 2018; Wenk et al., 2020)) but we are not aware of any identifiability guarantees.

heading: 3 GRAPHICAL MODELLING IN CONTINUOUS-TIME
text: We consider the underlying structure of an evolving process to be described by a multivariate dynamical system of d distinct stochastic processes x = (x1, . . . , xd) : [0, T ] → X d with each instantiation in time xj(t) for j = 1, . . . , d and t > 0 defined in a bounded open set X ⊂ R.
Definition 1 (Neural Dynamic Structural Model (NDSM)). We say that x = (x1, . . . , xd) : [0, T ]→ X d follows a Neural Dynamic Structural Model if there exist functions f1, . . . , fd ∈ F such that fj : X d → R and,
dxj(t) = fj(x(t))dt+ dwj(t), x(t0) = x0, t ∈ [0, T ], (3)
with F defined as the space of analytic feed-forward neural networks with sets of parameters θ ∈ Θ defined in bounded, real-valued intervals and wj(t) is standard Brownian motion independently generated across processes j1.
We will write fθ0 = (f1, . . . , fd) for the true underlying vector field, parameterized by a set of parameter values θ0. It will be useful to define each layer of each network precisely. Let A j 1 ∈ Rd×h denote the d × h weight matrix (we omit biases for clarity) in the input layer of fj , j = 1, . . . , d. Let Ajm ∈ Rh×h, for m = 2, . . . ,M − 1, denote the weight matrix of each hidden layer, and let AjM ∈ Rh×1 be the h× 1 dimensional output layers of each sub-network such that,
fj(X) := φ(· · ·φ(φ(XAj1)A j 2) · · · )A j M , j = 1, . . . , d, (4)
where φ(·) is an analytic activation function (e.g. tanh, sigmoid, arctan, softplus, etc.) and X ∈ Rn×d 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 (t1, . . . , tn) sampled from a temporal point process with positive intensity such that,
(x1, . . . ,xn) ∼ N (µ,Σn), (5)
with a dependency structure encoded in Σn ∈ Rn×n. The closer in time two observations are the more closely correlated we can expect them to be. We assume the data to be normalized, i.e. diagonal elements of Σn to be equal to 1. µ are the instantiations of the mean process that can be described by an system of ordinary differential equations dx(t) = f(x(t))dt,x(0) = x0, t ∈ [0, T ]. Time points at which observations are made are thus themselves assumed stochastic, driven by an independent temporal point process with intensity limdt→0 Pr(Observation in [t, t + dt]|Ht) > 0 for any t > 0 with respect to a filtrationHt that denotes sigma algebras generated by any sequence of prior observations. Perfectly homogeneous and systematic subsampling has measure zero under this probability model. This is important because it will enable, in principle, to infer local conditional independencies (defined below) arbitrarily well with increasing sample size.

heading: 3.1 GRAPHICAL PRESENTATION
text: 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 x is locally independent of y given z if, for each time point t, the past up until time t of z gives us the same predictable information about E(xt|Ht(y, z)) as the past of x and y until time t, whereHt(y, z) is the filtration generated by y and z up to time t. 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 a directed edge (xk → xj) ∈ E if and only if there exist no conditioning subset of V such xk ∈ V is locally conditionally independent of xj ∈ V. Lemma 1 (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 xk appears in the differential equation of xj , i.e. ||∂kfj ||L2 6= 0. Moreover, for any f ′ such that ||∂kf ′j ||L2 = 0 there exists an equivalent vector field f such that the euclidian norm of its column vectors ||[Aj1]·k||2 = 0. Proof. All proofs are given in Appendix B.
1For analytic function spacesF , the vector field is locally Lipschitz, i.e., the system is stable and the diffusion process has a unique stationary measure that is Gaussian. We assume unique solutions also as T →∞.
This Lemma specifies an equivalence relation between functional dependence graphs given by the underlying dynamical system and local independence graphs2. It is clear that enforcing ||[Aj1]·k||2 = 0 will remove the local dependence of the j-th stochastic process on the k-th stochastic process but it is not the case that all functions f with this particular local independence (||∂kfj ||L2 = 0) necessarily have zero-valued k-th column in its parameters Aj1. This proposition shows that in such cases there exists an equivalent vector field (i.e. that defines the exact same input-output map) that does have ||[Aj1]·k||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. a simple rearrangement of the nodes in the same hidden layer leads to a new configuration that produces the same mapping as the generating network). We define next local consistency as a desirable and target property for estimators in practice.
Definition 3 (Local consistency). An estimator fθ = (f1, . . . , fd) is locally consistent if for any δ > 0, there exists Nδ and Tδ such that for n > Nδ and T > Tδ, we have ||∂kfj ||L2 6= 03, for all k, j ∈ {1, . . . , d} such that xk is locally significant for xj , and have ||∂kfj ||L2 = 0 otherwise, with probability at least 1− δ.

heading: 3.2 FINITE-SAMPLE IDENTIFIABILITY
text: WriteRn(fθ) = 1n ∑n i=1 ||x(ti)− x̂(ti))||22 (the dependence on fθ is implicit in x̂) andR(fθ) for its population counterpart. The difficulty arises from the geometry of the loss function around the set of loss minimizers,
Θ? = {θ ∈ Θ : R(fθ) = R(fθ0)}, (6)
where Θ is the parameter space. The set Θ? (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 θ is not. It is sufficient for identifiablity of the local independence model to recover the group structure of input layer parameters [Aj1] exactly. In the context of analytic vector fields f , we define the adjacency matrix G ∈ {0, 1}d×d associated with G such that Gkj 6= 0 if and only if ||∂kfj ||L2 6= 0. If this pattern can be recovered exactly, then the inferred G corresponds to the functional structure. Optimization with this group structure desiderata has been often considered before as a penalized optimization problem,
arg min fθ Rn(fθ), subject to dx(t) = fθ(x(t))dt and ρ(fθ) ≤ η, (7)
where we have suppressed the dependence of ρ 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,
ρGL(fθ) := λGL d∑ k,j=1 ||[Aj1]·k||2, ρAGL(fθ) := λAGL d∑ k,j=1
1
||[Âj1]·k|| γ 2
||[Aj1]·k||2,
respectively, where Âj1 is the GL estimate to problem (7), λGL, λAGL determine the regularization strength and may vary with n and T , γ > 0 and || · ||2 is the Euclidian norm. As with other adaptive lasso estimators, AGL uses its base estimator to provide a rough data-dependent estimate to shrink groups of parameters with different regularization strengths. As n and T grow, the weights for non-significant features get inflated while the weights for significant ones remain bounded, allowing AGL to exactly identify significant parameters.
2Note, however, that this is not true in general for a marginalization of the local independence graph (that we do not consider), i.e. in 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.
3∂kfj denotes the partial derivative with respect the k-th argument of fj , || · ||L2 is the functional L2 norm.
The following generalization bound will be useful to define convergence rates for penalized solutions.
Lemma 2 (Generalization bound). Assume Σn to be invertible and let α = (α1, . . . , αn) such that α1 > · · · > αn > 0 are its eigenvalues. For any δ > 0, there exists Cδ > 0 such that,
|Rn(fθ)−R(fθ)| ≤ Cδ ( ||α||2 n )√ log ( n ||α||2 ) , (8)
with probability at least 1− δ. Remark on interpretation. Note that ||α||2 ≤ ||α||1 = n (the data is assumed to be scaled to have variance 1 for all observations) and that larger values of ||α||2 occur with a greater difference in magnitude in the entries of α. 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 α1 (the first component of α) as more variance in the data is explained by a single direction of variation and thus decreases the effective sample size – making ||α||2 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 an 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 to Dinh & Ho (2020) with the difference that the convergence speeds differ due to sample dependency.
Lemma 3 (Convergence of Adaptive Group Lasso). Let θ̃n ∈ Θ be the parameter solution of (7) with adaptive group lasso constraint. For any δ > 0, assuming that λAGL → 0 there exists v > 0, Cδ > 0, Nδ > 0 and Tδ > 0 such that,
min θ∈Θ∗ ||θ̃n − θ|| ≤ Cδ
( λAGL + ( ||α||2 n )√ log ( n ||α||2 )) 1ν , (9)
with probability at least 1− δ. Lemma 4 (Local consistency of Adaptive Group Lasso). Let γ > 0, > 0, ν > 0, λAGL = Ω(( n||α||2 ) −γ/ν+ ), and λAGL = Ω(λ γ+ GL ), 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 an 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 defined such that [G]jk = 0 if and only if ||[Ãj1]·k||2 = 0 is the local independence graph with high probability. The group lasso will generally not be locally consistent because it forces all parameters to 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,
|A1|min > Cδ ( λAGL + ( ||α||2 n )√ log ( n ||α||2 )) 1ν , (10)
where we have defined |A1|min := min{‖[Aj1]·k‖2 : j, k = 1, . . . , d, ||∂kfj ||L2 6= 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 it are essentially necessary in the context of structure learning via parameter estimation (Van de Geer et al., 2011; Zhao & Yu, 2006), and it allows us to describe next a guarantee on the finite sample consistency of structure learning with the Adaptive Group Lasso.
Lemma 5 (Finite sample local consistency of Adaptive Group Lasso). Under the conditions of Lemma 4 with the additional minimum restricted strength assumption in (10) on the problem design for particular values of n and α, the Adaptive Group Lasso recovers the structure G exactly with high probability.

heading: 3.3 ALGORITHM: NEURAL GRAPHICAL MODELLING
text: 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 a stochastic differential equation explicitly to be the solution to an ODE initial-value problem in which fθ is a free parameter specified by a neural network. For each estimate of fθ, the forward trajectory can be computed using any numerical ODE solver:
x̂(t1), . . . , x̂(tn) = ODESolve(fθ,x(t0), t1, . . . , tn), (11)
and thus implicitly enforcing the constraint dx(t) = fθ(x(t))dt while optimizing for Rn(fθ) and thus solving for (2). Gradients with respect to θ may be computed with adjoint sensitivities and a gradient descent algorithm can be used to backpropagate through the ODE solver and the continuous state dynamics to update the parameters of fθ, as shown by Chen et al. (2018).
Remark on optimization. The adaptive group lasso constraint is not differentiable and non-separable which precludes applying coordinate optimization algorithms. However, a 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 to exact zeros in the columns of the input matrices without having to use a cut-off value for selection. The proximal step for the group lasso penalty is given by a group soft-thresholding 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).

heading: 4 EXPERIMENTS
text: This section makes performance comparisons on controlled experiments designed to analyzed 4 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 a system.
We benchmark NGM against a variety of algorithms, namely: Three representative vector autoregression 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 et 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 representative two-stage collocation method we call Dynamic Causal Modelling (DCM) in which derivatives are first estimated on interpolations of the data and a 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 to 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 F1 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-offs 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.

heading: 4.1 LORENZ’S CHAOTIC MODEL
text: 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 (Heltberg et al., 2019), etc. The continuous dynamics in a d-dimensional Lorenz model are,
d dt xi(t) = (xi+1(t)− xi−2(t)) · xi−1(t)− xi(t) + F + σdwi(t), i = 1, . . . , d
where x−1(t) := xd−1(t), x0(t) := xp(t), xd+1(t) := x1(t), F is a forcing constant which determines the level of non-linearity and chaos in the series and wi(t) is standard independent Brownian motion across i = 1, . . . , d. The initial state of each variable is sampled from a standard Gaussian distribution, d is set to 10, F to 10 and σ to 0.5. We illustrate sample trajectories of this system in Figure 1.
• Irregularly sampled data is generated by removing randomly a percentage of the regularly sampled data (with a 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 are 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= 0. Thus looking at both measures together, NGM significantly improves performance over competing approaches. And moreover, NGM’s performance is 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 on 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 to error 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.

heading: 4.2 RÖSSLER’S HYPERCHAOTIC MODEL
text: Next, we consider data with non-linear dynamics, said to exhibit hyperchaotic behaviour, to demonstrate the flexibility of learned functional relationships.
Chaotic systems, such as Lorenz’s model, are characterized by one direction of exponential spreading. If the number of directions of spreading is greater than one the behavior of the system is hyperchaotic
(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 we consider a generalization of this model to arbitrary dimensions and non-linear vector fields as in (Meyer et al., 1997). The d−dimensional generalized Rössler model is given by,
d dt x1(t) = ax1(t)− x2(t) + σdw1(t),
d dt xi(t) = sin(xi−1(t))− sin(xi+2(t)) + σdwi(t), i = 2, . . . , d− 1,
d dt xd(t) = + bxd(t) · (xd−1(t)− q) + σdwd(t).
We use typical parameters for hyperchaotic behaviour: a = 0, = 0.1, b = 4, q = 2, as in (Meyer et al., 1997). This system is observed over a sequence of 1000 time points with a 0.1 time unit interval after randomly initializing each variable to a sample from a standard Gaussian distribution and d is set to 10 and 50 to evaluate also performance in higher dimensions. σ is set to 0.1.
Results. Performance results are given in Table 1 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 high-dimensional regimes. The strongest baseline is DCM which also models the underlying graph in continuous-time.

heading: 4.3 MODEL OF OSCILLATIONS IN YEAST GLYCOLYSIS
text: We conclude with an experiment from the pharmacology literature which makes extensive use of dynamical models to determine the interaction patterns of drugs in the body. In these systems, it is the difference in scale of different biochemicals that makes structure recovery difficult.
The glycolytic oscillator model is a standard benchmark for this kind of systems. It simulates the cycles of 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 assumptions, a feature of NGM is that it explicitly discovers f and a consistent graphical structure such that the model can be used to simulate the expected effect of interventions by modifying the weights of a trained network which may have a large impact on the design of (laboratory) experiments. We show, for illustration, NGM’s simulation of the mean behaviour of the glycolytic oscillator in Figure 3, which successfully recovers the true mean dynamics.
Results. Table 1 shows that performance on this data is more heterogeneous. Although NGM outperforms or is competitive with other methods, FDR is high. Better understanding how to model
structure in systems with different scales (see this explicitly in Figure 3 with numerous variables with values equal to a small fraction of the largest observations) thus remains an important challenge.

heading: 5 DISCUSSION
text: 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.) on an increasing resolution but with heterogeneous observation patterns. A graphical model that is stable over a large number of variables, over a large number of sampling frequencies and observation patterns has attractive properties to scientists in all these domains. In this paper, we have discussed graphical modelling from a continuous-time 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 a 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 irregularly-sampled multivariate time series. We conclude with some additional remarks.
• 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, a form of violation of this assumption was included in the Lorenz system (perturbed by a constant factor F which is not modelled) for which NGM achieves almost perfect structure recovery. Some robustness to violations of the assumption of no unobserved processes may thus be expected although we did not perform a formal investigation of this phenomenon.
• Switching dynamical systems. In the context of a Neural ODE, f was defined using a continuous neural network and x(t) is always continuous in t. Trajectories modeled by an ODE can thus have limited representation capabilities when discontinuities in the state occur. Examples are switching dynamical systems which are hybrid discrete and continuous systems that choose which dynamics to follow based on a discrete switch, and are popular in 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 in which case, to infer the causal graph, partial derivatives ∂kfj would be defined piece-wise.

heading: ACKNOWLEDGEMENTS
text: 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.

heading: Appendix
text: This appendix is outlined as follows:
• Section A discusses additional related work. – Section A.1 contrasts discrete-time and continuous-time structural models considering in-
consistencies in structure learning by modelling the underlying continuous-time process in discrete-time.
– Section A.2 discusses other work modelling differential equations. – Section A.3 comments on the philosophical debate around the nature of causality in dynamical
systems and its ramifications in machine learning. • Section B 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 as a 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 a 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.

heading: A RELATED WORK

heading: A.1 GRAPHICAL MODELLING IN DISCRETE-TIME
text: 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 xi is said to Granger-cause xj if omitting the past of a time series xi in a time series model including xj’s own and other covariates’ past increases the prediction error of the next time step of xj . 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 et al., 2018; Pamfil et al., 2020)), for example in the linear case assuming,
x(t+ ∆t) = B1x(t) + · · ·+Bkx(t− k∆t) + (t), (t) ∼ P, (12)
for some distribution P . Each one of the matrices B1, . . . , Bk then describe lagged causal relationships with different lags. Subsampling occurs when the frequency of observation is lower than ∆t 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-correlations 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 discretetime 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.

heading: A.1.1 EXAMPLE
text: Assume the true dynamics of two processes x(t) = (x1(t), x2(t))ᵀ to be given by,
d dt x(t) = Ax(t), A = ( −1 2 −4 −0.5 ) , (13)
The underlying dynamics are given by A: x2 causing an increase in the rate of change of x1 (with value 2) while x2 having a large negative causal effect on the rate of change of x1 (with value −4).
The corresponding discrete-time model (with one time lag for simplicity) may be defined as,
x(t+ ∆t) = B∆t · x(t), (14)
where ∆t is the time interval of observation. The two models are connected by a simple relation, given by B∆t = exp{A · ∆t}, that can be used to uniquely compute the off-diagonal entries [B∆t]12 and [B∆t]21 indicating the strength of "causal" effects x2 → x1 and x1 → x2 respectively, under the Granger-causality paradigm (Granger, 1969).
As a first observation, irrespective of the time interval ∆t note that a single discrete-time model B∆t may describe multiple underlying mechanisms A as matrix logarithms logB∆t are not identified if B∆t 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. B∆t) may change dramatically as a function of the time interval ∆t at which the process is observed, as shown in Figure 4. Describing causality within dynamical systems in discrete-time is inherently an ill-posed problem. In fact, only for dynamical systems that are bivariate, stable, and non-oscillating can we expect consistent conclusions at all measurement intervals.
Proposition 1 (Causal inconsistency in discrete time (Kuiper & Ryan, 2018)). The sign of offdiagonal entries of B∆t and A agree for all time intervals ∆t > 0 if A defines a bivariate, stable, and non-oscillating system of differential equations.
Proof. (Given here for completeness – can also be found in (Kuiper & Ryan, 2018)) Using the relationship B∆t = exp{A ·∆t} we can write,
B∆t = V −1 exp{DA·∆t}V (15)
where DA·∆t represents the diagonal matrix containing the scalar exponentials of the eigenvalues of A multiplied by the scalar ∆t and V is the matrix of eigenvectors. This can be used to relate two estimated discrete-time matrix coefficients with each other,
B∆2t = B (∆2t/∆1t) ∆1t
(16)
Consider now a bivariate process with eigenvalues of an estimated B∆1t denoted by λ1 and λ2, and let ∆2t = n ·∆1t. Then, misleading causal conclusions are obtained if the sign of the entries of B∆2t differ from B∆1t. In this case they may be computed explicitly. It holds that the sign of off-diagonal entries are equal if and only if (λ1 − λ2) (determining the sign of entries of B∆1t) and (λn1 − λn2 ) (determining the sign of entries of B∆2t) are equal. This is the case if both λ1 and λ2 lie between 0 and 1 i.e., a 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).

heading: A.2 RELATED WORK ON MODELLING DIFFERENTIAL EQUATIONS
text: Two-stage collocation methods were first proposed by Varah (1982). The authors proposed to fit a smoothing estimate x̂(·;h, θ) to the observations y1, . . . , yn with a smoothing parameter h, and using it and its derivative with respect to t in order to estimate the vector field f in a system of differential equations,
min θ
1
n n∑ i=1 ||dx̂(ti;h)− fθ(x̂(ti;h)))||22, (17)
where x̂(·;h) = arg min x(·;h)∈H 1 n
∑n i=1 ||yi − x(ti;h)||22 is the smooth interpolation function and most
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 θ of fθ rather than of the parameters of the graphical structure induced by fθ. 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 as multiple alternatives define the same input-output map of the vector field fθ, even in the infinite-sample regime.
Broader literature. Modelling time series is a wide and varied research topic. One can train feed-forward or recurrent neural networks to approximate a differential equation (Chen et 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.e., using discrete measurements of trajectories to infer the underlying structural model.

heading: A.3 PHILOSOPHICAL DEBATE ABOUT THE NATURE OF CAUSALITY
text: A mechanistic interpretation of causality views causal claims as claims about the existence of a mechanism or process that mediates events at one time with events at a later time and has been formalized mathematically as a system of structural differential equations. Proponents of this view hold that the metaphysical connection between mechanisms and causality is 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 a mechanism as a 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 a conserved quantity. In the latter, an interaction between two processes is causal if there is an exchange of a conserved quantity between them. Mechanistic theories of causality, from a philosophical standpoint, are normally contrasted with probabilistic (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 to be causal if and only if a change in one makes a difference to the other. This distinction is 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 a 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-independent) 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 in this paper as long as we assume that the causal system of interest may be represented as a Neural Dynamic Structural model. In such a model, causation across time in dynamical systems is 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 .

heading: B PROOFS
text: We will useRn(fθ) andRn(θ) interchangeably. With the notation introduced in the main body of this paper, we will use the following Lemma that extends the standard Taylor’s inequality around a local optimum to unbounded sets θ∗.
Lemma 7 (Lemma 3.2. (Dinh & Ho, 2020)) There exist c2, ν > 0 and such that |R(θ)−R(θ0)| ≥ c2d(θ,Θ
∗)ν for all θ ∈ Θ. d(θ,Θ∗) is the minimum Euclidian norm between θ and any element of Θ∗.

heading: B.1 PROOF OF LEMMA 1
text: The uniqueness of local independence graphs in systems of fully observed stochastic processes was given in Proposition 3.6 by Mogensen et al. (2020). We restate the Lemma for convenience.
Lemma 1 (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 xk appears in the differential equation of xj , i.e. ||∂kfj ||L2 6= 0. Moreover, for any f ′ such that ||∂kf ′j ||L2 = 0 there exists an equivalent vector field f such that the euclidian norm of its column vectors ||[Aj1]·k||2 = 0. Proof. Consider the class of NDSMs with F such that fj is independent of process xk, that is ||∂kfj ||L2 = 0, and the class of NDSMs F0 with fj parameterized such that ||[A j 1]·k||2 = 0. We will how that F = F0. It is clear that the class of NDSMs with F0 is contained in the class of NDSMs with F as a process xk interacts with the vector field fj of xj only if the corresponding entries in the first layer of the analytic neural network are non-zero: F0 ⊂ F . For the converse consider a NSDM whose vector field fj of xj is independent of the k-th process xk and consider a set of processes x̃ and x such that x̃ = x except for the k-th entry of x̃ which is set to the zero function x̃k : [0, T ]→ 0. Because of independence, and because the two processes differ only in their k-th entry: fj(x̃(t)) = fj(x(t)).
Now define Ã1 be the matrix such that [Ã1]ik = 0 and [Ã1]ik′ = [A1]ik for all k′ 6= k. Then it holds that Ã1x(t) = A1x̃(t). And therefore also fj(x(t);A1) = fj(x̃(t);A1) = fj(x(t); Ã1). By definition of F0, fj(·; Ã1) ∈ F0 and thus F ⊂ F0.

heading: B.2 PROOF LEMMA 2
text: We restate the Lemma for convenience.
Lemma 2 (Generalization bound). Assume Σn to be invertible and let α = (α1, . . . , αn) such that α1 > · · · > αn > 0 be its eigenvalues. For any δ > 0, n||α||2 > 3, there exists a C > 0 such that,
|Rn(fθ)−R(fθ)| ≤ ( ||α||2 n )√ C log ( n ||α||2 ) , (18)
with probability at least 1− δ. Proof. Without loss of generality, we consider each differential equation separately. nRn(fθ) =: XTX is a sum of dependent squared normal random variables. With this notation X = (X1, . . . , Xn) ∈ Rn where Xi = (Yij − x̂j(ti)) ∈ R is a univariate random variable and Yij ∈ R is the j-th random variable at time ti given of the observation model and underlying dynamical system. X = (X1, . . . , Xn) has a joint distribution defined by a mean µ and covariance matrix Σn, as defined in the main body of this paper. We may write Z = Σ−1/2n (X − µ) and,
XTX = (Z + Σ−1/2n µ)Σn(Z + Σ −1/2 n µ). (19)
Let Σn = V TAV be the eigendecomposition of Σn where V is an orthogonal basis of eigenvectors and A is a diagonal matrix of positive eigenvalues α = (α1, . . . , αn). U = V Z then is also
multivariate normal, with expectation zero and identity covariance matrix (since V TV = V V T = In). For u = V Σ−1/2n µ, rewriting the decomposition above in terms of U and u
XTX = (U + u)TA(U + u) = n∑ i=1 αi(Ui + ui) 2. (20)
The distribution of nRn(fθ) = XTX is thus a weighted non-central χ2 random variable.
Write W = V Σ−1/2n ∈ Rn×n. By applying the concentration results in Theorem 6 and 7 in (Zhang & Zhou, 2020) we have,
P(|Rn(θ)−R(θ)| > s) ≤ exp  −C1n 2s2
||α||22 + 2 ∑n i=1 (∑n j=1Wij(f(x(tj); θ)− f(x(tj); θ0)) )2 
(21)
≤ exp { −C2 ( ns
||α||2
)2} , (22)
for all s such that,
0 < s < ||α||2 n||α||∞
+ 2 ∑n i=1 (∑n j=1Wij(f(x(tj); θ)− f(x(tj); θ0)) )2 n . (23)
C1, C2 > 0 are two scalars not depending on n or α. The remaining of the proof follows the argument of (Dinh & Ho, 2020). We define events,
A(θ, s) = {|Rn(θ)−R(θ)| > s}, (24) B(θ, s) = {∃θ′ ∈ Θ such that ||θ′ − θ||2 ≤ s
4Mδ and |Rn(θ′)−R(θ′)| > s} (25)
C = {|Rn(θ)−Rn(θ′)| ≤Mδ||θ − θ′||2,∀θ, θ′ ∈ Θ}, (26)
where the last event C is defined with respect to the Lipschitz constant Mδ of Rn (assumed to be Lipschitz with probability at least 1 − δ, that is P(C) ≥ 1 − δ). Let m = dim(Θ), there exist C3(m) ≥ 1 and a finite setH ⊂ Θ such that,
Θ ⊂ ⋃ θ∈H V(θ, ), |H| ≤ C3/ m, (27)
where we choose = s/(4Mδ). V(θ, ) denotes the open ball centered at θ with radius , and |H| denotes the cardinality of H. In other words, H -covers Θ and the inequality involving the cardinality ofH follows because Θ is a bounded subset of Euclidian space, see e.g. section 27.1 in (Shalev-Shwartz & Ben-David, 2014). By a union bound over all elements inH,
P (∃θ ∈ H : |Rn(θ)−R(θ)| > s) ≤ C3(4Mδ)ms−m exp { −C2 ( n
||α||2
)2 s2 } . (28)
Since B(θ, s) ∩ C ⊂ A(θ, s) andH ⊂ Θ we have,
P(∃θ ∈ Θ : |Rn(θ)−R(θ)| > s) ≤ C3(4Mδ)ms−m exp { −C2 ( n
||α||2
)2 s2 } + δ. (29)
Now let k = n||α||2 and let s = √ C log(k)
k for notational simplicity. To complete the proof we need to choose C such that,
C3 ( 4Mδk√ C log(k) )m exp {−C2 · C · log(k)} ≤ δ. (30)
This inequality holds if,
C3(4Mδ) m · km−C2·C ≤ δ, (31)
since (C log(k))m/2 ≥ 1, which can be obtained if m−C2 ·C > 0 and C3(4Mδ)m · 3m−C2·C ≤ δ, since k > 3 by assumption, so that,
C ≥ 1 C2 log(3) (log(C3) +m log(4Mδ) + log(1/δ)) . (32)

heading: B.3 PROOF LEMMA 3
text: To traverse this result, we will start by considering the convergence of the group lasso.
Lemma 7 (Convergence of Group Lasso). For any δ > 0, assuming that λGL → 0 there exists v > 0, Cδ > 0, Nδ > 0 and Tδ > 0 such that,
min θ∈Θ∗ ||θ̂n − θ|| ≤ Cδ
( λ ν/ν−1 GL + ( ||α||2 n )√ log ( n ||α||2 )) 1ν , (33)
with probability at least 1− δ. Proof. Recall the group lasso and adaptive group lasso penalty terms,
ρGL(θ) := λGL d∑ k,j=1 ||[Aj1]·k|| and ρAGL(θ) := λAGL d∑ k,j=1
1
||[Âj1]·k||γ ||[Aj1]·k||.
By definition, we have,
Rn(θ̂n) + ρGL(θ̂n) ≤ Rn(θ0) + ρGL(θ0), (34)
where θ̂n = arg min θ∈Θ Rn(θ) + ρGL(θ) is the parameter solution to the group lasso.
It holds then that,
min θ∈Θ∗ c2||θ̂n − θ||ν2 ≤ R(θ̂n)−R(θ0) (35)
≤ |R(θ̂n)−Rn(θ̂n)|+ |R(θ0)−Rn(θ0)|+ |Rn(θ̂n)−Rn(θ0)| (36) ≤ 2 ( ||α||2 n )√ C log ( n ||α||2 ) + |ρGL(θ0)− ρGL(θ̂n)| (37)
≤ 2 ( ||α||2 n )√ C log ( n ||α||2 ) + λGL ·K · ||θ0 − θ̂n||2, (38)
where the first inequality is due to Lemma 6 (for some c2, ν > 0), the second inequality is due to the triangle inequality, the third inequality is due to equation (34) and Lemma 1, and the fourth inequality comes from the Lipschitzness of ρGL. We have used K > 0 to denote the Lipschitz constant of ρGL. The last step is given by Youngs’s inequality, e.g. as stated in section 5.1 (Dinh & Ho, 2020), to conclude that,
min θ∈Θ∗ ||θ̂n − θ||ν2 ≤ Cδ (( ||α||2 n )√ log ( n ||α||2 ) + λ ν ν−1 GL ) . (39)
for some constant Cδ .
We will now state and prove Lemma 3.
Lemma 3 (Convergence of Adaptive Group Lasso). For any δ > 0, assuming that λAGL → 0 there exists v > 0, Cδ > 0, Nδ > 0 and Tδ > 0 such that,
min θ∈Θ∗ ||θ̃n − θ|| ≤ Cδ
( λAGL + ( ||α||2 n )√ log ( n ||α||2 )) 1ν , (40)
with probability at least 1− δ.
Proof. By the convergence of the group lasso ||[Âj1]·k||2 is bounded away from zero for any process k that causally significant for process j, k, j ∈ {1, . . . , d}. Let the set of causally significant pairs (k, j) be denoted S. Then we can define,
M(θ) = ∑
(k,j)∈S
1
||[Âj1]·k|| γ 2
||[Aj1]·k||2 <∞. (41)
since ||[Âj1]·k||2 > 0.
Let θ̃n = arg min θ∈Θ Rn(θ) + ρAGL(θ) be the solution parameters of the adaptive group lasso problem. By a similar derivation to that used in (35),
R(θ̃n)−R(θ0) ≤ 2C ( ||α||2 n )√ log ( n ||α||2 ) + λAGL · (M(θ̃n)−M(θ0)) (42)
≤ 2C ( ||α||2 n )√ log ( n ||α||2 ) + λAGLM(θ̃n). (43)
And sinceM(θ̃) is a finite positive scalar, again by a similar derivation to that used in (35),
min θ∈Θ∗ ||θ̃n − θ||2 ≤ Cδ (( ||α||2 n )√ log ( n ||α||2 ) + λAGL )1/ν . (44)
for some constant Cδ > 0.

heading: B.4 PROOF OF LEMMA 4
text: We restate the Lemma for convenience.
Lemma 4 (Local consistency of Adaptive Group Lasso). Let γ > 0, > 0, ν > 0, λAGL = Ω(( n||α||2 ) −γ/ν+ ), and λAGL = Ω(λ γ+ GL ), then the adaptive group lasso is locally consistent.
Proof. By the convergence of the Group Lasso, for any pair (j, k) of non-significant processes,
||[Âj1]·k||2 ≤ Cδ ( λ ν/ν−1 GL + ( ||α||2 n )√ log ( n ||α||2 )) 1ν , (45)
with probability at least 1− δ. It holds therefore that,
lim n ||α||2 →∞
1
||[Âj1]·k|| γ 2
≥ ∞. (46)
Now assume for contradiction that there exists a pair (k, j) of non-locally significant processes (that is, xk is not locally significant for xj) such that ||[Ãj1]·k||2 6= 0 (the ∼ notation above the matrix Aj1 denotes estimation with the adaptive group lasso) and define φ(θ̃n) to be equal to θ̃n except that non-significant parameters are set to zero. By the definition of θ̃n as minimizing the empirical risk regularized by the adaptive group lasso constraint,
Rn(θ̃n) + λAGL 1
||[Âj1]·k|| γ 2
||[Ãj1]·k||2 ≤ Rn(φ(θ̃n)). (47)
Since the regularization term on the right-hand side is zero by the definition of φ(θ̃). Then by Assumption 3,
λAGL 1
||[Âj1]·k|| γ 2
||[Ãj1]·k||2 ≤ Rn(φ(θ̃n))−Rn(θ̃n) (48)
≤Mδ||φ(θ̃n)− θ̃n||2 (49) = Mδ||[Ãj1]·k||2. (50)
But since we have assumed ||[Ãj1]·k||2 6= 0 it follows from above that λAGL 1||[Âj1]·k||γ2 ≤ Mδ which is a contradiction of Lemma 7 that proved the convergence of the group lasso and in particular that lim n
||α||2 →∞
1
||[Âj1]·k|| γ 2
=∞ for non-locally significant pairs of processes (k, j).

heading: B.5 PROOF OF LEMMA 5
text: We restate the Lemma for convenience.
Lemma 5 (Finite sample local consistency of Adaptive Group Lasso). Under the conditions of Lemma 4 with the additional minimum restricted strength assumption in (10) on the problem design for particular values of n and α, the Adaptive Group Lasso recovers the structure G exactly with high probability.
Proof. First, note that for the set of loss minimizers Θ? defined in eq. (6) and by using the fact that neural networks are analytic, it does hold that for any two locally dependent processes xk and xj the first layer parameters of any model with minimum loss are bounded away from zero, i.e. ||[Aj1]·k||2 ≥ c for some c > 0.
To see this, assume for a contradiction that no such c exists, and therefore that there exists [Ãj1]·k ∈ Θ? such that ||[Aj1]·k||2 = 0 since neural networks are analytic and each one of the parameters is defined in bounded intervals. This would imply that there exists a neural network with the same input-output relationship as the true model fj that does not depend on its k-th input, which is a contradiction because ||∂kfj ||L2 6= 0. Next, given the minimum strength condition on the column norms of first layer parameters related to locally dependent processes,
|A1|min > Cδ ( λAGL + ( ||α||2 n )√ log ( n ||α||2 )) 1ν (51)
where recall that |A1|min := min{‖[Aj1]·k‖2 : j, k = 1, . . . , d, ||∂kfj ||L2 6= 0}, we have that by Lemma 3 that estimated parameters ||Ãj1]·k||2 are bounded away from zero for specific values of α and n since,
min θ∈Θ∗
||Ãj1]·k − [A j 1]·k||2 ≤ Cδ ( λAGL + ( ||α||2 n )√ log ( n ||α||2 ))1/ν , (52)
with high probability.

heading: C EXPERIMENTAL DETAILS

heading: C.1 RESULTS AS A FUNCTION OF FEATURE DIMENSIONALITY, NEURAL NETWORK PARAMETERIZATION AND RUN TIME COMPARISONS
text: 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 a function of more variables. We also include run-time comparisons and performance comparisons as a 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.
Note that the parameter of interest is a norm over the columns of parameter matrices in the first NN layer which we found to be sufficiently coarse not be sensitive to small variations in architecture choices. Our default parameterization worked well across all datasets we analyzed. In Table 4, we show additional experiments that cover 6 possible configurations that span reasonable choices that a practitioner may make. We find that performance results are largely invariant to differences in the number of layers and layer size.

heading: C.2 RESULTS USING THE AREA UNDER THE ROC CURVE
text: In this section we 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 5 contains comparisons for the Rössler and Glycolytic experiments.

heading: C.3 PURELY SYNTHETIC EXPERIMENT
text: This experiment investigates the behaviour of different regularization schemes with a purely synthetic data generating mechanism. We generate 1000 observations of 50 variables X ∈ R1000×50 using regular evaluations of the process defined by,
d dt xj(t) = gj(x(t)) + dwj(t), x(0) = x0, (53)
for j = 1, . . . , 50, where each vector field component gj : R50 → R is parameterized by a neural network with three hidden layers of 10 nodes, such that gj depends only on 5 (locally significant) of its 50 arguments. The initial state x0 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 at 0.1 units.
The results are given in Table 6. We write NGMGL for NGM with group lasso regularization, NGMAGL for NGM with adaptive group lasso regularization and NGML 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 enforce sparsity strictly. As a consequence, we may need cut-off values to interpret NGMGL locally which are difficult to specify in practice, while NGMAGL sets most non-local processes to zero exactly without further post-processing.

heading: C.4 ALGORITHM IMPLEMENTATION

heading: C.4.1 NEURAL GRAPHICAL MODELLING (NGM)
text: Proximal gradient descent. The proximal step for the group lasso penalty is given by a group soft-thresholding operation on the input weights.
In each iteration, proximal-gradient steps make two update computations to the relevant parameters, denoted here θ ∈ {[Aj1]·k, j = 1, . . . , d}. They are,
1. θ ← θ − α∇L(θ) 2. θ ← argminw∈Rd ||w − θ||2 + αλGL,n ∑d k,j=1 ||[A j 1]·k||
where the second part is the proximal operator with respect to the constraint, which for each θ ∈ [Aj1]·k is equivalent to a soft-threshold group-wise update,
[Aj1]·k ← [Aj1]·k
||[Aj1]·k|| max{0, ||[Aj1]·k|| − αλGL,n} (54)
Regularizing constants are chosen from the set {0.001, 0.01, 0.05, 0.1, 0.5, 1, 2} with γ = 2 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 converting the weights to the presence / absence of edges in the graph. A non-zero estimate of ||[Aj1]·k|| is considered as presence of an edge in the underlying graph and a zero estimate of ||[Aj1]·k|| is considered as absence of an edge.
Architecture. The integrand fθ was taken to be a feed-forward neural network as described with a single hidden layer of size 10 and elu activation functions after each layer except after the output layer. In each case we used the Adam optimiser as implemented by PyTorch. Starting learning rates varied between experiments (with values between 0.001 and 0.01) before being reduced by half if metrics failed to improve for a 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 are of interest which we found not to be sensitive to the exact epoch choice. The same architecture was used for all experiments.

heading: C.4.2 NEURAL GRANGER CAUSALITY
text: We implement Neural Granger Causality (Tank et al., 2018) with the code provided by the authors at https://github.com/iancovert/Neural-GC.
We considered two architectures: an MLP to fit the lagged time series explicitly to its next value in time, and a LSTM to model the hidden state capturing the relevant history information. We follow the author’s implementation and use a 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 / absence of edges in the underlying graph can simply be read off as the non-zero parameters.

heading: C.4.3 DYNAMIC CAUSAL MODELLING
text: 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) ∈ Rd (in our case computed separately with smoothing spline approximations).
Architecture. In our implementation, we parameterize f with neural networks. We use a single hidden layer of size 10 and elu activation functions after each layer except after the output layer. In each case we used the Adam optimiser as implemented by PyTorch.
The variant of Dynamic Causal Modelling we used, with neural networks to approximate the mapping between variable states and their derivatives, is our own implementation. The architecture is similar
to that used in NGM with the exception that derivatives are approximated a 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 to preserve the trajectory of the curves in each data separately.
The optimization problem thus consisted in fitting the vector field fθ explicitly such as to fit the approximated derivatives dxdt (t) evaluated at a given time t as well as possible. Derivatives are computed from a cubic spline interpolation of the time series. We manually tune interpolation hyperparameters to obtain a visually faithful approximation of the observed trajectory for each dataset.
The optimization objective is given as,
arg min fθ∈F
1
n n∑ i=1 ∣∣∣∣∣∣∣∣dXdt (ti)− fθ(x(ti)) ∣∣∣∣∣∣∣∣2 2 + ρAGL(fθ). (55)
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 to determine the presence / absence of edges in the underlying graph. Presence / absence of edges is defined by non-zero parameter values.

heading: C.4.4 PCMCI
text: PCMCI (Runge et al., 2017; Runge, 2018; Runge et al., 2019) is a discrete-time two-step approach that uses a version of the PC-algorithm with the momentary conditional independence test to account for autocorrelation in the time series.
PCMCI was implemented with the python package provided by the authors at https://jakobrunge.github.io/tigramite/.
Threshold selection. We chose to adjust for multiple testing with Benjamini-Hochberg’s procedure and considered associations significant, determining presence / absence of edges, with p-values below 0.00001 (chosen here because it gave a good trade-off between TPR and FDR i.e., the values with maximum F1 score).
C.4.5 SVAM
SVAM (Hyvärinen et al., 2010) was implemented with the implementation provided by the authors at https://github.com/cdt15/lingam with the BIC model selection criterion and a single lagged variable. Including more lagged variables did not alter our results much but adds an additional choice as to how to define causality as 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.

heading: C.5 DETAILS ON GLYCOLYSIS EXPERIMENT
text: Pharmacology is a branch of medicine that makes extensive use of 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 nonlinearities. The glycolytic oscillator model is 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.
The dynamics, defined by 7 biochemical components denoted x1, x2, . . . , x7, are given by the following system of equations (Daniels & Nemenman, 2015),
d dt x1(t) = 2.5− 100 · x1(t) · x6(t) 1 + (x6(t)/0.52)4 + 0.01dw1(t)
d dt x2(t) = 2 · 100 · x1(t) · x6(t) 1 + (x6(t)/0.52)4 − 6x2(t) · (1− x5(t))− 12x2(t)x5(t) + 0.01dw2(t)
d dt x3(t) = 6x2(t) · (1− x5(t))− 16 · x3(t) · (4− x6(t)) + 0.01dw3(t)
d dt x4(t) = 16 · x3(t) · (4− x6(t))− 100 · x4(t) · x5(t)− 13 · (x4(t)− x7(t)) + 0.01dw4(t)
d dt x5(t) = 6x2(t) · (1− x5(t))− 100 · x4(t) · x5(t)− 12x2(t)x5(t) + 0.01dw5(t)
d dt x6(t) = −2 · 100 · x1(t) · x6(t) 1 + (x6(t)/0.52)4 + 32 · x3(t) · (4− x6(t))− 1.28 · x6(t) + 0.01dw6(t)
d dt x7(t) = 1.3 · (x4(t)− x7(t))− 1.8 · x7(t) + 0.01dw7(t)
A sample of the trajectories is given in Figure 6, where after the first few time units, the system settles down onto a 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 1 of (Daniels & Nemenman, 2015). The data is stacked into two matrices for X ∈ RT×7 (and dX ∈ RT×7 for methods using approximated derivatives) where each row of X is a snapshot of the state x in time.

Review summaries:
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. 

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.