GRU-ODE-Bayes: Continuous modeling of sporadically-observed time series

Introduction

Multivariate time series are ubiquitous in various domains of science, such as healthcare (Jensen et al., 2014), astronomy (Scargle, 1982), or climate science (Schneider, 2001). Much of the methodology for time-series analysis assumes that signals are measured systematically at fixed time intervals. However, much real-world data can be sporadic (i.e., the signals are sampled irregularly and not all signals are measured each time). A typical example is patient measurements, which are taken when the patient comes for a visit (e.g., sometimes skipping an appointment) and where not every measurement is taken at every visit. Modeling then becomes challenging as such data violates the main assumptions underlying traditional machine learning methods (such as recurrent neural networks).

Recently, the Neural Ordinary Differential Equation (ODE) model (Chen et al., 2018) opened the way for a novel, continuous representation of neural networks. As time is intrinsically continuous, this framework is particularly attractive for time-series analysis. It opens the perspective of tackling the issue of irregular sampling in a natural fashion, by integrating the dynamics over whatever time interval needed. Up to now however, such ODE dynamics have been limited to the continuous generation of observations (e.g., decoders in variational auto-encoders (VAEs) (Kingma & Welling, 2013) or normalizing flows (Rezende et al., 2014)).

Instead of the encoder-decoder architecture where the ODE part is decoupled from the input processing, we introduce a tight integration by interleaving the ODE and the input processing steps. Conceptually, this allows us to drive the dynamics of the ODE directly by the incoming sporadic inputs. To this end, we propose (1) a continuous time version of the Gated Recurrent Unit and (2) a Bayesian update network that processes the sporadic observations. We combine these two ideas to form the GRU-ODE-Bayes method.

The tight coupling between observation processing and ODE dynamics allows the proposed method to model fine-grained nonlinear dynamical interactions between the variables. As illustrated in Figure 1, GRU-ODE-Bayes can (1) quickly infer the unknown parameters of the underlying stochastic process and (2) learn the correlation between its variables (red arrows in Figure 1). In contrast, the encoder-decoder based method NeuralODE-VAE proposed by Chen et al. (2018) captures the general structure of the process without being able to recover detailed interactions between the variables (see Section 4 for detailed comparison).

Our model enjoys important theoretical properties. We frame our analysis in a general way by considering that observations follow the dynamics driven by a stochastic differential equation (SDE). In Section 4 and Appendix I, we show that GRU-ODE-Bayes can exactly represent the corresponding Fokker-Planck dynamics in the special case of the Ornstein-Uhlenbeck process, as well as in generalized versions of it. We further perform an empirical evaluation and show that our method outperforms the state of the art on healthcare and climate data (Section 5).

We assume that observations $\mathbf{y}_{i}$ are sampled from the realizations of a $D$-dimensional stochastic process $\mathbf{Y}(t)$ whose dynamics is driven by an unknown SDE:

where $d\mathbf{W}(t)$ is a Wiener process. The distribution of $\mathbf{Y}(t)$ then evolves according to the celebrated Fokker-Planck equation (Risken, 1996). We refer to the mean and covariance parameters of its probability density function (PDF) as $\mu_{\mathbf{Y}}(t)$ and $\Sigma_{\mathbf{Y}}(t)$.

Our goal will be to model the unknown temporal functions $\mu_{\mathbf{Y}}(t)$ and $\Sigma_{\mathbf{Y}}(t)$ from the sporadic measurements $\mathbf{y}_{i}$. These are obtained by sampling the random vectors $\mathbf{Y}(t)$ at times $\mathbf{t}_{i}$ with some observation noise $\bm{\epsilon}$. Not all dimensions are sampled each time, resulting in missing values in $\mathbf{y}_{i}$. In contrast to classical SDE inference (Särkkä & Solin, 2019), we consider the functions $\mu_{\mathbf{Y}}(t)$ and $\Sigma_{\mathbf{Y}}(t)$ are parametrized by neural networks.

This SDE formulation is general. It embodies the natural assumption that seemingly identical processes can evolve differently because of unobserved information. In the case of intensive care, as developed in Section 5, it reflects the evolving uncertainty regarding the patient’s future condition.

Proposed method

At a high level, we propose a dual mode system consisting of (1) a GRU-inspired continuous-time state evolution (GRU-ODE) that propagates in time the hidden state $\mathbf{h}$ of the system between observations and (2) a network that updates the current hidden state to incorporate the incoming observations (GRU-Bayes). The system switches from propagation to update and back whenever a new observation becomes available.

