Year: 2020

Abstract:
The marriage of recurrent neural networks and neural ordinary differential networks (ODE-RNN) is effective in modeling irregularly sampled sequences. While ODE produces the smooth hidden states between observation intervals, the RNN will trigger a hidden state jump when a new observation arrives and thus cause the interpolation discontinuity problem. To address this issue, we propose the cubic spline smoothing compensation, which is a stand-alone module upon either the output or the hidden state of ODE-RNN and can be trained end-to-end. We derive its analytical solution and provide its theoretical interpolation error bound. Extensive experiments indicate its merits over both ODE-RNN and cubic spline interpolation.

text: The marriage of recurrent neural networks and neural ordinary differential networks (ODE-RNN) is effective in modeling irregularly sampled sequences. While ODE produces the smooth hidden states between observation intervals, the RNN will trigger a hidden state jump when a new observation arrives and thus cause the interpolation discontinuity problem. To address this issue, we propose the cubic spline smoothing compensation, which is a stand-alone module upon either the output or the hidden state of ODE-RNN and can be trained end-to-end. We derive its analytical solution and provide its theoretical interpolation error bound. Extensive experiments indicate its merits over both ODE-RNN and cubic spline interpolation.

heading: 1 INTRODUCTION
text: Recurrent neural networks (RNNs) are commonly used for modeling regularly sampled sequences (Cho et al., 2014). However, the standard RNN can only process discrete series without considering the unequal temporal intervals between sample points, making it fail to model irregularly sampled time series commonly seen in domains, e.g., healthcare (Rajkomar et al., 2018) and finance (Fagereng & Halvorsen, 2017). While some works adapt RNNs to handle such irregular scenarios, they often assume an exponential decay (either at the output or the hidden state) during the time interval between observations (Che et al., 2018; Cao et al., 2018), which may not always hold.
To remove the exponential decay assumption and better model the underlying dynamics, Chen et al. (2018) proposed to use the neural ordinary differential equation (ODE) to model the continuous dynamics of hidden states during the observation intervals. Leveraging a learnable ODE parametrized by a neural network, their method renders higher modeling capability and flexibility.
However, an ODE determines the trajectory by its initial state, and it fails to adjust the trajectory according to subsequent observations. A popular way to leverage the subsequent observations is ODE-RNN (Rubanova et al., 2019; De Brouwer et al., 2019), which updates the hidden state upon observations using an RNN, and evolves the hidden state using an ODE between observation intervals. While ODE produces smooth hidden states between observation intervals, the RNN will trigger a hidden state jump at the observation point. This inconsistency (discontinuity) is hard to reconcile, thus jeopardizing continuous time series modeling, especially for interpolation tasks (Fig. 1 top-left).
We propose a Cubic Spline Smoothing Compensation (CSSC) module to tackle the challenging discontinuity problem, and it is especially suitable for continuous time series interpolation. Our CSSC employs the cubic spline as a means of compensation for the ODE-RNN to eliminate the jump, as illustrated in Fig. 1 top-right. While the latent ODE (Rubanova et al., 2019) with an encoder-decoder structure can also produce continuous interpolation, CSSC can further ensure the interpolated curve pass strictly through the observation points. Importantly, we can derive the closed-form solution for CSSC and obtain its interpolation error bound. The error bound suggests two key factors for a good interpolation: the time interval between observations and the performance of ODE-RNN. Furthermore, we propose the hidden CSSC that aims to compensate for the hidden state of ODE-RNN (Fig. 1 bottom), which not only assuage the discontinuity problem but is more efficient when the observations are high-dimensional and only have continuity on the semantic level. We conduct extensive experiments and ablation studies to demonstrate the effectiveness of CSSC and hidden CSSC, and both of them outperform other comparison methods.

heading: 2 RELATED WORK
text: Spline interpolation is a practical way to construct smooth curves between a number of points (De Boor et al., 1978), even for unequally spaced points. Cubic spline interpolation leverages the piecewise third order polynomials to avoid the Runge’s phenomenon (Runge, 1901) and is applied as a classical way to impute missing data (Che et al., 2018).
Recent literature focuses on adapting RNNs to model the irregularly sampled time series, given their strong modeling ability. Since standard RNNs can only process discrete series without considering the unequal temporal intervals between sample points, different improvements were proposed. One solution is to augment the input with the observation mask or concatenate it with the time lag ∆t and expect the network to use interval information ∆t in an unconstrained manner (Lipton et al., 2016; Mozer et al., 2017). While such a flexible structure can achieve good performance under some circumstances (Mozer et al., 2017), a more popular way is to use prior knowledge for missing data imputation. GRU-D (Che et al., 2018) imputes missing values with the weighted sum of exponential decay of the previous observation and the empirical mean. Shukla & Marlin (2019) employs the radial basis function kernel to construct an interpolation network. Cao et al. (2018) let hidden state exponentially decay for non-observed time points and use bi-directional RNN for temporal modeling.
Another track is the probabilistic generative model. Due to the ability to model the missing data’s uncertainty, Gaussian processes (GPs) are adopted for missing data imputing (Futoma et al., 2017; Tan et al., 2020; Moor et al., 2019). However, this approach introduced several hyperparameters, such as the covariance function, making it hard to fine-tune in practice. Neural processes (Garnelo et al., 2018) eliminate such constraints by introducing a global latent variable that represents the whole process. Generative adversarial networks are also adopted for imputing (Luo et al., 2018).
Recently, neural ODEs (Chen et al., 2018) utilize a continuous state transfer function parameterized by a neural network to learn the temporal dynamics. Rubanova et al. (2019) combine the RNN and ODE to reconcile both the new observation and latent state evolution between observations. De Brouwer et al. (2019) update the ODE with the GRU structure with Bayesian inference at observations.
While the ODE produces the smooth hidden states between observation intervals, the RNN will trigger a jump of the hidden state at the observation point, leading to a discontinuous hidden state along the trajectory. This inconsistency (discontinuity) is hard to reconcile, thus jeopardizing the modeling of continuous time series, especially for interpolation tasks. The neural CDE (Kidger
et al., 2020) directly apply cubic splines interpolation at the input sequence to make the sparse input continuous and thus produce continuous output. On the contrary, our method tackles this jumping problem by introducing the cubic spline as a compensation for the vanilla ODE-RNN, at either the output or hidden space.

