Neural Differential Recurrent Neural Network with Adaptive Time Steps
The neural Ordinary Differential Equation (ODE) model has shown success in
learning complex continuous-time processes from observations on discrete time
stamps. In this work, we consider the modeling and forecasting of time series
data that are non-stationary and may have sharp changes like spikes. We propose
an RNN-based model, called RNN-ODE-Adap, that uses a neural ODE to represent
the time development of the hidden states, and we adaptively select time steps
based on the steepness of changes of the data over time so as to train the
model more efficiently for the "spike-like" time series. Theoretically,
RNN-ODE-Adap yields provably a consistent estimation of the intensity function
for the Hawkes-type time series data. We also provide an approximation analysis
of the RNN-ODE model showing the benefit of adaptive steps. The proposed model
is demonstrated to achieve higher prediction accuracy with reduced
computational cost on simulated dynamic system data and point process data and
on a real electrocardiography dataset.
Yixuan Tan, Liyan Xie, Xiuyuan Cheng
# Neural Differential Recurrent Neural Network with Adaptive Time Steps

###### Abstract

The neural Ordinary Differential Equation (ODE) model has shown success in learning complex continuous-time processes from observations on discrete time stamps. In this work, we consider the modeling and forecasting of time series data that are non-stationary and may have sharp changes like spikes. We propose an RNN-based model, called _RNN-ODE-Adap_, that uses a neural ODE to represent the time development of the hidden states, and we adaptively select time steps based on the steepness of changes of the data over time so as to train the model more efficiently for the "spike-like" time series. Theoretically, _RNN-ODE-Adap_ yields provably a consistent estimation of the intensity function for the Hawkes-type time series data. We also provide an approximation analysis of the RNN-ODE model showing the benefit of adaptive steps. The proposed model is demonstrated to achieve higher prediction accuracy with reduced computational cost on simulated dynamic system data and point process data and on a real electrocardiography dataset.

## 1 Introduction

We consider the modeling and forecasting of time series characterized by _irregular_ time steps and _non-stationary_ patterns, which are commonly observed in various applications, such as finance [10] and healthcare [3]. We treat the data as a sequence of ordered observations from an unknown underlying continuous-time process, sampled at discrete time points. Recurrent Neural Network (RNN) [29] is frequently employed to model such sequential data. This work proposes to use an RNN-based model with a neural Ordinary Differential Equation (ODE) to fit time series data.

The classical neural ODE type approaches typically assume _regular_ time grids in the data sequences, or the _same_ irregular time grids across different sequences [27]. When dealing with highly non-stationary data, such as a time series with sudden spikes, it becomes imperative to select sufficiently small time steps to ensure accurate modeling of regions with steep changes. These regions require more refined time steps, especially the ones with abrupt spikes. However, it is common for most of the time horizon to observe a slow-varying and less steep time series, which is more "flat" over time; thus, the refined time steps would result in unnecessarily high computational costs. Some examples of time series with abrupt spikes are shown in Figure 1, illustrating a spectrum ranging from continuous time series to discontinuous time series such as counting processes.

To train the neural ODE for data with non-stationary patterns more effectively, we propose an approach that employs _adaptive time steps_ in the neural ODE model, which we refer to as _RNN-ODE-Adap_. The model adaptively selects the time steps based on the local variation of the time series,
enabling it to capture underlying trends with potentially fewer steps. Our numerical experiments showed that, compared to other baseline models that use regular time steps, _RNN-ODE-Adap_ could achieve higher prediction accuracy with similar or lower time complexity. The contribution of the work is as follows.

* Based on a neural ODE model characterizing the dynamics of hidden states, we propose an algorithm to construct adaptive time steps, which assigns refined time steps to data around "spikes" while using rough time steps for data in "flat" segments. This can significantly reduce the computational cost in the training process with little impact on the modeling performance.
* We provide theoretical insights into the consistency of the model using the example of Hawkes process type data, and the approximation guarantee of the RNN-ODE model that illuminates the benefits of adaptive time steps.
* We conduct numerical experiments on both synthetic data and a real-world time-series data set to demonstrate the advantage of the proposed algorithm in terms of both modeling accuracy and computational efficiency.

### Related Works

Neural ODE.Our work is closely related to the neural ODE [2] model, which parameterizes the derivative of the hidden state using a neural network. In [2], a generative time-series model was proposed, which takes the neural ODE as the decoder. Furthermore, [27] proposed a non-generative model with continuous-time hidden dynamics to handle irregularly sampled data based on [2]. Compared with existing works related to neural ODE [27, 35, 37, 5, 22, 18, 23, 12, 11], we model the ODE that determines the progression of hidden states by including the data itself in the derivative of the hidden state. In contrast to existing works on non-stationary environments such as the piecewise-constant ODE [11], our work proposes to use adaptive time steps to automatically adapt to sparse spikes in the time series, without pre-defining the time period for each piece of ODE.

Neural CDE.We note that the Neural Controlled Differential Equation (CDE) [19] also incorporates the observations into the model continuously. Specifically, the hidden states in [19] follow the CDE \(h(t)=h(t_{0})+\int_{t_{0}}^{t}f_{\theta}(h(s))\mathrm{d}X_{s}\), where the integral is a Riemann-Stieltjes integral. We would

Figure 1: Illustration of spike-like time series. The crosses denote the discretely sampled time steps, which can be irregular. In the left panel, the subsequences enclosed with the orange, yellow, and green brackets represent the (training or testing) windows generated from this sequence.
like to emphasize some key differences between model (2) and Neural CDE. The \(X_{s}\) in Neural CDE is the natural cubic spline of \(\{(x(t_{i}),t_{i})\}_{i}\), and \(f_{\theta}:\mathbb{R}^{d_{h}}\to\mathbb{R}^{d_{h}\times(D+1)}\), where \(d_{h}\) is the number of hidden units and \(D\) is the data dimension. Thus, for the same number of hidden units, Neural CDE requires a more complex parameterized \(f_{\theta}\) to model \(h(t)\). Moreover, since \(X_{s}\) is obtained by cubic spline, it is less naturally adapted to the prediction task that requires extrapolation to the time stamps not seen when computing the spline. Therefore, it is hard to evaluate the prediction performance of Neural CDE and thus we defer the evaluation under the Neural CDE setting for future work.