To derive the GRU-based ODE, we first show that the GRU proposed by Cho et al. (2014) can be written as a difference equation. First, let $\mathbf{r}_{t}$, $\mathbf{z}_{t}$, and $\mathbf{g}_{t}$ be the reset gate, update gate, and update vector of the GRU:

where $\odot$ is the elementwise product. Then the standard update for the hidden state $\mathbf{h}$ of the GRU is

We can also write this as $\mathbf{h}_{t}=\operatorname{\mathbf{GRU}}(\mathbf{h}_{t-1},\mathbf{x}_{t})$. By subtracting $\mathbf{h}_{t-1}$ from this state update equation and factoring out $(1-\mathbf{z}_{t})$, we obtain a difference equation

This difference equation naturally leads to the following ODE for $\mathbf{h}(t)$:

where $\mathbf{z}$, $\mathbf{g}$, $\mathbf{r}$ and $\mathbf{x}$ are the continuous counterpart of Eq. 2. See Appendix A for the explicit form.

We name the resulting system GRU-ODE. Similarly, we derive the minimal GRU-ODE, a variant based on the minimal GRU (Zhou et al., 2016), described in appendix G.

In case continuous observations or control signals are available, they can be naturally fed to the GRU-ODE input $\mathbf{x}(t)$. For example, in the case of clinical trials, the administered daily doses of the drug under study can be used to define a continuous input signal. If no continuous input is available, then nothing is fed as $\mathbf{x}(t)$ and the resulting ODE in Eq. 3 is autonomous, with $\mathbf{g}(t)$ and $\mathbf{z}(t)$ only depending on $\mathbf{h}(t)$.

2 General properties of GRU-ODE

GRU-ODE enjoys several useful properties:

Boundedness. First, the hidden state $\mathbf{h}(t)$ stays within the $$rangeWe use the notation$$ to also mean multi-dimensional range (i.e., all elements are within $$).. This restriction is crucial for the compatibility with the GRU-Bayes model and comes from the negative feedback term in Eq. 3, which stabilizes the resulting system. In detail, if the$j$-th dimension of the starting state$\mathbf{h}(0)$is within$$, then$\mathbf{h}(t)_{j}$will always stay within$$ because

This can be derived from the ranges of $\mathbf{z}$ and $\mathbf{g}$ in Eq. 2. Moreover, would $\mathbf{h}(0)$ start outside of the $$region, the negative feedback would quickly push$\mathbf{h}(t)$ into this region, making the system also robust to numerical errors.

Continuity. Second, GRU-ODE is Lipschitz continuous with constant $K=2$. Importantly, this means that GRU-ODE encodes a continuity prior for the latent process $\mathbf{h}(t)$. This is in line with the assumption of a continuous hidden process generating observations (Eq. 1). In Section 5.5, we demonstrate empirically the importance of this prior in the small-sample regime.

General numerical integration. As a parametrized ODE, GRU-ODE can be integrated with any numerical solver. In particular, adaptive step size solvers can be used. Our model can then afford large time steps when the internal dynamics is slow, taking advantage of the continuous time formulation of Eq. 3. It can also be made faster with sophisticated ODE integration methods. We implemented the following methods: Euler, explicit midpoint, and Dormand-Prince (an adaptive step size method). Appendix C illustrates that the Dormand-Prince method requires fewer time steps.

3 GRU-Bayes

GRU-Bayes is the module that processes the sporadically incoming observations to update the hidden vectors, and hence the estimated PDF of $\mathbf{Y}(t)$. This module is based on a standard GRU and thus operates in the region $$that is required by GRU-ODE. In particular, GRU-Bayes is able to update$\mathbf{h}(t)$ to any point in this region. Any adaptation is then within reach with a single observation.

where $\mathbf{h}(t_{-})$ and $\mathbf{h}(t_{+})$ denote the hidden representation before and after the jump from GRU-Bayes update. We also investigate an alternative option where the $\mathbf{h}(t)$ is updated by each observed dimension sequentially. We call this variant GRU-ODE-Bayes-seq (see Appendix F for more details). In Appendix H, we run an ablation study of the proposed GRU-Bayes architecture by replacing it with a MLP and show that the aforementioned properties are crucial for good performance.

4 GRU-ODE-Bayes

The proposed GRU-ODE-Bayes combines GRU-ODE and GRU-Bayes. The GRU-ODE is used to evolve the hidden state $\mathbf{h}(t)$ in continuous time between the observations and GRU-Bayes transforms the hidden state, based on the observation $\mathbf{y}$, from $\mathbf{h}(t_{-})$ to $\mathbf{h}(t_{+})$. As best illustrated in Figure 2, the alternation between GRU-ODE and GRU-Bayes results in an ODE with jumps, where the jumps are at the locations of the observations.

