Learning Long-Term Dependencies in Irregularly-Sampled Time Series

Introduction

Irregularly-sampled time series, routine data streams in medical and business settings, can be modeled effectively by a time-continuous version of recurrent neural networks (RNNs). These class of RNNs whose hidden states are identified by ordinary differential equations, termed an ODE-RNN , provably suffer from the vanishing and exploding gradient problem (see Figure 1, the first two models), when trained by reverse-mode automatic differentiation .

An elegant solution to the vanishing gradient phenomenon , which results in difficulties in learning long-term dependencies in RNNs, is the long short term memory networks (LSTM) . LSTMs enforce a constant error propagation through the hidden states, learn to forget, and disentangle the hidden states (memory) from their output states. Despite becoming the standard choice in modeling regularly-sampled temporal dynamics, LSTMs similar to other discretized RNN models, face difficulties when the time-gap between the observations are irregular.

In this paper, we propose a compromise to design a novel recurrent neural network algorithm that simultaneously enjoys the approximation capability of ODE-RNNs in modeling irregularly-sampled time series and capability of learning long-term dependencies of the LSTMs’ computational graph.

To perform this, we let an LSTM cell compute its implicit memory mechanism by their typical (input, forget, and output) gates while receiving their feedback inputs from a time-continuous output state representation. This way, we incorporate a continuous-time dynamical flow within the LSTM network, enabling cells to respond to data arriving at arbitrary time-lags, while avoiding the vanishing gradient problem, a model we call ODE-LSTMs (See Figure 1, the last model).

We compare ODE-LSTMs to standard and advanced continuous-time RNN variants, on a set of synthetic and real-world sparse time-series tasks, and discover consistently better performance.

To put this in context, we first theoretically prove that the class of ODE-RNNs suffers from the exploding and vanishing gradient problem, making them unable to learn long-term dependencies efficiently. We show that learning ODE-RNNs by the adjoint method does not help with this problem. As a solution, we propose ODE-LSTMs, a continuous-time RNN model capable of learning long-term dependencies of irregularly-sampled time-series.

ODE-RNNs Instead of explicitly defining a state update function, ODE-RNNs identify an ordinary differential equations in the following form :

where $x_{t}$ is the input sequence, $h_{t}$ is an RNN’s hidden state, and $\tau$ is a dampening factor. The time-lag $T$ specifies at what times the inputs $x_{t}$ have been sampled.

ODE-RNNs were recently rediscovered and have shown promise in approximating irregularly-sampled data, thanks to the implicit definition of time in their resulting dynamical systems. ODE-RNNs can be trained by backpropagation through time (BPTT) through ODE solvers, or by treating the solver as a black-box and apply the adjoint method to gain memory efficiency . In Section 3, we show this family of recurrent networks faces difficulty to learn long-term dependencies.

Long Short-term Memory LSTMs express their discretized hidden states as a pair $(c_{t},h_{t})$ and its update function, $f_{\theta}(x_{t+1},(c_{t},h_{t}),1)\mapsto(c_{t+1},h_{t+1})$ is defined as follows:

where $\sigma$ is the sigmoid function $x\mapsto 1/(1+\exp(-x))$, the matrices $W_{x}$, $R_{x}$, and vectors $b_{x}$ for $x\in\{z,i,f,o\}$ are the weights of the RNN. The formulation shown in Equations (2-7) extends the original LSTM graph by a biased forget gate (as implemented in PyTorch and TensorFlow ). LSTMs demonstrate great performance on learning equidistant streams of data , however similar to other discrete-state RNNs, they are puzzled with the events arriving in-between observations. In Section 4, we introduce a continuous-time long short-term memory algorithm to tackle this.

In this section, we show that ODE-RNNs trained via backpropagation through time (BPTT) are susceptible to vanishing and exploding gradients. We also illustrate that the adjoint method is not immune to these gradient issues. We first formally define the gradient problems of the RNNs, and progressively construct Theorem 1.

Gradient propagation in recurrent networks Hochreiter discovered that the error-flow in the BPTT algorithm realizes a power series that determines the effectiveness of the learning process . In particular, the state-previous state Jacobian of an RNN:

governs whether the propagated error exponentially grows (explodes), exponentially vanishes, or stays constant. Formally:

