Graph Neural Flows for Unveiling Systemic Interactions Among Irregularly Sampled Time Series
Authors: Giangiacomo Mercatali
Andre Freitas
Jie Chen
Venue: NeurIPS 2024 poster (2024)
Keywords: Graph learning
neural flows
time series
Primary area: machine_learning_for_other_sciences_and_fields
TL;DR: A graph-based continuous-time model is proposed to unveil systemic interactions and improve time series tasks such as classification and forecasting.

Abstract: Interacting systems are prevalent in nature. It is challenging to accurately predict the dynamics of the system if its constituent components are analyzed independently. We develop a graph-based model that unveils the systemic interactions of time series observed at irregular time points, by using a directed acyclic graph to model the conditional dependencies (a form of causal notation) of the system components and learning this graph in tandem with a continuous-time model that parameterizes the solution curves of ordinary differential equations (ODEs). Our technique, a graph neural flow, leads to substantial enhancements over non-graph-based methods, as well as graph-based methods without the modeling of conditional dependencies. We validate our approach on several tasks, including time series classification and forecasting, to demonstrate its efficacy.

Graph Neural Flows for Unveiling Systemic Interactions Among Irregularly Sampled Time Series
Giangiacomo Mercatali
Andre Freitas
HES-SO Genève
Idiap Research Institute University of Manchester NBC, CRUK Manchester Institute andre frei tasdiiap..... ch
University Manchester giangiacomo. mat.eaaaaommmacaaallaalaaaaaaaae www.eraaaaaaaaaeeeeeeemmm Chesge.ch
Jie Chen
MIT-IBM Watson Al Lab IBM Research
chenjie@us, ieeus ibm. com
Abstract
Interacting systems prevalent in nature. is challenging to accurately predict the dynamics of the system if its constituent components are analyzed independently. We develop graph-based model that unveils the systemic interactions of time series observed at irregular time points, by using a directed acyclic graph to model the conditional dependencies (a form of causal notation) of the system components and learning this graph in tandem with a continuous-time model that parameterizes the solution curves of ordinary differential equations (ODEs). Our technique, a graph neural flow, leads to substantial enhancements over non-graph-based methods, as well as graph-bace based methods without the modeling of conditional dependencies. We validate our approach on several tasks, including time series classification and forecasting, demonstrate its efficacy.
Introduction
Real-life dynamical systems consist of a group of components interacting in a complex manner. With time series data for each component, predicting the system dynamics remains challenging, because modeling each component independently is straightforward while accounting for their interactions is hard without priori knowledge. Sometimes, these interactions are causal. For example, "phantom jams" in which a small disturbance (e.g., a driver hitting the brake too hard) in a heavy traffic can be amplified over large area of the transportation network (16). While traffic congestion often enjoys spatial proximity, outage of a power network can be propagated non-locally over the grid; i.e., sequence blackouts jumps across hundreds of kilometers (2I). We pose this question: What time series models best capture the interactive nature of different system dynamics?
modern tools to enhance time series models when multiple time series are interconnected by a given graph 33 In these models, time series are encoded by using recurrent neural network (RNN); at every time step, GNN is used to aggregate features over the graph. When the graph is unknown, graph-structure learning approaches have been proposed: some approaches learn single interaction graph over time 33 41 [1] while others infer
Graph neural networks (GNNs) 49
46
a different one at each time step (00). These models are in general discrete-time models, suitable for regularly spaced time points 39 40
WWrk done while the University of Manchester
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