Continuous-Time RNNs.Our model belongs to the extensive family of continuous-time RNNs, originating from [26]. Several existing studies explore various RNN architectures, such as [1, 17, 6, 30, 16] These RNN models leverage their structures to address the exploding and vanishing gradient problem. Our model also adopts a continuous-time ODE framework for time series data, and the proposed adaptive time stamp selection method can be viewed as effectively reducing the length of the discrete sequence when a significant part of the process is changing slowly. Meanwhile, our approach can also be used concurrently with the methodologies such as in [6]. As the focus of our work is to model the "spike-like" time series data, the combination of our model and the existing continuous-time RNN models can further improve the efficiency when applied to such data.

Time Adaptivity.Previous studies have investigated the incorporation of time adaptivity in continuous-time RNNs, such as GACTRNN [13], TARN [16], and LEM [31]. In these works, time adaptivity was incorporated by multiplying the ODE with an adaptively learned time modulator, usually parametrized by another sub-network. In contrast, our method adaptively selects time steps during the preprocessing phase, where the selection process only utilizes the steepness of change of the time series data. Therefore, the proposed model does not involve the training of a sub-network for the time modulator as in the previous models, which may incur an increase in model size and additional computational costs.

## 2 Problem Setup

### Training Data and Prediction Task

Consider a random continuous time series \(x(t)\in\mathbb{R}^{D}\) over the time horizon \([0,T]\) for some \(T\in\mathbb{R}^{+}\). We observe multiple independent and identically distributed samples of the continuous process \(x(t)\), where each sample is sampled at discrete time stamps, which can vary across different samples. We split the observed sequences into training and testing sequences. From the training sequences, we generate a total of \(K^{(\text{Tr})}\) training windows, denoted as \(\{\mathbf{x}^{(\text{Tr},k)}\}_{k=1}^{K^{(\text{Tr})}}\), each of window length \(N\) as our training data (see the left panel of Figure 2 for an illustrative example). Here \(\mathbf{x}^{(\text{Tr},k)}\coloneqq\{x^{(\text{Tr},k)}(t_{1}^{(\text{Tr},k)}), \ldots,x^{(\text{Tr},k)}(t_{N}^{(\text{Tr},k)})\}\) and \(0<t_{1}^{(\text{Tr},k)}<\cdots<t_{N}^{(\text{Tr},k)}\leq T\) are the corresponding time stamps for the \(k\)-th training window. Similarly, we create \(K^{(\text{Te})}\) testing windows of length \(N\) from the testing sequences, denoted as \(\{\mathbf{x}^{(\text{Te},k)}\}_{k=1}^{K^{(\text{Te})}}\). To simplify the notation, we may drop the superscripts and write \(\{x(t_{1}),x(t_{2}),\ldots,x(t_{N})\}\) for a given training window if it does not cause confusion.

Our goal is to make predictions based on historical data. Given a historical series \(x(t_{1}),\ldots,x(t_{n})\), we aim to perform either one-step or multi-step predictions. The one-step prediction involves predicting \(x(t_{n+1})\) at a single future time \(t_{n+1}\) based on \(\{x(t_{1})\ldots,x(t_{n})\}\), while the \(m\)-step prediction includes forecasting \(\{x(t_{n+1}),\ldots,x(t_{n+m})\}\) at future times \(t_{n+1}<\cdots<t_{n+m}\). The detailed formulas for measuring the prediction accuracy are provided in Appendix B.1.
### Training Objective

Given the time horizon \([0,T]\), recall that we have a collection of training windows with the \(k\)-th one denoted as \(\mathbf{x}^{(\mathrm{Tr},k)}=\{x^{(\mathrm{Tr},k)}(t_{1}^{(\mathrm{Tr},k)}),\ldots,x^ {(\mathrm{Tr},k)}(t_{N}^{(\mathrm{Tr},k)})\}\). Here, the time steps are allowed to be _heterogeneous_ for different training windows. We train the model parametrized by neural networks (see Section 3.1 for the neural ODE model adopted in this paper) with trainable parameters \(\Theta\) using the _mean-squared regression loss_ function

\[\mathcal{L}(\Theta;\{\mathbf{x}^{(\mathrm{Tr},k)}\}_{k=1}^{K^{(\mathrm{Tr})}})= \sum_{k=1}^{K^{(\mathrm{Tr})}}\sum_{i=1}^{N}\|\hat{x}^{(\mathrm{Tr},k)}(t_{i}^ {(\mathrm{Tr},k)})-x^{(\mathrm{Tr},k)}(t_{i}^{(\mathrm{Tr},k)})\|^{2}|t_{i}^{( \mathrm{Tr},k)}-t_{i-1}^{(\mathrm{Tr},k)}|, \tag{1}\]

where \(\Theta\) are the network parameters, \(\hat{x}^{(k)}(t)\) is the output of the neural ODE model under parameter \(\Theta\) conditioned on all past observation, and \(t_{0}^{(\mathrm{Tr},k)}\) is the added initial time stamp for each window.

The time difference term \(|t_{i}^{(\mathrm{Tr},k)}-t_{i-1}^{(\mathrm{Tr},k)}|\) ensures that the empirical mean-squared error loss (1) matches the \(\ell_{2}\) loss for function estimation. This term will be important to balance the fitting errors among time intervals with different time steps. This term would be necessary for the proposed scheme with adaptive (non-uniform) time steps. In our numerical examples, we also performed an ablation study regarding this term to demonstrate its necessity; see Figure A9 in Appendix B.5 for an example.

## 3 Method

We state the neural ODE model used for the hidden dynamic in Section 3.1. The algorithm for adaptive steps is introduced in Section 3.2, and the computational complexity is explained in Section 3.3. More implementation details, such as evaluation metrics and choice of thresholds, are given in Appendix B.1.

### Neural ODE for RNN model

To be able to model the observation \(x(t)\) as a function of a hidden value \(h(t)\), we follow the previous continuous-time RNN neural-ODE approach [1, 6] to model the hidden dynamics of \(h(t)\) as

\[h^{\prime}(t)=f(h(t),x(t);\theta_{h}), \tag{2}\]

where \(f\) is a neural network parameterized by \(\theta_{h}\). If one directly adopts the neural ODE model [2] to the hidden state \(h(t)\), the ODE model would be \(h^{\prime}(t)=f(h(t),t;\theta)\) without the observed time series data \(x(t)\). In contrast, the model (2) incorporates the observed incoming time series data \(x(t)\) as an input to \(f\), which is important for modeling the time series data especially when the underlying dynamics is non-stationary. The time evolution of the observed series \(x(t)\) is modeled by an output neural network \(g\) that maps the hidden value \(h(t)\) to \(x(t)\) as

