neuralode
Jupyter notebook with Pytorch implementation of Neural Ordinary Differential Equations
view repo
We introduce a new family of deep neural network models. Instead of specifying a discrete sequence of hidden layers, we parameterize the derivative of the hidden state using a neural network. The output of the network is computed using a blackbox differential equation solver. These continuousdepth models have constant memory cost, adapt their evaluation strategy to each input, and can explicitly trade numerical precision for speed. We demonstrate these properties in continuousdepth residual networks and continuoustime latent variable models. We also construct continuous normalizing flows, a generative model that can train by maximum likelihood, without partitioning or ordering the data dimensions. For training, we show how to scalably backpropagate through any ODE solver, without access to its internal operations. This allows endtoend training of ODEs within larger models.
READ FULL TEXT VIEW PDFJupyter notebook with Pytorch implementation of Neural Ordinary Differential Equations
Sample implementation of Neural Ordinary Differential Equations
Neural Ordinary Differential Equation
[IJCAI'19, NeurIPS'19] Anode: Unconditionally Accurate MemoryEfficient Gradients for Neural ODEs
Code for the NeurIPS 2019 paper: "Learning Dynamics of Attention: Human Prior for Interpretable Machine Reasoning"
Models such as residual networks, recurrent neural network decoders, and normalizing flows build complicated transformations by composing a sequence of transformations to a hidden state:
(1) 
where and . These iterative updates can be seen as an Euler discretization of a continuous transformation (Lu et al., 2017; Haber and Ruthotto, 2017; Ruthotto and Haber, 2018).
What happens as we add more layers and take smaller steps? In the limit, we parameterize the continuous dynamics of hidden units using an ordinary differential equation (ODE) specified by a neural network:
(2) 
Starting from the input layer , we can define the output layer to be the solution to this ODE initial value problem at some time . This value can be computed by a blackbox differential equation solver, which evaluates the hidden unit dynamics wherever necessary to determine the solution with the desired accuracy. Figure 1 contrasts these two approaches.
Defining and evaluating models using ODE solvers has several benefits:
In Section 2, we show how to compute gradients of a scalarvalued loss with respect to all inputs of any ODE solver, without backpropagating through the operations of the solver. Not storing any intermediate quantities of the forward pass allows us to train our models with constant memory cost as a function of depth, a major bottleneck of training deep models.
Euler’s method is perhaps the simplest method for solving ODEs. There have since been more than 120 years of development of efficient and accurate ODE solvers (Runge, 1895; Kutta, 1901; Hairer et al., 1987). Modern ODE solvers provide guarantees about the growth of approximation error, monitor the level of error, and adapt their evaluation strategy on the fly to achieve the requested level of accuracy. This allows the cost of evaluating a model to scale with problem complexity. After training, accuracy can be reduced for realtime or lowpower applications.
When the hidden unit dynamics are parameterized as a continuous function of time, the parameters of nearby “layers” are automatically tied together. In Section 3
, we show that this reduces the number of parameters required on a supervised learning task.
An unexpected sidebenefit of continuous transformations is that the change of variables formula becomes easier to compute. In Section 4, we derive this result and use it to construct a new class of invertible density models that avoids the singleunit bottleneck of normalizing flows, and can be trained directly by maximum likelihood.
Unlike recurrent neural networks, which require discretizing observation and emission intervals, continuouslydefined dynamics can naturally incorporate data which arrives at arbitrary times. In Section 5, we construct and demonstrate such a model.
The main technical difficulty in training continuousdepth networks is performing reversemode differentiation (also known as backpropagation) through the ODE solver. Differentiating through the operations of the forward pass is straightforward, but incurs a high memory cost and introduces additional numerical error.
We treat the ODE solver as a black box, and compute gradients using the adjoint sensitivity method (Pontryagin et al., 1962). This approach computes gradients by solving a second, augmented ODE backwards in time, and is applicable to all ODE solvers. This approach scales linearly with problem size, has low memory cost, and explicitly controls numerical error.
To optimize , we require gradients with respect to . The first step is to determining how the gradient of the loss depends on the hidden state at each instant. This quantity is called the adjoint
. Its dynamics are given by another ODE, which can be thought of as the instantaneous analog of the chain rule:
(4) 
We can compute by another call to an ODE solver. This solver must run backwards, starting from the initial value of . One complication is that solving this ODE requires the knowing value of along its entire trajectory. However, we can simply recompute backwards in time together with the adjoint, starting from its final value .
Computing the gradients with respect to the parameters requires evaluating a third integral, which depends on both and :
(5) 
The vectorJacobian products and in (4) and (5) can be efficiently evaluated by automatic differentiation, at a time cost similar to that of evaluating . All integrals for solving , and can be computed in a single call to an ODE solver, which concatenates the original state, the adjoint, and the other partial derivatives into a single vector. Algorithm 1 shows how to construct the necessary dynamics, and call an ODE solver to compute all gradients at once.
Most ODE solvers have the option to output the state at multiple times. When the loss depends on these intermediate states, the reversemode derivative must be broken into a sequence of separate solves, one between each consecutive pair of output times (Figure 2). At each observation, the adjoint must be adjusted in the direction of the corresponding partial derivative .
The results above extend those of Stapor et al. (2018, section 2.4.2). An extended version of Algorithm 1 including derivatives w.r.t. and can be found in Appendix C. Detailed derivations are provided in Appendix B. Appendix D provides Python code which computes all derivatives for scipy.integrate.odeint by extending the autograd
automatic differentiation package. This code also supports all higherorder derivatives. We have since released a PyTorch
(Paszke et al., 2017) implementation, including GPUbased implementations of several standard ODE solvers at github.com/rtqichen/torchdiffeq.In this section, we experimentally investigate the training of neural ODEs for supervised learning.
To solve ODE initial value problems numerically, we use the implicit Adams method implemented in LSODE and VODE and interfaced through the scipy.integrate package. Being an implicit method, it has better guarantees than explicit methods such as RungeKutta but requires solving a nonlinear optimization problem at every step. This setup makes direct backpropagation through the integrator difficult. We implement the adjoint sensitivity method in Python’s autograd framework (Maclaurin et al., 2015)
. For the experiments in this section, we evaluated the hidden state dynamics and their derivatives on the GPU using Tensorflow, which were then called from the Fortran ODE solvers, which were called from Python
autograd code.Test Error  # Params  Memory  Time  