GRU-ODE-Bayes is best understood as a filtering approach. Based on previous observations (until time $t_{k}$), it can estimate the probability of future observations. Like the celebrated Kalman filter, it alternates between a prediction (GRU-ODE) and a filtering (GRU-Bayes) phase. Future values of the time series are predicted by integrating the hidden process $\mathbf{h}(t)$ in time, as shown on the green solid line in Figure 2. The update step discretely updates the hidden state when a new measurement becomes available (dotted blue line). Let’s note that unlike the Kalman filter, our approach is able to learn complex dynamics for the hidden process.

In this way, we force our model to learn to mimic a Bayesian update.

5 Implementation

The pseudocode of GRU-ODE-Bayes is depicted in Algorithm 3, where a forward pass is shown for a single time series Code is available in the following anonymous repository : https://github.com/edebrouwer/gru_ode_bayes . For mini-batching several time series we sort the observation times across all time series and for each unique time point $\mathbf{t}[k]$, we create a list of the time series that have observations. The main loop of the algorithm iterates over this set of unique time points. In the GRU-ODE step, we propagate all hidden states jointly. The GRU-Bayes update and the loss calculation are only executed on the time series that have observation at that particular time point. The complexity of our approach then scales linearly with the number of observations and quadratically with the dimension of the observations. When memory cost is a bottleneck, the gradient can be computed using the adjoint method, without backpropagating through the solver operations (Chen et al., 2018).

Related research

Machine learning has a long history in time series modelling (Mitchell, 1999; Gers et al., 2000; Wang et al., 2006; Chung et al., 2014). However, recent massive real-world data collection, such as electronic health records (EHR), increase the need for models capable of handling such complex data (Lee et al., 2017). As stated in the introduction, their sporadic nature is the main difficulty.

To address the nonconstant sampling, a popular approach is to recast observations into fixed duration time bins. However, this representation results in missing observation both in time and across features dimensions. This makes the direct usage of neural network architectures tricky. To overcome this issue, the main approach consists in some form of data imputation and jointly feeding the observation mask and times of observations to the recurrent network (Che et al., 2018; Choi et al., 2016a; Lipton et al., 2016; Du et al., 2016; Choi et al., 2016b; Cao et al., 2018). This approach strongly relies on the assumption that the network will learn to process true and imputed samples differently. Despite some promising experimental results, there is no guarantee that it will do so. Some researchers have tried to alleviate this limitation by introducing more meaningful data representation for sporadic time series (Rajkomar et al., 2018; Razavian & Sontag, 2015; Ghassemi et al., 2015), like tensors (De Brouwer et al., 2018; Simm et al., 2017).

Others have addressed the missing data problem with generative probabilistic models. Among those, (multitask) Gaussian processes (GP) are the most popular by far (Bonilla et al., 2008). They have been used for smart imputation before a RNN or CNN architecture (Futoma et al., 2017; Moor et al., 2019), for modelling a hidden process in joint models (Soleimani et al., 2018), or to derive informative representations of time series (Ghassemi et al., 2015). GPs have also been used for direct forecasting (Cheng et al., 2017). However, they usually suffer from high uncertainty outside the observation support, are computationally intensive (Quiñonero-Candela & Rasmussen, 2005), and learning the optimal kernel is tricky. Neural Processes, a neural version of GPs, have also been introduced by Garnelo et al. (2018). In contrast with our work that focuses on continuous-time real-valued time series, continuous time modelling of time-to-events has been addressed with point processes (Mei & Eisner, 2017) and continuous time Bayesian networks (Nodelman et al., 2002). Yet, our continuous modelling of the latent process allows us to straightforwardly model a continuous intensity function and thus handle both real-valued and event type of data. This extension was left for future work.

Most recently, the seminal work of Chen et al. (2018) suggested a continuous version of neural networks that overcomes the limits imposed by discrete-time recurrent neural networks. Coupled with a variational auto-encoder architecture (Kingma & Welling, 2013), it proposed a natural way of generating irregularly sampled data. However, it transferred the difficult task of processing sporadic data to the encoder, which is a discrete-time RNN. In a work submitted concomitantly to ours (Rubanova et al., 2019), the authors proposed a convincing new VAE architecture that uses a Neural-ODE architecture for both encoding and decoding the data.

Related auto-encoder approaches with sequential latents operating in discrete time have also been proposed (Krishnan et al., 2015, 2017). These models rely on classical RNN architectures in their inference networks, hence not addressing the sporadic nature of the data. What is more, if they have been shown useful for smoothing and counterfactual inference, their formulation is less suited for forecasting. Our method also has connections to the Extended Kalman Filter (EKF) that models the dynamics of the distribution of processes in continuous time. However, the practical applicability of the EKF is limited because of the linearization of the state update and the difficulties involved in identifying its parameters. Importantly, the ability of the GRU to learn long-term dependencies is a significant advantage.

