Graph Neural Flows for Unveiling Systemic Interactions Among Irregularly Sampled Time Series

Authors: giangiacomo  mercatali
andre  freitas
jie  chen
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 08:39:41

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.

Main paper: 
 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 a 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 a 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., a sequence of blackouts jumps across hundreds of kilometers [21] . We pose this question: What time series models best capture the interactive nature of different system dynamics? Graph neural networks (GNNs) [49, 46] are 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 a recurrent neural network (RNN); at every time step, a GNN is used to aggregate features over the graph. When the graph is unknown, graph-structure learning approaches have been proposed; some approaches learn a single interaction graph over time [30, 41, 11] while others infer a different one at each time step [20] . These models are in general discrete-time models, suitable for regularly spaced time points [39, 40] .To handle irregular time points, we consider continuous-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 graphbased modeling [38, 24, 25, 26, 10, 27, 3] , where the equation accounts for multiple time series and the right-hand side uses a 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 a 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 acyclic graph (DAG), where a node is conditionally independent of its non-descendents given its parents [34, 12, 43, 22] . Such a 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 is 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 a 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 a 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 (e.g., GDE [38] , LG-ODE [24] , CG-ODE [25] , CF-GODE [26] , STG-NCDE [10] , MTGODE [27] , RiTINI [3] ) in two aspects. First, through learning, the graph reveals the conditional dependencies of the time series, offering a more intuitive structure for analysis. Second, it 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.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 are summarized in the introduction section and elaborated in the following sections. Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in 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. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