1Layer MLP  1.60%  0.24 M     
ResNet  0.41%  0.60 M  
RKNet  0.47%  0.22 M  
ODENet  0.42%  0.22 M 
We experiment with a small residual network which downsamples the input twice then applies 6 standard residual blocks He et al. (2016b), which are replaced by an ODESolve module in the ODENet variant. We also test a network with the same architecture but where gradients are backpropagated directly through a RungeKutta integrator, referred to as RKNet. Table 1 shows test error, number of parameters, and memory cost. denotes the number of layers in the ResNet, and is the number of function evaluations that the ODE solver requests in a single forward pass, which can be interpreted as an implicit number of layers.
We find that ODENets and RKNets can achieve around the same performance as the ResNet, while using fewer parameters. For reference, a neural net with a single hidden layer of 300 units has around the same number of parameters as the ODENet and RKNet architecture that we tested.
ODE solvers can approximately ensure that the output is within a given tolerance of the true solution. Changing this tolerance changes the behavior of the network. We first verify that error can indeed be controlled in Figure 3a. The time spent by the forward call is proportional to the number of function evaluations (Figure 3b), so tuning the tolerance gives us a tradeoff between accuracy and computational cost. One could train with high accuracy, but switch to a lower accuracy at test time.
Figure 3c) shows a surprising result: the number of evaluations in the backward pass is roughly half of the forward pass. This suggests that the adjoint sensitivity method is not only more memory efficient, but also more computationally efficient than directly backpropagating through the integrator, because the latter approach will need to backprop through each function evaluation in the forward pass.
It’s not clear how to define the ‘depth‘ of an ODE solution. A related quantity is the number of evaluations of the hidden state dynamics required, a detail delegated to the ODE solver and dependent on the initial state or input. Figure 3d shows that he number of function evaluations increases throughout training, presumably adapting to increasing complexity of the model.
The discretized equation (1) also appears in normalizing flows (Rezende and Mohamed, 2015) and the NICE framework (Dinh et al., 2014)
. These methods use the change of variables theorem to compute exact changes in probability if samples are transformed through a bijective function
:(6) 
An example is the planar normalizing flow (Rezende and Mohamed, 2015):
(7) 
Generally, the main bottleneck to using the change of variables formula is computing of the determinant of the Jacobian , which has a cubic cost in either the dimension of , or the number of hidden units. Recent work explores the tradeoff between the expressiveness of normalizing flow layers and computational cost (Kingma et al., 2016; Tomczak and Welling, 2016; Berg et al., 2018).
Surprisingly, moving from a discrete set of layers to a continuous transformation simplifies the computation of the change in normalizing constant:
Let
be a finite continuous random variable with probability
dependent on time. Let be a differential equation describing a continuousintime transformation of . Assuming that is uniformly Lipschitz continuous in and continuous in , then the change in log probability also follows a differential equation,(8) 
Proof in Appendix A. Instead of the log determinant in (6), we now only require a trace operation. Also unlike standard finite flows, the differential equation does not need to be bijective, since if uniqueness is satisfied, then the entire transformation is automatically bijective.
As an example application of the instantaneous change of variables, we can examine the continuous analog of the planar flow, and its change in normalization constant:
(9) 
Given an initial distribution , we can sample from and evaluate its density by solving this combined ODE.
While is not a linear function, the trace function is, which implies . Thus if our dynamics is given by a sum of functions then the differential equation for the log density is also a sum:
(10) 
This means we can cheaply evaluate flow models having many hidden units, with a cost only linear in the number of hidden units . Evaluating such ‘wide’ flow layers using standard normalizing flows costs , meaning that standard NF architectures use many layers of only a single hidden unit.
We can specify the parameters of a flow as a function of , making the differential equation change with . This is parameterization is a kind of hypernetwork (Ha et al., 2016). We also introduce a gating mechanism for each hidden unit, where is a neural network that learns when the dynamic should be applied. We call these models continuous normalizing flows (CNF).