Finally, other works have investigated the relationship between deep neural networks and partial differential equations. An interesting line of research has focused on deriving better deep architectures motivated by the stability of the corresponding patial differential equations (PDE) (Haber & Ruthotto, 2017; Chang et al., 2019). Despite their PDE motivation, those approaches eventually designed new discrete architectures and didn’t explore the application on continuous inputs and time.

Application to synthetic SDEs

Figure 1 illustrates the capabilities of our approach compared to NeuralODE-VAE on data generated from a process driven by a multivariate Ornstein-Uhlenbeck (OU) SDE with random parameters. Compared to NeuralODE-VAE, which retrieves the average dynamics of the samples, our approach detects the correlation between both features and updates its predictions more finely as new observations arrive. In particular, note that GRU-ODE-Bayes updates its prediction and confidence on a feature even when only the other one is observed, taking advantage from the fact that they are correlated. This can be seen on the left pane of Figure 1 where at time $t=3$, Dimension 1 (blue) is updated because of the observation of Dimension 2 (green).

By directly feeding sporadic inputs into the ODE, GRU-ODE-Bayes sequentially filters the hidden state and thus estimates the PDF of the future observations. This is the core strength of the proposed method, allowing it to perform long-term predictions.

In Appendix I, we further show that our model can exactly represent the dynamics of multivariate OU process with random variables. Our model can also handle nonlinear SDEs as shown in Appendix J where we present an example inspired by the Brusselator (Prigogine, 1982), a chaotic ODE.

Empirical evaluation

We evaluated our model on two data sets from different application areas: healthcare and climate forecasting. In both applications, we assume the data consists of noisy observations from an underlying unobserved latent process as in Eq. 1. We focused on the general task of forecasting the time series at future time points. Models are trained to minimize negative log-likelihood.

We used a comprehensive set of state-of-the-art baselines to compare the performance of our method. All models use the same hidden size representation and comparable number of parameters.

NeuralODE-VAE (Chen et al., 2018). We model the time derivative of the hidden representation as a 2-layer MLP. To take missingness across features into account, we add a mechanism to feed an observation mask.

Imputation Methods. We implemented two imputation methods as described in Che et al. (2018): GRU-Simple and GRU-D.

Sequential VAEs (Krishnan et al., 2015, 2017). We extended the deep Kalman filter architecture by feeding an observation mask and updating the loss function accordingly.

T-LSTM (Baytas et al., 2017). We reused the proposed time-aware LSTM cell to design a forecasting RNN with observation mask.

2 Electronic health records

Electronic Health Records (EHR) analysis is crucial to achieve data-driven personalized medicine (Lee et al., 2017; Goldstein et al., 2017; Esteva et al., 2019). However, efficient modeling of this type of data remains challenging. Indeed, it consists of sporadically observed longitudinal data with the extra hurdle that there is no standard way to align patients trajectories (e.g., at hospital admission, patients might be in very different state of progression of their condition). Those difficulties make EHR analysis well suited for GRU-ODE-Bayes.

We use the publicly available MIMIC-III clinical database (Johnson et al., 2016), which contains EHR for more than 60,000 critical care patients. We select a subset of 21,250 patients with sufficient observations and extract 96 different longitudinal real-valued measurements over a period of 48 hours after patient admission. We refer the reader to Appendix K for further details on the cohort selection. We focus on the predictions of the next 3 measurements after a 36-hour observation window.

3 Climate forecast

From short-term weather forecast to long-range prediction or assessment of systemic changes, such as global warming, climatic data has always been a popular application for time-series analysis. This data is often considered to be regularly sampled over long periods of time, which facilitates their statistical analysis. Yet, this assumption does not usually hold in practice. Missing data are a problem that is repeatedly encountered in climate research because of, among others, measurement errors, sensor failure, or faulty data acquisition. The actual data is then sporadic and researchers usually resort to imputation before statistical analysis (Junninen et al., 2004; Schneider, 2001).

We use the publicly available United State Historical Climatology Network (USHCN) daily data set (Menne et al., ), which contains measurements of 5 climate variables (daily temperatures, precipitation, and snow) over 150 years for 1,218 meteorological stations scattered over the United States. We selected a subset of 1,114 stations and an observation window of 4 years (between 1996 and 2000). To make the time series sporadic, we subsample the data such that each station has an average of around 60 observations over those 4 years. Appendix L contains additional details regarding this procedure. The task is then to predict the next 3 measurements after the first 3 years of observation.