Let $h_{t+T}=f(x_{t+T},h_{t},T)$ be a recurrent neural network, then we say unit $i$ of the network $f$ suffers from a vanishing gradient if for some small $\varepsilon>0$ it hold that

where $N$ is the dimension of the hidden state $h_{t}$ and super-script $v^{i}$ denotes the $i$-th entry of the vector $v$. We say unit $i$ of the network $f$ suffers from an exploding gradient if it holds that

We say the whole network $f$ suffers from a vanishing or respectively exploding gradient problem if the above condition hold for some of its units.

The factor $\varepsilon$ in Eq. 9 is essential as Gers et al. observed that a learnable vanishing factor in the form of a forget-gate significantly benefits the learning capabilities of RNNs, i.e., the network can learn to forget. Note that a RNN can simultaneously suffer from a vanishing and an exploding gradient by the definition above.

Now, consider an ODE-RNN given by Eq. 1 is implemented either by an Explicit Euler discretization or by a Runge-Kutta method . We can formulate their state-previous state Jacobian in the following two lemmas:

Let $\dot{h}=f_{\theta}(x,h,T)-h\tau$ be an ODE-RNN. Then state-previous state Jacobian of the explicit Euler is given by the following equation:

Let $\dot{h}=f_{\theta}(x,h,T)-h\tau$ be an ODE-RNN. Then state-previous state Jacobian of the Runge-Kutta method is given by

, where $\sum_{j=1}^{M}b_{i}=1$ and some $K_{i}$.

The proofs for Lemma 1 and Lemma 2 is provided in the supplements. Consequently, we have:

Let $\dot{h}=f_{\theta}(x,h,T)-h\tau$, and $h_{t}$ the RNN obtained by simulating the ODE by a solver based on the explicit Euler or Runge-Kutta method. Then the RNN suffers from a vanishing and exploding gradient problem, except for parameter configurations which give the non-trainable constant dynamics $f_{\theta}(h,x)=0$, and cases where $f_{\theta}(h,x)$ is constant, for a particular input sequence $x$ and $\theta$.

The proof is given in full in the supplementary materials. A brief outline of the proof is as follow: First, we look a the special cases of $\frac{\partial f}{\partial h}-\tau=0$. While such $f$ would enforce a constant error propagation by making the Jacobians equal to the identity, it also removes all dynamics from the ODE state. In other words, it would operate the ODE as a memory element. Intuitively, any interesting function $f_{\theta}$ pushes the Jacobians away from the identity matrix, creating a vanishing or exploding gradient depending on $f_{\theta}$.

(ODE-RNNs suffer from a vanish or exploding gradient regardless of the choice of ODE-solver) Let $\dot{h}=f_{\theta}(x,h,T)-h\tau$, with $f_{\theta}$ being uniformly Lipschitz continuous. Moreover, let $h(t)$ be the solution of the initial value problem with initial state $h_{0}$. Then, the gradients $\frac{\partial h(T)}{\partial h_{0}}$, i.e, the Jacobian of the ODE state at time $T$ with respect to the initial state $h_{0}$, can vanish and explode, except for parameter configurations which give rise to the non-trainable constant dynamics $f_{\theta}(h,x)=0$, and cases where $f_{\theta}(h,x)$ constant, for a particular input sequence $x$ and parameters $\theta$.

The proof is given in full in the supplementary materials. A brief outline of the proof is as follow: We start by approximating the initial-value problem by an explicit Euler method with a uniform step-size. We then let the step-size approach zero which due to the Picard–Lindelöf theorem, makes the series converge to the true solution of the ODE. Based on bounds on $\frac{\partial f}{\partial h}$, we can obtain bounds of the gradients in the limit, which can vanish or explode depending on $f_{\theta}$.

Does the adjoint method solve the vanishing gradient problem? Adjoint sensitivity method allows for performing memory-efficient reverse-mode automatic differentiation for training neural networks with their hidden states defined by ODEs . The method, however, possesses lossy reverse-mode integration steps, as it forgets the computed steps during the forward-pass. Consequently, at each reverse-mode step, the backward gradient pass diverges from the true forward pass . This is because the auxiliary differential equation in the adjoint sensitivity method, $a(t)$, still contains state-dependent components at each reverse-step, which depends on the historical values of the hidden states’ gradient. Therefore, both vanilla BPTT and the adjoint method face difficulties for learning long-term dependencies. In the next section, we propose a solution.