heading: 3 METHODS
text: In this section, we first formalize the irregularly sampled time series interpolation problem (Sec. 3.1), then introduce the background of ODE-RNN (Sec. 3.2). Based upon ODE-RNN, we present CSSC and its closed-form solution (Sec. 3.3), and illustrate the inference and training procedure (Sec. 3.4). Finally, we provide the interpolation error bound of CSSC (Sec. 3.5) and describe an useful extension of CSSC (Sec. 3.6).

heading: 3.1 PROBLEM DEFINITION
text: We focus on the interpolation task. Given an unknown underlying function x(t) : R→ Rd, t ∈ [a, b], and a set of n + 1 observations {xk|xk = x(tk)}nk=0 ∈ Rd sampled from x(t) at the irregularly spaced time points Π : a = t0 < t1 < ... < tn = b, the goal is to learn a function F (t) : R→ Rd to approximate x, such that F (tk) = xk.

heading: 3.2 BACKGROUND OF ODE-RNN
text: ODE-RNN (Rubanova et al., 2019) achieves the interpolation by applying ODE and RNN interchangeably through a time series, illustrated in top-left of Fig. 1. The function F on time interval t ∈ [tk, tk+1) is described by a neural ODE with the initial hidden state h(tk):
ḣ(t) = f(h(t)); (1) o(t) = g(h(t)), (2)
where the h ∈ Rm is the hidden embedding of the data, ḣ = dhdt is the temporal derivative of the hidden state, o ∈ Rd is the interpolation output of F (t). Here, f : Rm → Rm and g : Rm → Rd are the transfer function and the output function parameterized by two neural networks, respectively. At the observation time t = tk, the hidden state will be updated by an RNN as:
h(tk) = RNNCell(h(t − k ),xk); (3) o(tk) = g(h(tk)), (4)
where the input x ∈ Rd, t−k and t + k are the left- and right-hand limits of tk. The above formulation has two downsides. The first is the discontinuity problem: while the function described by ODE is right continuous o(tk) = o(t+k ), the RNN cell in Eq. (3) renders the hidden state discontinuity h(t−k ) 6= h(t + k ) and therefore output discontinuity o(t − k ) 6= o(t + k ). The second is that the model cannot guarantee o(tk) = xk without explicit constraints.

heading: 3.3 CUBIC SPLINE SMOOTHING COMPENSATION
text: To remedy the two downsides, we propose the module Cubic Spline Smoothing Compensation (CSSC), manifested in the top-right of Fig. 1. It computes a compensated output ô(t) as:
ô(t) = c(t) + o(t), (5)
where o(t) is the ODE-RNN output, and the c(t) is a compensation composed of piecewise continuous functions. Our key insight is that adding another continuous function to the already piecewise continuous o(t) will ensure the global continuity. For simplicity, we set c(t) as a piecewise polynomials function and then narrow it to a piecewise cubic function since it is the most commonly used polynomials for interpolation (Burden & Faires, 1997). As the cubic spline is computed for each dimension of c individually, w.l.o.g., we will discuss one dimension of the o, c, ô,x and thus denote them as o, c, ô, x, respectively. c(t) is composed with pieces as c(t) = ∑n−1 k=0 ck(t) with each piece ck defined at domain [tk, tk+1). To guarantee the smoothness, we propose four constraints to ô(t):
1. ô(t−k ) = ô(t + k ) = xk, k = 1, ..., n− 1, ô(t0) = x0, ô(tn) = xn (output continuity);
2. ˙̂o(t−k ) = ˙̂o(t + k ), k = 1, ..., n− 1 (first order output continuity);
3. ¨̂o(t−k ) = ¨̂o(t + k ), k = 1, ..., n− 1 (second order output continuity);
4. ¨̂o(t0) = ¨̂o(tn) = 0 (natural boundary condition).
The constraint 1 ensures the interpolation curves continuously pass through the observations. Constraint 2 and 3 enforce the first and second-order continuity at the observation points, which usually holds when the underline curve x is smooth. And constraint 4 specifies the natural boundary condition owing to the lack of information of the endpoints (Burden & Faires, 1997).
Given o(t) and such four constraints, c(t) has unique analytical solution expressed in Theorem 1.
Theorem 1. Given the first order and second order jump difference of ODE-RNN as
ṙk = ȯ(t + k )− ȯ(t − k ); (6) r̈k = ö(t + k )− ö(t − k ). (7)
where the analytical expression of ȯ and ö can be obtained as
ȯ = ∂g
∂h f ; ö = fᵀ
∂2g ∂h2 f + ∂g ∂h ᵀ ∂f ∂h f, (8)
and the error defined as
+k = xk − o(t + k ); (9) −k = xk − o(t − k ), (10)
then ck can be uniquely determined as
ck(t) = Mk+1 + r̈k+1 −Mk
6τk (t− tk)3 + Mk 2 (t− tk)2+
( −k+1 − + k
τk − τk(Mk+1 + r̈k+1 + 2Mk) 6 )(t− tk) + +k , (11)
where Mk is obtained as M = A−1d, (12)
A =  2 λ1 µ2 2 λ2 . . . . . .
. . . µn−2 2 λn−2
µn−1 2
 ,M =  M1 M2