4 Results

We report the performance using 5-fold cross-validation. Hyperparameters (dropout and weight decay) are chosen using an inner holdout validation set (20%) and performance are assessed on a left-out test set (10%). Those folds are reused for each model we evaluated for sake of reproducibility and fair comparison (More details in Appendix O). Performance metrics for both tasks (NegLL and MSE) are reported in Table 1. GRU-ODE-Bayes handles the sporadic data more naturally and can more finely model the dynamics and correlations between the observed features, which results in higher performance than other methods for both data sets. In particular, GRU-ODE-Bayes unequivocally outperforms all other methods on both data sets.

5 Impact of continuity prior

To illustrate the capabilities of the derived GRU-ODE cell presented in Section 2.1, we consider the case of time series forecasting with low sample size. In the realm of EHR prediction, this could be framed as a rare disease setup, where data is available for few patients only. In this case of scarce number of samples, the continuity prior embedded in GRU-ODE is crucially important as it provides important prior information about the underlying process.

To highlight the importance of the GRU-ODE cell, we compare two versions of our model : the classical GRU-ODE-Bayes and one where the GRU-ODE cell is replaced by a discretized autonomous GRU. We call the latter GRU-Discretized-Bayes. Table 2 shows the results for MIMIC-III with varying number of patients in the training set. While our discretized version matches the continuous one on the full data set, GRU-ODE cell achieves higher accuracy when the number of samples is low, highlighting the importance of the continuity prior. Log-likelihood results are given in Appendix M.

Conclusion and future work

We proposed a model combining two novel techniques, GRU-ODE and GRU-Bayes, which allows feeding sporadic observations into a continuous ODE dynamics describing the evolution of the probability distribution of the data. Additionally, we showed that this filtering approach enjoys attractive representation capabilities. Finally, we demonstrated the value of GRU-ODE-Bayes on both synthetic and real-world data. Moreover, while a discretized version of our model performed well on the full MIMIC-III data set, the continuity prior of our ODE formulation proves particularly important in the small-sample regime, which is particularly relevant for real-world clinical data where many data sets remain relatively modest in size.

In this work, we focused on time-series data with Gaussian observations. However, GRU-ODE-Bayes can also be extended to binomial and multinomial observations since the respective NegLL and KL-divergence are analytically tractable. This allows the modeling of sporadic observations of both discrete and continuous variables.

Acknowledgements

Edward De Brouwer is funded by a FWO-SB grant. Yves Moreau is funded by (1) Research Council KU Leuven: C14/18/092 SymBioSys3; CELSA-HIDUCTION, (2) Innovative Medicines Initiative: MELLODY, (3) Flemish Government (ELIXIR Belgium, IWT, FWO 06260) and (4) Impulsfonds AI: VR 2019 2203 DOC.0318/1QUATER Kenniscentrum Data en Maatschappij. Computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation - Flanders (FWO) and the Flemish Government – department EWI. We also gratefully acknowledge the support of NVIDIA Corporation with the donation of the Titan Xp GPU used for this research.

References

Appendix A Full formulation of the GRU-ODE cell

The full ODE equation for GRU-ODE is the following:

Appendix B Lipschitz continuity of GRU-ODE

As $\mathbf{h}$ is differentiable and continous on $t$, we know from the mean value theorem that for any $t_{a},t_{b}\in t$, there exists $t^{*}\in(t_{a},t_{b})$ such that

Taking the euclidean norm of the previous expression, we find

Furthermore, we showed that $\mathbf{h}$ is bounded on $$. Hence, because of the bounded functions appearing in the ODE (sigmoids and hyperbolic tangents), the derivative of$\mathbf{h}$is itself bounded by$$. We conclude that$\mathbf{h}(t)$is Lipschitz continuous with constant$K=2$.

Appendix C Comparison of numerical integration methods

We implemented three numerical integration methods, among which the classical Euler method and the Dormand-Prince method (DOPRI). DOPRI is a popular adaptive step size numerical integration method relying on 2 Runge-Kutta solvers of order 4 and 5. The advantage of adaptive step size methods is that they can tune automatically the number of steps to integrate the ODE until the desired point.

Figure 4 illustrates the number of steps taken by both solvers when given the same data and same ODE. We observe that using an adaptive step size results in half as many time steps. More steps are taken near the observations and as the underlying process becomes smoother, the step size increase, as observed on the right side of the figure. However, each time step requires significantly fewer computations for Euler than for DOPRI, so that Euler’s method appears more than competitive on the data and simulations we have considered so far. Nevertheless, DOPRI might still be preferred as default method because of its better numerical stability.

Appendix D Mapping to deal with missingness across features