The RNN state of a standard LSTM network, is represented by a pair $(c_{t},h_{t})$, where $c_{t}$ is the memory cell and $h_{t}$ the output state, i.e., see Equations (2- 7). The memory $c_{t}$ ensures a constant error propagation and the output state $h_{t}$ enables the LSTM to learn non-linear dynamics. We modify the way the output state $h_{t}$ is computed while preserving its gating mechanisms and memory cell.

To perform this, we declare the output dynamics of a cell by a continuous-time representation, which realizes an ODE-RNN. This way, the output state depend on the elapsed time when processing irregularly sampled time-series. Nonetheless, as the LSTM gates receive feedback connections from the cells’ outputs, the gating dynamics become dependent on the time-lag as well. The resulting architecture termed an ODE-LSTM is shown in algorithm 1.

Experimental evaluation

We constructed quantitative settings with synthetic and real-world benchmarks. We assessed the generalization performance of time-continuous RNN architectures on datasets that are deliberately created to express long-term dependencies and are of irregularly-sampled nature. All code and data is available at https://github.com/mlech26l/ode-lstms.

Baselines. We compare ODE-LSTM to a large variety of continuous-time RNNs introduced to model irregularly-sampled data. This set includes RNNs with continuous-state dynamics such as ODE-RNN and CT-RNNs , state-decay mechanisms such as CT-GRU , RNN Decay , CT-LSTM , and GRU-D , in addition to oscillatory models such as Phased-LSTM .

Furthermore, we tested ODE-LSTMs against intuitive time-gap modeling approaches we built here, termed an augmented LSTM topology as well as bi-directional RNNs . Experimental settings are given in the supplements.

We formulated a modified time-series variant of the XOR problem . In particular, the model observes a block of binary data in the form of a bit-after-bit time-series. The objective is then to learn an XOR function of the incoming bit-stream. This setup is equivalent to the binary-classification of the input sequence, where the labels are obtained by applying an XOR function to the inputs.

While any non-linear recurrent neural network architecture can learn the correct function, training the network to do so is non-trivial. For the model to make an accurate prediction, all bits in an upcoming chunk are required to be taken into account. However, the error signal is only provided after the last bit is observed. Consequently, during learning, the prediction error needs to be propagated to the first input time-step to precisely capture the dependencies, (see Figure 2).

We designed two modes, a dense encoding mode in which the input sequence is represented as a regular, periodically sampled time-series, and an event-based mode which compresses the data into irregularly sampled bit-streams, e.g., $1,1,1,1$ is encoded as $(1,t=4)$. (See Table 2). We observed that a considerable number of RNN variants faced difficulties in modeling these tasks, even in the dense-encoding model.

In particular, ODE-RNNs, CT-RNNs, RNN-Decay, Phased-LSTM, and GRU-ODE could not solve the XOR problem in the first mode. Phased-LSTM and RNN-Decay improved their performance in the second modality, whereas ODE-RNNs, CT-RNNs, and GRU-ODE still could not solve the task. The core reason for their mediocre performance is the exploitation of the vanishing gradient problem during training. The rest of the RNN variants (except CT-GRU) were successful in solving the task in both modes, with ODE-LSTM outperforming others in an event-based encoding scenario.

2 Person activity recognition with irregularly sampled time-series

We consider the person activity recognition dataset from the UCI repository . This task’s objective is to classify the current activity of a person, from four inertial measurement sensors worn on the person’s arms and feet. Even though the four sensors are measured at a fixed period of 211ms, the random phase-shifts between them creates an irregularly sampled time-series. Rubanova et al. showed that ODE-based RNN architectures perform remarkably well on this dataset. Here, we benchmarked the performance of the ODE-LSTM model against other variants.

This setting realizes a per-time-step classification problem. That is a new error signal is presented to the network at every time-step which makes the vanishing gradient less of an issue here. The results in Table 3 shows that the ODE-LSTM outperforms other RNN models on this dataset. While the significance of an evaluation on a single dataset is limited, it demonstrates that the supreme generalization ability of ODE-LSTM architecture.