... Mn−2 Mn−1
 ,d =  d1 d2 ...
dn−2 dn−1
 , (13)
τk = tk+1 − tk, µk = τk−1τk−1+τk , λk = τk τk−1+τk , dk = 6
[t+k ,t − k+1]− [t + k−1,t − k ]
τk−1+τk + 6ṙk−2r̈kτk−1−r̈k+1τkτk−1+τk ,
[t+k , t − k+1] =
−k+1− + k
τk , M0 = Mn = 0.
The proof for Theorem. 1 is in Appx. A. The c(t) is obtained by computing each c(t) individually according to Theorem. 1.
Computational Complexity. The major cost is the inverse of A, a tridiagonal matrix, whose inverse can be efficiently computed inO(n) complexity with the tridiagonal matrix algorithm (implementation detailed in Appx. C.1). Another concern is that ȯ and ö needs to compute Jacobian and Hessian in Eq. (8). We can circumvent this computing cost by computing the numerical derivative or an empirical substitution, detailed in Appx. C.2.
Model Reduction. Our CSSC can reduce to cubic spline interpolation if setting o in Eq. (5) as zero. In light of this, we further analyze our model with techniques used for cubic spline interpolation and experimentally show our advantages against it in Sec. 4.2.

heading: 3.4 INFERENCE AND TRAINING
text: For inference, firstly compute the predicted value o from ODE-RNN (Eq. (1-4)), then calculate the compensation c with CSSC(Eq. (11)); thus yielding smoothed output ô (Eq. (5)). For training, the CSSC is a standalone nonparametric module (since we have its analytical solution) on top of the ODE-RNN that allows the end-to-end training for ODE-RNN parameters. We employ Mean Squire Error (MSE) loss to supervise ô. In addition, we expect the compensation c to be small to push ODE-RNN to the leading role for interpolation and take full advantage of its model capacity. Therefore a 2-norm penalty for c is added to construct the final loss:
L = 1 N N∑ i=1 (||x(ti)− ô(ti)||2 + α||c(ti)||2) (14)
The ablation study (Sec. 4.5) shows that the balance weight α can effectively arrange the contribution of ODE-RNN and CSSC.
Gradient flow. Although c is non-parametric module, but the gradient can flow from it into the ODE-RNN because c(t) depends on the left and right limit of o(tk), ȯ(tk), ö(tk). We further analyze that ȯ(tk) plays a more important role than ö(tk) in the contribution to c(t), elaborated in Appx. C.3.

heading: 3.5 INTERPOLATION ERROR BOUND
text: With CSSC, we can even derive an interpolation error bound, which is hard to obtain for ODE-RNN. Without loss of generality, we analyze one dimension of ô, which is scalable to all dimensions. Theorem 2. Given the CSSC at the output space as Eq. (5), if x ∈ C4[a, b], f ∈ C3, g ∈ C4, then the error and the first order error are bounded as
||(x− ô)(r)||∞ ≤ Cr||(x− o)(4)||∞τ4−r, (r = 0, 1), (15) where || · ||∞ is uniform norm, (·)r is r-th derivative, C0 = 5384 , C1 = 1 24 , τ is the maximum interval over Π.
The proof of Theorem 2 is in Appx. B. The error bound guarantee the error can converge to zero if τ → 0 or ||(x − o)(4)|| → 0. This suggests that a better interpolation can come from a denser observation or a better ODE-RNN output. Interestingly, Eq. (15) can reduce to the error bound for cubic spline interpolation (Hall & Meyer, 1976) if o is set zero. Compared with ODE-RNN, which lacks the convergence guarantee for τ , our model more effectively mitigates the error for the densely sampled curve at complexity O(τ4) ; compared with the cubic spline interpolation, our error bound has an adjustable o that can leads smaller ||(x− o)(4)|| than ||x(4)||. An implicit assumption for this error bound is that the x should be 4-th order derivable; hence this model is not suitable for sharply changing signals.

