1 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 timeseries analysis assumes that signals are measured systematically at fixed time intervals. However, much realworld 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 timeseries 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 autoencoders (VAEs) (Kingma & Welling, 2013) or normalizing flows (Rezende et al., 2014)).
Instead of the encoderdecoder 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 GRUODEBayes method.
The tight coupling between observation processing and ODE dynamics allows the proposed method to model finegrained nonlinear dynamical interactions between the variables. As illustrated in Figure 1, GRUODEBayes 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 encoderdecoder based method NeuralODEVAE 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 H, we show that GRUODEBayes can exactly represent the corresponding FokkerPlanck dynamics in the special case of the OrnsteinUhlenbeck 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).
1.1 Problem statement
We consider the general problem of forecasting on sporadically observed dimensional time series. For example, data from patients where clinical longitudinal variables can potentially be measured. Each time series is measured at
time points specified by a vector of observation times
. The values of these observations are specified by a matrix of observations and an observation mask (to indicate which of the variables are measured at each time point).We assume that observations are sampled from the realizations of a dimensional stochastic process whose dynamics is driven by an unknown SDE:
(1) 
where is a Wiener process. The distribution of then evolves according to the celebrated FokkerPlanck equation (Risken, 1996)
. We refer to the mean and covariance parameters of its probability density function (PDF) as
and .Our goal will be to model the unknown temporal functions and from the sporadic measurements . These are obtained by sampling the random vectors at times with some observation noise . Not all dimensions are sampled each time, resulting in missing values in .
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.
2 Proposed method
At a high level, we propose a dual mode system consisting of (1) a GRUinspired continuoustime state evolution (GRUODE) that propagates in time the hidden state of the system between observations and (2) a network that updates the current hidden state to incorporate the incoming observations (GRUBayes). The system switches from propagation to update and back whenever a new observation becomes available.
We also introduce an observation model mapping
to the estimated parameters of the observations distribution
and (details in Appendix E). GRUODE then explicitly learns the FokkerPlanck dynamics of Eq. 1. This procedure allows endtoend training of the system to minimize the loss with respect to the sporadically sampled observations .2.1 GRUODE derivation
To derive the GRUbased ODE, we first show that the GRU proposed by Cho et al. (2014) can be written as a difference equation. First, let , , and be the reset gate, update gate, and update vector of the GRU:
(2)  
where is the elementwise product. Then the standard update for the hidden state of the GRU is
We can also write this as . By subtracting from this state update equation and factoring out , we obtain a difference equation
This difference equation naturally leads to the following ODE for :
(3) 
We name the resulting system GRUODE. Similarly, we derive the minimal GRUODE, 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 GRUODE input . 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 and the resulting ODE in Eq. 3 is autonomous, with and only depending on .
2.2 General properties of GRUODE
GRUODE enjoys several useful properties:
Boundedness. First, the hidden state stays within the range^{1}^{1}1We use the notation to also mean multidimensional range (i.e., all elements are within ).. This restriction is crucial for the compatibility with the GRUBayes model and comes from the negative feedback term in Eq. 3, which stabilizes the resulting system. In detail, if the th dimension of the starting state is within , then will always stay within because
This can be derived from the ranges of and in Eq. 2. Moreover, would start outside of the region, the negative feedback would quickly push into this region, making the system also robust to numerical errors.
Continuity. Second, GRUODE is Lipschitz continuous with constant . Importantly, this means that GRUODE encodes a continuity prior for the latent process . 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 smallsample regime.
General numerical integration. As a parametrized ODE, GRUODE 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 DormandPrince (an adaptive step size method). Appendix C illustrates that the DormandPrince method requires fewer time steps.
2.3 GRUBayes
GRUBayes is the module that processes the sporadically incoming observations to update the hidden vectors, and hence the estimated PDF of . This module is based on a standard GRU and thus operates in the region that is required by GRUODE. In particular, GRUBayes is able to update to any point in this region. Any adaptation is then within reach with a single observation.
To feed the GRU unit inside GRUBayes with a nonfullyobserved vector, we first preprocess it with an observation mask using , as described in Appendix D. For a given time series, the resulting update for its th observation at time with mask and hidden vector is
(4) 
where and
denote the hidden representation before and after the jump from GRUBayes update. We also investigate an alternative option where the
is updated by each observed dimension sequentially. We call this variant GRUODEBayesseq (see Appendix F for more details).2.4 GRUODEBayes
The proposed GRUODEBayes combines GRUODE and GRUBayes. The GRUODE is used to evolve the hidden state in continuous time between the observations and GRUBayes transforms the hidden state, based on the observation , from to . As best illustrated in Figure 2, the alternation between GRUODE and GRUBayes results in an ODE with jumps, where the jumps are at the locations of the observations.
Objective function
To train the model using sporadicallyobserved samples, we introduce two losses. The first loss, , is computed before the observation update and is the negative loglikelihood (NegLL) of the observations. For the observation of a single sample, we have (for readability we drop the time indexing):
where is the observation mask and are the parameters of the distribution before the update, for dimension . Thus, the error is only computed on the observed values of .
For the second loss, let denote the predicted distributions (from ) before GRUBayes. With , the PDF of given the noisy observation (with noise vector ), we first compute the analogue of the Bayesian update:
Let denote the predicted distribution (from ) after applying GRUBayes. We then define the postjump loss as the KLdivergence between and :
In this way, we force our model to learn to mimic a Bayesian update.
Similarly to the prejump loss, is computed only for the observed dimensions. The total loss is then obtained by adding both losses with a weighting parameter .
For binomial and Gaussian distributions, computing
can be done analytically. In the case of Gaussian distribution we can compute the Bayesian updated mean and variance aswhere for readability we dropped the dimension subindex. In many realworld cases, the observation noise , in which case is just the observation distribution: and .
2.5 Implementation
The pseudocode of GRUODEBayes is depicted in Algorithm 3, where a forward pass is shown for a single time series ^{2}^{2}2Code is available in the following anonymous repository : https://github.com/edebrouwer/gru_ode_bayes . For minibatching several time series we sort the observation times across all time series and for each unique time point
, 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 GRUODE step, we propagate all hidden states jointly. The GRUBayes 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).3 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 realworld 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).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ñoneroCandela & 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).
Most recently, the seminal work of Chen et al. (2018) suggested a continuous version of neural networks that overcomes the limits imposed by discretetime recurrent neural networks. Coupled with a variational autoencoder 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 discretetime RNN. Related autoencoder 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 longterm dependencies is a significant advantage.
4 Application to synthetic SDEs
Figure 1 illustrates the capabilities of our approach compared to NeuralODEVAE on data generated from a process driven by a multivariate OrnsteinUhlenbeck (OU) SDE with random parameters. Compared to NeuralODEVAE, 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 GRUODEBayes 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 , Dimension 1 (blue) is updated because of the observation of Dimension 2 (green).
By directly feeding sporadic inputs into the ODE, GRUODEBayes 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 longterm predictions.
In Appendix H
, 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
I where we present an example inspired by the Brusselator (Prigogine, 1982), a chaotic ODE.5 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 loglikelihood.
5.1 Baselines
We used a comprehensive set of stateoftheart baselines to compare the performance of our method. All models use the same hidden size representation and comparable number of parameters.
NeuralODEVAE (Chen et al., 2018). We model the time derivative of the hidden representation as a 2layer 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): GRUSimple and GRUD.
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.
TLSTM (Baytas et al., 2017). We reused the proposed timeaware LSTM cell to design a forecasting RNN with observation mask.
5.2 Electronic health records
Electronic Health Records (EHR) analysis is crucial to achieve datadriven 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 GRUODEBayes.
We use the publicly available MIMICIII 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 realvalued measurements over a period of 48 hours after patient admission. We refer the reader to Appendix J for further details on the cohort selection. We focus on the predictions of the next 3 measurements after a 36hour observation window.
5.3 Climate forecast
From shortterm weather forecast to longrange prediction or assessment of systemic changes, such as global warming, climatic data has always been a popular application for timeseries 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 K contains additional details regarding this procedure. The task is then to predict the next 3 measurements after the first 3 years of observation.
5.4 Results
We report the performance using 5fold crossvalidation. Hyperparameters (dropout and weight decay) are chosen using an inner holdout validation set (20%) (More details in Appendix
N). Performance metrics for both tasks (NegLL and MSE) are reported in Table 1. GRUODEBayes 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, GRUODEBayes unequivocally outperforms all other methods on both data sets.USHCNDaily  MIMICIII  
Model  MSE  NegLL  MSE  NegLL 
NeuralODEVAE  
NeuralODEVAEMask  
Sequential VAE  
GRUSimple  
GRUD  
TLSTM  
GRUODEBayes  