3 Event-based sequential MNIST

We determined a challenging sequence classification task by designing an event-based version for the sequential-MNIST dataset. For doing this we followed the procedure described below:

Sequentialization + encoding long-term dependencies transform the 28-by-28 image into a time-series of length 784

Compression + non-uniform sampling encode binary time-series in a event-based format, to get rid of consecutive occurrences of the same binary value, e.g., $1,1,1,1$ is transformed to $(1,t=4)$. (Read more about this experiment in supplements)

Using this sequentialization mechansim, we compress the sequences from 784 to padded sequences of 256 irregularly-sampled datapoints. To perform well on this task, RNNs must learn to store some information up to 256 time-steps, while taking the time-lags between them into account. Since an error signal is issued at the end of the sequence, only an RNN model immune to vanishing gradients can achieve high-degrees of accuracy.

Table 4 demonstrates that ODE-based RNN architectures, such as the ODE-RNN, CT-RNN, and the GRU-ODE struggle to learn a high-fidelity model of this dataset. On the other hand, RNNs built based on a memory mechanism, such as the Bi-directional RNN and GRU-D perform reasonably well, while the performance of ODE-LSTM surpasses that of other models.

4 Walker2d kinematic simulation

In this experiment, we evaluated how well ODE-LSTM can model a physical dynamical system. To this end, we collected simulation data of the Walker2d-v2 OpenAI gym environment using a pre-trained policy (see Figure 3). The objective of the model was to learn the kinematic simulation of the MuJoCo physics engine in an auto-regressive fashion and a supervised learning modality. We increased the complexity of this task by using the pre-trained policy at different training stages (between 500 to 1200 Proximal Policy Optimization (PPO) iterations ) and overwrote 1% of all actions by random actions. Moreover, we simulated frame-skips by removing 10% of the time-steps. Consequently, the dataset is irregularly-sampled. The results, shown in Table 5, indicate that ODE-LSTM can capture the kinematic dynamics of the physics engine better than other algorithms with a high margin.

Discussions, Scope and Limitations

What if we feed in samples’ time-lag as an additional input feature to network? The Augmented LSTM architecture we benchmarked against realizes this concept, which is a simplistic approach to making LSTMs compatible with irregularly sampled data. The RNN could then learn to make sense of the time input, for instance, by making its change proportional to the elapsed-time.

Nonetheless, the time characteristic of an augmented RNN depends purely on its learning process. Consequently, we can only hope that the augmented RNN generalize to unseen time-lags. Our experiments showed that an augmented LSTM performs reasonably well while being outperformed by models that explicitly declare their state by a continuous-time modality, such as ODE-LSTMs.

Difference between bidirectional RNNs and ODE-LSTM? A bi-directional architecture consists of two different types of RNNs reciprocally linked together in an auto-regressive fashion . In our context, the first RNN could be designed to handle irregularly-sample time series while the second one is capable of learning long-term dependencies . For example, an LSTM bidirectionally coupled with an ODE-RNN could, in principle, overcome both challenges. However, the use of heterogeneous RNN architectures might limit the learning process. In particular, due to different learning speeds, the LSTM could already be overfitting long before the ODE-RNN has learned useful dynamics.

Contrarily, our ODE-LSTM interlinks LSTMs and ODE-RNNs not in an autoregressive fashion, but at an architectural level, avoiding the problem of learning at different speeds. Our experiments showed that ODE-LSTMs consistently outperform a bi-directional LSTM-ODE-RNN architecture.

Related Works

Time-continuous RNNs The notion of CT-RNNs was introduced around three decades ago. It is identical to the ODE-RNN architecture with an additional dampening factor $\tau$. In our experiments, however, we observed a competitive performance to our ODE-LSTMs achieved by the GRU-D architecture . GRU-D encodes the dependence on the time-lags by a trainable decaying mechanism, similar to RNN-decay . While this mechanism enables modeling irregularly sampled time-series, it also introduces a vanishing gradient factor to the backpropagation path.