heading: 3.6 EXTEND TO INTERPOLATE HIDDEN STATE
text: Although CSSC has theoretical advantages from the error bound, it is still confronted with two challenges. Consider the example of video frames: each frame is high-dimensional data, and each pixel is not continuous through the time, but the spatial movement of the content (semantic manifold) is continuous. So the first challenge is that CSSC has linear complexity w.r.t. data dimension, which is still computation demanding when the data dimension becomes very high. The second challenge is that CSSC assumes the underlying function x is continuous, and it cannot handle the discontinuous data that is continuous in its semantic manifold (e.g. video).
To further tackle these two challenges, we propose a variant of CSSC that is applied to the hidden states, named as hidden CSSC, illustrated in bottom of Fig. 1. As we only compute the compensation to hidden state, which can keep a fixed dimension regardless the how large the data dimension is, the computational complexity of hidden CSSC is weakly related to data dimension. Also, the hidden state typically encodes a meaningful low-dimensional manifold of the data. Hence smoothing the hidden state is equivalent to smoothing the semantic manifold to a certain degree.
Hence, the Eq. (5) is modified to ĥ(t) = c(t) + h(t), and the output (Eq. (2),(4)) to o(t) = g(ĥ(t)) However, since there is no groundtruth for the hidden state like the constraint 1, we assume the h(t+k )
are the knots passed through by the compensated curve, rendering the constraints 1 as ĥ(t−k ) = ĥ(t+k ) = h(t + k ). The rationality for the knots is that h(t + k ) is updated by x(tk) and thus contain more information than h(t−k ). Given the above redefined variable, the closed-form solution of c is provided in Theorem. 3. Theorem 3. Given the second order jump at hidden state as
ṙk = ḣ(t + k )− ḣ(t − k ); r̈k = ḧ(t + k )− ḧ(t − k ), (16)
where ḣ = f , ḧ = ∂f∂hf , and the error defined as
+k = h(t + k )− h(t + k ) = 0;
− k = h(t + k )− h(t − k ), (17)
the ck is uniquely determined as Eq. (11).
Theorem 3 suggests another prominent advantage: hidden CSSC can be more efficiently implemented because its computation does not involve Hessian matrix.

heading: 4 EXPERIMENTS

heading: 4.1 BASELINE
text: We compare our method with baselines: (1) Cubic Spline Interpolation (Cubic Spline) (2) A classic RNN where ∆t is concatenated to the input (RNN-∆t) (3) GRU-D (Che et al., 2018) (4) ODERNN (Rubanova et al., 2019) which is what our compensation adds upon (5) Latent-ODE (Rubanova et al., 2019) which is a VAE structure employing ODE-RNN as the encoder and ODE as decoder (6) GRU-ODE-Bayes (De Brouwer et al., 2019) which extends the ODE as a continuous GRU. Our proposed methods are denoted as CSSC for data space compensation and hidden CSSC for hidden space compensation. Implementation Details is in Appx. D.

heading: 4.2 TOY SINUOUS WAVE DATASET
text: The toy dataset is composed of 1,000 periodic trajectories with variant frequency and amplitude. Following the setting of Rubanova et al. (2019), each trajectory contains 100 irregularly-sampled time points with the initial point sampled from a standard Gaussian distribution. Fixed percentages of observations are randomly selected with the first and last time points included. The goal is to interpolate the full set of 100 points.
The interpolation error on testing data is shown in Table 1, where our method outperforms all baselines in different observation percentages. Figure 2 (a) illustrates the benefit of the CSSC over cubic spline interpolation and ODE-RNN. Cubic spline interpolation cannot interpolate the curve when observations are sparse, without learning from the dataset. ODE-RNN performs poorly when a jump occurs at the observation. However, CSSC can help eliminate such jump and guarantee smoothness at the observation time, thus yielding good interpolation. In Figure 2 (b), we visualized the ODE-RNN output o and compensation c, and the CSSC output ô in detail to demonstrate how the the bad o becomes a good ô by adding c. Finally, Fig. 2 (c) demonstrates the first dimension of hidden states before and after smoothing, where the smoothed hidden state leads to a smoothed output.

heading: 4.3 MUJOCO PHYSICS SIMULATION
text: We test our proposed method with the "hopper" model provided by DeepMind Control Suite (Tassa et al., 2018) based on MuJoCo physics engine. To increase the trajectory’s complexity, the hopper is thrown up, then rotates and freely falls to the ground (Fig. 3). We will interpolate the 7-dimensional state that describes the position of the hopper. Both our hidden CSSC and CSSC achieve improved performance from ODE-RNN, especially when observations become sparser, as Tab. 1 indicates. The visual result is shown in Fig. 3, where CSSC resemble the GT trajectory most.