5.5 Impact of continuity prior
To illustrate the capabilities of the derived GRUODE 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 GRUODE is crucially important as it provides important prior information about the underlying process.
To highlight the importance of the GRUODE cell, we compare two versions of our model : the classical GRUODEBayes and one where the GRUODE cell is replaced by a discretized autonomous GRU. We call the latter GRUDiscretizedBayes. Table 2 shows the results for MIMICIII with varying number of patients in the training set. While our discretized version matches the continuous one on the full data set, GRUODE cell achieves higher accuracy when the number of samples is low, highlighting the importance of the continuity prior. Loglikelihood results are given in Appendix L.
Model  1,000 patients  2,000 patients  Full 

NeuralODEVAEMask 

GRUDiscretizedBayes  
GRUODEBayes 
^{2}^{2}2Statistically not different from best (pvalue with paired ttest). 


6 Conclusion and future work
We proposed a model combining two novel techniques, GRUODE and GRUBayes, 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 GRUODEBayes on both synthetic and realworld data. Moreover, while a discretized version of our model performed well on the full MIMICIII data set, the continuity prior of our ODE formulation proves particularly important in the smallsample regime, which is particularly relevant for realworld clinical data where many data sets remain relatively modest in size.
In this work, we focused on timeseries data with Gaussian observations. However, GRUODEBayes can also be extended to binomial and multinomial observations since the respective NegLL and KLdivergence are analytically tractable. This allows the modeling of sporadic observations of both discrete and continuous variables.
Acknowledgements
Edward De Brouwer is funded by a FWOSB grant. Yves Moreau is funded by (1) Research Council KU Leuven: C14/18/092 SymBioSys3; CELSAHIDUCTION, (2) Innovative Medicines Initiative: MELLODY, and (3) Flemish Government (ELIXIR Belgium, IWT, FWO 06260). 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
 Baytas et al. (2017) Baytas, I. M., Xiao, C., Zhang, X., Wang, F., Jain, A. K., and Zhou, J. Patient subtyping via timeaware lstm networks. In Proceedings of the 23rd ACM SIGKDD international conference on knowledge discovery and data mining, pp. 65–74. ACM, 2017.
 Bonilla et al. (2008) Bonilla, E. V., Chai, K. M., and Williams, C. Multitask gaussian process prediction. In Advances in neural information processing systems, pp. 153–160, 2008.
 Cao et al. (2018) Cao, W., Wang, D., Li, J., Zhou, H., Li, L., and Li, Y. Brits: Bidirectional recurrent imputation for time series. In Bengio, S., Wallach, H., Larochelle, H., Grauman, K., CesaBianchi, N., and Garnett, R. (eds.), Advances in Neural Information Processing Systems 31, pp. 6776–6786. Curran Associates, Inc., 2018.
 Che et al. (2018) Che, Z., Purushotham, S., Cho, K., Sontag, D., and Liu, Y. Recurrent neural networks for multivariate time series with missing values. Scientific reports, 8(1):6085, 2018.
 Chen et al. (2018) Chen, T. Q., Rubanova, Y., Bettencourt, J., and Duvenaud, D. Neural ordinary differential equations. In Advances in Neural Information Processing Systems, 2018, 2018.
 Cheng et al. (2017) Cheng, L.F., Darnell, G., Chivers, C., Draugelis, M. E., Li, K., and Engelhardt, B. E. Sparse multioutput gaussian processes for medical time series prediction. arXiv preprint arXiv:1703.09112, 2017.
 Cho et al. (2014) Cho, K., van Merrienboer, B., Gülçehre, Ç., Bougares, F., Schwenk, H., and Bengio, Y. Learning phrase representations using RNN encoderdecoder for statistical machine translation. CoRR, abs/1406.1078, 2014. URL http://arxiv.org/abs/1406.1078.
 Choi et al. (2016a) Choi, E., Bahadori, M. T., Schuetz, A., Stewart, W. F., and Sun, J. Doctor ai: Predicting clinical events via recurrent neural networks. In Machine Learning for Healthcare Conference, pp. 301–318, 2016a.
 Choi et al. (2016b) Choi, E., Bahadori, M. T., Sun, J., Kulas, J., Schuetz, A., and Stewart, W. Retain: An interpretable predictive model for healthcare using reverse time attention mechanism. In Advances in Neural Information Processing Systems, pp. 3504–3512, 2016b.
 Chung et al. (2014) Chung, J., Gulcehre, C., Cho, K., and Bengio, Y. Empirical evaluation of gated recurrent neural networks on sequence modeling. arXiv preprint arXiv:1412.3555, 2014.
 Du et al. (2016) Du, N., Dai, H., Trivedi, R., Upadhyay, U., GomezRodriguez, M., and Song, L. Recurrent marked temporal point processes: Embedding event history to vector. In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pp. 1555–1564. ACM, 2016.