Similarly, CT-GRU adds multiple decay factors in the form of extra dimensions to the RNN state. An attention mechanism inside the CT-GRU then selects which entry along the decay dimension to use for computing the next state update. The CT-GRU aims to avoid vanishing gradients by including a decay rate of 0, i.e., no decay at all. This mechanism nevertheless, fails as illustrated in Table 2.

Phased-LSTM adds a learnable oscillator to LSTM. The oscillator modulates LSTM to create dependencies on the elapsed-time, but also introduces a vanishing factor in its gradients.

GRU-ODE modifies the GRU topology by incorporating a continuous dynamical system. First, GRU is expressed as a discrete difference equation and then transformed into a continuous ODE. This process makes the error-propagation time-dependent, i.e., the near-constant error propagation property of GRU is abolished.

CT-LSTM combines the LSTM architecture with continuous-time neural Hawkes processes. At each time-step, the RNN computes two alternative next state options of its hidden state. The actual hidden state is then computed by interpolating between these two hidden states depending on the elapsed time.

Learning Irregularly-Sampled Data Statistical and functional analysis tools have long been studying non-uniformly-spaced data. An alternative and a natural fit for this problem is the use of time-continuous recurrent networks . We showed that although ODE-RNNs are performant models in these domains, their performance tremendously drops when the incoming samples have long-range dependencies. We solved this shortcoming by introducing ODE-LSTMs.

Learning Long-term Dependencies The notorious question of vanishing/exploding gradient was identified as the core reason for RNNs’ lack of generalizability when trained by gradient descent . Recent studies used state-regularization and long memory stochastic processes to analyze long-range dependencies. Apart from the original LSTM model and its variants that solve the problem in the context of RNNs, very few alternative researches exist .

As the class of CT RNNs become steadily popularized , it is important to characterize them better and understand their applicability and limitations . In this paper, we proposed a method to enable ODE-based RNNs to learn long-term dependencies.

Conclusion

We proposed a solution to learn long-term dependencies in irregularly-sampled input data streams. To perform this, we designed a novel long short term memory network, that possesses a continuous-time output state, and consequently modifies its internal dynamical flow to a continuous-time model. ODE-LSTMs resolve the vanishing and exploding of the gradient problem of the class of ODE-RNNs while demonstrating an attractive performance in learning long-term dependencies on data arriving at non-uniform intervals.

Broader Impact

Who will benefit from this research? Time series data with missing values and non-uniform intervals are the routine settings in many safety-critical application domains, such as medical, business, social, and the automation of industries.

The results of this paper enable users to construct learning systems that not only help handle irregularly sampled data efficiently but also to learn long-term dependencies that might be vital to their application.

For instance, consider the decision-critical domain of surgical processes or the treatment of patients in intensive care units (ICU) in which the medical team has to have access to the process actively, and the steps are taken throughout a surgical procedure, to make/take a current decision/action. An intelligent agent in use as an assistant during surgery must be able to do the same and carefully assign credits to the actions taken in the past (long-term dependencies) to output an accurate decision. This example simultaneously consists of irregularly-sampled inputs and long-term dependencies. Our proposed method enables these modalities.

Preventing failure of the system Like any other intelligent system, our proposed algorithm has to go through robustness analysis (perturbations, noise, and adversarial attack), before being deployed in high-stakes decision-making applications. This process would dramatically reduce the chance of failure of intelligent systems such as ours.

Whether the method leverages biases in the data The mechanisms of "learning to forget" and "learning long-term dependencies" are encoded in our proposed method. Both processes can be used as the controller of biases in data, and help us design fair machine learning systems.

Acknowledgments and Disclosure of Funding

M.L. is supported in parts by the Austrian Science Fund (FWF) under grant Z211-N23 (Wittgenstein Award). R.H. is partially supported by the Horizon-2020 ECSEL Project grant No. 783163 (iDev40), and Boeing.

References

S1 Proofs

Derivation of the Euler’s method Jacobian Let $\dot{h}=f_{\theta}(x,h,T)-h\tau$ be an ODE-RNN. Then the explicit Euler’s method with step-size $T$ is defined as the discretization

Therefore, state-previous state Jacobian is given by

Derivation of the Runge-Kutta Jacobian Let $\dot{h}=f_{\theta}(x,h,T)-h\tau$ be an ODE-RNN. Then the Runge-Kutta method with step-size $T$ is defined as the discretization