heading: 4.4 MOVING MNIST
text: In addition to low-dimensional data, we further evaluate our method on high-dimensional image interpolation. Moving MNIST consists of 20-frame video sequences where two handwritten digits are drawn from MNIST and move with arbitrary velocity and direction within the 64×64 patches, with potential overlapping and bounce at the boundaries. As a matter of expediency, we use a subset of 10k videos and resize the frames into 32×32. 4 (20%) and 6 (30%) frames out of 20 are randomly observed, including the starting and ending frames. We encode the image with 2 ResBlock (He et al., 2016) into 32-d hidden vector and decode it to pixel space with a stack of transpose convolution layers. Since the pixels are only continuous at the semantic level, only hidden CSSC is evaluated with comparison methods. As shown in Tab. 1, the hidden CSSC can further improve ODE-RNN’s result, and spline interpolation behaves the worse since it can only interpolate at the pixel space, which is discontinuous through time. The visual result (Fig. 4) shows that the performance gain comes from the smoother movement and the clearer overlapping.

heading: 4.5 ABLATION STUDY
text: The effect of end-to-end training. Apart from our standard end-to-end training of CSSC, two alternative training strategies are pre-hoc and post-hoc CSSN. Pre-hoc CSSC is to train a standard CSSC but only use the ODE-RNN part when inference. On the contrary, post-hoc CSSC is to train an ODE-RNN without CSSC, but apply CSSC upon the output of ODE-RNN when inference. The
comparison of pre-hoc, post-hoc, and standard CSSC is presented in Tab. 2, where standard CSSC behaves the best. The pre-hoc CSSC is the worst because training with CSSC can tolerate the error of the ODE-RNN; thus, inference without CSSC exposes the error of ODE-RNN, and it even performs worse than standard ODE-RNN. The post-hoc CSSC can increase the performance of ODE-RNN with simple post-processing when inference. However, such performance gain is not guaranteed when the observation is sparse. For example, the post-hoc CSSC even decreases the performance of ODE-RNN in 10% observation setting. Standard CSSC has higher performance than ODE-RNN in all our experiments, indicating the importance of end-to-end training.
The effect of α. We study the effect of α ranging from 0 to 10000 given different percentages of observation on MuJoCo dataset. The performance of CSSC is quite robust to the choice of α, shown in Tab. 4 (in Appendix), especially the MSE only fluctuated from 0.000375 to 0.000463 as α ranges from 1 to 10000 in 30% observation setting. Interestingly, Fig. 5 (in Appendix) visually compares the interpolation for o, c, and ô under variant α and indicates higher α contributes to lower c, thus o is more dominant of smoothed output ô. On the other hand, the smaller α will make o less correlated with the ground truth, but the CSSC c can always make ô a well-interpolated curve.

heading: 5 DISCUSSION
text: Limitations. While the CSSC model interpolates the trajectory that can strictly cross the observation points, such interpolation is not suitable for noisy data whose observations are inaccurate. Moreover, interpolation error bound (Eq. (15)) requires the underlying data is fourth-order continuous, which indicates that CSSC is not suitable to interpolate sharply changed data, e.g., step signals.
Future work. The CSSC can not only be applied to ODE-RNN but can smooth any piecewise continuous function. Applying CSSC to other more general models is a desirable future work. Also, while interpolation for noisy data is beyond this paper’s scope, but hidden CSSC shows the potential to tolerate data noise by capturing the data continuity at the semantic space rather than the observation space, which can be a future direction.

heading: 6 CONCLUSION
text: We introduce the CSSC that can address the discontinuity issue for ODE-RNN. We have derived the analytical solution for the CSSC and even proved its error bound for the interpolation task, which is hard to obtain in pure neural network models. The CSSC combines the modeling ability of deep neural networks (ODE-RNN) and the smoothness advantage of cubic spline interpolation. Our experiments have shown the benefit of such a combination. The hidden CSSC extends the smoothness from output space to the hidden semantic space, enabling a more general format of continuous signals.

heading: A PROOF OF INTERPOLATION
text: Substitute ô with Eq. (5) we will have
ck(t + k ) + o(t + k ) = xk (18)
ck(t − k+1) + o(t − k+1) = xk+1 (19)
ċk−1(t − k ) + ȯ(t − k ) = ċk(t + k ) + ȯ(t + k ) (20) c̈k−1(t − k ) + ö(t − k ) = c̈k(t + k ) + ö(t + k ) (21)
c̈0(t + 0 ) + ö(t + 0 ) = 0 (22)
c̈n−1(t − n ) + ö(t − n ) = 0 (23)
And we let c̈k(t+k ) = Mk, k = 0, 1, ..., n − 1, and c̈n−1(t−n ) = Mn. We define the first order and second order jump difference of ODE-RNN as
ṙk = ȯ(t + k )− ȯ(t − k ); (24) r̈k = ö(t + k )− ö(t − k ). (25)
With Eq. (21), we have M̂k+1 = c̈k(t − k+1) = Mk+1 + r̈(tk+1). (26)
Using constraint Eq. (18)(19), we denote
ck(tk) = xk − o(t+k ) = + k (27) ck(tk+1) = xk+1 − o(t−t+1) = − k+1. (28)
Also denote the step size τk = tk+1 − tk. Then applying constraint Eq. (18)(19)(21) we have the piece cubic function expressed as
ck(t) = M̂k+1 −Mk
6τk (t−tk)3+ Mk 2
(t−tk)2+( −k+1 − + k
τk −τk(M̂k+1 + 2Mk) 6 )(t−tk)+ +k . (29)
Next, we try to solve all Mk. We firstly express ċk−1(t−k ) and ċk(t + k ) as
ċk−1(t − k ) = M̂k −Mk−1 2 τk − 1 +Mk−1τk−1 + −k − + k−1 τk−1 − M̂k + 2Mk−1 6 τk−1, (30)
ċk(t + k ) =
−k+1 − + k
τk − M̂k+1 + 2Mk 6 τk (31)
Applying Eq. (20) we have
2Mk+ τk−1
τk−1 + τk Mk−1+ τk τk−1 + τk
Mk+1 = 6 [t+k , t − k+1]− [t + k−1, t − k ]
τk−1 + τk + 6ṙk − 2r̈kτk−1 − r̈k+1τk τk−1 + τk
(32)
where [t+k , t − k+1] =
−k+1− + k
τk .
µk = τk−1
τk−1 + τk (33)
λk = τk
τk−1 + τk (34)
dk = 6 [t+k , t − k+1]− [t + k−1, t − k ]
τk−1 + τk + 6ṙk − 2r̈kτk−1 − r̈k+1τk τk−1 + τk . (35)
Then Mk can be obtained by solving by system of linear equations: AM = d, (36)
where
A =  2 λ1 µ2 2 λ2
. . . . . . . . . µn−2 2 λn−2
µn−1 2
 ,M =  M1 M2