\[\hat{x}(t)=g(h(t);\theta_{d}), \tag{3}\]

where \(g\) is called the output function parameterized by \(\theta_{d}\).

Given the neural network functions \(f\) and \(g\) (which generally can adopt any architecture) and the observed time series data \(x(t)\), from any initial input \(h(0)\), we can numerically solve the RNN neural ODE model (2) to obtain the \(h\) values at any time \(t\in(0,T)\) as \(h(t)=h(0)+\int_{0}^{t}f(h(s),x(s);\theta_{h})ds\), and then predict the value of \(x(t)\) by \(\hat{x}(t)=g(h(t);\theta_{d})\). The neural ODE integration can be solved
by existing first-order or higher-order schemes, and the back-propagation can be computed by the adjoint method [24, 2]. If one uses the forward Euler scheme, the discrete-time dynamic of \(h(t)\) (after incorporating the time step into the network function \(f\)) becomes \(h_{i+1}=h_{i}+f(h_{i},\theta_{i})\), which recovers the structure of Residual networks [28, 14]. In this work, we adopt the forward Euler scheme in experiments due to its better stability than higher-order schemes when the dynamic has steep changes. Our methodology of adaptive time grids can potentially be extended to higher-order differential schemes.

### Adaptive Time Steps

We propose to learn a neural-ODE RNN model using adaptive time stamps, and thus the method is called _RNN-ODE-Adap_. The construction of adaptive time steps is summarized in Algorithm 1. The intuition behind the proposed algorithm is to assign longer (rough) time intervals during time regions where the time series is slowly time-varying (such as "flat" curves), while assigning shorter (fine) time intervals during those regions with "spikes" (highly non-stationary and fast time-varying regimes). For constructing the adaptive time stamps, we assume the initial time grid is sufficiently fine and adopt a dyadic-partition type algorithm to be detailed as follows.

Given a raw (discrete-time) training window \(x(t_{0}),x(t_{1}),\ldots,x(t_{N})\) sampled at the finest level of the time stamps \(0\leq t_{0}<\cdots<t_{N}\leq T\). For simplicity, below we write it as \(x_{0},x_{1},x_{2},\ldots,x_{N}\). Without loss of generality, we assume \(N\) is a power of two. We first define a monitor function \(M(\cdot)\) that measures the variation of the sub-sequence \(\{x_{i},\ldots,x_{j}\}\), \(i<j\). In this paper, we mainly adopt the _maximum variation_ defined as

\[M(\{x_{i},\ldots,x_{j}\}):=\max_{i+1\leq k\leq j}\frac{\|x_{k}-x_{k-1}\|_{2}} {|t_{k}-t_{k-1}|}, \tag{4}\]

which captures the maximum variation among any two adjacent time stamps. Here we may also choose \(\ell_{p}\) norms for any \(p\geq 1\).

We then screen from the finest level of time grids and adaptively merge neighboring time grids if their maximum variation is below a pre-specified threshold \(\epsilon>0\). In detail, for the first level \(l=1\), we group the original \(N\) time intervals into \(N/2\) sub-intervals (as demonstrated in Figure 2) and each

Figure 2: Illustration of adaptive time steps resulted from Algorithm 1. In this example, \(N=8\) and \(L=2\); three samples are removed in phase \(l=1\), and one sample is removed in phase \(l=2\).
sub-interval contains three consecutive time stamps: \(\{x_{0},x_{1},x_{2}\},\{x_{2},x_{3},x_{4}\},\ldots,\{x_{N-2},x_{N-1},x_{N}\}\). Then we calculate the maximum variation \(M(x_{0},x_{1},x_{2}),\ldots,M(x_{N-2},x_{N-1},x_{N})\) for each sub-interval. We then merge the two consecutive time intervals into one, i.e., remove the middle time stamp \(x_{2n+1}\), for \(n=0,1,\ldots,N/2-1\), if

\[M(x_{2n},x_{2n+1},x_{2n+2})<\epsilon.\]

In other words, we only keep the time stamps on which the maximum variation exceeds \(\epsilon\).

The above selection procedure is repeated similarly for \(l=2,3,\ldots\) until a pre-specified maximum integer \(L\). The value \(L\) corresponds to the roughest time interval. The algorithm is detailed in Algorithm 1, in which we maintain a set \(\mathcal{D}\) that characterizes which time stamps to be removed. Meanwhile, we also keep a Flag vector in each round as an indicator of whether the midpoint time stamp was removed in the _last_ round and Flag\({}_{\texttt{new}}\) indicating whether the middle time stamp will be removed in the _current_ round. The elements in the Flag vector equals \(0\) for intervals with slight variation (\(M(\cdot)\leq\epsilon\)) and \(1\) otherwise. The primary usage of such Flag vector is that for two consecutive intervals in round \(l^{\prime}\), we only merge the intervals if both of them are slow-varying intervals (i.e., merged from a previous round \(l<l^{\prime}\)).

```
1:Input: Data series \(\{x_{0},x_{1},x_{2},\ldots,x_{N}\}\); threshold \(\epsilon>0\); \(L\in\mathbb{Z}_{+}\).
2:Initialize: \(\mathcal{D}=\emptyset\). A flag vector Flag = \(\{0,0,\ldots,0\}\) of length \(N\).
3:for\(l=1\)to\(L\)do
4: Define a new flag: Flag\({}_{\texttt{new}}=\{0,0,\ldots,0\}\) of length \(\lfloor N/2^{l}\rfloor\).
5:for\(i=1\)to\(\lfloor N/2^{l}\rfloor\)do
6:if Flag\([2(i-1)+1]=\texttt{Flag}[2i]=0\)then
7: Compute the monitoring function \(M(\{x_{2^{l}(i-1)},x_{2^{l}(i-1)+2^{l-1}},x_{2^{l}i}\})\).
8:if\(M<\epsilon\)then
9:\(\mathcal{D}=\mathcal{D}\cup(2^{l}(i-1)+2^{l-1})\).
10:else
11: Mark Flag\({}_{\texttt{new}}[i]=1\).
12:endif
13:else
14: Mark Flag\({}_{\texttt{new}}[i]=1\).
15:endif
16:endfor
17: Update Flag = Flag\({}_{\texttt{new}}\).
18:endfor
19:Output: Indexes of removed time steps \(\mathcal{D}\).
```