Appendix E Observation model mapping

The mapping from hidden $\mathbf{h}$ to the parameters of the distribution $\mu_{\mathbf{Y}(t)}$ and $\log(\Sigma_{\mathbf{Y}(t)})$. For this purpose we use a classical multi-layer perceptron architecture with a 25 dimensional hidden layer. Note that me map to the log of the variance in order to keep it positive.

Appendix F GRU-ODE-Bayes-seq

On top of the architecture described in the main bulk of this paper, we also propose a variant which process the sporadic inputs sequentially. In other words, GRU-Bayes will update its prediction on the hidden $\mathbf{h}$ for one input dimension after the other rather than jointly. We call this approach GRU-ODE-Bayes-seq.

In this sequential approach for GRU-Bayes, we process one-by-one all dimensions of $\mathbf{y}[k]$ that were observed at time $t[k]$ by first applying the preprocessing to each and then sending them to the GRU unit. The preprocessing steps are the same as in the nonsequential scheme (Appendix D) but without concatenation at the end because only one dimension is processed at a time. Note that the preprocessing of dimensions cannot be done in parallel as the hidden state $\mathbf{h}$ changes after each dimension is processed, which affects the computed $\theta_{d}$ and thus the resulting vector $\mathbf{q}_{d}$.

Appendix G Minimal GRU-ODE

Following the same reasoning as for the full GRU cell, we also derived the minimal GRU-ODE cell, based on the minimal GRU cell. The minimal GRU writes :

This can be rewritten as the following difference equation :

Appendix H Ablation study of GRU-Bayes

In order to demonstrate the fitness of the GRU-Bayes module for our architecture, we ran an ablation study where we replaced the GRU-Bayes with a 2 layers multi-layer perceptron. We used a tanh activation function for the hidden units and a linear activation for the output layer. We evaluate the performance of this modified architecture on the MIMIC dataset for the forecasting task. The results are presented in table 3. The proposed architecture outperforms a simple MLP module, due to the properties described in sections 2.2 and 2.3.

Appendix I Application to Ornstein-Uhlenbeck SDEs

We demonstrate the capabilities of our approach on data generated from a process driven by an SDE as in Eq. 1. In particular, we focus on extensions of the multidimensional correlated Ornstein-Uhlenbeck (OU) process with varying parameters. For a particular sample $i$, the dynamics is given by the following SDE:

where $\mathbf{W}(t)$ is a $D$-dimensional correlated Wiener process, $\mathbf{r}_{i}$ is the vector of targets, and $\theta_{i}$ is the reverting strength constant. For simplicity, we consider $\theta_{i}$ and $\sigma_{i}$ parameters as scalars. Each sample $\mathbf{y}_{i}$ is then obtained via the realization of process (5) with sample-specific parameters.

We now show that our model exactly captures the dynamics of the distribution of $\mathbf{Y}(t)$ as defined in Eq. 5. The evolution of the PDF of a diffusion process is given by the corresponding Fokker-Planck equation. For the OU process, this PDF is Gaussian with time-dependent mean and covariance. Conditioned on a previous observation at time $t^{*}$, this gives

Correlation of $\mathbf{Y}(t)$ is constant and equal to $\rho$, the correlation of the Wiener processes. The dynamics of the mean and variance parameters can be better expressed in the following ODE form:

With initial conditions $\mu_{\mathbf{Y}}(0,t^{*})=\mathbf{Y}(t^{*})$ and $\sigma_{\mathbf{Y}}^{2}(0,t^{*})=0$. We next investigate how specific versions of this ODE can be represented by our GRU-ODE-Bayes.

In standard OU, the parameters $\mathbf{r}_{i}$, $\sigma_{i}$, and $\theta_{i}$ are fixed and identical for all samples. The ODE (6) is linear and can then be represented directly with GRU-ODE by storing $\mu_{\mathbf{Y}}(t,t^{*})$ and $\sigma_{\mathbf{Y}}^{2}(t,t^{*})$ in the hidden state $\mathbf{h}(t)$ and matching the Equations (3) and (6). The OU parameters $\mathbf{r}_{i}$, $\sigma_{i}$ and $\theta_{i}$ are learned and encoded in the weights of GRU-ODE. GRU-Bayes then updates the hidden state and stores $\mu_{\mathbf{Y}}(t,t^{*})$ and $\sigma_{\mathbf{Y}}^{2}(t,t^{*})$.

I.1.2 Generalized Ornstein-Uhlenbeck processes