where the coefficients $b_{i}$ and the values $K_{i}$ are taken according to the Butcher tableau with $\sum_{j=1}^{M}b_{i}=1$ and $K_{1}=h_{t}$.

Then state-previous state Jacobian of the Runge-Kutta method is given by the following equation :

Note that the explicit Euler method is an instance of the Runge-Kutta method with $M=1$ and $b_{1}=1$.

Proof of ODE-RNN suffering from vanishing or exploding gradients Let $\dot{h}=f_{\theta}(x,h,T)-h\tau$ be an ODE-RNN with latent dimension $N$. Without loss of generality let $h_{0}$ be the initial state at $t=0$ and $h_{T}$ denote the ODE state which should be computed by a numerical ODE-solver. Then ODE-solvers, including fixed-step methods and variable-step methods such as the Dormand-Prince method , discretize the interval $[0,T]$ by a series $t_{0},t_{1},\dots t_{n}$, where $t_{0}=0$ and $t_{n}=T$ and each $h_{t_{i}}$ is computed by a single-step explicit Euler or Runge-Kutta method from $h_{t_{i-1}}$.

Our proof closely aligns with the analysis in Hochreiter and Schmidhuber . We refer the reader to for a rigorous discussion on the vanishing and exploding gradients.

We first prove the theorem for a scalar RNN, i.e., $n=1$, and then extend the discussion to the general case. The error-flow per RNN step between $t=0$ and $t=T$ is given by

which realizes a power series depending on the value

Obviously, the condition that this term is equal to 1 is not enforced during training and violated for any non-trivial $f_{\theta}$, such as $f_{\theta}(h,x)=\sigma(W_{h}h+W_{x}x+\hat{b})$ with $\sigma$ being a sigmoidal or rectified-linear activation function. The exact magnitude depends on the weights $W_{h}$, as

A non-zero time-constant $\tau$ pushes the gradient toward a vanishing region.

Note that the Equation (S6) only becomes equal to 1, if $\sum_{j=1}^{M}b_{i}\frac{\partial f}{\partial h}\Big|_{h=K_{m_{i}}}=\tau$. This would imply that $\frac{\partial h_{t_{m}}}{h_{t_{m}-1}}=0$, i.e., when the change in ODE-state between two time-points is zero. A variable that does not change over time is a memory element. Thus the only solution of enforcing a constant-error propagation is to include an explicit memory element in the architecture which does not change its value between two arbitrary time-points $t_{m}$ and $t_{m-1}$.

For the general case $n\geq 1$, the error-flow per RNN step between $t=0$ and $t=T$ is given by

As $h$ is a vector, we need to consider all possible error-propagation paths. The error-flow from unit $u$ to unit $v$ is then given by summing all $N^{n-1}$ possible paths between $u$ to $v$,

The arguments of the scalar case hold for every individual path in Equation (S9). The only difference between the the scalar case and the individual paths in the vectored version is the non-diagonal connections in the general case do not include the constant 1 and $\tau$. The error-propagation magnitude between $u$ and $v$ with $u\neq v$ is given by

Again, for $f_{\theta}(h,x)=\sigma(W_{h}h+W_{x}x+\hat{b})$ we obtain an error-flow that depends on the weights $W_{h}$ and can be either vanishing or exploding, depending on its magnitude.

Proof that even gradients of the ODE solution can vanish or explode Let $\dot{h}=f_{\theta}(x,h,T)-h\tau$ be an ODE-RNN with latent dimension $N$, with $f_{\theta}$ being uniformly Lipschitz continuous. Without loss of generality let $h_{0}$ be the initial state at $t=0$ and $h_{T}$ denote the ODE state which should be computed by a numerical ODE-solver. We approximate the interval $[0,T]$ by a uniform discretization grid, i.e. $t_{i}-t_{i-1}=t_{j}-t_{j-1}=T/n$ for all $i,j$ $t_{0},t_{1},\dots t_{n}$, where $t_{0}=0$ and $t_{n}=T$ and each $h_{t_{i}}$ is computed by a single-step explicit Euler from $h_{t_{i-1}}$.