Esteva et al. (2019)
Esteva, A., Robicquet, A., Ramsundar, B., Kuleshov, V., DePristo, M., Chou, K.,
Cui, C., Corrado, G., Thrun, S., and Dean, J.
A guide to deep learning in healthcare.
Nature medicine, 25(1):24, 2019.  Futoma et al. (2017) Futoma, J., Hariharan, S., and Heller, K. Learning to detect sepsis with a multitask gaussian process rnn classifier. arXiv preprint arXiv:1706.04152, 2017.
 Garnelo et al. (2018) Garnelo, M., Schwarz, J., Rosenbaum, D., Viola, F., Rezende, D. J., Eslami, S., and Teh, Y. W. Neural processes. arXiv preprint arXiv:1807.01622, 2018.
 Gers et al. (2000) Gers, F. A., Schmidhuber, J., and Cummins, F. Learning to forget: Continual prediction with lstm. Neural Computation, 12(10):2451–2471, 2000.
 Ghassemi et al. (2015) Ghassemi, M., Pimentel, M. A., Naumann, T., Brennan, T., Clifton, D. A., Szolovits, P., and Feng, M. A multivariate timeseries modeling approach to severity of illness assessment and forecasting in icu with sparse, heterogeneous clinical data. In AAAI, pp. 446–453, 2015.
 Goldstein et al. (2017) Goldstein, B. A., Navar, A. M., Pencina, M. J., and Ioannidis, J. Opportunities and challenges in developing risk prediction models with electronic health records data: a systematic review. Journal of the American Medical Informatics Association, 24(1):198–208, 2017.
 Jensen et al. (2014) Jensen, A. B., Moseley, P. L., Oprea, T. I., Ellesøe, S. G., Eriksson, R., Schmock, H., Jensen, P. B., Jensen, L. J., and Brunak, S. Temporal disease trajectories condensed from populationwide registry data covering 6.2 million patients. Nature communications, 5:4022, 2014.
 Johnson et al. (2016) Johnson, A. E., Pollard, T. J., Shen, L., Liwei, H. L., Feng, M., Ghassemi, M., Moody, B., Szolovits, P., Celi, L. A., and Mark, R. G. Mimiciii, a freely accessible critical care database. Scientific data, 3:160035, 2016.
 Junninen et al. (2004) Junninen, H., Niska, H., Tuppurainen, K., Ruuskanen, J., and Kolehmainen, M. Methods for imputation of missing values in air quality data sets. Atmospheric Environment, 38(18):2895–2907, 2004.
 Kingma & Welling (2013) Kingma, D. P. and Welling, M. Autoencoding variational bayes. arXiv preprint arXiv:1312.6114, 2013.
 Krishnan et al. (2015) Krishnan, R. G., Shalit, U., and Sontag, D. Deep kalman filters. arXiv preprint arXiv:1511.05121, 2015.