**Algorithm 1** A dyadic algorithm for selecting adaptive time steps.

The output of Algorithm 1 is the set \(\mathcal{D}\) of time stamps to be removed. The model is trained on the remaining time steps only. We provide an illustration in Figure 2 of the algorithm for selecting adaptive time steps. From the final results in Figure 2, it can be seen that the output of the adaptive time steps uses longer time steps to model stationary periods (from \(i=0\) to \(4\) and from \(i=6\) to \(8\)), and uses _shorter_ time steps to model _spikes_ (from \(i=4\) to \(6\)) in the sequence.

### Computational Complexity

The computational complexity of applying Algorithm 1 in the preprocessing stage to \(K^{\rm(Tr)}\) training windows is \(O(K^{\rm(Tr)}ND)\), where \(D\) is the data dimension. For the neural ODE model described as in
(2)-(3), when \(f\) possesses the same network structure as a vanilla RNN with \(d_{h}\) hidden units and \(g\) is a one-layer fully connected network, the complexity in the training process is \(O(n_{e}K^{\rm(Tr)}\bar{N}_{a}d_{h}(d_{h}+D))\), where \(n_{e}\) and \(\bar{N}_{a}\) represent the number of training epochs and the average length of the adaptive windows, respectively.

Since the computational cost in the training process usually dominates that in the preprocessing step (which happens as long as \(n_{e}d_{h}\geq 2^{L}\)), the overall complexity of the _RNN-ODE-Adap_ model is \(O(n_{e}K^{\rm(Tr)}\bar{N}_{a}d_{h}(d_{h}+D))\). This is of the same order as the complexity of training a vanilla RNN (we refer to Appendix B.2 for the specific structure) with \(d_{h}\) hidden units in \(n_{e}\) epochs, using \(K^{\rm(Tr)}\) training windows with the same length \(\bar{N}_{a}\). Therefore, compared with the complexity when training with the original finest \(N\) time grids, the complexity associated with the adaptive method will be reduced by a factor of \(\bar{N}_{a}/N\). The smallest achievable complexity will be reduced by a factor of \(1/2^{L}\) when choosing a sufficiently large threshold \(\epsilon\).

## 4 Theory

In this section, we provide the recovery consistency of the training objective (Section 4.1) and approximation error guarantee of the RNN-ODE model revealing the benefit of adaptive step size (Section 4.2). All proofs can be found in Appendix A.

### Function Estimation for Event-type Data

We present the theoretical analysis for function estimation based on the proposed model under counting-type time series. It is worthwhile noting that counting-type time series represent a special class of continuous-time models since they exemplify the _extreme_ case of "spike-like" data, where we have discontinuities from zero to one, as shown at the right end of Figure 1.

For event-type sequences, the raw data contains a list of event times \(0<t_{1}<t_{2}<\ldots<t_{n}<T\) on the time horizon \([0,T]\). Each timestamp is the time when an event happens. In practice, the estimation is performed on _discrete-time grids_. Define the counting process \(N(t):=\sum_{i=1}^{n}\mathbf{1}(t_{i}\leq t)\) as the total number of events happened before time \(t\). We convert such continuous-time data into discrete observations by discretizing the time interval \([0,T]\) into \(M\) intervals of equal length \(\Delta t=T/M\), and then let \(x_{m}=N(m\Delta t)-N((m-1)\Delta t)\), \(m=0,1,\ldots,M\) (by convention \(x_{0}=0\)). When \(\Delta t\) is chosen sufficiently small, it becomes the Bernoulli process where \(x_{i}\in\{0,1\}\).

We consider the temporal Hawkes processes [25], in which the values \(x_{i}\) are mostly zero under mild assumptions, corresponding to sparse "spikes". Such temporal Hawkes processes can be characterized by its _conditional intensity_ function defined as

\[\lambda^{*}(t)=\lim_{\Delta\to 0}\Delta^{-1}\mathbb{E}[N(t+\Delta)-N(t)| \mathcal{F}_{t}],\]

where the filtration \(\mathcal{F}_{t}\) stands for the information available up to time \(t\). In the case of Hawkes processes, \(\lambda(t)=\mu+\alpha\int_{0}^{t}\phi(t-s)dN(s)\) is simply a linear function of past jumps of the process, where \(\phi(\cdot)\) is the influence kernel. For example, under the special case of exponential kernels, the intensity function becomes \(\lambda^{*}(t)=\mu+\alpha\beta\int_{0}^{t}e^{-\beta(t-\tau)}dN(\tau)\).

The intensity function recovery consistency by minimizing least-square population loss is proved in Theorem 4.2 under a memory constraint. We parameterize the function by a neural network (NN) based structure characterized as in (2)-(3). We define the prototypical network architecture below.

**Definition 4.1**.: _Define the function class NN-ODE\((d_{\rm out},L_{h},p_{h},L_{d},p_{d})\) as_

\[\begin{split}\text{NN-ODE}(d_{\rm out},L_{h},p_{h},L_{d},p_{d})& \coloneqq\{F:\mathbb{R}\mapsto\mathbb{R}^{d_{\rm out}}|F(t)=g(h(t)),h^{ \prime}(t)=f(h(t),x(t)),\\ g&\text{ is NN with $L_{d}$ layers and max-width $p_{d},h$ is NN with $L_{h}$ layers and max width $p_{h}.$}\}\end{split} \tag{5}\]
**Theorem 4.2**.: _Assume there exist \(d\) buffer time steps with samples \(x_{-d},\ldots,x_{-1}\) prior to the Hawkes count data \(\{x_{0},\ldots,x_{M}\}\) and each time step has duration \(\Delta t=T/M\). We further assume NN-ODE\((d_{\mathrm{out}},L_{h},p_{h},L_{d},p_{d})\) is rich enough to model the true intensity function. Then the minimizer \(F^{*}\) to the population loss function_

\[\Psi(F)\coloneqq\sum_{m=1}^{M}\mathbb{E}\big{[}(x_{m}-F(m\Delta t)\Delta t)^{2 }|x_{m-d}\ldots x_{m-1}\big{]},\]