Even when making the discretization grid $t_{0},t_{1},\dots t_{n}$ finer and finer, the gradient propagation issue is not resolved. Let $h_{i}$ denote the intermediate values computed by the Picard-iteration, i.e., the explicit Euler. By the Picard–Lindelöf theorem, we know that $h_{T}$ converges to the true solution $h(T)$.

First, we assume there exists a $\xi>0$ such that $\xi\leq\frac{\partial f}{\partial h}\Big|_{h=h_{m}}-\tau\text{ for all }m$. Note that this situation can naturally occur if we have a $f_{\theta}(h,x)=\sigma(W_{h}h+W_{x}x+\hat{b})$. In the limit $n\rightarrow\infty$ we get

Conversely, lets assume there exists a $\xi<0$ such that $\xi\geq\frac{\partial f}{\partial h}\Big|_{h=h_{m}}-\tau\text{ for all }m$. Note that this situation can also naturally occur, for instance if $\tau>0$ and regions where $f^{\prime}$ is small. In the limit $n\rightarrow\infty$ we get

Similar to the argument in the proof above, we can extend the scalar case to the general case. However, summing over all possible path might not be trivial, as the number of possible path also growths to infinity.

Instead, we assume $u=v=l_{1}=\dots l_{n}-1$, i.e., we only look at the error-propagation through the diagonal element $u$.

which is equivalent to the scalar case. For an interesting $f$ such as $f_{\theta}(h,x)=\sigma(W_{h}h+W_{x}x+\hat{b})$, the term $\frac{f}{h}$ depends on the value $W^{u,u}_{h}$. By assuming $W^{w,z}$ for any $(w,z)\neq(u,u)$ is neglectable small, we can infer that the effects of the gradient by any other path in Equation (S11) is neglectable small. Thus the global error flow depends on $W^{u,u}_{h}$, which can make the error-flow either explode or vanish depending on its value.

Note that this argument is similar to arguing that as the multi-dimensional case properly contains the scalar case, the multi-dimensional case can express an exploding or vanishing gradient too.

Proof that the ODE-LSTM does not suffer from a vanishing or exploding gradient

Recall that we assume that $R_{z},R_{i},R_{f},W_{f}$ and $b_{f}$ are initialized close to 0 and that we are at the beginning of the training process, i.e., we assume the weights do not differ significantly from their initialized values.

For the derivative of the input update activation we can simply apply the chain-rule and get

where $\tanh^{\prime}$ denotes the functional derivative of the hyperbolic tangent. As $0\leq\tanh^{\prime}\leq 1$, $0\leq o_{t}^{u}\leq 1$ and most importantly $R_{z}$ is initialized close to 0, we can safely assume that

Similar argument holds for the input and forget gate derivatives, where we assumed that $R_{i}$ and $R_{f}$ are initialized close to 0. Therefore

where $\sigma^{\prime}$ denotes the functional derivatives of the sigmoid function.

Consequently, with a proper weight initialization, the Jacobian simplifies to

We assumed that $W_{f}$ and $b_{f}$ are initialized close to 0. Hence,

, which is less than 1 (no exploding) but much greater than 0 (no vanishing) and ensures a near-constant error propagation at the beginning of the training process.

As already mentioned in the paper, the exact value of the error flow can be controlled by changing the forget gate bias from its default value of 1. If the underlying data distribution contains dependencies with a very long time-lag, we can bring the error flow factor closer to 1 by increasing forget gate bias. Thus enabling the ODE-LSTM to learn even very long-term dependencies in the data.

S2 Experimental evaluation

For models containing differential equations, we used the ODE-solvers as listed in Table S1. Hyperparameter settings used for our evaluation is shown in Table S2.

Batching Sequences of our event-based bit-stream classification task and event-based seqMNIST can have different lengths. To allow an arbitrary batching of several sequences, we pad all sequences to equal length and apply a binary mask during training and evaluation.

The individual datasets are created as follows:

Bit-stream XOR dataset Every data point is a block of 32 random bits. The binary labels are created by applying an XOR function on the bit block, i.e., class A if the number of 1s in the bit-stream are even, class B if the number of 1s in the bit-stream is odd. For training, a cross-entropy loss on these two classes is used. The training set consists of 100,000 samples, which are less than $0.0024\%$ of all possible bit-streams that can occur. The test set consists of 10,000 samples.