... Mn−2 Mn−1
 ,d =  d1 d2 ...
dn−2 dn−1  (37) And A is non-singular since it is strict diagonally dominant matrix, hence guarantee single solution for Mk. Hence M = A−1d. (38) Now we still need to calculate ȯ(t) and ö(t).
ȯ(t) = ∂g
∂h
ᵀ ∂h
∂t = ∂g ∂h
ᵀ
f (39)
ḧ(t) = df(h(t)) dt = ∂f ∂h
ᵀ dh
dt = ∂f ∂h
ᵀ
f, (40)
ö(t) = d ∂g∂h dt
ᵀ
ḣ + ∂g
∂h
ᵀ
ḧ = ( ∂2g
∂h2 f)ᵀf +
∂g
∂h
ᵀ ∂f
∂h f (41)

heading: B PROOF FOR ERROR BOUND
text: The following proof are based on the default notation and setting for scalar time series interpolation: given a set of n+1 points {x(ti)}ni=0 irregularly spaced at time points Π : a = t0 < t1 < ... < tn = b, the goal is to approximate the underlying ground truth function x(t), t ∈ Ω = [a, b]. Let o(t) as the ODE-RNN prediction, and ô(t) = c(t) + o(t) as our smoothed output where c(t) is the compensation defined by Eq. (11). To simplify the notation, we drop the argument t for a function, e.g. x(t)→ x. We firstly introduce several lemmas, then come the the main proof for the interpolation error bound. Lemma 4. (Hall & Meyer, 1976) Let e be any function inC2(Π) ⋂n k=1 C
4(tk−1, tk) with constraints x(tk) = 0, k = 0, ..., n, and natural boundary e′′(a) = e′′(b) = 0, then
|e(r)(tk)| ≤ ρr||e(4)||∞τ4−r (k = 0, ..., n; r = 1, 2), (42) where ρ1 = 124 , ρ2 = 1 4 .
Lemma 5. Let e = x− ô, if x ∈ C4(Ω), f ∈ C3, g ∈ C4, then e ∈ C2(Π) ⋂n k=1 C 4(tk−1, tk)
Proof. Since ô ∈ C2(Π), and c ∈ ⋂n k=1 C
4(tk−1, tk), hence we only have to prove that o ∈⋂n k=1 C
4(tk−1, tk). Since o(t) in each time period (tk−1, tk) is calculated by ODE (Eq. (1), and (2)), we can express different orders of the derivative of o as:
ȯ = ∂g
∂h
ᵀ
f, (43)
ö = ( ∂2g
∂h2 f)ᵀf +
∂g
∂h
ᵀ ∂f
∂h f, (44)
o(3) = fᵀ( ∂f
∂h
∂2g ∂h2 + ∂3g ∂h3 ◦ f + 2 ∂ 2g ∂h2 ∂f ∂h )f + ∂g ∂h
ᵀ
( ∂2f ∂h2 ◦ f + ∂f ∂h ∂f ∂h )f, (45)
o(4) =fᵀ ( ∂4g ∂h4 ◦ f ◦ f + ∂ 3g ∂h3 ◦ (∂f ∂h ∂f ∂h ) + 3 ∂3g ∂h3 ◦ f ∂f ∂h + 2 ∂f ∂h ∂3g ∂h3 ◦ f + 3 ∂ 2g ∂h2 ∂2f ∂h2 ◦ f
+ 3 ∂2g ∂h2 ∂f ∂h ∂f ∂h + ∂2f ∂h2 ◦ f ∂ 2g ∂h2 + ∂f ∂h ∂f ∂h ∂2g ∂h2 + 3 ∂f ∂h ∂2g ∂h2 ∂f ∂h
) f
+ ∂g
∂h ᵀ(∂3f ∂h3 ◦ f ◦ f + ∂ 2f ∂h2 ◦ (∂f ∂h f) + 2 ∂2f ∂h2 ◦ f ∂f ∂h + ∂f ∂h ∂2f ∂h2 ◦ f + ∂f ∂h ∂f ∂h ∂f ∂h ) f,
(46)
where ∂f∂h is the Jacobian matrix since f is a multi-valued function, ∂2g ∂h2 is the Hessian matrix since o is a scalar and g is single-valued function. The higher order derivative than second order (Hessian) becomes multi-dimensional matrix which can not be mathmatically expressed for standard matrix production, so we indicate the product of a multi-dimensional matrix and a vector as ◦, which has higher computing priority then normal matrix production. From Eq. (46), o(4) depends on ∂
4g ∂h4 and
∂3f ∂h3 . Hence given g ∈ C 4 and f ∈ C3, we obtain o ∈ ⋂n k=1 C 4(tk−1, tk).
Lemma 6. (Birkhoff & Priver, 1967) Given any function v ∈ C4(tk, tk + 1), let u be the cubic Hermite interpolation matching v, we have
|(v(t)− u(t))(r)| ≤ Ar(t)||v(4)||∞τ4−r (r = 0, 1), (47)
with
Ar(t) = τ ri [(t− ti)(ti+1 − t)]2−r
r!(4− 2r)!τ4−r if t ∈ [ti, tj). (48)
This is the error bound for cubic Hermite interpolation.
Given the above lemmas, we are ready to the formal proof for Theorem 2. Proof of Theorem 2:
Proof. According to Lemma 6, we let v = x−o (since x−o ∈ C4(tk, tk+1) according to Lemma 5), let u be the cubic Hermite interpolation matching v, then we have
|(x(t)− o(t)− u(t))(r)| ≤ Ar(t)||(x− o)(4)||∞τ4−r (r = 0, 1), (49)
with
Ar(t) = τ ri [(t− ti)(ti+1 − t)]2−r
r!(4− 2r)!τ4−r if t ∈ [ti, tj). (50)
Let e = x− ô, ε = u− c, then ε is cubic Hermite interpolation implying
ε(tk) = x(tk)− o(tk)− c(tk) = x(tk)− ô(tk) = 0, (51)
and ε̇(tk) = ẋ(tk)− ȯ(tk)− ċ(tk) = ẋ(tk)− ˙̂o(tk) = ė(tk). (52)
Therefore, ε(t) in [tk, tk+1) is a cubic function with endpoints satisfying can be constructed for each interval [tk, tk+1) as Eq. (51) and (52), which can be reconstructed as
ε(t) = ė(tk)H1(t) + ė(tk+1)H2(t), (53)
where
H1(t) = (t− tk)(t− tk+1)2/τ2k , (54) H2(t) = (t− tk+1)(t− tk)2/τ2k . (55)
From Lemma 4 and Lemma 5, e is bounded as Eq. (42). Combining Eq. (42) and (53) yields: |ε(r)(t)| ≤ ρ1||e(4)||∞ ( |H(r)1 (t)|+ |H (r) 2 (t)| ) τ3 (r = 0, 1), (56)
which is rewritten as
|(x(t)− ô(t))(r)| ≤ Br(x)||e(4)||∞τ4−r = Br(x)||(x− o)(4)||∞τ4−r, (57)
with Br(t) = ρ1 ( |H(r)1 (t)|+ |H (r) 2 (t)| ) τ r−1.
Using Triangle inequality, Lemma 6, and Eq. 57, it issues
|x(r)(t)− ô(r)(t)| = |x(r)(t)− o(r)(t)− u(r)(t) + u(r)(t)− c(r)(t)| ≤ |x(r)(t)− o(r)(t)− u(r)(t)|+ |u(r)(t)− c(r)(t)| ≤ (Ar(x) +Br(x))||(x− o)(4)||∞τ4−r. (58)
Let Cr(x) = Ar(x) +Br(x), and an analysis of the optimality (Hall & Meyer, 1976) yields
C0(x) ≤ 5
384 ; C1(x) ≤
1
12 . (59)