Krishnan et al. (2017)
Krishnan, R. G., Shalit, U., and Sontag, D.
Structured inference networks for nonlinear state space models.
In
ThirtyFirst AAAI Conference on Artificial Intelligence
, 2017.  Lee et al. (2017) Lee, C., Luo, Z., Ngiam, K. Y., Zhang, M., Zheng, K., Chen, G., Ooi, B. C., and Yip, W. L. J. Big healthcare data analytics: Challenges and applications. In Handbook of LargeScale Distributed Computing in Smart Healthcare, pp. 11–41. Springer, 2017.
 Lipton et al. (2016) Lipton, Z. C., Kale, D., and Wetzel, R. Directly modeling missing data in sequences with rnns: Improved classification of clinical time series. In Machine Learning for Healthcare Conference, pp. 253–270, 2016.
 (26) Menne, M., Williams Jr, C., and Vose, R. Longterm daily climate records from stations across the contiguous united states.
 Mitchell (1999) Mitchell, T. M. Machine learning and data mining. Communications of the ACM, 42(11):30–36, 1999.
 Moor et al. (2019) Moor, M., Horn, M., Rieck, B., Roqueiro, D., and Borgwardt, K. Temporal convolutional networks and dynamic time warping can drastically improve the early prediction of sepsis. arXiv preprint arXiv:1902.01659, 2019.
 Prigogine (1982) Prigogine, I. From being to becoming. Freeman, 1982.
 QuiñoneroCandela & Rasmussen (2005) QuiñoneroCandela, J. and Rasmussen, C. E. A unifying view of sparse approximate gaussian process regression. Journal of Machine Learning Research, 6(Dec):1939–1959, 2005.
 Rajkomar et al. (2018) Rajkomar, A., Oren, E., Chen, K., Dai, A. M., Hajaj, N., Hardt, M., Liu, P. J., Liu, X., Marcus, J., Sun, M., et al. Scalable and accurate deep learning with electronic health records. npj Digital Medicine, 1(1):18, 2018.
 Razavian & Sontag (2015) Razavian, N. and Sontag, D. Temporal convolutional neural networks for diagnosis from lab tests. arXiv preprint arXiv:1511.07938, 2015.
 Rezende et al. (2014) Rezende, D. J., Mohamed, S., and Wierstra, D. Stochastic backpropagation and approximate inference in deep generative models. In International Conference on Machine Learning, pp. 1278–1286, 2014.
 Risken (1996) Risken, H. Fokkerplanck equation. In The FokkerPlanck Equation, pp. 63–95. Springer, 1996.
 Scargle (1982) Scargle, J. D. Studies in astronomical time series analysis. iistatistical aspects of spectral analysis of unevenly spaced data. The Astrophysical Journal, 263:835–853, 1982.
 Schneider (2001) Schneider, T. Analysis of incomplete climate data: Estimation of mean values and covariance matrices and imputation of missing values. Journal of climate, 14(5):853–871, 2001.
 Soleimani et al. (2018) Soleimani, H., Hensman, J., and Saria, S. Scalable joint models for reliable uncertaintyaware event prediction. IEEE transactions on pattern analysis and machine intelligence, 40(8):1948–1963, 2018.
 Wang et al. (2006) Wang, J., Hertzmann, A., and Fleet, D. J. Gaussian process dynamical models. In Advances in neural information processing systems, pp. 1441–1448, 2006.
 Zhou et al. (2016) Zhou, G.B., Wu, J., Zhang, C.L., and Zhou, Z.H. Minimal gated unit for recurrent neural networks. International Journal of Automation and Computing, 13(3):226–234, 2016.