_optimized within the neural network class \(F\in\text{NN-ODE}(D,L_{h},p_{h},L_{d},p_{d})\), satisfies \(F^{*}(m\Delta t)=\tilde{\lambda}(m):=\frac{1}{\Delta t}\int_{(m-1)\Delta t}^{m \Delta t}\lambda^{*}(t)dt\), which is the discretized intensity._

_Remark 4.3_.: We have the recovered intensity function \(F^{*}(m\Delta t)=\tilde{\lambda}(m)\) and is extendable to the entire time horizon as \(F^{*}(t)=F^{*}(m\Delta t)\mathbf{1}\{(m-1)\Delta t<t\leq m\Delta t\}\) for any \(t\in[0,T]\). In Appendix A.1 we extend the analysis to show that under the asymptotic scenario when \(M\to\infty\), we have that \(\int_{0}^{T}|F^{*}(t)-\lambda^{*}(t)|dt\to 0\).

The above argument is made under the population loss, showing that using the least-square loss function can recover the actual intensity function for the Hawkes process. It is mainly due to the generality of the ODE model (2) and (3), which is consistent with Hawkes process and most time series models. The argument may be extended to empirical processes by utilizing the empirical concentration of the process.

### Approximation Analysis of RNN-ODE-Adap

For theoretical generality, in this subsection, we consider the continuous-time process \(y(t)\in\mathbb{R}^{D^{\prime}}\) satisfying

\[h^{\prime}(t)=f(h(t),x(t)),\quad y(t)=g(h(t)),\quad h(0)=h_{0},\quad t\in[0,T], \tag{6}\]

where \(x(t)\in\mathbb{R}^{D}\) is the observable input data and \(h(t)\in\mathbb{R}^{d_{h}}\) is the underlying hidden process from some initial value \(h_{0}\). Taking \(y(t)\) to be \(x(t)\) reduces the model to the case (2)(3) considered in the other parts of the work. We provide two theorems: Theorem 4.6 proves the uniform approximation to \(y(t)\) by continuous-time RNN-ODE model without time discretization; Theorem 4.9 further takes into account the discrete-time scheme and obtains the approximation on a time grid.

Approximation of the continuous-time model.We will use neural network functions \(f_{\theta}\) and \(g_{\phi}\) to approximate the functions \(f\) and \(g\), respectively, see Lemma 4.5. Given \(x(t)\) on \([0,T]\), let \(h_{\mathrm{NN}}(t)\) be the solution to the hidden-process ODE \(h^{\prime}_{\mathrm{NN}}(t)=f_{\theta}(h_{\mathrm{NN}}(t),x(t))\) from \(h_{\mathrm{NN}}(0)=h_{0}\). This leads to the output process \(y_{\mathrm{NN}}(t)\) defined by

\[h^{\prime}_{\mathrm{NN}}(t)=f_{\theta}(h_{\mathrm{NN}}(t),x(t)),\quad y_{ \mathrm{NN}}(t)=g_{\phi}(h_{\mathrm{NN}}(t)),\quad h_{\mathrm{NN}}(0)=h_{0}, \quad t\in[0,T]. \tag{7}\]

The approximation of \(y_{\mathrm{NN}}(t)\) to \(y(t)\) will be based on the approximation of \(f_{\theta}\) and \(g_{\phi}\), which calls for the regularity condition of the system (6). We take the following technical conditions.

**Assumption 4.4**.:
1. _The observed process_ \(x:[0,T]\to[-1,1]^{D}\) _and is Lipschitz continuous over_ \(t\)_; The hidden process_ \(h:[0,T]\to[-1,1]^{d_{h}}\)_._
2. \(f:[-1.1,1.1]^{d_{h}}\times[-1,1]^{D}\to\mathbb{R}^{d_{h}},(\eta,x)\mapsto f( \eta,x)\)_, and is Lipschitz continuous with respect to both_ \(\eta\) _and_ \(x\)_._
3. \(g:[-1.1,1.1]^{d_{h}}\to[-1,1]^{D^{\prime}},\eta\mapsto g(\eta)\) _is Lipschitz continuous._
We let \(L_{g}\) denote the global Lipschitz constant of \(g\) on \([-1.1,1.1]^{d_{h}}\). For \(f\), both global and local Lipschitz constants on the domain \([-1.1,1.1]^{d_{h}}\times[-1,1]^{D}\) are used. More detailed definitions of these constants will be introduced in Lemma 4.5 (for the global constant) and Theorem 4.6 (for the local constant).

The next lemma directly follows by applying [36] to the case where \(f\) and \(g\) have 1st-order regularity (Lipschitz continuity). The proof is given in appendix A.2.

**Lemma 4.5**.: _For any \(\epsilon_{f},\epsilon_{g}>0\), there exist neural networks \(f_{\theta},g_{\phi}\) such that_

\[\max_{\eta\in[-1.1,1.1]^{d_{h}},x\in[-1,1]^{D}}\left\|f(\eta,x)-f_{\theta}( \eta,x)\right\|_{2}<\epsilon_{f},\max_{\eta\in[-1.1,1.1]^{d_{h}}}\left\|g( \eta)-g_{\phi}(\eta)\right\|_{2}<\epsilon_{g}, \tag{8}\]

_and_

* \(f_{\theta}\) _has_ \(O(\ln\frac{C_{f}}{\epsilon_{f}}+\ln d_{h}+1)\) _layers and_ \(O((C_{f}/\epsilon_{f})^{d_{h}+D}(\ln\frac{C_{f}}{\epsilon_{f}}+\ln d_{h}+1))\) _trainable parameters._
* \(g_{\phi}\) _has_ \(O(\ln\frac{C_{g}}{\epsilon_{g}}+\ln D^{\prime}+1)\) _layers and_ \(O((C_{g}/\epsilon_{g})^{d_{h}}(\ln\frac{C_{g}}{\epsilon_{g}}+\ln D^{\prime}+1))\) _trainable parameters._

_The constants in big-\(O\) may depend on \(D,D^{\prime}\), and \(d_{h}\). Here \(C_{f}\coloneqq\max\{L^{f,h},L^{f,x},M_{f}\}\), where \(M_{f}=\sup_{(\eta,x)\in[-1.1,1.1]^{d_{h}}\times[-1,1]^{D}}\left\|f(\eta,x)\right\|\) and \(L^{f,h},L^{f,x}\) are denote the Lipschitz constant of \(f\) on \([-1.1,1.1]^{d_{h}}\times[-1,1]^{D}\) (see formal definitions in (A1) in the proof of Lemma 4.5 in Appendix A.2). \(C_{g}\coloneqq\max\{L_{g},M_{g}\}\), and \(M_{g}=\sup_{\eta\in[-1.1,1.1]^{d_{h}}}\left\|g(\eta)\right\|\)._