We first compare continuous and discrete planar flows at learning to sample from a known distribution. We show that a planar CNF with hidden units can be at least as expressive as a planar NF with layers, and sometimes much more expressive.
We configure the CNF as described above, and train for 10,000 iterations using Adam (Kingma and Ba, 2014)
. In contrast, the NF is trained for 500,000 iterations using RMSprop
(Hinton et al., 2012), as suggested by Rezende and Mohamed (2015). For this task, we minimize as the loss function where is the flow model and the target density can be evaluated. Figure 4 shows that CNF generally achieves lower loss.


Continuoustime normalizing flows are reversible, so we can train on a density estimation task and still be able to sample from the learned density efficiently.
A useful property of continuoustime normalizing flows is that we can compute the reverse transformation for about the same cost as the forward pass, which cannot be said for normalizing flows. This lets us train the flow on a density estimation task by performing maximum likelihood estimation, which maximizes where is computed using the appropriate change of variables theorem, then afterwards reverse the CNF to generate random samples from .
For this task, we use 64 hidden units for CNF, and 64 stacked onehiddenunit layers for NF. Figure 5
shows the learned dynamics. Instead of showing the initial Gaussian distribution, we display the transformed distribution after a small amount of time which shows the locations of the initial planar flows. Interestingly, to fit the Two Circles distribution, the CNF rotates the planar flows so that the particles can be evenly spread into circles. While the CNF transformations are smooth and interpretable, we find that NF transformations are very unintuitive and this model has difficulty fitting the two moons dataset in Figure
4(b).Applying neural networks to irregularlysampled data such as medical records, network traffic, or neural spiking data is difficult. Typically, observations are put into bins of fixed duration, and the latent dynamics are discretized in the same way. This leads to difficulties with missing data and illdefined latent variables. Missing data can be addressed using generative timeseries models (Álvarez and Lawrence, 2011; Futoma et al., 2017; Mei and Eisner, 2017; Soleimani et al., 2017a)
or data imputation
(Che et al., 2018). Another approach concatenates timestamp information to the input of an RNN (Choi et al., 2016; Lipton et al., 2016; Du et al., 2016; Li, 2017).We present a continuoustime, generative approach to modeling time series. Our model represents each time series by a latent trajectory. Each trajectory is determined from a local initial state, , and a global set of latent dynamics shared across all time series. Given observation times and an initial state , an ODE solver produces , which describe the latent state at each observation.We define this generative model formally through a sampling procedure:
(11)  
(12)  
(13) 
Function is a timeinvariant function that takes the value at the current time step and outputs the gradient: . We parametrize this function using a neural net. Because is timeinvariant, given any latent state , the entire latent trajectory is uniquely defined. Extrapolating this latent trajectory lets us make predictions arbitrarily far forwards or backwards in time.
We can train this latentvariable model as a variational autoencoder
(Kingma and Welling, 2014; Rezende et al., 2014), with sequencevalued observations. Our recognition net is an RNN, which consumes the data sequentially backwards in time, and outputs . A detailed algorithm can be found in Appendix E. Using ODEs as a generative model allows us to make predictions for arbitrary time points … on a continuous timeline.The fact that an observation occurred often tells us something about the latent state. For example, a patient may be more likely to take a medical test if they are sick. The rate of events can be parameterized by a function of the latent state: . Given this rate function, the likelihood of a set of independent observation times in the interval is given by an inhomogeneous Poisson process (Palm, 1943):
We can parameterize using another neural network. Conveniently, we can evaluate both the latent trajectory and the Poisson process likelihood together in a single call to an ODE solver. Figure 7 shows the event rate learned by such a model on a toy dataset.
A Poisson process likelihood on observation times can be combined with a data likelihood to jointly model all observations and the times at which they were made.
We investigate the ability of the latent ODE model to fit and extrapolate time series. The recognition network is an RNN with 25 hidden units. We use a 4dimensional latent space. We parameterize the dynamics function with a onehiddenlayer network with 20 hidden units. The decoder computing is another neural network with one hidden layer with 20 hidden units. Our baseline was a recurrent neural net with 25 hidden units trained to minimize negative Gaussian loglikelihood. We trained a second version of this RNN whose inputs were concatenated with the time difference to the next observation to aid RNN with irregular observations.
# Observations  30/100  50/100  100/100 