heading: C COMPUTATIONAL COMPLEXITY
text: C.1 IMPLEMENTATION OF THE INVERSE OF MATRIX
For sake of simplicity, our Pytorch implementation adopts torch.inverse to compute A−1 in Eq. (12), which is actually the implementation of the LU composition using partial pivoting with best complexity O(n2). Its complexity is higher than the complexity of tridiagonal matrix algorithm O(n), whose implementation will be left for future work.
C.2 COMPUTATION OF ȯ AND ö
As Eq. 8 indicates, computing ȯ and ö requires the Hessian and Jacobian of g and the Jacobian of f . These Jacobians and Hessians will first participate in the inference stage and then are involved in the gradient backpropagation. However, the latest Pytorch1 does not support the computing of Jacobian and Hessian in batch, so the training process will be very slow in practice, rendering the computing cost prohibitively high. Therefore, we propose two ways to circumvent such issues from both numerical and analytical views.
1Pytorch Version 1.6.0, updated in July 2020.
Numerical Differentiation. The first solution is to use a numerical derivative. We approximate the left limitation and right limitation of ȯ(t) and ö(t) as
ȯ(t−) = o(t−)− o(t−∆t)
∆ (60)
ö(t−) = o(t−)− 2o(t−∆) + o(t− 2∆)
∆2 (61)
ȯ(t+) = o(t+ ∆t)− o(t+)
∆ (62)
ö(t+) = o(t+ 2∆t)− 2o(t+ ∆t) + o(t+)
∆2 (63)
where ∆ = 0.001. In this way, we can avoid computing Jacobian or Hessian, and the computational complexity almost remains the same as ODE-RNN. The last row of Table. 3 shows that the numerical differentiation can maintain the same performance with the analytical solution with 30% and 50% observation. However, when observations become sparser, e.g., 10%, the analytical differentiation will gain a better performance.
Analytical Approximation. For analytical solution, the major computing burden is the Hessian matrix; thus, we approximate ö with Eq. (65) by simply dropping the Hessian term. The motivation and rationality is detailed in Appx. C.3.
C.3 COMPUTATION REDUCTION FOR ANALYTICAL DERIVATIVE
To further reduce the computation for the analytical derivative of ȯ and ö, we investigate whether blocking the gradient of ȯ or ö will affect the interpolation performance. From first three rows in Table. 3, we can see that performance of blocking the gradient of ȯ is worse than that of blocking ö, indicating ȯ is more important than ö in terms of computing the compensation c. In light of this, we further drop ö in the Eq. (7), meaning the second order jump difference r̈k is zero. According to the performance shown in the fourth row of Table. 3, ö has a minor impact on the compensation c when the observation is dense. We investigate the reason by check how ö impact the computation of c and find that they are correlated by dk in Eq. (12). We write dk again here for better clarification:
dk = 6 [t+k , t − k+1]− [t + k−1, t − k ]
τk−1 + τk + 6ṙk − 2r̈kτk−1 − r̈k+1τk τk−1 + τk . (64)
From (Eq. (6, 7)) and second term of above equation (Eq. (64)), we noticed that ṙk, r̈k serve as the bridge between ȯ and ö and dk, indicating that the importance of ȯ and ö in generating compensation c can be examined by estimating the relative significance of ṙk, r̈k appear in dk. We can denote the relative significance as fraction of terms included r̈k and ṙk as s =
|2r̈kτk−1+r̈k+1τk| |6ṙk| ≤
3max{|r̈kτk−1|,|r̈k+1τk|} 6|ṙk| . Without loss of generality, we let max{|r̈kτk−1|, |r̈k+1τk|} = |r̈kτk−1|, then we have s ≤ 3max{|r̈kτk−1|,|r̈k+1τk|}6|ṙk| = |r̈kτk−1| 2|ṙk| . It shows that r̈k has a coefficient of τk−1 on the numerator, which is the time interval between two adjacent observations. In our experiment, 100 samples is uniformly sampled in 5s, and a certain percent of the samples will be selected as observations. In this setting, if we have 50% samples as observations, the average τk = 5/50 = 0.1; thus s can be informally estimated as s ≤ |r̈k|20|ṙk| , indicating ȯ is more important than ö in our experiments. As the observation ratio is higher, the τk becomes smaller, hence the relative significance of r̈k and ṙk will become even larger. Armed with the above intuition that ö is less important and the fact that he Hessian contribute the major complexity O(W 2), we drop the term with Hessian and find a better approximation of ö, leading the final approximation of ö as:
ö = ∂g
∂h
ᵀ ∂f
∂h f. (65)
Such approximation yields descent performance in practice. In addition, because o is multi-variable and g is a multi-valued function in practice, and such Jacobian needs to run in batch. We implement our Jacobian operation to tackle these difficulties and can run fast.
D IMPLEMENTATION DETAILS
The neural ODE state h has size 15. f is a 5-layer MLP with hidden state size 300. g is a 2-layer MLP with hidden state size 300. RNNCell has hidden state size 100. To obey the condition f ∈ C3 and g ∈ C4 required by the error bound in Theorem 2, we select all the nonlinear function as tanh function. The α is selected as 1000 based on the ablation study of α in Sec. 4.5. The network is optimized by AdaMax (Kingma & Ba, 2014) with β1 = 0.9, β2 = 0.999, and learning rate 0.02.

Review summaries:
The paper builds on ODE-RNN model that allows to represent a time series as a continuous trajectory. The authors address the limitation of the ODE-RNN model that the trajectory is continuous everywhere except the observation points. They introduce a compensation term based on cubic splines that transforms the output trajectory into a continuous one. This approach is also applied to correct a hidden state trajectory. The authors demonstrate better interpolation properties on sparse time series data.

This work addresses the discontinuity issues caused by jumps in hidden state/output at the arrival of new observations in a ODE-RNN. This problem is tackled by adding a cubic spline smoothing component on top of ODE-RNN to produce smooth and continuous hidden state/outputs. They derive a closed form solution for the cubic spline component based on the output of ODE-RNN and obtain an error bound for 4th order derivable inputs. Although the cubic spline component has no trainable component, this work shows that the gradient can flow through it to ODE-RNN and perform end-to-end training.