To handle irregular time points, we consider continuous-time time models. The celebrated neural ordinary differential equation (ODE) technique 7 models a time series as the solution of an unknown ODE and optimizes the ODE parameters by using gradients computed through the adjoint. Neural ODEs were later adopted for latent variable modeling [39 5 which introduced discontinuities at observed time points for reducing prediction errors and variances. Neural ODEs were also adopted for graph- based modeling [38] 24 25 26
10, 27 3 where the equation accounts for multiple time series and the right-hand side uses graph to associate the different series. These models construct a graph, which could be time-dependent, in various manners, such as based on node features, co-observations within a sliding window, latent representations, or attentions.
In this work, we propose learning graph that reveals the dependency structure of the time series. To this end, we consider a form of causal notation the Bayesian network [36] 37 which is a directed acyclie graph (DAG), where a node conditionally independent of its non-descendents
given its parents [34 12 43 22) Such conditional dependence structure specifies how component
dynamics depends on each other. This model has a potential for causal discovery when one interprets the learned graph unknown a priori (e.g., how blackouts cascade over the power grid, whose known topology differs from the unknown influence graph (21)). More importantly, when modeled properly, the graph can improve the performance of downstream tasks because of the capturing of systemic interactions.
Our proposed model is a graph neural flow (GNeuralFlow). A neural flow 4 is the learned solution of an unknown ODE based on irregularly sampled time series; it advantageous over the neural ODE technique in that it models directly the ODE solution rather than the right-hand side, thus avoiding repeating calls of numerical solver, whose cost could be expensive. We condition multiple neural flows, one for each time series, on the DAG, and we instantiate their interactions as GNN; e.g., a graph convolutional network, GCN (31]. The graph convolution therein augments the parameterization of the ODE solution by aggregating the information of the neighboring time series at each time point, fitting a graph-conditioned ODE that models the interacting system.
Thus, GNeuralFlow is advantageous over prior graph ODE approaches (c.g.. GDE 38 LG-ODE 24]. CG-ÔDE (25), CF-GODE (26). STG-NCDE [10], MTGODE (27). RiTINI
two aspects. First, through learning, the graph reveals the conditional dependencies of the time series, offering a more intuitive structure for analysis. Second, removes the reliance on numerical ODE solvers and gains computational efficiency. We demonstrate empirical evidence to show that GNeuralFlow outperforms graph ODE approaches in several downstream tasks, on both synthetic and real-life data.
We highlight the following contributions of this work:
We propose a novel graph-based continuous-time model GNeuralFlow for learning systemic interactions. The interactions are modeled as a Bayesian network, which can be learned in tandem with other model parameters.
We design model parameterizations by leveraging GNNs to encode the systemic interactions. These parameterizations can additionally be used in latent variable modeling.
We demonstrate the use of GNeuralFlow in regression problems and latent variable modeling and show the performance improvement in several time series classification and forecasting benchmarks.
2 Background: Neural ODE and Neural Flows
Denote by x(t e Rd the solution of an ODE
X = ftt,x)
(I)
under well-behaving conditions (e.g., a specified initial condition of x and Lipschitz continuity of f). Neural ODE (7) is modeling technique that allows uncovering the trajectory x(t) without known right-hand side f. The technique parameterizes f by a neural network with parameters O such that the trajectory x is function of Ó. Through matching the trajectory with observed data at a few (possibly irregular) time points by using a loss function L, the vector field f is unveiled. Given the initial condition x(ln) = Xo. we write x(lo).... x x(IN) = ODESolve(f. Xo, (lo..... IN)), where the solutions at times 4, ..,.... In are obtained by invoking any blackbox numerical ODE solver (such as Kutta (7).. The training of the model parameters 0 requires the gradient VeL, which can
2

---

be economically computed by using the adjoint == Vx rather than expensively pooopprpgggting through the ODE solver.
Neural ODE has two far-reaching impacts. First, it is a continuous-time technique, which is a better alternative to discrete-ting timeete--mmeeeeee techniques (such RNNs) for modeling irregularly sampled time series [39] Second, it leads to a continuous version of the normalizing flow [35].
Other than directly modeling the observed data Chen et al. 7 proposed to use neural ODEs as latent variable models, which model the latents z instead; that is, - f(t, z). A straightforward idea is build variational autoencoder (VAE) [29]. where the encoder is RNN that evolves the hidden state h(t) over training time points and concludes a latent variable Zu as the ODE initial condition. A drawback is that this approach still uses discrete-time model (RNN) to handle irregularly sampled observations. Two approaches mitigating this drawback are ODE-RNN [39 and GRU-ODE-Bayes (5] Both approaches demonstrate smaller prediction errors and variances compared with the vanilla neural ODE + VAE approach.
Another drawback of neural ODE is that it invokes numerical solver, often multiple times in the adjoint computation because of multiple time intervals, which can be rather time consuming. A neural flow 4 is an alternative to neural ODE as it models the solution of ( directly:
x = Ft,,
by using a parameterized function F' that depends on the initial condition X Optimizing the parameters of F can be more efficient because ODE solvers no longer needed. A neural flow is not to be confused with a normalizing flow. Neural flows can replace the use of neural ODEs in latent variable models ODE-RNN and GRU-ODE-Bayes.
3 DAG-Based ODE for Modeling Systemic Interactions
In this section, we motivate the form of ODE considered in this paper based on DAG modeling.
3.1 DAG Model for Systemic Interactions
In probabilistic graphical models, the conditional dependence structure a principled framework for modeling systemic interactions. Therein, Bayesian network (36) of n random variables ......, y .....''''
a DAG with these variables as the nodes. Let A E R Rx be the (weighted) adjacency matrix of the DAG, where dij # 0 means that y is parent of y². A Bayesian network describes the conditional dependence structure of the variables; namely, a node is conditionally independent of its non-descendents given its parents. Therefore, the joint probability p(y'.....
....,... .......")) be factorized v"") = [T p(y" ppayy")) ))), where pa(v') = {v² aij * 0}
into a much simpler form: p(y'. ".y',...,,
denotes the parent set of y². The conditional dependence is necessary condition for the causal relationship between parent o and child y² [37].
The (linear) structural equation model, SEM 45
is commonly used tool to further quantify the conditional probabilities. Without loss of generality, assume that the random variables are topologically sorted according to the partial ordering x, where i < iff there exists an edge from to j. Then, the DAG adjacency matrix
|a:j] is strictly upper triangular. The SEM model reads
y'
y² d12y²
+ E2
y = d1ay" + A23yy²
63
8" ainy' 1 + agn.y² +... an-1nn Ens
where the residuals c....., En are (possibly correlated) Gaussian noises. In this model, y² depends on y² only when i j. Importantly, some of the above a,y can be zero. Then, yd independent of such y' given the rest.

---

ron-interacinng
interacting
non-interacting. new
Interacting. new Inc.
02
200
42
Figure 1: Left two: Trajectories of a non-interacting system and an interacting system (using interaction matrix A), under the same initial conditions. Right two: Replica of the left two systems but the initial conditions are changed. Trajectories change on the rightmost plot.
3.2 DAG-Based ODE as Continuous- Time Models
Consider an autonomous ODE Bx with the initial condition x(0) Xo E TR2, This ODE describes 2D vector field Bx; each solution curve x(t) expmmBttxx given the initial point x is a streamline instantancously tangential to the vector field. Throughout the paper, we use expmn to denote matrix exponential for matrix arguments and exp to denote element-wiss exponential.
Now consider n trajectories x (tt)....,x"(tt Inspired by SEM, we model that (i) the vector fields that generate the n trajectories follow the same conditional dependence structure governed by A and (ii) the residual Bx E/ aiy ax' gives the velocity x for each Mathematically,
Bx
x
Bx B12xx
+x²
Bx a13xx Bx
+ x
(2)
Bx" = ain din Bx 02n Bx ... I,, Bx +x"
In other words, the solution curve x for each j and any initial point xy is streamline instantaneously tangential to the residual field.
The behaviors of non-interacting and interacting systems are fundamentally different. Figure | illustrates an example. The first plot shows three independent trajectories satisfying x* Bxx with initial conditions x for i 1, 3. The second plot shows three conditionally dependent trajectories, under the same initial conditions, but interacting through the DAG adjacency matrix
0
0.5
0
0
0
0.7
0
0
Because x is independent of the rest, it is the same in both plots: but in the second plot, because x² depends on x! and x depends on x2, these two trajectories are different from their counterparts in the first plot. Moreover, in the third and fourth plots, change the initial conditions. As long as the new initial points are along the original trajectories, the new trajectories still follow the old ones when A does not exist; however, when A exists, x² and x3 deviate from the original trajectories, because of the conditional dependence.
3.3 From Linear Dependence to General
In general, the dependence among the time series may not be linear, even though structurally is governed by the matrix A. For example, nonlinear SEM may inspire the ODE system Bx x!. Bx² Bx1 + x², Bx Bx @ Bx + x3, where the dependence of x on x and x" is not linear. Thus, consider general ODE systems of the form
ftt {x') U pa(x))),
1= 1.
where recall that pa(x)) denotes the set of parents of x5. Equivalently, we write X = (t, X, A) in
the matrix form. Note that 2 is special case of (3), which can be written X = (I - A T
Note also that (3) is permutation equivariant and A can be any DAG matrix, not necessarily upper triangular ones.

---

GNeuralFlow: An ODE-Solver-Free Method
4.1 Problem Sctup and Model Framework
Problem 1. Let A E R be the weighted adjacency matrix of a DAG and let X(t): R + Raxd be the solution curve the initial-value ODE system
X f(t, with X(0) = Xo,
where the right-hand side f unknown P Given data ((t)..... ((,)) at irregular time points, develop model that predicts (t) for any t 2 to as well as A. To account for practical use, at some time point ',, some rows X (1,) may be missing.
growing body of research addresses the problem when X has a single row (i.e., n = 1; hence, A is irrelevant), notably through using neural network to parameterize f and using numerical solver to evaluate X at to..... iN 7 39 5]. When n > 1 and the rows of the system are independent (i.e., f identically the same function for each row), these methods straightforwardly apply through batch training. However, when the rows of the system are not independent, the problem becomes rather challenging. One may flatten 4 into nd-dimensional problem, but such high dimension renders approaches using numerical solvers too costly.
Instead, we use neural network to parameterize the solution of (4) directly; i.e..
(() = F(t, Xo. A);
(5)
where the solution a function but depends on the initial Xo as well as A. The neural network parameterization cannot be entirely free. First, the solution F should satisfy the initial condition F(0, X., A) = Xo. Second, the fundamental theorem on flows [32] Theorem 9.12] asserts that every smooth f with an initial condition determines unique Ftt, Xo. A) and for any t, F(t, is diffeomorphism. Therefore, we formulate our model for Problem I Tin the following.
Solution framework. The model for X() is a neural network Ftt, X, A) that satisfies:
1. F(O, Xo, A) Xo:
2. F(t, X. A) is invertible in X for any t and A; equivalently, the streamline Ft, A) given any
Xo and A self-intersecting.
4.2 Graph Encoder
We will make heavy use of GNN to encode the DAG adjacency matrix A for parameterizing F., For simplicity, we employ the seminal architecture GCN. A (popularly used) two-layer GCN reads
X GCN(G, X) Á Rel.UAXX)J,,
where is the input node feature matrix, X contains the transformed features, and W and U are parameters. GCN defines Å as symmetric normalization of A, but many alternatives are viable, such as simple scaling of A. We also consider
Â = I-I-A AT/7, where 7 = max, { EEtz BByl) and B-A+AT.
This definition motivated by SEM, where I A the operator (Section 3.3
Here, A is not symmetrized because doing cannot distinguish edge directions. A benefit of taking (7) is that the scaling factor 7 leads to a bounded spectral norm, which is an ingredient of invertibility required by the Solution framework.
Theorem 1. For any DAG adjacency matrix A, the matrix À defined in 7 admits Å 1/2 S 2.
"I practice, the initial time point to may not be zero, in which case one may shift the ODE along the temporal dimension by to.
5

---

4.3 Parameterization of F
The neural network F (5) can be defined in several ways by incorporating the graph encoder while satisfying the Solution framework.
ResNet flow. The first design is the ResNet architecture
F(I, X + (1() A,),
which building block of invertible networks M.. Here, (p(t) satisfies (o(0) such that the
requirement F(0, Xo. A) = Xo is met. Additionally, if (p-) E [0, 1] and g(t, A) is a contractive mapping, then Ftt,,, A) is invertible.
We let be the tanh function and parameterize by using two MLPs together with a GCN:
g(t, X) MLP (X|XX(I (X
X -CC(AA,
where || denotes concatenation row-wise and each MLP acts on the input matrix row-wise indepen- dently. The neural network is generally not contractive, but bounding the spectral norm of each linear layer can theoretically guarantee contraction of MLP [9]] Moreover, Theorem iindicates that bounding the spectral norm of the GCN parameters can guarantee contraction of GCN as well (because A ll2 is bounded). Thus, in theory, g can be made contractive. In practice, regularization is used to encourage a small Lipschitz constant of g [19].
GRU flow. The second design mimics the GRU 9
F(,,X,A) A) X + ()
X =
(10)
where 1tht k = 1, is computed by
r* (1,X) 3 si (1, (I,)
(1,1,) = (anh(S (1, (1,,,()) X)),
21,,, X) = a sigmoid (1)), (1, =1, @ (1, X) X).
The base form X + (t)) (tt, X), analogous to ResNet, comes from (5) who derived an ODE with the right-hand side being h through algebraic manipulation of the GRU. We extend the base form by including an analogous term 1/2 that incorporates the graph encoder. It can be shown that F is invertible under deliberate choice of a and when the MLPs ft. f f and the GCN are contractive. As discussed earlier, Theorem indicates that GCN can be made contractive similarly as the MLPs through bounding the spectral norm of their parameters.
Theorem 2. If f (t, -), (t.-). -). f² (t, -), and GCN(A, contractive, the function F(t, A) defined in (TOp is invertible whenever (53 + 6) 2.
Coupling flow. For the third design, use normalizing flows, because they invertible by definition. An example of the normalizing flow is the coupling flow [13). Let {U, V be a partitioning of the column indices I...., With the graph encoder, we extend usual coupling flow block the following:
F(I, X, A) X O exp (p.(() u(1, Xv, Xv)) + (v.(()) (00))-(( o(1, Xv, Xv)
(1I)
F(1,X,A)y= =X,,
Xv,
Xy = GCN(A,,
where
u(t, Xv, x,) MLP3 (MLP' II MLP MLP2(Xv(1)) (xvltt)
(t, Xv. Xv) = MLP (MLP'
|| MLP² (vvl())
Here, Pu and P are two functions that have a range [0, 1] and attain at the origin. The input X split into X and Xy and the second part goes through the GCN encoder, producing Xv. Note that one cannot apply the GCN encoder on the entire X; otherwise, the U block will have dependency on itself in the scaling and the shift. The scaling network u and the shift network u are essentially MLPs that share initial layers; they take both Xy and Xy inputs.

---

4.4 Learning the Graph
GNeuralFlow contains two sets of parameters: the DAG matrix A and other parameters of (call them 0), including the flow parameters, the graph encoder parameters, and possibly other parameters (e.g.. in latent variable modeling. Appendix F). Let us use [(A, 0) denote the training loss (which could be the quadratic loss in regression models, or likelihood/El .BO loss in latent variable models), making an explicit distinction between and Then, the learning problem is:
min
((A, 0) s.t. A corresponds to DAG. A,
(12)
The DAG constraint is combinatorial, which makes the problem NP-hard [8] Fortunately, it is known that A is a DAG matrix iff tr(expm(A @ A)) = or tr((1 - aA @ A)") n for any a * 48 47). Hence, (12) becomes an equality-constraine problem over continuous variables, to which the augmented Lagrangian method [2] is an effective solution. We discuss the optimization details in
Appendix(D D
5 Experiments
We conduct a comprehensive set of experiments to demonstrate that the proposed graph-based approach effectively improves the performance of time series tasks. The experiments are done on four synthetically generated interacting systems and four real-life datasets. Details of these datasets and their tasks given in Appendix G In all experiments, test sets. We train with early stopping by using Adam and report the results on the test set. Standard errors obtained by performing five repetitive runs. Hyperparmeter details are given in Appendix H All experiments are conducted on a machine with an Nvidia A A100 GPU, 8 CPU cores, and 80GB main memory.
split the data into train, validation, and
5.1 Synthetic Systems
We generate four synthetic systems with the graph size varying from to 30. These systems follow
the graph-based equation (4) solution (5). They are named "Sink," "Trianllee" "Sawtooth," and "Square," " following those defined in 4 but we add a graph to make the system SEM-like. For
example, Triangle is generated by following Ft,, X, A) = (I - ) + S sign(sin(w)) du); see Appendix G
for other systems and details.
1.0 x100>
Sink
x10-> Trianglo
x10 1SaWttotot
x10-> Square
1.5
1.5
Neural Flow
NSE
0.5
1.0
1.0
GNeura/Flow graph)
0.5
0.5
GNeuralFlow (true graph)
0.0
20
20
20
30
20 #groph nodes
#graph nodes
#graph nodes
#groph nodes
Figure 2: Comparison with neural flow for forecasting on synthetic systems. (ResNet flow).
The forecast results are reported in Figure Two observations follow. First, across datasets and across system sizes, GNeuralFlow lowers the MSE of neural flows. This should not be surprising, since the systems are generated to interacting and thus standard neural flows that treat the trajectories independently struggle to capture the influence of the graph. Second, GNeuralFlow performs similarly when the graph is either learned or supplied by the ground truth. This observation indicates that the learned graph sufficiently encodes the interacting nature of the trajectories.
We further compare GNeuralFlow with various baselines, including neural ODE, neural flows, graph ODEs (GDE [8) LG-ODE (24). and CF-GODE (26)). graph learning methods (NRI [31 and dNRI (20) ) and non-graph GRU variant for time series (GRU-D (6). For graph ODEs, the ground- truth graph is used. For neural flows and GNeuralFlow, all three flow designs are experimented with. T that across all datasets, GNeuralFlow significantly outperforms the baselines; moreover, GNeuralFlon also significantly outperforms neural flow for each flow design. These findings indicate that our graph encoder is rather effective and the modeling of conditional dependencies (a DAG structure) is advantageous over that of other graph structures.
7

---

Table 1: Comparison with non non-graph neural flows/ODE, graph ODE, and other time series methods on synthetic systems (5-node graphs). Best is boldfaced and second-best highlighted in gray
Sink
Triangle
Sawtooth
Square MSE (x 10
MSE (x 10-4)
MSE (x 1010 1x10->)
MSE (x (x10-))
Neural ODE
(:: 0.03)
8.32 (:: 0.24)
9.32 (#: 0.36)
16.8 (#: 0.39)
No Graph
Neural flow (ResNet)
8,41 (:: 0.05)
4.01 (#: 0.52)
4.73 (#: 0.06)
9.61 (#: 0.02)
Neural flow (GRU)
10.9 (:: 0.43)
10.3 (4: 0.45)
16.1 0.41)
17.2 (#: 0.51)
Neural flow (Coupling)
9.31 (#: 0.23)
12.2 (: 0.4))
14.2 0:244)
13.0 (# 0.63)
GRU-D
12.3 (# 0.23)
11.3 (# 0.32)
17.6 (=
18.7 (# 0.31)
Graph
ODE
GDE
10.4 (1 0.20)
3.99 (1 0.05)
7.65 (1: 0.03)
15.89 0.81)
LG-ODE
8.57 (1. 0.06)
3.58 (1 0.21)
7.07 (1 0.02)
13.99 0.73)
CF-GODE
8.60 (d: 0.14)
7.19 (d. 0.02)
8.19 (1: 0.03)
13.53 (d: 0.11)
Graph
Method Leam
NRI
5.25 (#: 0.02)
3.96 (#: 0.16)
4.99 (# 0.12)
9.39 (# 0.45)
dNRT
5,40 (d: 0.04)
3.39 (d: 0.09)
4.97 (#: 0.21)
9.78 (#: 0.21)
GNowralFlow (ResNet)
3,95 (d: 0.32)
2.32 (d: 0.13)
3.84 (#: 0.06)
8.24 (E: 0.64)
Our
GoewralFlow (GRU)
6,83 (::: 0.23)
5,41
(#: 0.23)
5.11 (# 0.13)
9.14 (# 0.61)
GNeuralFlow (Coupling)
4.45 (: 0.51)
3.21 (#: 0.34)
4.25 0.09)
8.33
0.23)
Truth
Noise
0.01
Noise 0.05
Noise 0.1
Noise 0.3
0.75
0.50
0.25
0.00
Data
MSE 1.31e-02
MSE 5.63e-03
7.95e-03
2. 2.76c-03
1.0
1.0
Triangle
1.0
1.0
Square
DAG metric values
0.8
0.9
0.8
0.8
TPR
0.6
0.6
0.6
0.6
FDR
FPR
190 SHD
0.2
0.2
0.2
0.2
0.0
0.0
0.2
0.4
0.0
0.0
0.2
0.4
0.0
0.2 noise
0.4
0.0.0
0.2
0.4
Figure 3: Graph learning quality and forecast quality. Top two rows: Sink (20 nodes); bottom row:
all four datasets (20 nodes).
With the above encouraging results, we investigate the quality of graph learning. While the learning approach discussed Section[4.4] works for a generally initialized A (as demonstrated by the previous plots), we perform a more in-depth investigation by initializing A through perturbing each entry of the ground truth with zero-mean Gaussian. evaluate graph quality by using the metrics proposed by [8) TPR (true positive rate), FDR (false positive rate), FPR (false prediction rate), and SHD (structural Hamming distance).
Figure shows that the standard deviation of the Gaussian increases, the learned DAG is more and more different from the ground truth, with the TPR decreasing and the FDR, FPR, and SHD increasing. However, the forecast MSE remains relatively flat (in fact, the MSE from the ground truth slightly higher). Note that with the nonzeros of the ground truth lying between and a Gaussian with standard deviation 0.3 as the initial guess barely carries the signal of the original graph. We have not been able to establish identifiability conditions for A from the ODE as data generation model; and suspect that the conditions, if at all exist, may be unrealistically restrictive, given the

---

empirical findings that a better downstream performance can be achieved by a DAG significantly different from the ground truth. However, even though the ground truth is not recovered, better downstream performance showcases the robust advantage of a graph-based model that intends to capture the complex interplay inside a system.
We also compare the time costs of neural ODE, neural flows, and GNeuralFlow. Table indicates that Flow is more expensive than the corresponding neural flow, because of the additional modeling of the graph. However, GNeuralFlow is more economic than neural ODE, as expected, because does not run numerical solver.
Table 2: Time comparison (in seconds) with neural ODE and neural flows on synthetic systems.
Sink
Triangle
Sawtooth
Square
Neural ODE
1.529
1.527
1.742
2.206
Neural flow (ResNet)
1.022
1.013
1.021
1.020
Neural flow (GRU)
0.251
0.249
0.247
0.247
Neural flow (Coupling)
0.136
0.137
0.136
0.133
GNeuralFlow (ResNet)
1.521
1.521
1.534
1,533
GNeuralFlow (GRU)
0.275
0.283
0.286
0,284
GNeuralFlow (Coupling)
1.215
1.214
1.212
1,213
5.2 Latent Variable Modeling: Smoothing
Just like neural ODE and neural flows, GNeuralFlow can be used for latent variable modeling (see details in Appendix F). To illustrate the effectiveness of GNeuralFlow for this application, we first perform experiments with the smoothing approach in this subsection, by using real-life datasets Activity, Physionet, and MuJoCo [ For Activity, we treat each sensor as a graph node; while for the other two, we treat each feature as a node. The tasks are reconstruct the time series, to classify the activity at each time step (Activity), and to predict the mortality of patients based on the entire time series (Physionet). We again compare our methods with baselines including neural ODE, neural flows, and graph ODEs. For graph ODEs, we construct dynamic graph at each time step by utilizing the covariance matrix the time series data within each batch.
Table Comparison with non-graph neural flows/ODE and graph ODE methods for the smoothing approach. Left two: classification task; right three: reconstruction.
Activity Accuracy
Physionet AUC
Activity MSE (x10
Physionet MSE (x10- -3)
MujoCo MSE (x10 -3)
ODE-RNN
0.785 0.003)
0.781
0.004)
6.050 (#: 0.10)
4.52 (+: 0.03)
2.540 0.12)
No Graph
Neural flow (ResNet)
0.760 (= 0.004)
0.784 (= 0.010)
6.279 (#: 0.09)
4,90 (:: 0.12)
8.403 (#: 0.14)
Neural flow (GRU)
0.783 (= 0.008)
0.788
(= 0.008)
5.837
(:: 0.m)
5.04 (1: 0.13)
4,249 (d: 0.07)
Neural flow (Coupling)
0.752 (= 0.012)
0.788 (= 0.004)
6.579
(d: 0.04)
4.86 (4: 0.07
4.217 (#: 0.14)
Graph ODE
GDE
0.721 (= 0.014)
0.757 (= 0.010)
6.491
(#: 0.011)
4.83 (:: 0.38)
5,220 (:: 0.42)
LG-ODE
0.743 0.023)
0.748
0.018)
5.738
(#: 0.089)
4.87 (4: 0.27)
6.699 (4: 0.83)
CG-ODE
0.768 0.048)
0.783
0.082)
6.241
(# 0.012)
4.73 (+: 0.07)
4.312 (#: 0.17)
Our Method
GNeuralFlow (ResNet)
0.786
0.009)
0,800
0.009)
5.947
(#: 0.03)
4.31 ( 0.06)
2.916 (#: 0.21)
GNeuralFlow (GRU)
0.804 0.003)
0.812
0.001)
5.169 0.05)
4.23
(1: 0.15)
4.112 0.13)
GNeuralFlow (Coupling)
0.808 0.005)
0.808
0.008)
5.431 (1: 0.10)
4.59 (1: 0.23)
3.849 (1: 0.07)
From Table we see that GNeuralFlow generally performs the best compared with the various same These findings are consistent with those the synthetic data case, demonstrating a good utility of our method for real-life applications.
baselines. Moreover, by using the flow design, GNeuralFlow is always better than neural flows.

---

5.3 Latent Variable Modeling:
Filtering
We also perform experiments with the filtering approach. on the MIMIC-IV dataset 4 We treat each feature as a node to set up a graph of 97 longitudinal features, including lab tests, outputs, and prescriptions in clinical events. This graph is the largest among all experiments in this paper.
Table 2 reports the forecast error on the next three time points and the estimated likelihood of the time series. We that GNeuralFlow performs the best with the GRU flow design, while some graph ODE approaches come second. Moreover, GNeuralFlow consistently outperforms neural flows, across flow architecture designs and evaluation metrics. The improvement over corresponding neural flow versions is at least the sum of the standard errors of the two compared methods, which are demonstratively significant.
Table 4: Comparison with non-graph neural flows/ODE and graph ODE methods for the filtering approach. For both metrics, the lower the better.
MSE
NLL
GRU-ODE-Bayes
0.379 (d. 0.005)
0.748 (.L. 0.045)
No Graph
Neural flow (ResNet)
0.379 (1: 0.005)
0.774 (L 0.059)
Neural flow (GRU)
0.364 0.008)
0.734 (:: 0.054)
Neural flow (Coupling)
0,366 0.002)
0.675 (:: 0.003)
Graph ODE
GDE
0.342 0.001)
0.657 (d: 0.007)
L.GODDE
0.349 (#: 0.002)
0.649 (1: 0.005)
CG-ODE
0,372 (#: 0.011)
0,825 (:: 0.018)
Our Method
GNeuralFlow (ResNet)
0.356 (# 0.0007)
0.663 (:: 0.008)
GNeuralFlow (GRU)
0.335 (#: 0.003)
0.606 (:: 0.001)
GNeuralFlow (Coupling)
0.350 0.004)
0.662 (:: 0.008)
Conclusions and Discussions
In this work, we address the challenge of learning the systemic interactions of time series, by proposing a graph-based model GNeuralFlow and learning the graph structure in tandem with the system dynamics. GNeuralFlow is a continuous-time model, which can be used for irregularly sampled time series. Moreover, the systemic interactions are modeled by a conditional dependence structure. We apply GNeuralFlow to latent variable modeling and demonstrate that incorporating the DAG structure improves time series classification and forecasting noticeably.
Several time-series approaches do not learn a static graph but a dynamic one (0)). Such graph is interpreted as a latent structure, which varies depending on past data. In contrast, our approach learns an explicit structure that governs the dynamics over time. Nevertheless, our mathematical framework can be straightforwardly adapted to learning time-varying latent graph, if desired. To achieve so, we reuse the loss calculation in [12 remove the constraint, and parameterize A as a function of X(t). This way, we sacrifice the DAG interpretation of the interactions but gain a time-dependent graph.
A limitation of the proposed model is that the number of parameters on the A part grows quadratically with the number of time series (nodes). Hence, this part of the computational cost be cubic, because the evaluation of the DAG constraint and the gradient involves the computation of the matrix exponential. Such scalability challenge a common problem for DAG structure learning. While past research showeased the feasibility of learning a graph with a few hundred nodes [47), going beyond is generally believed require either computational technique or a new modeling approach. One potential direction to introduce structures into A (such as low-rankness [5)), which admit faster matrix evaluation.
Acknowledgments and Disclosure of Funding
GM acknowledges support from the Engineering and Physical Sciences Research Council (EPSRC) and the BBC under iCASE. AF is partially funded by the CRUK National Biomarker Centre, by the Manchester Experimental Cancer Medicine Centre and the NIHR Manchester Biomedical Research Centre. JC is supported by the MIT-IBM Watson AI Lab.
10

---

References
[1] Jens Behrmann, Will Grathwohl, Ricky T. Q. Chen, David Duvenaud, and Jörn-Henrik Jacobsen. Invertible residual networks. In ICML, 2019.
[2] Dimitri P. Bertsekas. Nonlinear Programming. Athena Scientific, 2nd edition, 1999.
[3] Dhananjay Bhaskar, Sumner Magruder, Edward De Brouwer, Aarthi Venkat, Frederik Wenkel, Guy Wolf, and Smita Krishnaswamy. Inferring dynamic regulatory interaction graphs from time series data with perturbations. In LoG, 2024.
[4] Marin Bilos, Johanna Sommer, Syama Sundar Rangapuram, Tim Januschowski, and Stephan Günnemann. Neural flows: Efficient alternative to neural ODEs. In NeurIPS, 2021.
[5] Edward De Brouwer, Jaak Simm, Adam Arany, and Yves Moreau. BRU-ODE-Bayes Continu- ous modeling of sporadically
time series. NeurlPS, 2019.
[6] Zhengping Che, Sanjay Purushotham, Kyunghyun Cho, David Sontag, and Yan Liu. Recurrent neural networks for multivariate time series with missing values. Scientific Reports, 8(6085), 2018.
[7] Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. In NewrIPS, 2018.
[8] David Maxwell Chickering, David Heckerman, and Christopher Meek. Large-sample learning of Bayesian networks is NP-hard. Journal of Machine Learning Research, 5:1287-1330, 2004.
[9] Kyunghyun Cho, Bart Van Merriënboer, Dzmitry Bahdanau, and Yoshua Bengio. On the properties of neural machine translation: Encoder-decoder approaches. arXiv preprint 1409 1259, 2014.
[10] Jeongwhan Choi, Hwangyong Choi, Jeehyun Hwang, and Noseong Park. Graph neural con- trolled differential equations for traffic forecasting. In AAAI, 2022.
[11] Enyan Dai and Jie Chen. Graph-augmented normalizing flows for anomaly detection of multiple time series. In ICLR, 2022.
[12] Tristan Deleu, Mizu Nishikawa- Toowa-- Jithendaraa Subramanian, Nikolay Malkin, Laurent Charlin, and Yoshua Bengio. Joint Bayesian inference of graphical structure and parameters with single generative flow network. In NeurIPS, 2023.
[13] Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio. Density estimation using real NVP. In ICLR, 2017.
[14] O. Duncan. Introduction to Structural Equation Models. Academic Press, New York, 1975.
[15] Zhuangyan Fang, Shengyu Zhu, Jiji Zhang, Yue Liu, Zhitang Chen, and Yangbo He. On low rank directed acyclic graphs and causal structure learning. TNNLS, 35(4):42244493, 2024.
[16] Morris R. Flynn, Aslan R. Kasimov, Chai-Christophe Nave, Rodolfo R. Rosales, and Benjamin Seibold. Self-sustained nonlinear waves in traffic flow. Physical Review E, 79(5):056113, 2009.
[17] George E. Forsythe, Michael A. Malcolm, and Cleve B. Moler. Computer Methods for Mathe- matical Computations. Prentice Hall, first edition edition, 1977.
[18] Ary Goldberger, Luis A. N. Amaral, Leon Glass, Jeffrey M. Hausdorff, Plamen Ch. Ivanov, Roger G. Mark, George B. Moody Joseph E. Mietus, Chung- Kang Peng, and H. Eugene Stanley, PhysioBank, Toosio, and PhysioNet: components of new research resource for complex physiologic signals. Circulation, (01(2))e211-e222 2000.
[19] Henry Gouk, Eibe Frank, Bernhard Pfahringer, and Michael J. Cree. Regularisation of neural networks by enforcing Lipschitz continuity. Machine Learning, 110:393-416, 2021.
[20] Colin Graber and Alexander G. Schwing. Dynamic neural relational inference. In CVPR, 2020.

---

[21] Paul D. H. Hines, Ian Dobson, and Pooya Rezaei. Cascading power outages propagate locally an influence graph that not the actual grid topology. IEEE Transactions on Power Systems, 32(2):958-966 2017.
[22] Sujai Hiremath, Jacqueline R.M.A. Maasch, Mengxiao Gao, Promit Ghosal, and Kyra Gan. Hybrid top-down global causal discovery with local search for linear and nonlinear additive noise models. Preprint arXiv:2405. 14496, 2024.
[23] Carlos Hoppen, Juan Monsalve, and Vilmar Trevisan. Spectral norm of oriented graphs. Linear Algebra and Applications, 574:167-181, 2019.
[24] Zijie Huang, Yizhou Sun, and Wei Wang. Learning continuous system dynamics from pregularly-smpplee partial observations. In NeurlPS, 2020.
[25] Zijie Huang, Yizhou Sun, and Wei Wang. Coupled graph ODE for learning interacting system dynamics. In SIGKDD, pages 705-715, 2021. 2021.
[26] Song Jiang, Zijie Huang, Xiao Luo, and Yizhou Sun. CF-GODE: Contimuous-time causal inference for multi-agent dynamical systems. In SIGKDD, 2023.
[27] Ming Jin, Yu Zheng, Yuan-Fang Li, Siheng Chen, Bin Yang, and Shirui Pan. Multivariate time series forecasting with dynamic graph neural ODEs. TKDE, 35(9)9168-91989 2023.
[28] Alistair Johnson, Lucas Bulgarelli, Tom Pollard, Steven Horng, Leo Anthony Celi, and Roger Mark. MIMIC-IV (version 1.0). PhysioNet, 2021.
[29] Diederik Kingma and Max Welling, Auto-encoding variational Bayes. ICLR, 2014.
[30] Thomas Kipf, Ethan Fetaya, Kuan-Chieh Wang, Max Welling, and Richard Zemel. Neural relational inference for interacting systems. In ICML, 2018
[31] Thomas N. Kipf and Max Welling, Semi-supervised classification with graph convolutional networks. ICLR, 2017.
[32] John M. Lee. Introduction to Smooth Manifolds. Springer, 2nd edition edition, 2012.
[33] Yaguang Li, Rose Yu, Cyrus Shahabi, and Yan Liu. Diffusion convolutional recurrent neural network: Data-driven traffic forecasting. In ICLR, 2018.
[34] Mizu Nishikawa- Tristan Deleu, Jithendaraa Subramanian, Yoshua Bengio, and Laurent Charlin. Bayesian learning of causal structure and mechanisms with GFlowNets and variational Bayes. Preprint arXiv:221 1.02763, 2022.
[35] George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshminarayanan. Normalizing flows for probabilistic modeling and inference. JMLR, 22(57): 1-64, 2021.
[36] Judea Pearl. Bayesian networks: A model of self-activated memory for evidential reasoning. In Proceedings of the 7th Conference of the Cognitive Science Society, 1985.
[37] Judea Pearl. Causality: Models, Reasoning, and Inference. Cambridge University Press, 2000.
[38] Michael Poli, Stefano Massaroli, Junyoung Park, Atsushi Yamashita, Hajime Asama, and Jinkyoo Park. Graph neural ordinary differential equations. Preprint arXiv: 19171113122,,,,111111111111 2019.
[39] Yulia Rubanova, Ricky T. Q. Chen, and David Duvenaud. Latent ODEs for tiregularlyssamplee time series. In NewrIPS, 2019.
[40] Mona Schirmer, Mazin Eltayeb, Stefan Lessmann, and Maja Rudolph. Modeling irregular time series with continuous recurrent units. In ICML, 2022.
[41] Chao Shang, Jie Chen, and Jinbo Bi. Discrete graph structure learning for forecasting multiple time series. In ICLR, 2021.
12

---

[42] Ikaro Silva, George Moody. Daniel J Scott, Leo A Celi, and Roger G Mark. Predicting in- hospital mortality of ICU patients: The physionet/computing cardiology challenge 2012. Comput Cardiol, 2012(39):245-225 248, 2010.
[43] Stephen Smith and Qing Zhou. Coordinated multi-neighborhhodd learning directed acyclic graph. Preprint arXiv:2405.15356 2024.
[44] Yuval Tassa, Yotam Doron, Alistair Muldal, Tom Erez, Yazhe Li, Diego de Las Casas, David Budden, Abbas Abdolmaleki, Josh Merel, Andrew Lefrancq. Timothy Lillicrap, and Martin Riedmiller. DeepMind control suite. Preprint arXiv: 1111180006000 2018.
[45] S. Wright. Correlation and causation. Journal of Agricultural Research, 20:557-585, 1921.
[46] Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and Philip S. Yu. A comprehensive survey on graph neural networks. IEEE Transactions on Neural Networks and Learning Systems, 32(1):4-24 2021.
[47] Yue Yu, Jie Chen, Tian Gao, and Mo Yu. DAG-GNN: DAG structure learning with graph neural networks. ICML, 2019.
[48] Xun Zheng, Bryon Aragam, Pradeep Ravikumar, and Eric P. Xing. DAGs with NO TEARS: Continuous optimization for structure learning. In NearIPS, 2018.
[49] Jie Zhou, Ganqu Cui, Zhengyan Zhang, Cheng Yang, Zhiyuan Liu, Lifeng Wang, Changcheng Li, and Maosong Sun. Graph neural networks: A review of methods and applications. Al Open,
1 -157-81,20 2020.
13

---

A
Supporting Code
Code is available attcppe///ghabb..oo///aaae/aaele
B
Proofs
Proof of Theorem 7 When A DAG adjacency matrix, its symmetrization B is the adjacency
matrix of the corresponding undirected graph. We call the DAG an orientation of the undirected graph. The Gershgorin circle theorem asserts that the spectral radius of B, p(B). is bounded by ?
Meanwhile, 23 show that the spectral norm
| S (BB). Then, || 12 S 1 + 7/7 = 2. L
Proof of Theorem 2 We follow the proof 4 Theorem 1] (see A.3 the paper), which concludes
that [h(x;) hyy) S { )) y/- Since GCN is contractive, such inequality applies to both h h² and A²2 Then, applying Eqn (14) of the paper,
| (x) (w))
i (x))
Lip(h²)
| (x:) Lip(h") |
2a x-
<a( <a(2)+))
<a( (++))+))
Therefore, when a53 6) S the product h h/2 contractive and therefore F is invertible.
C Details of Example in Section 3
The matrix is
4 5
3 1
and the initial conditions are
0.6 0.5
0.7
0.2
x
x
x
0.1
0.3
D Training Method
In this section, briefly describe the augmented Lagrangian method for solving the equality- constrained problem
min A.
C(A, 0)
M(A) -,
where the constraint can either be ((A) tr(expm(A A -)) =((AAAAA))))
Define the augmented Lagrangian
- C(A, AM(A)+ + (1
(13)
where À and e denote the Lagrange multiplier and the penalty parameter, respectively. The general idea of the method is to gradually increase the penalty parameter to ensure that the constraint is eventually satisfied. Over iterations, as a dual variable will converge to the Lagrange multiplier of the original problem. The upate rule at the kth iteration reads
A 0h argmin Let. A,O
A k- = xt + (h
A- n
if (((AA*) (((A*=)))
ck
else,
where 77 (1, +00) and y E (0, 1) are hyperparameters to be tuned.
14

---

The subproblem of optimizing A and & can be solved by using the Adam optimizer. It requires the gradient of L. and hence of h. For (AA) tr(expm (A A)) n, it can be derived that Vh expm A @A) 2A, which can be obtained virtually for free after h has been evaluated. For h(A) = tr aA & A)") n, one may use automatic differentiation to obtain the gradient.
Algorithm T summarizes the training procedure. Note that an effective initialization of A would use an empty diagonal. Moreover, in every update of A, one may keep its diagonal zero throughout.
Algorithm 1 Training algorithm of GNeuralFlow
1: Initialize e <- and
2: for k= 0, 2, .......6 do
3: Compute A & and 0 minimizer of 13 by using the Adam optimizer
4: Update Lagrange multiplier À + X + AA(A
5:
i > and (M(A*) )) > 7 ( (-1) then
6:
C ne
7:
end if
8: if (T((A*)] threshold then break
9:
10: end if
11:
end for
E Handling Missing Data
At a particular time t, some rows of the observed data X(t) may be empty (i.e., measurements of some time series at time are missing). In this case, one evaluates the graph neural flow subgraph of present measurements. Rather than extracting this subgraph and the corresponding rows of X, the GCN encoder (6) offers convenient approach for evaluation: masking. This approach is particularly favorable in batching, because tensor dimensions do not change. In particular, we mask out (i.e., setting zero) the part of A corresponding to missing data. Then, for the output X, the part corresponding to present data is correctly calculated, while the part of missing data becomes zero and this condition is invariant across layers, Note that the GCN layers must not have bias terms to maintain this invariance.
F GNeuralFlow for Latent Variable Modeling
In addition to straightforwardly modeling the data space, neural flows find successful use in the latent space. Here, we discuss two popular approaches in latent variable modeling and how GNeuralFlow can be incorporated.
Of particular consideration is the role the graph. One may straightforwardly extend 4 by using the flow in the latent/hidden space; however, this method models the interactions among the latent/hidden dimensions, which are less interpretable than those among the time series. Hence, we propose to use the flow on augmented space, part of which carries the graph information, as a new design complementary those proposed in Section
On a high level, smoothing [39 and filtering 5 approaches use a neural ODE or neural flow to continuously evolve the hidden state from time f--1 (denoted as Htt)-1)) to time tú (denoted as Ht,t;)); and then use an RNN to introduce a jump (denoted as H )) on observing input data X(t,).
To model and learn the interaction graph in the data space, we use the graph encoder (6) to produce a transformed data X(t,) and use a second RNN to introduce the paired jump H (t,) given (ty)).
Then, the pair of hidden states, H(t,) and (ty)) are concatenated and a standard neural flow evolve the concatenated state to the next time point, the result of which then projected to the proper hidden dimension. The smoothing and filtering approaches differ in fine details, including different uses of the RNNs and hidden states. As a result, the loss function L in (12) also different. Details are presented in the following.
15

---

F1
Smoothing Approach
Given observation data ((to)....... X(tx) this approach produces latent quantity Zo by combined
use of LSTM and neural flow, and then traces out smooth curve Z(t) using another flow, taking ((tu)
Zu as the initial condition. Then, the observation data ((y)) recovered from Z(t)).
VAE used to set up the training loss. The decoder pl (lo), ((xx) (Zo) factorized as
N
p(X(t)). ((( = II p ((t)))
==
where each Z(t;) is computed by running standard neural flow: Z(ty) F(t,, Zo). The encoder, on the other hand, produces a latent Gaussian Zo with mean μ and diagonal covariance diag(o)); ); that is,
(Z ((to)
((tn)) = (Z μr, diag(or)) log oj = 0 ((H(t)))
where H(N is the hidden state to be elaborated soon and g is neural network projection The ELBO loss for the VAE is
LELOO = DxL (2 (Zoo)((t)). ((tx)))lp(Zu)
Ez x((.)....X((x))
..., (oo)
Our GNeuralFlow uses a pair of LSTMs together with another neural flow evolve the hidden state. Specifically, we maintain pair of states H(t) and Htt., the latter which includes the graph information. At time 11- we concatenate H and H(tz-1), run the neural flow F to evolve the concatenated state time t,, and apply projection 9proj so that the net result H' (t(y) remains in the same dimension as H(t
H' (4,))) Oproj( Fly,
Then, we run a pair of LSTMs to obtain the paired states at time ty:
H(tt) LSTM' (+), (t)),
x(t))), H(t,) LSTM² (H')), X(t)).
where the second LSTM applied to the transformed observation data x()) produced by the GCN encoder (6). By doing so, the graph models the interaction inside the data x(() rather than the hidden
states H
F.2 Filtering Approach
As opposed to the preceding approach, the filtering approach uses only a decoder. Each time, it first evolves the hidden state to Ht, and then runs GRU to update the hidden state H(t,). This approach maintains two Gaussians, the first one models the observation X(t,):
(tx(t,) Ma diag(o))s)) [Myes log o 9abs HH (t,)).
while the second one models the jump caused by the GRU:
N (H (iag()((a))) [Mppaa ppuaa] Apost (H())))
The training loss aims maximizing the observation data likelihood while minimizing the KL divergence of the two Gaussians:
- Σ (X Wheee diag(o () (NiMi ()) (M diag( oppoaa))
==1
Our GNeuralFlow uses pair of GRUs together with a standard neural flow to evolve the hidden state. Specifically, we maintain pair of states Htt) and Htt), the latter of which includes the graph
When the encoder is run backward in time, Htoo) instead of Htto).
16

---

information. At time 1--1, we concatenate H(t))1) and H(tz-1), run the neural flow F to evolve the concatenated state time t,, and apply projection 9pooj so that the result (ty) remains in the same dimension as H(
H (ty) = Oproj Flty,
Then, we run a pair of GRUs to obtain the paired states at time lj:
H(ty) GRU
GRU (H'(t)), 9pee Ht(t), H' (tz), Sprep (t()),, ((,,,, (t))) H(t,)
where the second GRU is applied to the transformed observation data X(t,) produced by the GCN encoder (6). By doing so, the graph models the interaction inside the data (() rather than the hidden states H t))
G Datasets and Tasks
Table 5: Experiment settings and datasets.
Dataset
Method
Tasks & Metrics
#Nodes (v)
#Times (N)
#Samples
Split
Synthetic
regression
forecast MSE graph metries
5-30
500
1000
60:20:20
Synthetic
regression
15
500
1000
60:20:20
Activity
smoothing
reconstruction MSE
4
50
6554
75:5:20
classification accuracy
Physionet
smoothing
reconstruction MSE
41
8000
60:20:20
classification AUC
MuJoCo
smoothing
reconstruction MSE
14
100
10000
60:20:20
MIMIC-IV
filtering
forecast MSE
97
19
17874
70:15:15
log-likelihood
The datasets used in this paper include four synthetic ODE systems and four real-life datasets. Table summarizes the basic information of these datasets, tasks, evaluation metrics, and learning methods.
G.1
Synthetic Datasets
We define four interacting systems based on either the ODE X = f(t,X,A) or the solution ((t) Ft,, Xo.
A):
Sink (2D): f(t, X, A)
A A)X where B [= 10
2
Triangle (ID): F(t, X, A) = (I - A AT) (X + f/¹ sign(sin(()) du)
Sawtooth (ID): (t,, X, X) A)=
-I (I- A T)(X xx [t])
Square (ID): F(t, X, A) = =1 A X sign(sin(t))))
For Sink, X Rax2 while for the other three systems, Rax1, where n is the number of trajectories (graph nodes) in the system. The initial condition X, is uniformly sampled from [0, 1|ax2 for Sink, and from [-2, 2 nx1 for the other three systems. The time interval is [0, 10] and the time points are uniformly random.
The DAG adjacency matrix is generated by using the following procedure:
1. Generate sparse matrix A with pre-defined density, where the nonzero locations are random and the nonzero values are uniformly random.
2. Keep only the strict upper triangular part of A (i.e., diagonal is zero).
3. Perform symmetric row/column permutation on A.
The task is to predict the trajectories X(t) given Xo.

---

G.2 Real-Life Datasets
We four real-life datasets preprocessed by [
Activity [39] contains time series recorded by four sensors, on individuals performing various activities: walking, sitting, lying, etc. The task is to classify the activities at each time point. Additionally, since the smoothing approach for latent variable modeling reconstructs the time series, also evaluate different models the reconstruction quality. We treat each sensor node.
graph
Physionet 42 contains time series of patients' measurements (37 variables in total) from the first 48 hours after being admitted to ICU. The task is to predict the mortality of the patients. Additionally, since the smoothing approach for latent variable modeling reconstructs the time series, also evaluate different models on the reconstruction quality. We treat each variable as one graph node.
MuJoCo 44 contains physics simulations by randomly sampling initial positions and velocities and letting the dynamics evolve deterministically in time. Each sequence includes 14 features. We treat each feature as one graph node. We evaluate different models on the reconstruction quality.
MIMIC-IV [18 28 contains time series of ICU patients' measurements, including their vital signs, laboratory test results, medication, and any output data during their ICU stay (97 variables in total). The task is predict the three measurements in the 12 hour interval after the observation window of 36 hours. We treat each variable as one graph node.
H Hyperparameter Details
Table 6: Graph learning hyperparameters.
Synthetic systems
(all architectures)
Real-life datascts
# points
ResNet
GRU
Coupling
3
3
0.3
"
?
n
Y
1
5
0.25
Activity
0.21
0.21
7
0.21
15
0.21
Physionet
10
0.5
10
0.5
20
7
0.19
MuJaCo
15
0.5
10
0.5
15
0.5
25
7
0.19
MIMIC-IV
10
0.15
10
0.15
10
0.15
30
0.16
We reuse the architecture parameters and training hyperparameters 4 and only tune the graph learning hyperparameters (sce Algorithm| I in Section D). Table 6 lists the tuned values.
18

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The contributions of the paper summarized in the introduction section and elaborated the following sections.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification:
The limitations of the work discussed in the concluding section.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
the paper. The authors should use their best
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
19

---

Justification:
Assumptions are given in the theorems and proofs are given in the appendix.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: Experiment information provided in part in the main text and in part in the appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm. (b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome deseribe the particular way they provide for reproducibility, In the case of closed-source models, it may be that access to the model is limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
20

---

Answer: [Yes]
Justification:
Code will be released upon publication of the paper.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips. cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
/
wwt.Suubmisiooooiiiooooooooo for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data,
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: Experiment settings and details are provided in part in the main text and in part in the appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: Error bars are obtained by performing five repetitive runs for each method and dataset.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
21

---

It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative еπхо rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: Compute information is given in the main text.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics
Answer: [Yes]
Justification: We confirm.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of work, none which we feel must be specifically highlighted here.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
22

---

The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: The paper poses no such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g.. code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: Codes and datasets used for experiments are publicly available under permis- sive licenses.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the Datasets.papersstt paperswi thcode. cww/ datasets package should be provided. For popular datasets,
has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
23

---

If this information is not available online, the authors are encouraged to reach out to the asset's 's creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: Code implementation of the proposed method will be released upon publication. Comprehensive documentation will be provided for reprodueibility and access.
Guidelines:
The answer NA means that the paper does release new assets.
Researchers should communicate the details the dataset/code/moodll as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors to adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.


Reviews:
[{"review_id": "YSTIgXGecq", "reviewer": "Reviewer_ufSP", "summary": "The paper proposes a continuous-time model to discover the causal structure from irregular multivariate time series data. The idea is to introduce the DAG (directed acyclic graph) to model the interaction of different time series at the same time step in the vector field that generates the multivariate time series（that is, the $n$ trajectories). Instead of learning the vector field, the paper proposes to learn the solution of the ODE system directly with neural flows. The structural interaction is implemented with graph convolution neural networks and three particular parameterization methods are proposed for the solution function $F$ to guarantee its invertibleness. The experiments on both synthetic and real-world irregular time series datasets are conducted to assess the efficacy of the proposed method for both prediction and classification tasks.", "questions": "See Weaknesses.", "limitations": "No potential negative societal impact.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "* The proposed method is well motivated and clearly illustrated, and the presentation is easy to follow. The reviewer appreciates the way to establish the solution framework presented in Section 4.1, which is very well explained.\n* It is interesting to introduce the DAG structure learning to irregular time series modeling and the proposed method seems reasonable to the reviewer. \n* The proposed method achieves better results than the existing method.", "weaknesses": "* In Section 3.2, it states \"(i) the vector fields $\\mathbf{B}\\mathbf{x}^1, \\ldots, \\mathbf{B}\\mathbf{x}^n$ follow the same conditional dependence structure governed by $\\mathbf{A}$\",  but according to Eq. (2), the vector fields seems to be $\\mathbf{B}\\mathbf{x}^j - \\sum_{i=1}^{j-1}a_{ij} \\mathbf{Bx}^{i}$ for $j=1,\\ldots,n$, so I would suggest changing the description to \"the vector fields that generate the $n$ trajectories follow the same conditional dependence structure governed by $\\mathbf{A}$\" to avoid confusion.\n* To learn the DAG matrix, the proposed method employs NOTEARS and DAG-GNN. How do these DAG learning algorithms impact the model performance?\n\n* The proposed method (graph neural flows) is dubbed GNeuralFlow in the main text but is referred to as GNF in the tables of experiments. It is better to unify the abbreviation. Besides, no code is released.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**RE: ..I suggest changing the description to \"the vector fields..**\n\nThanks for the feedback. We will update the text.\n\n\n**RE: .. How do these DAG learning algorithms impact the model performance**\n\nThanks for pointing out this part. We address it as follows.\n\nA limitation of the proposed model is that the number of parameters on the $A$ part grows quadratically with the number of time series (nodes). Such a scalability challenge is a common problem for DAG structure learning. While past research showcased the feasibility of learning a graph with a few hundred nodes (Yu et al 2019), going beyond is generally believed to require a new computational technique.\n\n\n**References**\n\nYu et al 2019. DAG-GNN: DAG structure learning with graph neural networks. ICML\n\n\n**RE: unify the abbreviation GNF**\n\nGNF is commonly an abbreviation for Graph Normalizing Flow. We will change the abbreviation.\n\n\n**Final comment**\n\nThank you for the comments. Please consider raising the score if we have addressed them appropriately.", "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": "", "extra_scores": {}}, {"review_id": "t4hbNbxssy", "reviewer": "Reviewer_4eJC", "summary": "This paper addresses the problem of multivariate time series prediction with feature interactions. It proposes learning a Directed Acyclic Graph (DAG) to model interactions, encoded by a Graph Neural Network (GNN), and using a neural flow to model the dynamics. The experiments demonstrate improvements over graph-free models and models with given graphs obtained from covariance.", "questions": "For the experiments on synthetic datasets (Figure 1), how are the graphs for graph ODEs obtained? Are all models using the ground truth graph? If so, how is the conclusion that “DAG is better than other graph structures at modeling dependencies” (line 252) drawn?\n\nDoes the performance gain over existing graph ODE models come from the learnable graph? Is there an ablation study for that?\n\nIs the model learning a constant graph independent of time? How does it perform when there is a significant change in data dependence?", "limitations": "The authors have adequately addressed the limitations.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The combination of graph learning and neural flow, along with DAG-conditioned optimization, is original. Since flows are popular in generative modeling, the idea of converting them to the graph setting seems promising. The paper also learns graph parameters using the matrix A to learn (or relearn the DAG structure), which can be useful in cases where the graph is not previously known with certainty.", "weaknesses": "Weaknesses:\nFlow matching inherently has the idea of using a linear interpolant (straight line) to encode optimal transport paths in time series data. In spaces where this is not true, or where the optimal transport paths are perhaps in a latent dimension the assumptions are violated. This is why Flow matching has been more useful for generative modeling than time series. This is not noticed by the authors and the restrictive experiments they perform are insufficient to test this. \n\nI don't see the need for a DAG here. See the RITINI model [Bhaskar et al, LOG 2023], when there is time series data there can be feedback loops in the variables rather than DAG-style causality. I would appreciate trying to loosen this requirement, which would enable cyclic dependencies. \n\nThe set of comparisons featured here is rather limited. Here are several works that could be compared against: \n\nKipf et al., Neural Relational Inference for Interacting Systems, ICML, 2018.\nGraber and Schwing. Dynamic Neural Relational Inference, CVPR, 2020 (dynamic graph inference for discrete-time systems)\nNishikawa-Toomey et al., Bayesian learning of Causal Structure and Mechanisms with GFlowNets and Variational Bayes (jointly learn the DAG and the parameters of a linear-Gaussian causal model)\nDeleu et al., Joint Bayesian Inference of Graphical Structure and Parameters with a Single Generative Flow Network, NeurIPS, 2023 (jointly learn the posterior over the structure of a Bayesian Network, and also the parameters of its conditional probability distributions)\nSmith and Zhou, Coordinated Multi-Neighborhood Learning on a Directed Acyclic Graph, arXiv:2405.15358 (constraint-based approach that exploits conditional independence relationships between variables to infer the underlying causal model)\nHiremath et al. Hybrid Global Causal Discovery with Local Search, arXiv:2405.14496 (global causal discovery by learning local causal substructures using topological sorting and pruning of spurious edges)\nBhaskar et al. Inferring dynamic regulatory interaction graphs from time series data with perturbations, LOG 2023", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**RE: flow matching (W1)**\n\nOur effort consists in extending Neural Flows (Bilos et al 2021), to integrate additional interdependency information in the form of a learned graph. Note that (Bilos et al 2021)  was proposed in the context of time series, differently from Flow Matching. We are not proposing a flow-matching approach. \n\n\n**RE: no need for a DAG here (W2)**\n\nThanks for pointing out this perspective. Our view on the matter is the following.\n\nA DAG provides a clear and interpretable representation of causal relationships between variables. In the context of time series data from multiple sensors, understanding these causal relationships can help identify how changes in one sensor's readings might influence others. This is crucial in medical monitoring, where understanding causality can aid in diagnosing conditions or understanding physiological responses.\n\nIn addition, DAG structure learning is a long-lasting and challenging problem in probabilistic graphical models (in particular, Bayesian networks). The structure of a Bayesian network takes the form of a DAG and plays a vital role in causal inference (Pearl, 1985). Learning the DAG structure is a combinatorial problem; it is not only theoretically intimidating (NP-hard; see Chickering et al. (2004)) but also practically challenging even when approximate solutions are sought.\n\n**References:**\n\nBilos et al (2021).  Neural flows: Efficient alternative to neural ODEs. Neurips.\n\nPearl (1985). Bayesian networks: A model of self-activated memory for evidential reasoning. Proceedings of the 7th Conference of the Cognitive Science Society.\n\nChickering et al. (2004). Large-sample learning of Bayesian networks is NP-hard. JMLR.\n\n\n\n**RE: related works (W3)**\n\nThanks for providing these related works, we will place them appropriately in the main text. \n\nHere we add the comparison with two of the requested baselines. Specifically, we provide a comparison with NRI (Kipf et al., 2018) and dNRI (Graber and Schwing 2020) methods on the synthetic datasets.    \n\n\n| model      | Sink            | Triangle        | Sawtooth        | Square          |\n| ---------- | --------------- | --------------- | --------------- | --------------- |\n| NRI        | 5.25(±0.02)     | 3.96(±0.16)     | 4.99(±0.12)     | 9.39 (±0.45)    |\n| dNRI       | 5.40(±0.04)     | 3.39(±0.09)     | 4.97(±0.21)     | 9.78(±0.21)     |\n| GNF-resnet | **3.95(±0.32)** | **2.32(±0.11)** | **3.84(±0.06)** | **8.24(±0.64)** |\n\nWe notice that three of the requested comparisons were posted after the conference deadline, namely, Nishikawa-Toomey, (2024), Smith and Zhou (2024), and Hiremath et al. (2024).\n\n\n  \n**References:**\n\nKipf et al. (2018). Neural Relational Inference for Interacting Systems, ICML.   \n  \nGraber and Schwing. (2020). Dynamic Neural Relational Inference, CVPR\n\nHiremath et al. (2024). Hybrid Global Causal Discovery with Local Search, arXiv\n\nSmith and Zhou (2024). Coordinated Multi-Neighborhood Learning on a Directed Acyclic Graph. Arxiv.\n\nNishikawa-Toomey et al. (2024). Bayesian learning of Causal Structure and Mechanisms with GFlowNets and Variational Bayes. Arxiv.\n\n\n**RE: how are the graphs for graph ODEs obtained, and why DAG is better than other graph structures at modeling dependencies (Q1)**\n\nThanks for raising this point. The graph ODEs are given the ground truth graph in input for the synthetic data, while they are given the covariance matrix in the real-world datasets. \n\nEmpirically, we find that when we employ our DAG method, we achieve improvements over graph-ODEs, even while they are given the ground truth graph. \n\n\n**RE: performance gain over existing graph ODE (Q2)**\n\nThanks for raising this point. The performance gain is due to the DAG component. For example, in Figure 2, we show that a NeuralFlow model without graph dependencies performs worse than both GNeuralFlow with a learned graph and with the ground truth graph.\n\n\n**RE: learning a constant graph independent of time (Q3)**\n\nWhile we learn a constant graph, and certain datasets may not work under such an assumption, in the experiments we achieved success also in real-life datasets. Moreover, our framework can support learning a dynamical graph by parameterizing it with time-dependent node embeddings.", "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": "", "extra_scores": {}}, {"review_id": "Vm4M6tfnsH", "reviewer": "Reviewer_3M7x", "summary": "The paper focuses on a graph-based model to capture dependencies in irregularly sampled time series data. The framework employs a causal prior—a directed acyclic graph—where nodes are conditionally independent of non-descendants given their parents, specifying component dynamics dependencies. The proposed model, termed a graph neural flow, learns the solution of an unknown ODE directly from irregularly sampled time series, which contrasts with neural ODEs by avoiding repeated calls to a costly numerical solver. Multiple neural flows, one for each time series, are conditioned on the DAG, with their interactions instantiated using a GNN, such as a graph convolutional network (GCN). The GCN enhances the ODE solution parameterization by aggregating neighboring time series information at each time point, effectively modeling a graph-conditioned ODE that captures system interactions.", "questions": "Based on the *weaknesses* above please focus on the following aspects:\n1. **Experiments:** Explain choices of datasets and baselines or extend the experimental evaluation (W1). Conduct qualitative evaluations of graph-based causal representations (W2).\n2. **Limitations:** Perform a comprehensive computational cost analysis comparing GNeuralFlow with simpler and more complex methods to highlight efficiency gains (W3).\n3. **Contribution:** Please better position the presented contribution among relevant works (W4).", "limitations": "It would be more complete to experimentally showcase the computational limitations, demonstrating in practice the scalability of the proposed method for real-world datasets.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The proposed model is nicely presented, along with precise and meaningful mathematical formulations for capturing interactions in irregularly sampled time series. GNeuralFlow shows significant performance in both classification and regression problems spanning several synthetic and real-world datasets. Hyperparameters and models’ configurations are presented enabling reproducibility and fair comparisons.", "weaknesses": "1. *W1:* Except for the standard experiments on synthetic data representing dynamical systems, the authors explore different tasks on various datasets, slightly deviating from some setups followed in existing related works. For instance, some studies include additional experiments on interpolation/extrapolation [1] or use additional datasets (e.g., classification on MIMIC-III). Moreover, incorporating additional standard baselines in addition to the presented ones could strengthen the evaluation (recurrent/attention-based ones, e.g., mTAND, GRU-D, and others).\n2. *W2:* The proposed GNeuralFlow is primarily based on the concept of Neural Flows for efficient computation of ODE system solutions. The significant methodological extensions in the paper include the incorporation of graph-based representations and causal formulation. However, these extensions are not qualitatively evaluated by analyzing causal relationships or visualizing interactions within the learned graphs. For instance, authors in [2] present studies/experiments to quantify causal discovery in the studied datasets.\t\t\n3. *W3:* Although the approach of replacing the ODE solver offers a computational advantage over classical Neural-ODE-based methods, a comprehensive computational cost analysis comparing different methods, including simpler (e.g., sequential, non-ODE-based) yet effective models, is crucial.\n4. *W4:* The presentation of the related work is quite brief, causing some confusion, particularly regarding closely related graph-based methods like CF-GODE and LG-ODE. Including comparisons and comments on additional relevant papers could be beneficial [3,4] (please explain if not relevant).\n\n[1] Schirmer, M., Eltayeb, M., Lessmann, S., & Rudolph, M. (2022, June). Modeling irregular time series with continuous recurrent units. In International conference on machine learning (pp. 19388-19405). PMLR.\n\n[2] Löwe, S., Madras, D., Zemel, R., & Welling, M. (2022, June). Amortized causal discovery: Learning to infer causal graphs from time-series data. In Conference on Causal Learning and Reasoning (pp. 509-525). PMLR.\n\n[3] Choi, J., Choi, H., Hwang, J., & Park, N. (2022, June). Graph neural controlled differential equations for traffic forecasting. In Proceedings of the AAAI conference on artificial intelligence (Vol. 36, No. 6, pp. 6367-6374).\n\n[4] Jin, M., Zheng, Y., Li, Y. F., Chen, S., Yang, B., & Pan, S. (2022). Multivariate time series forecasting with dynamic graph neural odes. IEEE Transactions on Knowledge and Data Engineering, 35(9), 9168-9180.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**RE: additional experiments (Q1)**\n\nThanks for your comment. \nHere we provide additional baselines. Specifically, we provide the comparison with GRU-D (Che et al 2016), NRI (Kipf et al., 2018), and dNRI (Graber and Schwing 2020) on the synthetic datasets.\n\nWe observe that our GNF-gru can improve over both Neural-flow with GRU and GRU-D. While GNF-gru is outperformed by NRI and dNRI on most datasets, our GNR-resnet achieves the best performance overall.\n\n| model             | Sink            | Triangle        | Sawtooth        | Square          |\n| ----------------- | --------------- | --------------- | --------------- | --------------- |\n| NRI               | 5.25(±0.02)     | 3.96(±0.16)     | 4.99(±0.12)     | 9.39 (±0.45)    |\n| dNRI              | 5.40(±0.04)     | 3.39(±0.09)     | 4.97(±0.21)     | 9.78(±0.21)     |\n| NF-GRU            | 10.9(±0.43)     | 10.3(±0.45)     | 16.1(±0.41)     | 17.2(±0.51)     |\n| GRU-D             | 12.31(±0.23)    | 11.25(±0.32)    | 17.55(±0.53)    | 18.73(±0.31)    |\n| GNF-resnet (ours) | **3.95(±0.32)** | **2.32(±0.11)** | **3.84(±0.06)** | **8.24(±0.64)** |\n| GNF-gru (ours)    | 6.83(±0.23)     | 5.41(±0.23)     | 5.11(±0.13)     | 9.14(±0.61)     |\n\n \n\n\n**RE: computational cost (Q2)**\n\nThanks for the feedback. We present the training cost as the mean training time for one epoch (in seconds) over synthetic datasets. As expected, our graph neural flow (GNF) models are more expensive than the neural flow corresponding ones (because of the additional modeling of the graph). However, the GNF are still cheaper than the Neural ODEs (because we do not need a numerical solver).\n\n\n| Model               | Sink  | Triangle | Sawtooth | Square |\n| ------------------- | ----- | -------- | -------- | ------ |\n| Neural-ODE          | 1.529 | 1.527    | 1.742    | 2.206  |\n| NF-resnet           | 1.022 | 1.013    | 1.021    | 1.020  |\n| NF-coupling         | 0.136 | 0.137    | 0.136    | 0.133  |\n| NF-gru              | 0.251 | 0.249    | 0.247    | 0.247  |\n| GNF-resnet (ours)   | 1.521 | 1.521    | 1.534    | 1.533  |\n| GNF-coupling (ours) | 1.215 | 1.214    | 1.212    | 1.213  |\n| GNF-gru (ours)      | 0.275 | 0.283    | 0.286    | 0.284  |\n\n\n**RE: position the presented contribution among relevant works (Q3)**\n\nThanks for pointing out the related works from Choi et al (2022) and Jin et al (2022). We will add them to the main text of our paper. There are some differences with our method.\n\nFirst, these papers take the Neural ODE approach, whereas our method follows neural flows. While both aim at solving an unknown ODE, neural flows are computationally more economical than Neural ODE, because they do not require the repeated calls of ODESolve. It is unclear if the referenced papers can straightforwardly swap Neural ODE with neural flows for this benefit, without changing other components of their models.\n\nSecondly, these papers consider time-varying graphs. Here, we are not to argue if a time-varying graph or a constant graph is superior; rather, they come from different modeling beliefs. Time-varying graph modeling typically constructs the graph from data (e.g., an affinity graph of given feature vectors at a time, or a co-occurrence graph of observations within a sliding time window) or parameterizes the graph based on node embeddings; whereas DAG structure learning treats the graph as a free parameter to learn. Our method models a constant graph that is assumed to generate the data over time.\n\n\n**References:** \n\nChoi et al (2022). Graph neural controlled differential equations for traffic forecasting. AAAI.\n\nJin et al (2022). Multivariate time series forecasting with dynamic graph neural odes. IEEE Transactions on Knowledge and Data Engineering.\n\nChe et al (2016). Recurrent Neural Networks for Multivariate Time Series with Missing Values.\n\nKipf et al. (2018). Neural Relational Inference for Interacting Systems, ICML.   \n  \nGraber and Schwing. (2020). Dynamic Neural Relational Inference, CVPR\n\n\n**Final comment**\n\nThank you for the feedback. Please consider raising the score if we have addressed your concerns.", "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": "", "extra_scores": {}}, {"review_id": "vQkHcfmkOI", "reviewer": "Reviewer_zrVL", "summary": "This paper proposes a novel graph-based continuous-time model GNeuralFlow for learning systemic interactions and interactions are modeled as a Bayesian network. Experimental results for regression problems show that the proposed GNF achieves state-of-the-art performances on time series classification and forecasting benchmarks.", "questions": "Please see comments in Weaknesses.", "limitations": "Not applicable.", "rating": 3, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "+ The paper is well-written and easy to follow. \n+ The introduction give a nice overview and motivation for the problem.\n+ The results are good and evaluation is reasonable.", "weaknesses": "- The proposed method lacks a theoretical analysis.\n- The authors can consider using state-of-the-art GNN as backbone for graph encoder.\n- Can the authors apply this method to real-world datasets?\n- Computational complexity and corresponding comparisons with baselines are missing.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**RE: theoretical analysis**\n\nOur work mainly focuses on modeling. We offered certain analyses to justify some modeling choices, including guaranteeing contractive mappings required by the neural flows (Theorem 1, Theorem 2, and other inline text within section 4.3).\n\n**RE: sota GNN**\n\nThank you for pointing out this possible direction. We did not explore further as the network we used was already able to outperform the considered baselines.\n\n\n**RE: real-word dataset**\n\nThanks for this request. We note that we already performed experiments over four real-life datasets, including Activity, Physionet, Mujoco, and Mimic-IV. See our paper in sections 5.2 and 5.3, and Appendix F for dataset description.\n\n\n\n**RE: computational complexity**\n\nThanks for this point. We have added a table with the wall-clock runtime. As expected, the graph-neural-flow (GNF) models are more expensive than the neural-flow corresponding ones (because of the additional modeling of the graph). However, the GNF are still cheaper than the Neural ODEs (because we do not need a numerical solver).\n\n\n| Model               | Sink  | Triangle | Sawtooth | Square |\n| ------------------- | ----- | -------- | -------- | ------ |\n| Neural-ODE          | 1.529 | 1.527    | 1.742    | 2.206  |\n| NF-resnet           | 1.022 | 1.013    | 1.021    | 1.020  |\n| NF-coupling         | 0.136 | 0.137    | 0.136    | 0.133  |\n| NF-gru              | 0.251 | 0.249    | 0.247    | 0.247  |\n| GNF-resnet (ours)   | 1.521 | 1.521    | 1.534    | 1.533  |\n| GNF-coupling (ours) | 1.215 | 1.214    | 1.212    | 1.213  |\n| GNF-gru (ours)      | 0.275 | 0.283    | 0.286    | 0.284  |\n\n\n**Final comment**\n\nThank you for your review. Please consider raising the score if we have addressed your points appropriately.", "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": "", "extra_scores": {}}, {"review_id": "dm2rfvTi8K", "reviewer": "Reviewer_2pqA", "summary": "Authors propose a novel idea to learn the interactions between time series which is crucial to make reliable forecasts in interacting systems. Paper proposes a new algorithm GNeuralFlow; a new graph-based continuous time model to learn systematic interaction. A GNN is used to learn the interactions between time series, where each time series is associated with a directed acyclic graph. GNeuralFlow learns the graph structure for irregularly sampled time series.  \n\nThe experiments performed on synthetic and real data shows they GNeuralFlow outperform the other models based regression and classification tasks. The experiments to model latent variable modelling shows GNeuralFlow's robustness in learning and utilizing the underlying graph structure for improved predictive performance.", "questions": "- The paper discusses about the increase in the complexity due to added Graphs in the model. Is there an estimate how much is the computation time increase (is it quadratic?) ? Does this time power increase (from quadratic to cubic) as # of interacting time series are increased in the experiment?", "limitations": "The paper is well structured but the biggest strength becomes its limitation as well. The addition of graph and learning the parameters allow the measurement of interactions but at the same time it increase the complexity of the model. This would be a the main challenge to use the model in real applications and might need approximate solutions to adapt this version of paper.", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "Originality:  Authors identify an important problem in non-uniform time samples time series to model the interactions between them. The paper introduces the original idea (as per my knowledge) of using Graph Neural Flow (GNeuralFlow) which leverages the power of graph based modelling and continuous time technique to capture the relationship between time series.   The method has wide application in climate, finance and more where a dynamic system is crucial. \n\nQuality:   The main development of paper GNeuralFlow is built on a strong mathematical foundation and authors provide empirical and experimental results to support their claims. \n\nClarity: The paper is well structured and describes the problem clearly and then proceeds to describe GNeuralFlow as a solution to it and advance the field. The text provides a clear explanation of the idea with a complete mathematical framework to support the same.  \n\nSignificance: The idea of capturing relationship between time series using GNeural is significant improvement in the methodology. Learning the graph as a solution to ODE gives the flexibility to use it for unequal samples of time series. This has several applications in time series domain including streaming data. \n\nThe idea is to use traditional GCN to learn the interaction between time series for continuous time GNeuralFlow can handle the time series data observed at non-uniform time stamps without the need for interpolation or computing expensive pre-processing.", "weaknesses": "Although the paper address a big problem faced in steaming data application with unequal time samples using novel approach the idea comes with additional complexity. The complex nature of the solution is computationally expensive and is not scalable after a certain size.  The fact has been identified by authors as well and documented in the limitations sections.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the thoughtful comments, we reply to your questions below.\n\n**RE: computational cost**\n\nThe computational cost can be cubic for the number of nodes because the evaluation of the DAG constraint and the gradient involves the computation of the matrix exponential. One possible mitigation is that rather than modeling the graph adjacency matrix as a free variable, we impose structures on the matrix (such as diagonal plus low rank). Then, the matrix exponential is less expensive to compute.", "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": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: The proposed approach is found to be well motivated and well presented whereas the results are also found to be convincing. The approach is further considered both original and to have wide applications.

During the rebuttal phase, the authors included additional experimentation and comparison to baselines as well as clarified computational aspects of their modeling framework. This was appreciated by the reviewers and found to improve the paper. Consequently, the reviewers in general agreed that the manuscript is acceptable for publication. However, it was also discussed that the authors should include code for reproducibility.

The authors are thus encouraged to update their manuscript according to the additional experimentations and clarifications in their final version. Furthermore, they are also most strongly encouraged to include code in their final manuscript for reproducibility that unfortunately was not included in their original submission.
Decision: Accept (poster)