RNN  0.3937  0.3202  0.1813 
Latent ODE  0.1642  0.1502  0.1346 
We generated a dataset of 1000 2dimensional spirals, each starting at a different point, sampled at 100 equallyspaced timesteps. The dataset contains two types of spirals: half are clockwise while the other half counterclockwise. To make the task more realistic, we add gaussian noise to the observations.
To generate irregular timestamps, we randomly sample points from each trajectory without replacement (). We report predictive rootmeansquared error (RMSE) on 100 time points extending beyond those that were used for training. Table 2 shows that the latent ODE has substantially lower predictive RMSE.
Figure 8 shows examples of spiral reconstructions with 30 subsampled points. Reconstructions from the latent ODE were obtained by sampling from the posterior over latent trajectories and decoding it to dataspace. Examples with varying number of time points are shown in Appendix F. We observed that reconstructions and extrapolations are consistent with the ground truth regardless of number of observed points and despite the noise.


Figure 7(c) shows latent trajectories projected onto the first two dimensions of the latent space. The trajectories form two separate clusters of trajectories, one decoding to clockwise spirals, the other to counterclockwise. Figure 9 shows that the latent trajectories change smoothly as a function of the initial point , switching from a clockwise to a counterclockwise spiral.
The use of minibatches is less straightforward than for standard neural networks. One can still batch together evaluations through the ODE solver by concatenating the states of each batch element together, creating a combined ODE with dimension . In some cases, controlling error on all batch elements together might require evaluating the combined system times more often than if each system was solved individually. However, in practice the number of evaluations did not increase substantially when using minibatches.
When do continuous dynamics have a unique solution? Picard’s existence theorem (Coddington and Levinson, 1955) states that the solution to an initial value problem exists and is unique if the differential equation is uniformly Lipschitz continuous in and continuous in . This theorem holds for our model if the neural network has finite weights and uses Lipshitz nonlinearities, such as tanh or relu.
Our framework allows the user to trade off speed for precision, but requires the user to choose an error tolerance on both the forward and reverse passes during training. For sequence modeling, the default value of 1.5e8 was used. In the classification and density estimation experiments, we were able to reduce the tolerance to 1e3 and 1e5, respectively, without degrading performance.
Reconstructing the state trajectory by running the dynamics backwards can introduce extra numerical error if the reconstructed trajectory diverges from the original. This problem can be addressed by checkpointing: storing intermediate values of on the forward pass, and reconstructing the exact forward trajectory by reintegrating from those points. We did not find this to be a practical problem, and we informally checked that reversing many layers of continuous normalizing flows with default tolerances recovered the initial states.
The use of the adjoint method for training continuoustime neural networks was previously proposed (LeCun et al., 1988; Pearlmutter, 1995), though was not demonstrated practically. The interpretation of residual networks He et al. (2016a) as approximate ODE solvers spurred research into exploiting reversibility and approximate computation in ResNets (Chang et al., 2017; Lu et al., 2017). We demonstrate these same properties in more generality by directly using an ODE solver.
One can adapt computation time by training secondary neural networks to choose the number of evaluations of recurrent or residual networks (Graves, 2016; Jernite et al., 2016; Figurnov et al., 2017; Chang et al., 2018). However, this introduces overhead both at training and test time, and extra parameters that need to be fit. In contrast, ODE solvers offer wellstudied, computationally cheap, and generalizable rules for adapting the amount of computation.
Recent work developed reversible versions of residual networks (Gomez et al., 2017; Haber and Ruthotto, 2017; Chang et al., 2017), which gives the same constant memory advantage as our approach. However, these methods require restricted architectures, which partition the hidden units. Our approach does not have these restrictions.
Much recent work has proposed learning differential equations from data. One can train feedforward or recurrent neural networks to approximate a differential equation (Raissi and Karniadakis, 2018; Raissi et al., 2018a; Long et al., 2017), with applications such as fluid simulation (Wiewel et al., 2018). There is also significant work on connecting Gaussian Processes (GPs) and ODE solvers (Schober et al., 2014). GPs have been adapted to fit differential equations (Raissi et al., 2018b) and can naturally model continuoustime effects and interventions (Soleimani et al., 2017b; Schulam and Saria, 2017). Ryder et al. (2018) use stochastic variational inference to recover the solution of a given stochastic differential equation.
The dolfin library (Farrell et al., 2013) implements adjoint computation for general ODE and PDE solutions, but only by backpropagating through the individual operations of the forward solver. The Stan library (Carpenter et al., 2015) implements gradient estimation through ODE solutions using forward sensitivity analysis. However, forward sensitivity analysis is quadratictime in the number of variables, whereas the adjoint sensitivity analysis is linear (Carpenter et al., 2015; Zhang and Sandu, 2014). Melicher et al. (2017) used the adjoint method to train bespoke latent dynamic models.
In contrast, by providing a generic vectorJacobian product, we allow an ODE solver to be trained endtoend with any other differentiable model components. While use of vectorJacobian products for solving the adjoint method has been explored in optimal control (Andersson, 2013; Andersson et al., In Press, 2018), we highlight the potential of a general integration of blackbox ODE solvers into automatic differentiation (Baydin et al., 2018)
for deep learning and generative modeling.
We investigated the use of blackbox ODE solvers as a model component, developing new models for timeseries modeling, supervised learning, and density estimation. These models are evaluated adaptively, and allow explicit control of the tradeoff between computation speed and accuracy. Finally, we derived an instantaneous version of the change of variables formula, and developed continuoustime normalizing flows, which can scale to large layer sizes.
We thank Wenyi Wang and Geoff Roeder for help with proofs, and Daniel Duckworth, Ethan Fetaya, Hossein Soleimani, Eldad Haber, Ken Caluwaerts, Daniel FlamShepherd, and Harry Braviner for feedback. We thank Chris Rackauckas, Dougal Maclaurin, and Matthew James Johnson for helpful discussions.
Journal of Machine Learning Research
, 12(May):1459–1500, 2011.Learning to Detect Sepsis with a Multitask Gaussian Process RNN Classifier.
ArXiv eprints, 2017.Proceedings of the IEEE conference on computer vision and pattern recognition
, pages 770–778, 2016a.Hidden physics models: Machine learning of nonlinear partial differential equations.
Journal of Computational Physics, pages 125–141, 2018.Let be a finite continuous random variable with probability dependent on time. Let be a differential equation describing a continuousintime transformation of . Assuming that is uniformly Lipschitz continuous in and continuous in , then the change in log probability also follows a differential equation:
To prove this theorem, we take the infinitesimal limit of finite changes of through time. First we denote the transformation of over an change in time as
(14) 
We assume that is Lipschitz continuous in and continuous in , so every initial value problem has a unique solution by Picard’s existence theorem. We also assume is bounded. These conditions imply that , , and are all bounded. In the following, we use these conditions to exchange limits and products.
We can write the differential equation using the discrete change of variables formula, and the definition of the derivative:
(15)  
(16)  
(17)  
(18)  
(19)  
(20) 
The derivative of the determinant can be expressed using Jacobi’s formula, which gives
(21)  
(22)  
(23) 
Substituting with its Taylor series expansion and taking the limit, we complete the proof.
(24)  
(25)  
(26)  
(27) 
∎
Let , then . Since the trace of an outer product is the inner product, we have
(28) 
This is the parameterization we use in all of our experiments.
The continuous analog of NICE (Dinh et al., 2014) is a Hamiltonian flow, which splits the data into two equal partitions and is a volumepreserving transformation, implying that . We can verify this. Let
(29) 
Then because the Jacobian is all zeros on its diagonal, the trace is zero. This is a volumepreserving flow.
The FokkerPlanck equation is a wellknown partial differential equation (PDE) that describes the probability density function of a stochastic differential equation as it changes with time. We relate the instantaneous change of variables to the special case of FokkerPlanck with zero diffusion, the Liouville equation.
As with the instantaneous change of variables, let evolve through time following . Then Liouville equation describes the change in density of –a fixed point in space–as a PDE,
(30) 
However, (30) cannot be easily used as it requires the partial derivatives of , which is typically approximated using finite difference. This type of PDE has its own literature on efficient and accurate simulation (Stam, 1999).
Instead of evaluating at a fixed point, if we follow the trajectory of a particle , we obtain
(31) 
We arrive at the instantaneous change of variables by taking the log,
(32) 
While still a PDE, (32) can be combined with to form an ODE of size ,
(33) 
Compared to the FokkerPlanck and Liouville equations, the instantaneous change of variables is of more practical impact as it can be numerically solved much more easily, requiring an extra state of for following the trajectory of . Whereas an approach based on finite difference approximation of the Liouville equation would require a grid size that is exponential in .
We present an alternative proof to the adjoint method (Pontryagin et al., 1962) that is short and easy to follow.
Let follow the differential equation , where are the parameters. We will prove that if we define an adjoint state
(34) 
then it follows the differential equation
(35) 
For ease of notation, we denote vectors as row vectors, whereas the main text uses column vectors.
The adjoint state is the gradient with respect to the hidden state at a specified time . In standard neural networks, the gradient of a hidden layer depends on the gradient from the next layer by chain rule
(36) 
With a continuous hidden state, we can write the transformation after an change in time as
(37) 
and chain rule can also be applied
(38) 
The proof of (35) follows from the definition of derivative:
(39)  
(40)  
(41)  
(42)  
(43)  
(44)  
(45) 
We pointed out the similarity between adjoint method and backpropagation (eq. 38). Similarly to backpropagation, ODE for the adjoint state needs to be solved backwards in time. We specify the constraint on the last time point, which is simply the gradient of the loss wrt the last time point, and can obtain the gradients with respect to the hidden state at any time, including the initial value.
(46) 
Here we assumed that loss function depends only on the last time point . If function depends also on intermediate time points , etc., we can repeat the adjoint step for each of the intervals , in the backward order and sum up the obtained gradients.
We can generalize (35) to obtain gradients with respect to –a constant wrt. –and and the initial and end times, and . We view and as states with constant differential equations and write
(47) 
We can then combine these with to form an augmented state^{1}^{1}1Note that we’ve overloaded to be both a part of the state and the (dummy) independent variable. The distinction is clear given context, so we keep as the independent variable for consistency with the rest of the text. with corresponding differential equation and adjoint state,
(48) 
Note this formulates the augmented ODE as an autonomous (timeinvariant) ODE, but the derivations in the previous section still hold as this is a special case of a timevariant ODE. The Jacobian of has the form
(49) 
where each is a matrix of zeros with the appropriate dimensions. We plug this into (35) to obtain
(50) 
The first element is the adjoint differential equation (35), as expected. The second element can be used to obtain the total gradient with respect to the parameters, by integrating over the full interval and setting .
(51) 
Finally, we also get gradients with respect to and , the start and end of the integration interval.
(52) 
Between (35), (46), (51), and (52) we have gradients for all possible inputs to an initial value problem solver.
This more detailed version of Algorithm 1 includes gradients with respect to the start and end times of integration.
To obtain the latent representation , we traverse the sequence using RNN and obtain parameters of distribution . The algorithm follows a standard VAE algorithm with an RNN variational posterior and an ODESolve model:
Run an RNN encoder through the time series and infer the parameters for a posterior over :