Appendix A Full formulation of the GRUODE cell
The full ODE equation for GRUODE is the following:
with
Matrices ,
and bias vectors
are the parameters of the cell. and are the dimension of the hidden process and of the inputs respectively.Appendix B Lipschitz continuity of GRUODE
As is differentiable and continous on , we know from the mean value theorem that for any , there exists such that
Taking the euclidean norm of the previous expression, we find
Furthermore, we showed that is bounded on . Hence, because of the bounded functions appearing in the ODE (sigmoids and hyperbolic tangents), the derivative of is itself bounded by . We conclude that is Lipschitz continuous with constant .
Appendix C Comparison of numerical integration methods
We implemented three numerical integration methods, among which the classical Euler method and the DormandPrince method (DOPRI). DOPRI is a popular adaptive step size numerical integration method relying on 2 RungeKutta 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
The preprocessing step for GRUBayes takes in the hidden state and computes the parameters for the observation PDFs . In the case of a Gaussian, contains the means and logvariances for dimension of . Then, we create a vector that concatenates with the observed value and the normalized error term, which for the Gaussian case is , where and
are the mean and standard deviation derived from
. We then multiply the vectors by a dimensionspecific weight matrixand apply a ReLU nonlinearity. Next, we zero all results that did not have an observation (by multiplying them with mask
). Finally, the concatenation of the results is fed into the GRU unit of GRUBayes.Appendix E Observation model mapping
The mapping from hidden to the parameters of the distribution and
. For this purpose we use a classical multilayer 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 GRUODEBayesseq
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, GRUBayes will update its prediction on the hidden for one input dimension after the other rather than jointly. We call this approach GRUODEBayesseq.
In this sequential approach for GRUBayes, we process onebyone all dimensions of that were observed at time 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 changes after each dimension is processed, which affects the computed and thus the resulting vector .
Appendix G Minimal GRUODE
Following the same reasoning as for the full GRU cell, we also derived the minimal GRUODE cell, based on the minimal GRU cell. The minimal GRU writes :
This can be rewritten as the following difference equation :
Which leads to the corresponding ODE :
Appendix H Application to OrnsteinUhlenbeck 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 OrnsteinUhlenbeck (OU) process with varying parameters. For a particular sample , the dynamics is given by the following SDE:
(5) 
where is a dimensional correlated Wiener process, is the vector of targets, and is the reverting strength constant. For simplicity, we consider and parameters as scalars. Each sample is then obtained via the realization of process (5) with samplespecific parameters.
h.1 Representation capabilities
We now show that our model exactly captures the dynamics of the distribution of as defined in Eq. 5. The evolution of the PDF of a diffusion process is given by the corresponding FokkerPlanck equation. For the OU process, this PDF is Gaussian with timedependent mean and covariance. Conditioned on a previous observation at time , this gives
Correlation of is constant and equal to , the correlation of the Wiener processes. The dynamics of the mean and variance parameters can be better expressed in the following ODE form:
(6) 
With initial conditions and . We next investigate how specific versions of this ODE can be represented by our GRUODEBayes.
h.1.1 Standard OrnsteinUhlenbeck process
In standard OU, the parameters , , and are fixed and identical for all samples. The ODE (6) is linear and can then be represented directly with GRUODE by storing and in the hidden state and matching the Equations (3) and (6). The OU parameters , and are learned and encoded in the weights of GRUODE. GRUBayes then updates the hidden state and stores and .
h.1.2 Generalized OrnsteinUhlenbeck processes
When parameters are allowed to vary over samples, these have to be encoded in the hidden state of GRUODEBayes, rather than in the fixed weights. For and , GRUBayes computes and stores their current estimates as the observations arrive. This is based on previous hidden and current observation as in Eq. 4. The GRUODE 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 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 OrnsteinUhlenbeck processes with sampledependent and .
Perfect representation for sample dependent requires the multiplication of inputs in Eq. 6, which GRUODE 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. GRUBayes is again able to update the hidden vector with the new estimates.
h.1.3 Nonaligned 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 GRUBayes module.
h.2 Case Study: 2D OrnsteinUhlenbeck Process
h.2.1 Setup
We evaluate our model on a 2dimensional 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, varies across samples as and . The correlation between the Wiener processes is set to . We also set and . 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
. The third setting is identical to the first but with (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 . The training set contains 10,000 samples with an average of 20 observations scattered over a 10second time interval. Models are trained with a negative loglikelihood objective function, but mean square errors (MSE) are also reported. We compare our methods to NeuralODEVAE (Chen et al., 2018). Additionally, we consider an extended version of this model where we also feed the observation mask, called NeuralODEVAEMask.
Negative LogLikelihood  MSE  
Model  Random  Random Lag  Random r  Random Lag  
NeuralODEVAEMask  
NeuralODEVAE  
GRUODEBayes  
GRUODEBayesminimal  
GRUODEBayesseq  
GRUODEBayesseqminimal  

h.2.2 Empirical evaluation
Figure 1 shows a comparison of predictions between NeuralODEVAE and GRUODEBayes for the same sample issued from the random setting. Compared to NeuralODEVAE, 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 GRUODEBayes 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 Dimension 1 (blue) is updated because of the observation of Dimension 2 (green).
By directly feeding sporadic inputs into the ODE, GRUODEBayes 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 longterm predictions. In contrast, NeuralODEVAE 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 3. 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., ), the observation of one dimension contains no information on the other and this results in lower performance.
Figure 6 illustrates how GRUODEBayes updates its prediction and confidence as more and more observations are processed. This example is for the first setup (randomized ). 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 . As more data is processed, the predictions converge to the asymptotic distribution of the underlying process.
Appendix I 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 2dimensional process characterized by the following equations:
Where and stand for the two dimensions of the process and and are parameters of the ODE. This system becomes unstable when . We add a stochastic component to this process to make it the following SDE, which we will model:
(7) 
Where and are correlated Brownian motions with correlation coefficient . We simulate 1,000 trajectories driven by the dynamics given in Eq. 7 with parameters and 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 J MIMICIII: preprocessing details
MIMICIII is a publicly available database containing deidentified healthrelated 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.
j.1 Admission/Patient cleanup
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.
j.2 Variables preprocessing
The subset of 96 variables that we use in our study are shown in Table 4
. 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 30minute 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
4, 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 K USHCNDaily: 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.essdive.lbl.gov/ftp/ushcn_daily/. All states files contain daily measurements for 5 variables: precipitation, snowfall, snow depth, maximum temperature and minimum temperature.
k.1 Cleaning and subsampling
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 L Smallsample 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 5, we present the negative loglikelihood results. They further illustrate that the continuity prior embedded in our GRUODEBayes strongly helps in the smallsample regime.
1,000 Patients  2,000 Patients  Full  
Model  NegLL  NegLL  NegLL 
NeuralODE 

GRUDiscBayes  
GRUODEBayes  

Appendix M 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 N Hyperparameters used
All methods were trained using the same dimension for the hidden , for sake of fairness. We tuned the following hyperparameters using a 20% left out validation set:
Dropout rate of , , and .
Weight decay: , , , , , and .
Learning rate : 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). We performed 5fold cross validation and present the test performance average and standard deviation in all tables.
Comments
There are no comments yet.