background: neural ode and neural flows: 
 Denote by x(t) ∈ R d the solution of an ODEEQUATIONunder well-behaving conditions (e.g., a specified initial condition of x and Lipschitz continuity of f ). Neural ODE [7] is a modeling technique that allows uncovering the trajectory x(t) without a known right-hand side f . The technique parameterizes f by a neural network with parameters θ such that the trajectory x is a 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(t 0 ) = x 0 , we write x(t 0 ), . . . , x(t N ) = ODESolve(f, x 0 , (t 0 , . . . , t N )), where the solutions at times t 1 , . . . , t N are obtained by invoking any blackbox numerical ODE solver (such as Runge-Kutta [17] ). The training of the model parameters θ requires the gradient ∇ θ L, which can be economically computed by using the adjoint a := ∇ x L rather than expensively back-propagating 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-time techniques (such as 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 x, 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 to build a variational autoencoder (VAE) [29] , where the encoder is an RNN that evolves the hidden state h(t) over training time points and concludes a latent variable z 0 as the ODE initial condition. A drawback is that this approach still uses a 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 a 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 (1) directly:x(t) = F (t, x 0 ),by using a parameterized function F that depends on the initial condition x 0 . Optimizing the parameters of F can be more efficient because ODE solvers are 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.Question: Does the paper discuss the limitations of the work performed by the authors? Answer: [Yes] Justification: The limitations of the work are 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 in 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 on 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 a few datasets or with a 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, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a 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 be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
dag-based ode for modeling systemic interactions: 
 In this section, we motivate the form of ODE considered in this paper based on DAG modeling.In probabilistic graphical models, the conditional dependence structure is a principled framework for modeling systemic interactions. Therein, a Bayesian network [36] of n random variables y 1 , . . . , y n is a DAG with these variables as the nodes. Let A ∈ R n×n be the (weighted) adjacency matrix of the DAG, where a ij ̸ = 0 means that y i is a parent of y j . 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 1 , . . . , y n ) can be factorized into a much simpler form: p(y 1 , . . . , y n ) = n j=1 p(y j | pa(y j )), where pa(y j ) = {y i : a ij ̸ = 0} denotes the parent set of y j . The conditional dependence is a necessary condition for the causal relationship between parent y i and child y j [37] .The (linear) structural equation model, SEM [45, 14] , is a 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 ≺, where i ≺ j iff there exists an edge from i to j. Then, the DAG adjacency matrix A = [a ij ] is strictly upper triangular. The SEM model readsy 1 = ϵ 1 y 2 = a 12 y 1 + ϵ 2 y 3 = a 13 y 1 + a 23 y 2 + ϵ 3 . . . y n = a 1n y 1 + a 2n y 2 + • • • + a n-1,n y n-1 + ϵ n ,where the residuals ϵ 1 , . . . , ϵ n are (possibly correlated) Gaussian noises. In this model, y j depends on y i only when i < j. Importantly, some of the above a ij 's can be zero. Then, y j is independent of such y i 's given the rest.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.Consider an autonomous ODE ẋ = Bx with the initial condition x(0) = x 0 ∈ R 2 . This ODE describes a 2D vector field Bx; each solution curve x(t) = expm(Bt)x 0 given the initial point x 0 is a streamline instantaneously tangential to the vector field. Throughout the paper, we use expm to denote matrix exponential for matrix arguments and exp to denote element-wise exponential. Now consider n trajectories x 1 (t), . . . , x n (t). 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 jj-1 i=1 a ij Bx i gives the velocity ẋj for each j. Mathematically,EQUATIONIn other words, the solution curve x j for each j and any initial point x j 0 is a streamline instantaneously tangential to the residual field.The behaviors of non-interacting and interacting systems are fundamentally different. Figure 1 illustrates an example. The first plot shows three independent trajectories satisfying ẋi = Bx i with initial conditions x i 0 for i = 1, 2, 3. The second plot shows three conditionally dependent trajectories, under the same initial conditions, but interacting through the DAG adjacency matrix A = 0 0.5 0 0 0 0.7 0 0 0 . Because x 1 is independent of the rest, it is the same in both plots; but in the second plot, because x 2 depends on x 1 and x 3 depends on x 2 , these two trajectories are different from their counterparts in the first plot. Moreover, in the third and fourth plots, we 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 2 and x 3 deviate from the original trajectories, because of the conditional dependence.In general, the dependence among the time series may not be linear, even though structurally it is governed by the matrix A. For example, a nonlinear SEM may inspire the ODE system Bx 1 = ẋ1 ,Bx 2 = Bx 1 + ẋ2 , Bx 3 = Bx 1 ⊙ Bx 2 + ẋ3, where the dependence of x 3 on x 1 and x 2 is not linear. Thus, we consider general ODE systems of the form ẋj = f (t, {x j } ∪ pa(x j )), j = 1, . . . , n,where recall that pa(x j ) denotes the set of parents of x j . Equivalently, we write Ẋ = f (t, X, A) in the matrix form. Note that (2) is a special case of (3), which can be written as Ẋ = (I -A ⊤ )XB ⊤ . Note also that (3) is permutation equivariant and A can be any DAG matrix, not necessarily upper triangular ones.4 GNeuralFlow: An ODE-Solver-Free Method EQUATIONwhere the right-hand side f is unknown. 2 Given data X(t 0 ), . . . , X(t N ) at irregular time points, develop a model that predicts X(t) for any t ≥ t 0 as well as A. To account for practical use, at some time point t j , some rows of X(t j ) may be missing.A growing body of research addresses the problem when X has a single row (i.e., n = 1; hence, A is irrelevant), notably through using a neural network to parameterize f and using a numerical solver to evaluate X at t 0 , . . . , t N [7, 39, 5] . When n > 1 and the rows of the system are independent (i.e., f is 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 an nd-dimensional problem, but such a high dimension renders approaches using numerical solvers too costly.Instead, we use a neural network to parameterize the solution of (4) directly; i.e.,EQUATIONwhere the solution F is a function of t but depends on the initial X 0 as well as A. The neural network parameterization cannot be entirely free. First, the solution F should satisfy the initial condition F (0, X 0 , A) = X 0 . Second, the fundamental theorem on flows [32, Theorem 9.12] asserts that every smooth f with an initial condition determines a unique F (t, X 0 , A) and for any t, F (t, •, A) is a diffeomorphism. Therefore, we formulate our model for Problem 1 in the following.Solution framework. The model for X(t) is a neural network F (t, X, A) that satisfies:1. F (0, X 0 , A) = X 0 ;2. F (t, X, A) is invertible in X for any t and A; equivalently, the streamline F (t, X 0 , A) given any X 0 and A is not self-intersecting.The task is to predict the trajectories X(t) given X 0 .Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Answer: [Yes]• It should be clear whether the error bar is the standard deviation or the standard error of the mean. Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or 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.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a 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 to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
graph encoder: 
 We will make heavy use of a 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 readsEQUATION)where X is the input node feature matrix, X contains the transformed features, and W and U are parameters. GCN defines A as a symmetric normalization of A, but many alternatives are viable, such as a simple scaling of A. We also considerEQUATION)This definition is motivated by SEM, where I -A ⊤ is the operator (Section 3.3). Here, A is not symmetrized because doing so cannot distinguish edge directions. A benefit of taking (7) is that the scaling factor γ 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 A defined in (7) admits ∥ A∥ 2 ≤ 2.The neural network F in ( 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 architectureEQUATIONwhich is a building block of invertible networks [1] . Here, φ(t) satisfies φ(0) = 0 such that the requirementF (0, X 0 , A) = X 0 is met. Additionally, if φ(•) ∈ [0, 1] and g(t, •, A) is a contractive mapping, then F (t, •, A) is invertible.We let φ be the tanh function and parameterize g by using two MLPs together with a GCN:EQUATIONwhere || denotes concatenation row-wise and each MLP acts on the input matrix row-wise independently. The neural network g is generally not contractive, but bounding the spectral norm of each linear layer can theoretically guarantee contraction of an MLP [19] . Moreover, Theorem 1 indicates that bounding the spectral norm of the GCN parameters can guarantee contraction of GCN as well (because ∥ A∥ 2 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] :EQUATIONwhereh k , k = 1, 2, is computed by r k (t, X) = β • sigmoid(f k r (t, X)), c k (t, X) = tanh(f k c (t, r k (t, X) ⊙ X)), z k (t, X) = α • sigmoid(f k z (t, X)), h k (t, X) = z k (t, X) ⊙ (c k (t, X) -X).The base form X + φ(t) • h(t, 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 h 2 that incorporates the graph encoder. It can be shown that F is invertible under a deliberate choice of α and β when the MLPs f k z , f k r , f k c and the GCN are contractive. As discussed earlier, Theorem 1 indicates that GCN can be made contractive similarly as the MLPs through bounding the spectral norm of their parameters.Theorem 2. If f k z (t, •), f k r (t, •), f k c (t,•), and GCN(A, •) are contractive, the function F (t, •, A) defined in (10) is invertible whenever α(5β + 6) ≤ 2.Coupling flow. For the third design, we use normalizing flows, because they are invertible by definition. An example of the normalizing flow is the coupling flow [13] . Let {U, V } be a partitioning of the column indices 1, . . . , d. With the graph encoder, we extend a usual coupling flow block to the following:EQUATIONwhereu(t, X V , X V ) = MLP 3 MLP 1 (X V ||t) || MLP 2 ( X V ||t) v(t, X V , X V ) = MLP 4 MLP 1 (X V ||t) || MLP 2 ( X V ||t) .Here, φ u and φ v are two functions that have a range [0, 1] and attain 0 at the origin. The input X is split into X U and X V and the second part goes through the GCN encoder, producing X V . Note that one cannot apply the GCN encoder on the entire X; otherwise, the U block will have a dependency on itself in the scaling and the shift. The scaling network u and the shift network v are essentially MLPs that share initial layers; they take both X V and X V as inputs.GNeuralFlow contains two sets of parameters: the DAG matrix A and other parameters of F (call them θ), including the flow parameters, the graph encoder parameters, and possibly other parameters (e.g., in latent variable modeling, Appendix F). Let us use L(A, θ) to denote the training loss (which could be the quadratic loss in regression models, or likelihood/ELBO loss in latent variable models), making an explicit distinction between A and θ. Then, the learning problem is:EQUATION)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)) = n or tr((I + αA ⊙ A) n ) = n for any α ̸ = 0 [48, 47] . Hence, ( 12) becomes an equality-constrained problem over continuous variables, to which the augmented Lagrangian method [2] is an effective solution. We discuss the optimization details in Appendix D. 
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 are given in Appendix G. In all experiments, we split the data into train, validation, and test sets. We train with early stopping by using Adam and report the results on the test set. Standard errors are obtained by performing five repetitive runs. Hyperparmeter details are given in Appendix H. All experiments are conducted on a machine with an Nvidia A100 GPU, 8 CPU cores, and 80GB main memory.We generate four synthetic systems with the graph size varying from 3 to 30. These systems follow the graph-based equation (4) or solution (5) . They are named "Sink," "Triangle," "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 F (t, X, A) = (I -A ⊤ )(X + t 0 sign(sin(u)) du); see Appendix G for other systems and details. The forecast results are reported in Figure 2 . 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 be 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 [38] , LG-ODE [24] , and CF-GODE [26] ), graph learning methods (NRI [31] and dNRI [20] ), and a non-graph GRU variant for time series (GRU-D [6] ). For graph ODEs, the groundtruth graph is used. For neural flows and GNeuralFlow, all three flow designs are experimented with. Table 1 shows that across all datasets, GNeuralFlow significantly outperforms the baselines; moreover, GNeuralFlow 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. With the above encouraging results, we investigate the quality of graph learning. While the learning approach discussed in 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 a zero-mean Gaussian. We evaluate graph quality by using the metrics proposed by [48] : TPR (true positive rate), FDR (false positive rate), FPR (false prediction rate), and SHD (structural Hamming distance).Figure 3 shows that as 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 is slightly higher). Note that with the nonzeros of the ground truth lying between 0 and 1, 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 (4) as a data generation model; and we 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, a 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 2 indicates that GNeuralFlow 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 it does not run a numerical solver. 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 [4] . 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 to 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 a dynamic graph at each time step by utilizing the covariance matrix of the time series data within each batch. From Table 3 , we see that GNeuralFlow generally performs the best compared with the various baselines. Moreover, by using the same flow design, GNeuralFlow is always better than neural flows. These findings are consistent with those in the synthetic data case, demonstrating a good utility of our method for real-life applications.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 4 reports the forecast error on the next three time points and the estimated likelihood of the time series. We see 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.  
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 [20] . Such a 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 a 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 can be cubic, because the evaluation of the DAG constraint and the gradient involves the computation of the matrix exponential. 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 [47] , going beyond is generally believed to require either a new computational technique or a new modeling approach. One potential direction is to introduce structures into A (such as low-rankness [15] ), which admit faster matrix evaluation. 
safeguards: 
 Question: 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)?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 a 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 permissive 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 or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.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 package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a 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. 
 
 Appendix: Code is available at https://github.com/gmerca/GNeuralFlow.Proof of Theorem 1. When A is 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, ρ(B), is bounded by γ.Meanwhile, [23] show that the spectral norm ∥A∥ 2 ≤ ρ(B). Then,∥ A∥ 2 ≤ 1 + γ/γ = 2.Proof of Theorem 2. We follow the proof of [4, Theorem 1] (see A.3 of the paper), which concludes that |h(x) -h(y)| ≤ α( 5 4 β + 3 2 )|x -y|. Since GCN is contractive, such an inequality applies to both h = h 1 and h = h 2 . Then, applying Eqn ( 14) of the paper,|h 1 (x)h 2 (x)-h 1 (y)h 2 (y)| < |h 1 (x)| <1 • Lip(h 2 ) <α( 5 4 β+ 3 2 ) + |h 2 (x)| <1 • Lip(h 1 ) <α( 5 4 β+ 3 2 ) |x-y| < 2α( 5 4 β+ 3 2 )|x-y|.Therefore, when α(5β + 6) ≤ 2, the product h 1 h 2 is contractive and therefore F is invertible.The matrix isB = -4 5 -3 1 ,and the initial conditions arex 1 0 = 0.6 0.5 , x 2 0 = 0.7 0.1 , x 3 0 = 0.2 0.3 .In this section, we briefly describe the augmented Lagrangian method for solving the equalityconstrained problemmin A,θ L(A, θ) s.t. h(A) = 0,where the constraint can either be h(A) = tr(expm(A⊙A))-n or h(A) = tr((I+αA⊙A) n )-n.EQUATION)where λ and c 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 readsA k , θ k = argmin A,θ L c k λ k+1 = λ k + c k h(A k ) c k+1 = ηc k if |h(A k )| > γ|h(A k-1 )| c k else,where η ∈ (1, +∞) and γ ∈ (0, 1) are hyperparameters to be tuned.The subproblem of optimizing A and θ can be solved by using the Adam optimizer. It requires the gradient of L c and hence of h. For h(A) = tr(expm(A ⊙ A)) -n, it can be derived that ∇h(A) = expm(A ⊙ A) ⊤ ⊙ 2A, which can be obtained virtually for free after h has been evaluated. For h(A) = tr((I + αA ⊙ A) n ) -n, one may use automatic differentiation to obtain the gradient.Algorithm 1 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. Compute A k and θ k as a minimizer of ( 13) by using the Adam optimizer 4:Update Lagrange multiplier λ ← λ + ch(A k )5: if k > 0 and |h(A k )| > γ|h(A k-1 )| then 6: c ← ηc 7:end if 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 t are missing). In this case, one evaluates the graph neural flow F on a subgraph of present measurements. Rather than extracting this subgraph and the corresponding rows of X, the GCN encoder (6) offers a 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.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 of 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 an augmented space, part of which carries the graph information, as a new design complementary to those proposed in Section 4.3.On a high level, smoothing [39] and filtering [5] approaches use a neural ODE or a neural flow to continuously evolve the hidden state from time t j-1 (denoted as H(t j-1 )) to time t j (denoted as H ′ (t j )); and then use an RNN to introduce a jump (denoted as H(t j )) on observing input data X(t j ).To model and learn the interaction graph in the data space, we use the graph encoder (6) to produce a transformed data X(t j ) and use a second RNN to introduce the paired jump H(t j ) given X(t j ). Then, the pair of hidden states, H(t j ) and H(t j ), are concatenated and a standard neural flow evolve the concatenated state to the next time point, the result of which is 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) is also different. Details are presented in the following.Given observation data X(t 0 ), . . . , X(t N ), this approach produces a latent quantity Z 0 by a combined use of LSTM and neural flow, and then traces out a smooth curve Z(t) using another flow, taking Z(t 0 ) = Z 0 as the initial condition. Then, the observation data X(t j ) is recovered from Z(t j ).A VAE is used to set up the training loss. The decoder p(X(t 0 ), . . . , X(t N )|Z 0 ) is factorized asp(X(t 0 ), . . . , X(t N )|Z 0 ) = N j=0 p(X(t j )|Z(t j )),where each Z(t j ) is computed by running a standard neural flow: Z(t j ) = F (t j , Z 0 ). The encoder, on the other hand, produces a latent Gaussian Z 0 with mean µ and diagonal covariance diag(σ); that is,q(Z 0 |X(t 0 ), . . . , X(t N )) = N (Z 0 | µ, diag(σ)), [µ, log σ] = g(H(t N )), where H(t N ) is the hidden state to be elaborated soon and g is a neural network projection. 3 The ELBO loss for the VAE isL ELBO = D KL q(Z 0 |X(t 0 ), . . . , X(t N )) || p(Z 0 ) -E Z0∼q(Z0|X(t0),...,X(t N )) log p(X(t 0 ), . . . , X(t N )|Z 0 ) .Our GNeuralFlow uses a pair of LSTMs together with another neural flow to evolve the hidden state. Specifically, we maintain a pair of states H(t) and H(t), the latter of which includes the graph information. At time t j-1 , we concatenate H(t j-1 ) and H(t j-1 ), run the neural flow F to evolve the concatenated state to time t j , and apply a projection g proj so that the net result H ′ (t j ) remains in the same dimension as H(t j-1 ):H ′ (t j ) = g proj ( F (t j , H(t j-1 )|| H(t j-1 )) ).Then, we run a pair of LSTMs to obtain the paired states at time t j :H(t j ) = LSTM 1 H ′ (t j ), X(t j ) , H(t j ) = LSTM 2 H ′ (t j ), X(t j ) ,where the second LSTM is applied to the transformed observation data X(t j ) produced by the GCN encoder (6) . By doing so, the graph models the interaction inside the data X(t) rather than the hidden states H(t).As opposed to the preceding approach, the filtering approach uses only a decoder. Each time, it first evolves the hidden state to H ′ (t j ) and then runs a GRU to update the hidden state to H(t j ). This approach maintains two Gaussians, the first one models the observation X(t j ):N (X(t j ) | µ j obs , diag(σ j obs )), [µ j obs , log σ j obs ] = g obs (H ′ (t j )), while the second one models the jump caused by the GRU:N (µ j post , diag(σ j post )), [µ j post , log σ j post ] = g post (H(t j )).The training loss aims at maximizing the observation data likelihood while minimizing the KL divergence of the two Gaussians:L = - N j=1 log N (X(t j ) | µ j obs , diag(σ j obs ))+λD KL N (µ j obs , diag(σ j obs )) || N (µ j post , diag(σ j post )) .Our GNeuralFlow uses a pair of GRUs together with a standard neural flow to evolve the hidden state. Specifically, we maintain a pair of states H(t) and H(t), the latter of which includes the graph information. At time t j-1 , we concatenate H(t j-1 ) and H(t j-1 ), run the neural flow F to evolve the concatenated state to time t j , and apply a projection g proj so that the net result H ′ (t j ) remains in the same dimension as H(t j-1 ):H ′ (t j ) = g proj ( F (t j , H(t j-1 )|| H(t j-1 )) ).Then, we run a pair of GRUs to obtain the paired states at time t j :H(t j ) = GRU 1 H ′ (t j ), g prep (X(t j ), H ′ (t j )) , H(t j ) = GRU 2 H ′ (t j ), g prep ( X(t j ), H ′ (t j )) ,where the second GRU is applied to the transformed observation data X(t j ) produced by the GCN encoder (6) . By doing so, the graph models the interaction inside the data X(t) rather than the hidden states H(t). The datasets used in this paper include four synthetic ODE systems and four real-life datasets. Table 5 summarizes the basic information of these datasets, tasks, evaluation metrics, and learning methods.We define four interacting systems based on either the ODE Ẋ = f (t, X, A) or the solution X(t) = F (t, X 0 , A):• Sink (2D): f (t, X, A) = (I -A ⊤ )XB ⊤ where B = -4 10 -3 2 • Triangle (1D): F (t, X, A) = (I -A ⊤ )(X + t 0 sign(sin(u)) du) • Sawtooth (1D): F (t, X, A) = (I -A ⊤ )(X + t -⌊t⌋) • Square (1D): F (t, X, A) = (I -A ⊤ )(X + sign(sin(t)))For Sink, X ∈ R n×2 ; while for the other three systems, X ∈ R n×1 , where n is the number of trajectories (graph nodes) in the system. The initial condition X 0 is uniformly sampled from [0, 1] n×2 for Sink, and from [-2, 2] n×1 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 a sparse n × n matrix A with a 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).We use four real-life datasets preprocessed by [4] .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, we also evaluate different models on the reconstruction quality. We treat each sensor as one graph node.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, we 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 to predict the next three measurements in the 12 hour interval after the observation window of 36 hours. We treat each variable as one graph node. We reuse the architecture parameters and training hyperparameters in [4] and only tune the graph learning hyperparameters (see Algorithm 1 in Section D). Table 6 lists the tuned values.In practice, the initial time point t0 may not be zero, in which case one may shift the ODE along the temporal dimension by t0.When the encoder is run backward in time, one uses H(t0) instead of H(t0). 