For the event-based encoding, we introduce a time-dimension. The time is normalized such that the complete sequence equals 1 unit of time, i.e., 32 bits corresponds to exactly 1 second. An illustration of the two different encodings is shown in Figure S1.

Person Activity We consider a variation of the "Human activity" dataset described in form the UCI machine learning repository . The dataset is comprised of 25 recordings of human participants performing different physical activities. The eleven possible activities are ”lying down”, ”lying”, ”sitting down”, ”sitting”, ”standing up from lying”, ”standing up from”, ”sitting”, ”standing up from sitting on the ground”, ”walking”, ”falling”, ”on all fours”, and ”sitting on the ground”. The objective of this task is to recognize the activity from inertial sensors worn by the participant, i.e., a per-time-step classification problem. We group the eleven activities listed above into seven different classes, as proposed by .

The input data consists of sensor readings from four inertial measurement units placed on the participant’s arms and feet. The sensors are read at a fixed period of 211 ms but have different phase-shifts in the 25 recordings. Therefore, we treat the data as irregularly sampled time-series.

The 25 recordings are split into partially overlapping sequences of length 32, to allow an efficient training of the machine learning models.

Our results are not directly comparable to the experiments in , as we use a different representation of the input features. While represents each input feature as a value-mask pair, i.e., 24 input features, we represent the data in the form of a 7-dimensional feature vector. The first four entries of the input indicate the senor ID, i.e., which arm or foot, whereas the remaining three entries contain the sensor reading.

Event-based seqMNIST The MNIST dataset consists of 70,000 data points split into 60,000 training and 10,000 test samples . Each sample is a 28-by-28 grayscale image, quantized with 8-bits and represents one out of 10 possible digits, i.e., a number from 0 to 10.

We pre-process each sample as follows: We first apply a threshold to transform the 8-bits pixel values into binary values. The threshold is 128, on a scale where 0 represents the lowest possible and 255 the larges possible pixel value. We further transform the 28-by-28 image into a time-series of length 784. Next, we encode binary time-series in a event-based format. Essentially, the encoding step gets rid of consecutive occurrences of the same binary value, i.e., $1,1,1,1$ is transformed into $(1,t=4)$. By introducing a time dimension, we can compress the sequences from 784 to an average of 53 time-steps.

To allow an efficient batching and training, we pad each sequence to a length of 256. Note that no information was lost during this process. We normalize the added time dimension such that 256 symbols correspond to 1 second or unit of time. The resulting task is a per-sequence classification problem of irregularly sampled time-series.

Walker2d kinematic modeling Here we create a dataset based on the Walker2d-v2 OpenAI gym environment and the MuJoCo physics engine . Our objective is to benchmark how well the RNN architecture can model kinematic dynamical systems in an irregularly sampled fashion. The learning setup is based on an auto-regressive supervised learning, i.e., the model predicts the next state of the Walker2d environment based on the current state.

In order to obtain interesting simulation rollouts, we trained a non-recurrent policy by Proximal Policy Optimization (PPO) using the Rllib reinforcement learning framework. We then collect the training data for our benchmark by performing rollouts on the Walker2d-v2 environment using our pre-trained policy. Note that because the policy is deterministic, there is no need to include the actions produced by the policy in the training data.

We introduce three sources of uncertainty to make this task more challenging. First of all, for each rollout we uniformly sample a checkpoint of policy at 562, 822, 923, or 1104 PPO iterations. Secondly, we overwrite 1% of all actions by random actions. Thirdly, we exclude 10% of the time-steps, i.e., we simulate frame-skips/frame-drops. Note that the last step transforms the rollouts into irregularly sampled time-series and introduces a time dimension.

In total, we collected 400 rollouts, i.e., 300 used for training, 40 for validation, and 60 for testing. For an efficient training, we align the rollouts into sequences of length 64. We use the mean-square-error as training loss and evaluation metric. We train each RNN for 200 epochs and log the validation error after each training epochs. At the end, we restore the weights that achieved the best (lowest) validation error and evaluate them on the test set.

Reproducibility statement We publish all code and data used in our experimental setup at this link https://github.com/mlech26l/ode-lstms.