For the spike-like data, the majority of the regions are slow-varying, with the spikes occupying only a minor part of the whole interval \([0,T]\). Thus, the whole interval \([0,T]\) may be partitioned into two disjoint sets \(\mathcal{D}_{1}\) and \(\mathcal{D}_{2}\), each of which consisting of unions of disjoint intervals in \([0,T]\). To characterize this partition more precisely, we define the constants related to an interval in \([0,T]\) as follows:

* For an interval \([s,t]\subset[0,T]\), we define the domains \(B^{h},B^{x}\) as \[B^{h}:=(h([s,t])+B_{r}^{d_{h}})\subset[-1.1,1.1]^{d_{h}},\quad B^{x}:=(x([s,t] )+B_{r}^{D})\cap[-1,1]^{D},\] (9) with \(r=0.1\), and \(B_{r}^{d_{h}},B_{r}^{D}\) represent balls with radius \(r\) in \(\mathbb{R}^{d_{h}},\mathbb{R}^{D}\) respectively (see Figure 3 for illustration). Here, \(h([s,t])+B_{r}^{d_{h}}\) means the Minkowski addition, namely \(\{h_{1}+h_{2},h_{1}\in h([s,t]),h_{2}\in B_{r}^{d_{h}}\}\), and \(x([s,t])+B_{r}^{D}\) is defined in the same way. Then, we denote \[L_{[s,t]}^{f,h} \coloneqq\sup_{x\in B^{x}}\sup_{\eta,\eta_{2}\in B^{h}}\frac{ \left\|f(\eta_{1},x)-f(\eta_{2},x)\right\|}{\|\eta_{1}-\eta_{2}\|},\] \[L_{[s,t]}^{f,x} \coloneqq\sup_{h\in B^{h}}\sup_{x_{1},x_{2}\in B^{x}}\frac{\left\| f(\eta,x_{1})-f(\eta,x_{2})\right\|}{\|x_{1}-x_{2}\|},\] (10) as the local Lipschitz constants of \(f\) within the domain \(B^{h}\times B^{x}\), and \[M_{[s,t]}^{f}\coloneqq\sup_{(\eta,x)\in B^{h}\times B^{x}}\|f(\eta,x)\|_{2}.\] (11)

With the local Lipschitz constant defined as above, we suppose that any time grid \([s_{1},t_{1}]\) in \(\mathcal{D}_{1}\) corresponds to a local Lipschitz constant \(L_{[s_{1},t_{1}]}^{f,h}\leq L_{\mathrm{low}}\). On contrast, if \([s_{2},t_{2}]\) belongs to \(\mathcal{D}_{2}\), the local Lipschitz constant \(L_{\mathrm{low}}<L_{[s_{2},t_{2}]}^{f,h}\leq L_{\mathrm{high}}(\leq L^{f,h})\). Here, \(\mathcal{D}_{1}\) is comprised of regions with slow
variations, while \(\mathcal{D}_{2}\) encompasses regions with sharp changes, as demonstrated in Figure 3. It may often be the case that \(|D_{1}|\) is greater than \(|D_{2}|\). Then, we define

\[L^{\rm(avg)}\coloneqq\frac{1}{T}(L_{\rm low}|D_{1}|+L_{\rm high}|D_{2}|). \tag{12}\]

Following Lemma 4.5 and the partition described above, Theorem 4.6 below provides the approximation results for the continuous-time process \(y(t)\) using (7).

**Theorem 4.6**.: _Under Assumption 4.4 and for \(L^{\rm(avg)}\) defined as in (12), suppose \(\epsilon_{f},\epsilon_{g}>0\) and \(\epsilon_{f}\) satisfies_

\[Te^{L^{\rm(avg)}T}\epsilon_{f}<0.1, \tag{13}\]

_and let \(f_{\theta},g_{\phi}\) be the neural networks satisfying (8) (the model complexity is bounded as in Lemma 4.5), then_

\[\max_{t\in[0,T]}\|y(t)-y_{\text{NN}}(t)\|<\epsilon_{g}+L_{g}Te^{L^{\rm(avg)}T} \epsilon_{f}. \tag{14}\]

_Remark 4.7_ (Interpretation of \(L^{\rm(avg)}\) and local Lipschitz constants).: (14) can provide an improved bound because when the data have sharp changes, \(L_{\rm high}\) (as the \(\infty\)-norm of the Lipschitz constant over time) can be large while \(L^{\rm(avg)}=(L_{\rm low}|D_{1}|+L_{\rm high}|D_{2}|)/T\) (as certain \(L^{1}\)-norm of the Liptshictz constant over time) may stay at a smaller value. The partition \(\mathcal{D}_{1}\cup\mathcal{D}_{2}\) reflects how adaptively choosing grids may help improve the theoretical results, and this will be further explored in the next subsection. Therein \(x(t)\) will only be observed at a discrete time grid, which can be adaptively chosen according to the local Lipschitz constants (see Theorem 4.9 for more details).

_Remark 4.8_ (Arbitrary desired accuracy in (14)).: For any \(\varepsilon>0\), we can choose

\[\epsilon_{f}<\frac{1}{T\exp(L^{\rm(avg)}T)}\min\{0.1,\frac{\varepsilon}{2L_{g }}\},\quad\epsilon_{g}<\frac{\varepsilon}{2},\]

then the right-hand side of (14) is bounded by \(\varepsilon\).

Approximation under time discretization.We assume that \(x(t)\) is only observed at discrete time grids \(\{t_{i}\}_{i=1}^{N}\) instead of on the whole interval \([0,T]\). Below, the time grids can be chosen adaptively, which will be detailed in Remark 4.10. Given the time grids \(\{t_{i}\}_{i=1}^{N}\), we define the following constants that will be used in the theorems later:

Figure 3: Demonstration of the domains \(B^{x}\) and \(B^{h}\) defined as in (9) for the time interval \([s,t]\) (here \(d_{h}=2,D=1\)). The domains \(\mathcal{D}_{1}\) and \(\mathcal{D}_{2}\) that correspond to slowly and fast varying regions are colored in orange and blue respectively.
* By (A2), for each \(i\), let \(L_{i}^{f,h},L_{i}^{f,x}\) and \(M_{i}^{f}\) be defined as in (10) and (11) respectively, where we take the interval \([s,t]\) as \([t_{i-1},t_{i}]\).
* By (A1), for each \(i\), let \(L_{i}^{x}\) be the Lipschitz constant of \(x(t)\) on \(t\in[t_{i-1},t_{i}]\). \(i=1,\ldots,N+1\) (we follow the convention that \(t_{0}=0,t_{N+1}=T\)).

For \(\Delta t_{i}\coloneqq t_{i}-t_{i-1}\) and \(\hat{h}_{\text{NN}}(0)=h_{0}\), the forward Euler scheme is applied on (7) as follows:

\[\hat{h}_{\text{NN}}(t_{i})=\hat{h}_{\text{NN}}(t_{i-1})+\Delta t_{i}f_{\theta} (\hat{h}_{\text{NN}}(t_{i-1}),x(t_{i-1})),\quad\hat{y}_{\text{NN}}(t_{i})=g_{ \phi}(\hat{h}_{\text{NN}}(t_{i})),\quad i=1,\ldots,N. \tag{15}\]

Compared to Theorem 4.6, Theorem 4.9 below additionally accounts for the discretization error from the numerical integration, providing an upper bound of the approximation error using \(x(t)\) observed at discrete time grids. Theorem 4.9 focuses on the forward Euler method, and we refer to Remark 4.11 for its extension to higher-order schemes.

**Theorem 4.9**.: _Under Assumption 4.4 and given a time grid \(\{t_{i}\}_{i=1}^{N}\) on \([0,T]\) at which \(x(t)\) is observed. Suppose \(\epsilon_{f},\epsilon_{g}>0\) and \(\epsilon_{f},\Delta t_{j}\) satisfy_

\[T\exp(\sum_{i=1}^{N}L_{i}^{f,h}\Delta t_{i})\left(\epsilon_{f}+\max_{j}\{\mu_{ j}\Delta t_{j}\}\right)<0.1, \tag{16}\]

_where_

\[\mu_{j}\coloneqq L_{j}^{f,h}M_{j}^{f}+L_{j}^{f,x}L_{j}^{x},\]

_and let \(f_{\theta},g_{\phi}\) be the neural networks satisfying (8) (the model complexity is bounded as in Lemma 4.5), then_

\[\max_{i}\|y(t_{i})-\hat{y}_{\text{NN}}(t_{i})\|\leq\epsilon_{g}+L_{g}T\exp( \sum_{i=1}^{N}L_{i}^{f,h}\Delta t_{i})\left(\epsilon_{f}+\max_{j}\{\mu_{j} \Delta t_{j}\}\right). \tag{17}\]

The condition (16) in Theorem 4.9 is imposed to guarantee that the numerically integrated hidden states \(\{\hat{h}_{\text{NN}}(t_{i})\}\) belong to \([-1.1,1.1]^{d_{h}}\), so that the approximation results in Lemma 4.5 are applicable.

_Remark 4.10_ (Arbitrary desired accuracy in (17)).: For any \(\varepsilon>0\), suppose the time grids satisfy that

\[\max_{j}\{\mu_{j}\Delta t_{j}\}\}<\frac{1}{T\exp(\sum_{i=1}^{N}L_{i}^{f,h} \Delta t_{i})}\min\{0.05,\ \frac{\varepsilon}{3L_{g}}\}, \tag{18}\]

then we can choose

\[\epsilon_{f}<\frac{1}{T\exp(\sum_{i=1}^{N}L_{i}^{f,h}\Delta t_{i})}\min\{0.05,\ \frac{\varepsilon}{3L_{g}}\},\quad\epsilon_{g}<\frac{\varepsilon}{3},\]

to make the right-hand side of (17) bounded by \(\varepsilon\).

_Remark 4.11_ (Extension to higher-order integration schemes).: The numerical integration scheme (15) can be extended to the multi-step explicit methods of higher orders (e.g. Runge-Kutta methods), given that the time grid selection appropriately fulfills the requirements of the integration scheme. For example, we may choose \(t_{i+1}-t_{i}=t_{i}-t_{i-1}\) for adjacent sub-intervals \([t_{i-1},t_{i}],[t_{i},t_{i+1}]\) to apply the commonly used RK4 method.

Theorem 4.9 provides insights into the utility of the adaptive steps for improving the model fitting performance, which is reflected in the last term, involving \(\max_{j}\{\mu_{j}\Delta t_{j}\}\), in Eq. (17). Specifically, time grids may be selected such that \(\Delta t_{i}\) is small if \(L_{i}^{x}\) is great, indicating a steep change in \(x(t)\) for \(t\in[t_{i-1},t_{i}]\). On the contrary, when the variation in \(x(t)\) is smaller, we employ larger \(\Delta t_{i}\) to reduce the total number of required time grids.
## 5 Numerical Experiments