When parameters are allowed to vary over samples, these have to be encoded in the hidden state of GRU-ODE-Bayes, rather than in the fixed weights. For $\mathbf{r}_{i}$ and $\sigma_{i}$, GRU-Bayes computes and stores their current estimates as the observations arrive. This is based on previous hidden and current observation as in Eq. 4. The GRU-ODE module then simply has to keep these estimates unchanged between observations:

This can be easily done by switching off the update gate (i.e., setting $\mathbf{z}(t)$ to 1 for these dimensions). These hidden states can then be used to output the mean and variance in Eq. 6, thus enabling the model to represent generalized Ornstein-Uhlenbeck processes with sample-dependent $\mathbf{r}_{i}$ and $\sigma_{i}$.

Perfect representation for sample dependent $\theta_{i}$ requires the multiplication of inputs in Eq. 6, which GRU-ODE is not able to perform exactly but should be able to approximate reasonably well. If an exact representation is required, the addition of a bilinear layer is sufficient.

Furthermore, the same reasoning applies when parameters are also allowed to change over time within the same sample. GRU-Bayes is again able to update the hidden vector with the new estimates.

I.1.3 Non-aligned time series

Our approach can also handle samples that would be dephased in time (i.e, the observation windows are not aligned on an intrinsic time scale). Longitudinal patient data recorded at different stages of the disease for each patient is one example, developed in Section 5. This setting is naturally handled by the GRU-Bayes module.

I.2 Case Study: 2D Ornstein-Uhlenbeck Process

We evaluate our model on a 2-dimensional OU process with correlated Brownian motion as defined in Eq. 5. For best illustration of its capabilities, we consider the three following cases.

In the first setting, $\mathbf{r}_{i}$ varies across samples as $\mathbf{r}_{i}^{1}\sim\mathcal{U}(0.5,1.5)$ and $\mathbf{r}_{i}^{2}\sim\mathcal{U}(-1.5,-0.5)$. The correlation between the Wiener processes $\rho$ is set to $0.99$. We also set $\sigma=0.1$ and $\theta=1$. The second case, which we call random lag is similar to the first one but adds an extra uniformly distributed random lag to each sample. Samples are then time shifted by some $\Delta_{T}\sim\mathcal{U}(0,0.5)$. The third setting is identical to the first but with $\rho=0$ (i.e., both dimensions are independent and no information is shared between them).

We evaluate all methods and settings on the forecast of samples after time $t=4$. The training set contains 10,000 samples with an average of 20 observations scattered over a 10-second time interval. Models are trained with a negative log-likelihood objective function, but mean square errors (MSE) are also reported. We compare our methods to NeuralODE-VAE (Chen et al., 2018). Additionally, we consider an extended version of this model where we also feed the observation mask, called NeuralODE-VAE-Mask.

I.2.2 Empirical evaluation

Figure 1 shows a comparison of predictions between NeuralODE-VAE and GRU-ODE-Bayes for the same sample issued from the random $\mathbf{r}_{i}$ setting. Compared to NeuralODE-VAE, which retrieves the average dynamics of the sample, our approach detects the correlation between both features and updates its predictions more finely as the observations arrive. In particular, note that GRU-ODE-Bayes updates its prediction and confidence on a feature even when only the other one is observed, taking advantage from the fact that they are correlated. This can be seen on the left pane of Figure 1 where at time $t=3$ Dimension 1 (blue) is updated because of the observation of Dimension 2 (green).

By directly feeding sporadic inputs into the ODE, GRU-ODE-Bayes sequentially filters the hidden state and thus estimates the PDF of the future observations. This is the core strength of the proposed method, allowing it to perform long-term predictions. In contrast, NeuralODE-VAE first stores the whole dynamics in a single vector and later maps it to the dynamics of the time series (illustrated in Figure 1).

This analysis is confirmed by the performance results presented in Table 4. Our approach performs better on all setups for both NegLL and MSE. What is more, the method deals correctly with lags (i.e., the second setup) as it results in only marginal degradation of NegLL and MSE. When there is no correlation between both dimensions (i.e., $\rho=0$), the observation of one dimension contains no information on the other and this results in lower performance.

Figure 6 illustrates how GRU-ODE-Bayes updates its prediction and confidence as more and more observations are processed. This example is for the first setup (randomized $\mathbf{r}_{i}$). Initially, the predictions have large confidence intervals and reflect the general statistics of the training data. Then, observations gradually reduce the variance estimate as the model refines its predictions of the parameter $\mathbf{r}_{i}$. As more data is processed, the predictions converge to the asymptotic distribution of the underlying process.

Appendix J Application to synthetic nonlinear SDE: the Brusselator

On top of the extended multivariate OU process, we also studied a nonlinear SDE. We derived it from the Brusselator ODE, which was proposed by Ilya Prigogine to model autocatalytic reactions (Prigogine, 1982). It is a 2-dimensional process characterized by the following equations:

Where $x$ and $y$ stand for the two dimensions of the process and $a$ and $b$ are parameters of the ODE. This system becomes unstable when $b>1+a$. We add a stochastic component to this process to make it the following SDE, which we will model:

Where $dW_{1}(t)$ and $dW_{2}(t)$ are correlated Brownian motions with correlation coefficient $\rho$. We simulate 1,000 trajectories driven by the dynamics given in Eq. 7 with parameters $a=0.3$ and $b=1.4$ such that the ODE is unstable. Figure 7 show some realization of this process. The data set we use for training consists in random samples from those trajectories of length 50. We sample sporadically with an average rate of 4 samples every 10 seconds.

Figures 8 show the predictions of the trained model on different samples of the proposed stochastic Brusselator process (newly generated samples). At each point in time are displayed the means and the standard deviation of the filtered process. We stress that it means that those predictions only use the observations prior to them. Red arrows show that information is shared between both dimensions of the process. The model is able to pick up the correlation between dimensions to update its belief about one dimension when only the other is observed. The model presented in these figures used 50 dimensional latents with DOPRI solver.

Appendix K MIMIC-III: preprocessing details

MIMIC-III is a publicly available database containing deidentified health-related data associated for about 60,000 admissions of patients who stayed in critical care units of the Beth Israel Deaconess Medical Center between 2001 and 2012. To use the database, researchers must formally request access to the data via http://mimic.physionet.org.

We only take a subset of admissions for our analysis. We select them on the following criteria:

Keep only patient who are in the metavision system.

Keep only patients with single admission.

Keep only patients whose admission is longer than 48 hours, but less than 30 days.

Remove patients younger than 15 years old at admission time.

Remove patients without chart events data.

Remove patients with fewer than 50 measurements over the 48 hours. (This corresponds to measuring only half of retained variable a single time in 48 hours.)

This process restricts the data set to 21,250 patients.

K.2 Variables preprocessing

The subset of 96 variables that we use in our study are shown in Table 5. For each of those, we harmonize the units and drop the uncertain occurrences. We also remove outliers by discarding the measurements outside the 5 standard deviations interval. For models requiring binning of the time series, we map the measurements in 30-minute time bins, which gives 97 bins for 48 hours. When two observations fall in the same bin, they are either averaged or summed depending on the nature of the observation. Using the same taxonomy as in Table 5, lab measurements are averaged, while inputs, outputs, and prescriptions are summed.

This gives a total of 3,082,224 unique measurements across all patients, or an average of 145 measurements per patient over 48 hours.

Appendix L USHCN-Daily: preprocessing details

The United States Historical Climatology Network (USHCN) data set contains data from 1,218 centers scattered across the US. The data is publicly available and can be downloaded at the following address: https://cdiac.ess-dive.lbl.gov/ftp/ushcn_daily/. All states files contain daily measurements for 5 variables: precipitation, snowfall, snow depth, maximum temperature and minimum temperature.

We first remove all observations with a bad quality flag, then remove all centers that do not have observation before 1970 and after 2001. We then only keep the observations between 1950 and 2000. We subsample the remaining observations to keep on average 5% of the observations of each center. Lastly, we select the last 4 years of the kept series to be used in the analysis.

This process leads to a data set with 1,114 centers, and a total of 386,068 unique observations (or an average of 346 observations per center, sporadically spread over 4 years).

Appendix M Small-sample regime: additional results

In the main text of the paper, we presented the results for the Mean Square Error (MSE) for the different data subsets of MIMIC. In Table 6, we present the negative log-likelihood results. They further illustrate that the continuity prior embedded in our GRU-ODE-Bayes strongly helps in the small-sample regime.

Appendix N Computing Infrastructure

All models were run using a NVIDIA P100 GPU with 16GB RAM and 9 CPU cores (Intel(R) Xeon(R) Gold 6140). Implementation was done in Python, using Pytorch as autodifferentitation package. Required packages are available in the code

Appendix O Hyper-parameters used

All methods were trained using the same dimension for the hidden $\mathbf{h}$, for sake of fairness. For each fold, we tuned the following hyper-parameters using a 20% left out validation set:

Dropout rate of , $0.1$, $0.2$ and $0.3$.

Weight decay: $0.1$, $0.03$, $0.01$, $0.003$, $0.001$, $0.0001$ and .

Best model was selected using early stopping and performance were assessed by applying the best model on a held out test set (10% of the total data). The different folds were reused for each compared model for sake of reproducibility and fair comparison. We performed 5-fold cross validation and present the test performance average and standard deviation in all tables.