Decision: 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.

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: * 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.
* To learn the DAG matrix, the proposed method employs NOTEARS and DAG-GNN. How do these DAG learning algorithms impact the model performance?

* 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.
questions: questions: See Weaknesses.
limitations: limitations: No potential negative societal impact.
paper_summary: paper_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.
strength_weakness: strength_weakness: Weaknesses:
Flow 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. 

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

The set of comparisons featured here is rather limited. Here are several works that could be compared against: 

Kipf et al., Neural Relational Inference for Interacting Systems, ICML, 2018.
Graber and Schwing. Dynamic Neural Relational Inference, CVPR, 2020 (dynamic graph inference for discrete-time systems)
Nishikawa-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)
Deleu 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)
Smith 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)
Hiremath 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)
Bhaskar et al. Inferring dynamic regulatory interaction graphs from time series data with perturbations, LOG 2023
questions: 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?

Does the performance gain over existing graph ODE models come from the learnable graph? Is there an ablation study for that?

Is the model learning a constant graph independent of time? How does it perform when there is a significant change in data dependence?
limitations: limitations: The authors have adequately addressed the limitations.
paper_summary: paper_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.
strength_weakness: strength_weakness: 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).
2. *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.		
3. *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.
4. *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).

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

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

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

[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.
questions: questions: Based on the *weaknesses* above please focus on the following aspects:
1. **Experiments:** Explain choices of datasets and baselines or extend the experimental evaluation (W1). Conduct qualitative evaluations of graph-based causal representations (W2).
2. **Limitations:** Perform a comprehensive computational cost analysis comparing GNeuralFlow with simpler and more complex methods to highlight efficiency gains (W3).
3. **Contribution:** Please better position the presented contribution among relevant works (W4).
limitations: 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.
paper_summary: paper_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.
strength_weakness: strength_weakness: - The proposed method lacks a theoretical analysis.
- The authors can consider using state-of-the-art GNN as backbone for graph encoder.
- Can the authors apply this method to real-world datasets?
- Computational complexity and corresponding comparisons with baselines are missing.
questions: questions: Please see comments in Weaknesses.
limitations: limitations: Not applicable.
paper_summary: paper_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.  

The 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.
strength_weakness: strength_weakness: 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.
questions: 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: 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.

Comments: title: Clarified neural flows
Thank you.
comment: Dear Reviewer zrVL,

Thank you for your valuable suggestions to enhance our work. In response to your request for further improvement, we have explored an additional GNN model. Our findings in the table below indicate that significant improvement can be achieved by combining the ResNet model with a Message Passing GNN (Veličković, 2022; Bronstein et al., 2021).


| model             | Sink            | Triangle        | Sawtooth        | Square          |
| ----------------- | --------------- | --------------- | --------------- | --------------- |
| GNF-resnet (ours) | 3.95(±0.32)     | 2.32(±0.11)     | 3.84(±0.06)     | 8.24(±0.64)     |
| GNF-MPGNN (ours)  | **3.45(±0.13)** | **2.00(±0.03)** | **2.79(±0.05)** | **4.32(±0.12)** |

We kindly ask you to consider raising the score if this update, or any of our previous revisions, addresses your concerns either fully or partially.


**References**

Veličković 2022. Message passing all the way up. ICLR

Bronstein et al 2021. Geometric deep learning. Arxiv.
Thanks for the clarification of the neural flows approach, in other words I see that you directly model the integral curve.  I think this is also interesting for a graph. I also acknowledge the effort of adding the baselines. The idea that you would use a DAG because others do is not a great justification, we know medicine is filled with circularly dependent variables. Therefore I will raise my score slightly to 5.
I appreciate the author for the detailed response. After carefully reading the rebuttal, I am retaining my score.
Thanks for further clarification. I will keep my score.
Thank you for the replies, which have partially addressed the issues raised by my side (e.g., baselines and computational costs). Considering other reviewers' comments as well as the impact of the contribution (when it comes to mathematical background, performance improvements and complexity), I would prefer to maintain my rating.