We validate the performance of the proposed method using three types of datasets: the simulated spiral series, the simulated event data, and a real ECG dataset. We present the complexity vs. accuracy tradeoff curves for different methods, and demonstrate the advantage of the _RNN-ODE-Adap_ method. The code can be found at [https://github.com/Yixuan-Tan/RNN_ODE_Adap](https://github.com/Yixuan-Tan/RNN_ODE_Adap).

### Trained Models

In this section, we examine and report the performance of two models, _RNN-ODE_ and _RNN-ODE-Adap_. Both models are trained to minimize the loss function as in (1), and the difference lies in the choice of the time grid. Specifically, _RNN-ODE_ is trained using regular (non-adaptive) time steps, and _RNN-ODE-Adap_ is trained with adaptively selected time steps by Algorithm 1. The architecture of both models is the same as vanilla RNN; see more details in Appendix B.2.

As explained in Section 3.3, the computational cost of training the models is proportional to the average length of the training windows. Hence, in this section, when we compare the performance under varying complexities, the "complexities" are discussed in terms of averaged "numbers of grids" of the training data. More details of the experiment settings are in Appendix B, with the boxplots of the error plots and additional results provided in Appendix B.5.

### Simulated Spiral Data

We first investigate the capability of our method to fit and capture the underlying dynamics of the simulated spiral data. For a given matrix \(A\in\mathbb{R}^{2\times 2}\), one spiral is generated by integrating the ODE

\[x^{\prime}(t)=f(x(t))=\|x(t)\|^{-2}Ax(t), \tag{19}\]

over the time span \([0,T]\), with the initial value \(x(0)=x_{0}\in\mathbb{R}^{2}\). The initial training and testing windows are of length \(N=64\), corresponding to the largest complexity shown in Figure 4.

We first compare the on-sample prediction performance with RNN and LSTM [15] under different computational complexities. For each testing window, we use the first half as available historical

Figure 4: Comparison of the MSE prediction errors on the simulated spiral data generated from Eq. (19) for RNN, LSTM, RNN-ODE, RNN-ODE-Adap. The \(x\)-axis represents the average length of the training windows, which reflects the complexity of the models (see Section 3.3).
data and perform predictions for the second half. Figure 4 shows the averaged Mean Squared Errors (MSEs) computed as in Eq. (A6) for the models, with varying complexities. From the trade-off curves in Figure 4, we observe that when the training cost is relatively low, namely the training windows are short, all the models underfit and do not learn the dynamics well, resulting in poor prediction performance. When the training cost increases and the training windows consist of more time steps, the models overfit, and prediction accuracy worsens.

Moreover, we observe from Figure 4 that for the same complexity, RNN-ODE significantly improves the forecasting performance compared to the vanilla RNN, especially when the number of grids is not too small so that the models begin to learn the dynamics well. Additionally, RNN-ODE-Adap further achieves smaller prediction errors than RNN-ODE since it selects data points more informatively with the same number of grids. Finally, we note that while LSTM performs best in most cases, it possesses a more complex network structure. We refer to Appendix B.5 for additional results about the LSTM and Lipschitz-RNN [6] variants of the adaptive model.

Figure 5 shows examples of spiral reconstructions using about 30% of the data. The time steps might be obtained by interpolation and regular (the upper half of Figure 5), or be chosen adaptively by Algorithm 1 and thus irregular (the lower half of Figure 5). We can see that there are mismatches between shapes reconstructed by RNN and LSTM and the ground truth spiral shape. In contrast, we note that RNN-ODE and RNN-ODE-Adap are consistent with the underlying spirals.

### Simulated Point-Process Data

We further apply our method to a simulated example of event times data generated from temporal Hawkes processes [25] as described in Section 4. We train the model (2)-(3) and estimate the true intensity function \(\lambda(t)\) using the output \(x(t)\). The mean squared loss \(L=\int_{0}^{T}(dN(t)/dt-\lambda(t))^{2}dt\) is used when fitting the neural ODE model.

The fitting errors of the four models versus the complexity are shown below in Figure 6 (left). It can be observed that for all the models, the fitting errors decrease as the training complexity

Figure 5: Comparison of RNN, LSTM, RNN-ODE, RNN-ODE-Adap on the reconstruction of simulated spiral data generated from Eq. (19), using regular (upper) v.s. irregular (lower) time series.
increases. Furthermore, RNN-ODE-Adap achieves the best fitting performance. The right panel shows the log-log plot of RNN-ODE and RNN-ODE-Adap, from which we can see more clearly that for fixed model complexity (network structure), the proposed model approaches the true intensity function.

Figure 6 (right) shows two examples of fitting performance. In this example, all models use 33 grids on average. Thus, the complexity is 50% of the largest one. It can be observed that RNN fails to capture the smooth decay of the kernel. Furthermore, we can see that with the adaptive choice of time steps which is more refined, RNN-ODE-Adap can learn the dynamics of the intensity function much better - it can estimate the "jump" in the intensity accurately.

### Real Data: ECG Time Series

We validate the proposed RNN-ODE-Adap on one public electrocardiography (ECG) dataset PTB-XL [34, 9]. We focus on learning the underlying dynamics of ECG signals using the RNN-ODE model and use adaptive time steps for "spikes" in data series. We remark that windows of the highest sampling rate are chosen to have \(N=96\) time grids, in which usually two cycles are contained. In this way, the prediction of the second half given the first half would be more meaningful.

Figure 7 (left) shows the on-sample prediction MSEs of the four methods for two different prediction lengths, 24 and 48, which are 1/4 and 1/2 of the whole window. Here, the prediction is performed with the original finest grids by integrating the ODE function. It can be seen that RNN-ODE has smaller prediction errors than RNN on average, and adding addictive steps helps achieve slightly better performance. Furthermore, LSTM still achieves the smallest error most of the time. The reason for this is explained in Section 5.2 and may be due to its more complex network structure.

Figure 7 (right) and Figure A11 in Appendix B.5 present examples of prediction on the testing windows for prediction lengths 48 and 24, respectively. These examples demonstrate that RNN-ODE-Adap models capture the cycles and trends of the ECG more effectively than RNN. The good performance implies that the proposed algorithm could be used to fit and predict the ECG-type signal well. This also implies potential future applications of RNN-ODE-Adap in real healthcare problems.

Figure 6: Left: Comparison of the fitting errors of the underlying intensity function of the simulated event-type data generated from the Hawkes process for RNN, LSTM, RNN-ODE, RNN-ODE-Adap. \(x\)-axis represents computational complexity, \(y\)-axis is the fitting error computed as in Eq. (A7). Right: Examples of fitted intensity function of the simulated event times data generated from the Hawkes process using RNN, LSTM, RNN-ODE, RNN-ODE-Adap.
## 6 Discussion

In this paper we propose a general framework for constructing adaptive time steps when using the neural ODE combining the observed data to model times series. We demonstrate that it tends to be more efficient for modeling "spike-like" time series. The proposed algorithm of adaptive time steps is widely applicable to other types of models, not limited to neural ODE and RNN models. Moreover, the selection of adaptive time steps can be generalized to a broad class of non-stationary time series with different kinds of non-stationarities.

This highlights the potential for further research in this filed, which may be approached from several angles. Firstly, there is a need for more theoretical analysis of the proposed RNN-ODE-Adap framework under a broad spectrum of non-stationary time series, ranging from continuous to discontinuous data sequences. Additionally, there is a need to investigate the more flexible adaptive scheme, which could freely add middle steps adaptively and merge time grids without being restricted to dyadic partitioning.

## Acknowledgement

Y.T. and X.C. are partially supported by Simons Foundation (ID 814643) and NSF (DMS-2007040).
