LONG EXPRESSIVE MEMORY FOR SEQUENCE MODELING
Year: 2022

Abstract:
We propose a novel method called Long Expressive Memory (LEM) for learning long-term sequential dependencies. LEM is gradient-based, it can efficiently process sequential tasks with very long-term dependencies, and it is sufficiently expressive to be able to learn complicated input-output maps. To derive LEM, we consider a system of multiscale ordinary differential equations, as well as a suitable time-discretization of this system. For LEM, we derive rigorous bounds to show the mitigation of the exploding and vanishing gradients problem, a wellknown challenge for gradient-based recurrent sequential learning methods. We also prove that LEM can approximate a large class of dynamical systems to high accuracy. Our empirical results, ranging from image and time-series classification through dynamical systems prediction to keyword spotting and language modeling, demonstrate that LEM outperforms state-of-the-art recurrent neural networks, gated recurrent units, and long short-term memory models.

heading: 1 INTRODUCTION
text: Learning tasks with sequential data as inputs (and possibly outputs) arise in a wide variety of contexts, including computer vision, text and speech recognition, natural language processing, and time series analysis in the sciences and engineering. While recurrent gradient-based models have been successfully used in processing sequential data sets, it is well-known that training these models to process (very) long sequential inputs is extremely challenging on account of the so-called exploding and vanishing gradients problem (Pascanu et al., 2013). This arises as calculating hidden state gradients entails the computation of an iterative product of gradients over a large number of steps. Consequently, this (long) product can easily grow or decay exponentially in the number of recurrent interactions.
Mitigation of the exploding and vanishing gradients problem has received considerable attention in the literature. A classical approach, used in Long Short-Term Memory (LSTM) (Hochreiter & Schmidhuber, 1997) and Gated Recurrent Units (GRUs) (Cho et al., 2014), relies on gating mechanisms and leverages the resulting additive structure to ensure that gradients do not vanish. However, gradients might still explode, and learning very long-term dependencies remains a challenge for these architectures (Li et al., 2018). An alternative approach imposes constraints on the structure of the hidden weight matrices of the underlying recurrent neural networks (RNNs), for instance by requiring these matrices to be unitary or orthogonal (Henaff et al., 2016; Arjovsky et al., 2016; Wisdom et al., 2016; Kerg et al., 2019). However, constraining the structure of these matrices might lead to significantly reduced expressivity, i.e., the ability of the model to learn complicated inputoutput maps. Yet another approach relies on enforcing the hidden weights to lie within pre-specified bounds, leading to control on gradient norms. Examples include Li et al. (2018), based on independent neurons in each layer, and Rusch & Mishra (2021a), based on a network of coupled oscillators. Imposing such restrictions on weights might be difficult to enforce, and weight clipping could reduce expressivity significantly.
This brief survey highlights the challenge of designing recurrent gradient-based methods for sequence modeling which can mitigate the exploding and vanishing gradients problem, while at the same time being sufficiently expressive and possessing the ability to learn complicated input-output maps efficiently. We seek to address this challenge by proposing a novel gradient-based method.
The starting point for our method is the observation that realistic sequential data sets often contain information arranged according to multiple (time, length, etc., depending on the data and task) scales. Indeed, if there were only one or two scales over which information correlated, then a simple model with a parameter chosen to correspond to that scale (or, e.g., scale difference) should be able to model the data well. Thus, it is reasonable to expect that a multiscale model should be considered to process efficiently such multiscale data. To this end, we propose a novel gradientbased architecture, Long Expressive Memory (LEM), that is based on a suitable time-discretization of a set of multiscale ordinary differential equations (ODEs). For this novel gradient-based method (proposed in Section 2):
• we derive bounds on the hidden state gradients to prove that LEM mitigates the exploding and vanishing gradients problem (Section 4);
• we rigorously prove that LEM can approximate a very large class of (multiscale) dynamical systems to arbitrary accuracy (Section 4); and
• we provide an extensive empirical evaluation of LEM on a wide variey of data sets, including image and sequence classification, dynamical systems prediction, keyword spotting, and language modeling, thereby demonstrating that LEM outperforms or is comparable to state-of-the-art RNNs, GRUs and LSTMs in each task (Section 5).
We also discuss a small portion of the large body of related work (Section 3), and we provide a brief discussion of our results in a broader context (Section 6). Much of the technical portion of our work is deferred to Supplementary Materials.

heading: 2 LONG EXPRESSIVE MEMORY
text: We start with the simplest example of a system of two-scale ODEs,
dy dt = τy (σ (Wyz + Vyu + by)− y) , dz dt = τz (σ (Wzy + Vzu + bz)− z) . (1)
Here, t ∈ [0, T ] is the continuous time, 0 < τy ≤ τz ≤ 1 are the two time scales, y(t) ∈ Rdy , z(t) ∈ Rdz are the vectors of slow and fast variables and u = u(t) ∈ Rm is the input signal. For simplicity, we set dy = dz = d. The dynamic interactions between the neurons are modulated by weight matrices Wy,z,Vy,z , bias vectors by,z and a nonlinear tanh activation function σ(u) = tanh(u). Note that refers to the componentwise product of vectors. However, two scales (one fast and one slow), may not suffice in representing a large number of scales that could be present in realistic sequential data sets. Hence, we need to generalize (1) to a multiscale version. One such generalization is provided by the following set of ODEs,
dy dt = σ̂ (W2y + V2u + b2) (σ (Wyz + Vyu + by)− y) ,
dz dt = σ̂ (W1y + V1u + b1) (σ (Wzy + Vzu + bz)− z) .
(2)
In addition to previously defined quantities, we need additional weight matrices W1,2,V1,2, bias vectors by,z and sigmoid activation function σ̂(u) = 0.5(1 + tanh(u/2)). As σ̂ is monotone, we can set W1,2 = V1,2 ≡ 0 and (b1)j = by, (b2)j = bz , for all 1 ≤ j ≤ d, with σ̂(by,z) = τy,z to observe that the two-scale system (1) is a special case of (2). One can readily generalize this construction to obtain many different scales in (2). Thus, we can interpret (τ z(y, t), τ y(y, t)) = (σ̂ (W1y + V1u + b1) , σ̂ (W2y + V2u + b2)) in (2) as input and state dependent gating functions, which endow ODE (2) with multiple time scales. These scales can be learned adaptively (with respect to states) and dynamically (in time). Moreover, it turns out that the multiscale ODE system (2) is of the same general form (see SM§C) as the well-known HodgkinHuxley equations modeling the dynamics of the action potential for voltage-gated ion-channels in biological neurons (Hodgkin & Huxley, 1952).
Next, we propose a time-discretization of the multiscale ODE system (2), providing a circuit to our sequential model architecture. As is common with numerical discretizations of ODEs, doing so properly is important to preserve desirable properties. To this end, we fix ∆t > 0, and we discretize (2) with the following implicit-explicit (IMEX) time-stepping scheme to arrive at LEM, written in compact form as,
∆tn = ∆tσ̂(W1yn−1 + V1un + b1),
∆tn = ∆tσ̂(W2yn−1 + V2un + b2),
zn = (1−∆tn) zn−1 + ∆tn σ(Wzyn−1 + Vzun + bz), yn = (1−∆tn) yn−1 + ∆tn σ(Wyzn + Vyun + by).
(3)
For steps 1 ≤ n ≤ N , the hidden states in LEM (3) are yn, zn ∈ Rd, with input state un ∈ Rm. The weight matrices are W1,2,z,y ∈ Rd×d and V1,2,z,y ∈ Rd×m and the bias vectors are b1,2,z,y ∈ Rd. We also augment LEM (3) with a linear output state ωn ∈ Ro with ωn = Wyyn, and Wy ∈ Ro×d.

heading: 3 RELATED WORK
text: We start by comparing our proposed model, LEM (3), to the widely used LSTM of Hochreiter & Schmidhuber (1997). Observe that ∆tn,∆tn in (3) are similar in form to the input, forget and output gates in an LSTM (see SM§D), and that LEM (3) has exactly the same number of parameters (weights and biases) as an LSTM, for the same number of hidden units. Moreover, as detailed in SM§D, we show that by choosing very specific values of the LSTM gates and the ∆tn,∆tn terms in LEM (3), the two models are equivalent. However, this analysis also reveals key differences between LEM (3) and LSTMs, as they are equivalent only under very stringent assumptions. In general, as the different gates in both LSTM and LEM (3) are learned from data, one can expect them to behave differently. Moreover in contrast to LSTM, LEM stems from a discretized ODE system (2), which endows it with (gradient) stable dynamics.
The use of multiscale neural network architectures in machine learning has a long history. An early example was provided in Hinton & Plaut (1987), who proposed a neural network with each connection having a fast changing weight for temporary memory and a slow changing weight for long-term learning. More recently, one can view convolutional neural networks as multiscale architectures for processing multiple spatial scales in data (Bai et al., 2020).
The use of ODE-based learning architectures has also received considerable attention in recent years with examples such as continuous-time neural ODEs (Chen et al., 2018; Queiruga et al., 2020; 2021) and their recurrent extensions ODE-RNNs (Rubanova et al., 2019), as well as RNNs based on discretizations of ODEs (Chang et al., 2018; Erichson et al., 2021; Chen et al., 2020; Lim et al., 2021; Rusch & Mishra, 2021a;b). In addition to the specific details of our archiecture, we differ from other discretized ODE-based RNNs in the explicit use of multiple (learned) scales in LEM.

heading: 4 RIGOROUS ANALYSIS OF LEM
text: Bounds on hidden states. The structure of LEM (3) allows us to prove (in SM§E.1) that its hidden states satisfy the following pointwise bound.
Proposition 4.1. Denote tn = n∆t and assume that ∆t ≤ 1. Further assume that the initial hidden states are z0 = y0 ≡ 0. Then, the hidden states zn,yn of LEM (3) are bounded pointwise as,
max 1≤i≤d
max{|zin|, |yin|} ≤ min ( 1,∆ √ tn ) , ∀1 ≤ n, with ∆ = 1 + ∆t√
2−∆t . (4)
On the exploding and vanishing gradient problem. For any 1 ≤ n ≤ N , let Xn ∈ R2d, denoted the combined hidden state, given by Xn = [ z1n,y 1 n, . . . . . . , z i n,y i n, . . . . . . , z d n,y d n ] . For simplicity of the exposition, we consider a loss function: En = 12‖yn − yn‖2, with yn being the underlying ground truth. The training of our proposed model (3) entails computing gradients of the above loss function with respect to its underlying weights and biases θ ∈ Θ = [W1,2,y,z,V1,2,y,z,b1,2,y,z],
at every step of the gradient descent procedure. Following Pascanu et al. (2013), one uses chain rule to show,
∂En ∂θ
= ∑
1≤k≤n
∂E (k) n
∂θ ,
∂E (k) n
∂θ = ∂En ∂Xn ∂Xn ∂Xk ∂+Xk ∂θ
(5)
In general, for recurrent models, the partial gradient ∂E (k) n
∂θ , which measures the contribution to the
hidden state gradient at step n arising from step k of the model, can behave as ∂E (k) n
∂θ ∼ γn−k, for some γ > 0 Pascanu et al. (2013). If γ > 1, then the partial gradient grows exponentially in sequence length, for long-term dependencies k << n, leading to the exploding gradient problem. On the other hand, if γ < 1, then partial gradients decays exponentially for k << n, leading to the vanishing gradient problem. Thus, mitigation of the exploding and vanishing gradient problem entails deriving bounds on the gradients. We start with the following upper bound (proved in SM§E.2), Proposition 4.2. Let zn,yn be the hidden states generated by LEM (3). We assume that ∆t << 1 is chosen to be sufficiently small. Then, the gradient of the loss function En with respect to any parameter θ ∈ Θ is bounded as ∣∣∣∣∂En∂θ
∣∣∣∣ ≤ (1 + Ŷ)tn + (1 + Ŷ)Γt2n, Ŷ = ‖yn‖∞, η = max{‖W1‖∞, ‖W2‖∞,‖Wz‖∞, ‖Wy‖∞}, Γ = 2 (1 + η) (1 + 3η)
(6)
If we choose the hyperparameter ∆t = O(n−1) (see SM (17) for the order-notation), then one readily observes from (6) that the gradient ∂θEn is uniformly bounded for any sequence length n and the exploding gradient problem is clearly mitigated for LEM (3). Even if one chooses ∆t = O(n−s), for some 0 ≤ s ≤ 1, we show in SM Remark E.1 that the gradient can only grow polynomially (e.g. as O(n) for s = 1/2), still mitigating the exploding gradient problem. Following Pascanu et al. (2013), one needs a more precise characterization of the partial gradient ∂θE (k) n , for long-term dependencies, i.e., k << n, to show mitigation of the vanishing gradient problem. In SM§E.3, we state and prove proposition E.2, which provides a precise formula for the asymptotics of the partial gradient. Here, we illustrate this formula in a special case as a corollary, Proposition 4.3. Let yn, zn be the hidden states generated by LEM (3) and the ground truth satisfy yn ∼ O(1). Then, for any k << n (long-term dependencies) we have,
∂E (k) n
∂θ = O
( ∆t 3 2 ) . (7)
Here, constants in O(∆t 32 ) depend on only on η (6) and η = ‖W2‖1 and are independent of n, k.
This formula (7) shows that although the partial gradient can be small, i.e., O(∆t 32 ), it is in fact independent of k, ensuring that long-term dependencies contribute to gradients at much later steps and mitigating the vanishing gradient problem.
Universal approximation of general dynamical systems. The above bounds on hidden state gradients show that the proposed model LEM (3) mitigates the exploding and vanishing gradients problem. However, this by itself, does not guarantee that it can learn complicated and realistic input-output maps between sequences. To investigate the expressivity of the proposed LEM, we will show in the following proposition that it can approximate any dynamical system, mapping an input sequence un to an output sequence on, of the (very) general form,
φn = f ( φn−1,un ) , on = o(φn), ∀ 1 ≤ n ≤ N, (8)
with φn ∈ Rdh ,on ∈ Rdo denoting the hidden and output states, respectively. The input signal is un ∈ Rdu and maps f : Rdh × Rdu 7→ Rdh and o : Rdh 7→ Rdo are Lipschitz continuous. For simplicity, we set the initial state φ0 = 0. Proposition 4.4. For all 1 ≤ n ≤ N , let φn,on be given by the dynamical system (8) with input signal un. Under the assumption that there exists a R > 0 such that max{‖φn‖, ‖un‖} < R, for all 1 ≤ n ≤ N , then for any given > 0 there exists a LEM of the form (3), with hidden states yn, zn ∈ Rdy and output state ωn = Wyyn ∈ Rdo , for some dy such that the following holds,
‖on − ωn‖ ≤ , ∀1 ≤ n ≤ N. (9)
From this proposition, proved in SM§E.4, we conclude that, in principle, the proposed LEM (3) can approximate a very large class of dynamical systems.
Universal approximation of multiscale dynamical systems. While expressing a general form of input-output maps between sequences, the dynamical system (8) does not explicitly model dynamics at multiple scales. Instead, here we consider the following two-scale fast-slow dynamical system of the general form,
φn = f(φn−1,ψn−1,un), ψn = τg(φn,ψn−1,un), on = o(ψn). (10)
Here, 0 < τ << 1 and 1 are the slow and fast time scales, respectively. The underlying maps (f ,g) : Rdh×dh×du 7→ Rdh are Lipschitz continuous. In the following proposition, proved in SM§E.5, we show that LEM (3) can approximate (10) to desired accuracy. Proposition 4.5. For any 0 < τ << 1, and for all 1 ≤ n ≤ N , let φn,ψn,on be given by the two-scale dynamical system (10) with input signal un. Under the assumption that there exists a R > 0 such that max{‖φn‖, ‖ψn‖, ‖un‖} < R, for all 1 ≤ n ≤ N , then for any given > 0, there exists a LEM of the form (3), with hidden states yn, zn ∈ Rdy and output state ωn ∈ Rdo with ωn = Wyn such that the following holds,
‖on − ωn‖ ≤ , ∀1 ≤ n ≤ N. (11) Moreover, the weights, biases and size (number of neurons) of the underlying LEM (3) are independent of the time-scale τ .
This argument can be readily generalized to more than two time scales (see SM Proposition E.4). Hence, we show that, in principle, the proposed model LEM (3) can approximate multiscale dynamical systems, with model size being independent of the underlying timescales. These theoretical results for LEM (3) point to the ability of this architecture to learn complicated multiscale inputoutput maps between sequences, while mitigating the exploding and vanishing gradients problem. Although useful prerequisities, these theoretical properties are certainly not sufficient to demonstrate that LEM (3) is efficient in practice. To do this, we perform several benchmark evaluations, and we report the results below.

heading: 5 EMPIRICAL RESULTS
text: We present a variety of experiments ranging from long-term dependency tasks to real-world applications as well as tasks which require high expressivity of the model. Details of the training procedure for each experiment can be found in SM§A. As competing models to LEM, we choose two different types of architectures—LSTMs and GRUs—as they are known to excel at expressive tasks such as language modeling and speech recognition, while not performing well on long-term dependency tasks, possibly due to the exploding and vanishing gradients problem. On the other hand, we choose state-of-the-art RNNs which are tailor-made to learn tasks with long-term dependencies. Our objective is to evaluate the performance of LEM and compare it with competing models. All code to reproduce our results can be found at https://github.com/tk-rusch/LEM.
Very long adding problem. We start with the well-known adding problem (Hochreiter & Schmidhuber, 1997), proposed to test the ability of a model to learn (very) long-term dependencies. The input is a two-dimensional sequence of length N , with the first dimension consisting of random numbers drawn from U([0, 1]) and with two non-zero entries (both set to 1) in the second dimension, chosen at random locations, but one each in both halves of the sequence. The output is the sum of two numbers of the first dimension at positions, corresponding to the two 1 entries in the second dimension. We consider three very challenging cases, namely input sequences with length N = 2000, 5000 and 10000. The results of LEM together with competing models including stateof-the-art RNNs, which are explicitly designed to solve long-term dependencies, are presented in Fig. 1. We observe in this figure that while baseline LSTM is not able to beat the baseline meansquare error of 0.167 (the variance of the baseline output 1) in any of the three cases, a proper weight initialization for LSTM, the so-called chrono-initialization of Tallec & Ollivier (2018) leads to much better performance in all cases. For N = 2000, all other architectures (except baseline LSTM) beat the baseline convincingly. However for N = 5000, only LEM, chrono-LSTM and coRNN are able
to beat the baseline. In the extreme case of N = 10000, only LEM and chrono-LSTM are able to beat the baseline. Nevertheless, LEM outperforms chrono-LSTM by converging faster (in terms of number of training steps) and attaining a lower test MSE than chrono-LSTM in all three cases.

heading: Sequence length N = 5000
text: Sequential image recognition. We consider three experiments based on two widely-used image recognition data sets, i.e., MNIST (LeCun et al., 1998) and CIFAR-10 (Krizhevsky et al., 2009), where the goal is to predict the correct label after reading in the whole sequence. The first two tasks are based on MNIST images, which are flattened along the rows to obtain sequences of length N = 784. In sequential MNIST (sMNIST), the sequences are fed to the model one pixel at a time in streamline order, while in permuted sequential MNIST (psMNIST), a fixed random permutation is applied to the sequences, resulting in much longer dependency than for sMNIST. We also consider the more challenging noisy CIFAR-10 (nCIFAR-10) experiment (Chang et al., 2018), where CIFAR-10 images are fed to the model row-wise and flattened along RGB channels, resulting in 96-dimensional sequences, each of length 32. Moreover, a random noise padding is applied after the first 32 inputs to produce sequences of length N = 1000. Hence, in addition to classifying the underlying image, a model has to store this result for a long time. In Table 1, we present the results for LEM on the three tasks together with other SOTA RNNs, which were explicitly designed to solve long-term dependency tasks, as well as LSTM and GRU baselines. We observe that LEM outperforms all other methods on sMNIST and nCIFAR-10. Additionally on psMNIST, LEM performs as well as coRNN, which has been SOTA among single-layer RNNs on this task.
EigenWorms: Very long sequences for genomics classification. The goal of this task (Bagnall et al., 2018) is to classify worms as belonging to either the wild-type or four different mutants, based on 259 very long sequences (length N = 17984) measuring the motion of a worm. In addition to the nominal length, it was empirically shown in Rusch & Mishra (2021b) that the EigenWorms sequences exhibit actual very long-term dependencies (i.e., longer than 10k).
Following Morrill et al. (2021) and Rusch & Mishra (2021b), we divide the data into a train, validation and test set according to a 70%, 15%, 15% ratio. In Table 2, we present results for LEM together with other models. As the validation and test sets, each consist of only 39 sequences, we report the mean (and standard deviation of) accuracy over 5 random initializations to rule out lucky outliers. We observe from this table that LEM outperforms all other methods, even the 2-layer UnICORNN architecture, which has been SOTA on this task.
Healthcare application: Heart-rate prediction. In this experiment, one predicts the heart rate from a time-series of measured PPG data, which is part of the TSR archive (Tan et al., 2020) and has been collected at the Beth Isreal Deaconess medical center. The data set, consisting of 7949 sequences, each of length N = 4000, is divided into a train, validation and test set according to a 70%,15%,15% ratio, (Morrill et al., 2021; Rusch & Mishra, 2021b). The results, presented in Table 3, show that LEM outperforms the other competing models, including having a test L2 error of 35% less than the SOTA UnICORNN.
Multiscale dynamical system prediction. The FitzHugh-Nagumo system (Fitzhugh, 1955)
v′ = v − v 3
3 − w + Iext, w′ = τ(v + a− bw), (12)
is a prototypical model for a two-scale fast-slow nonlinear dynamical system, with fast variable v and slow variable w and τ << 1 determining the slow-time scale. This relaxation-oscillator is an approximation to the Hodgkin-Huxley model (Hodgkin & Huxley, 1952) of neuronal actionpotentials under an external signal Iext ≥ 0. With τ = 0.02, Iext = 0.5, a = 0.7, b = 0.8 and initial data (v0, w0) = (c, 0), with c randomly drawn from U([−1, 1]), we numerically approximate (12) with the explicit Runge-Kutta method of order 5(4) in the interval [0, 400] and generate 128 training and validation and 1024 test sequences, each of length N = 1000, to complete the data set. The results, presented in Table 4, show that LEM not only outperforms LSTM by a factor of 6 but also all other methods including coRNN, which is tailormade for oscillatory time-series. This reinforces our theory by demonstrating efficient approximation of multiscale dynamical systems with LEM.
Google12 (V2) keyword spotting. The Google Speech Commands data set V2 (Warden, 2018) is a widely used benchmark for keyword spotting, consisting of 35 words, sampled at a rate of 16 kHz from 1 second utterances of 2618 speakers. We focus on the 12-label task (Google12) and follow the pre-defined splitting of the data set into train/validation/test sets and test different sequential models. In order to ensure comparability of different architectures, we do not use performanceenhancing tools such as convolutional filtering or multi-head attention. From Table 5, we observe that both LSTM and GRU, widely used models in this context, perform very well with a test accuracy of around 95%. Nevertheless, LEM is able to outperform both on this task and provides the best performance.
Language modeling: Penn Tree Bank corpus. Language modeling with the widely used small scale Penn Treebank (PTB) corpus (Marcus et al., 1993), preprocessed by Mikolov et al. (2010), has been identified as an excellent task for testing the expressivity of recurrent models (Kerg et al., 2019). To this end, in Table 6, we report the results of different architectures, with a similar number of hidden units, on the PTB char-level task and observe that RNNs, designed explicitly for learning long-term dependencies, perform significantly worse than LSTM and GRU. On the other hand, LEM is able to outperform both LSTM and GRU on this task by some margin (a test bpc of 1.25 in contrast with approximately a bpc of 1.36). In fact, LEM provides the smallest test bpc among all reported single-layer recurrent models on this task, to the best of our knowledge. This superior performance is further illustrated in Table 7, where the test perplexity for different models on the PTB word-level task is presented. We observe that not only does LEM significantly outperform (by around 40%) LSTM, but it also provides again the best performance among all single layer recurrent models, including the recently proposed TARNN (Kag & Saligrama, 2021). Moreover, the singlelayer results for LEM are better than reported results for multi-layer LSTM models, such as in Gal & Ghahramani (2016) (2-layer LSTM, 1500 units each: 75.2 test perplexity) or Bai et al. (2018) (3-layer LSTM, 700 units each: 78.93 test perplexity).

heading: 6 DISCUSSION
text: The design of a gradient-based model for processing sequential data that can learn tasks with longterm dependencies while retaining the ability to learn complicated sequential input-output maps is
very challenging. In this paper, we have proposed Long Expressive Memory (LEM), a novel recurrent architecture, with a suitable time-discretization of a specific multiscale system of ODEs (2) serving as the circuit to the model. By a combination of theoretical arguments and extensive empirical evaluations on a diverse set of learning tasks, we demonstrate that LEM is able to learn long-term dependencies while retaining sufficient expressivity for efficiently solving realistic learning tasks.
It is natural to ask why LEM performs so well. A part of the answer lies in the mitigation of the exploding and vanishing gradients problem. Proofs for gradient bounds (6),(7) reveal a key role played by the hyperparameter ∆t. We observe from SM Table 8 that small values of ∆t might be needed for problems with very long-term dependencies, such as the EigenWorms dataset. On the other hand, no tuning of the hyperparameter ∆t is necessary for several tasks such as language modeling, keyword spotting and dynamical systems prediction and a default value of ∆t = 1 yielded very good performance. The role and choice of the hyperparameter ∆t is investigated extensively in SM§B.1. However, mitigation of exploding and vanishing gradients problem alone does not explain high expressivity of LEM. In this context, we proved that LEMs can approximate a very large class of multiscale dynamical systems. Moreover, we provide experimental evidence in SM§B.2 to observe that LEM not only expresses a range of scales, as it is designed to do, but also these scales contribute proportionately to the resulting multiscale dynamics. Furthermore, empirical results presented in SM§B.2 show that this ability to represent multiple scales correlates with the high accuracy of LEM. We believe that this combination of gradient stable dynamics, specific model structure, and its multiscale resolution can explain the observed performance of LEM.
We conclude with a comparison of LEM and the widely-used gradient-based LSTM model. In addition to having exactly the same number of parameters for the same number of hidden units, our experiments show that LEMs are better than LSTMs on expressive tasks such as keyword spotting and language modeling, while also providing significantly better performance on long-term dependencies. This robustness of the performance of LEM with respect to sequence length paves the way for its application to learning many different sequential data sets where competing models might not perform satisfactorily.

heading: ACKNOWLEDGEMENTS.
text: The research of TKR and SM was performed under a project that has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 770880). NBE and MWM would like to acknowledge IARPA (contract W911NF20C0035), NSF, and ONR for providing partial support of this work. Our conclusions do not necessarily reflect the position or the policy of our sponsors, and no official endorsement should be inferred.
The authors thank Dr. Ivo Danihelka (DeepMind) for pointing out that the hidden states for LEM satisfy the maximum principles (19), (20).

heading: Supplementary Material for:
text: Long Expressive Memory for Sequence Modeling

heading: A TRAINING DETAILS
text: All experiments were run on CPU, namely Intel Xeon Gold 5118 and AMD EPYC 7H12, except for Google12, PTB character-level and PTB word-level, which were run on a GeForce RTX 2080 Ti GPU. All weights and biases of LEM (3) are initialized according to U(−1/ √ d, 1/ √ d), where d is the number of hidden units.
The hyperparameters are selected based on a random search algorithm, where we present the rounded hyperparameters for the best performing LEM model (based on a validation set) on each task in Table 8.
We base the training for the PTB experiments on the following language modelling code: https://github.com/deepmind/lamb, where we fine-tune, based on a random search algorithm, only the learning rate, input-, output- and state-dropout, L2-penalty term and the maximum gradient norm.
We train LEM for 100 epochs on sMNIST, psMNIST and nCIFAR-10, after which we decrease the learning rate by a factor of 10 and proceed training for 20 epochs. Moreover, we train LEM for 50, 60 as well as 400 epochs on EigenWorms, Google12 and FitzHugh-Nagumo. We decrease the learning rate by a factor of 10 after 50 epochs on Google12. On the Healthcare task, we train LEM for 250 epochs, after which we decrease the learning rate by a factor of 10 and proceed training for 250 epochs.

heading: B FURTHER EXPERIMENTAL RESULTS
text: B.1 ON THE CHOICE OF THE HYPERPARAMETER ∆t.
The hyperparameter ∆t in LEM (3) measures the maximum allowed (time) step in the discretization of the multi-scale ODE system (2). In propositions 4.1, 4.2 and 4.3, this hyperparameter ∆t plays a key role in the bounds on the hidden states (4) and their gradients (6). In particular, setting ∆t = O(N−1) will lead to hidden states and gradients, that are bounded uniformly with respect to the underlying sequence length N . However, these upper bounds on the hidden states and gradients account for worst-case scenarios and can be very pessimistic for the problem at hand.
Thus, in practice, we determine ∆t through a hyperparameter tuning procedure as described in section A. To this end, we perform a random search within ∆t < 2 and present the resulting optimal values of ∆t for each of the considered data sets in Table 8. From this table, we observe that for data sets such as PTB, FitzHugh-Nagumo and Google 12 we do not need any tuning of ∆t and a
default value of ∆t = 1 resulted in very good empiricial performance. On the other data sets such as sMNIST, nCIFAR-10 and the healthcare example, where the sequence length (N = O(103)) is larger, we observe that values of ∆t ≈ 0.1 yielded the best performance. The notable exception to this was for the EigenWorms data set, with a very long sequence length of N = 17984 as well as demonstrated very long range dependencies in the data, see Rusch & Mishra (2021b). Here, a value of ∆t = 1.6× 10−3 resulted in the best observed performance. To further investigate the role of the hyperparameter ∆t in the EigenWorms experiment, we perform a sensitivity study where the value of ∆t is varied and the corresponding accuracy of the trained LEM is observed. The results of this sensitivity study are presented in Fig. 2, where we plot the test accuracy (Y-axis) vs. the value of ∆t (X-axis). From this figure, we observe that the accuracy is rather poor for ∆t ≈ 1 but improves monotonically as ∆t is reduced till a value of approximately 10−2, after which it saturates. Thus, in this case, a value of ∆t = O(N− 12 ) (for sequence length N ) suffices to yield the best empirical performance.
Given this observation, we further test whether ∆t = O(N− 12 ) suffices for other problems with long-term dependencies. To this end, we consider the adding problem and vary the input sequence length by an order of magnitude, i.e., from N = 250 to N = 2000. The value of ∆t is now fixed at ∆t = 1√
N and the resulting test loss (Y-axis) vs the number of training steps (X-axis) is plotted in
Fig. 3. We see from this figure that this value of ∆t sufficed to yield very small average test errors for this problem for all considered sequence lengths N . Thus, empirically a value of ∆t in the range 1√ N ≤ ∆t ≤ 1 yields very good performance.
Even if we set ∆t = O( 1√ N ), it can happen for very long sequences N >> 1 that the gradient can be quite small from the gradient asymptotic formula (7). This might lead to saturation in training, resulting in long training times. However, we do not observe such long training times for very long sequence lengths in our experiment. To demonstrate this, we again consider Fig. 3 where the number of training steps (X-axis) is plotted for sequence lengths that vary an order of magnitude. The figure clearly shows that the approximately the same number of training steps are needed to attain a low test error, irrespective of the sequence length. This is further buttressed in Fig. 1, where similar number of training steps where needed for obtaining the same very low test error, even for long sequence lengths, with N up to 10000. Moreover, from section A, we see that the number of epochs for different data sets is independent of the sequence length. For instance, only 50 epochs were necessary for EigenWorms with a sequence length of N = 17984 and ∆t = 1.6 × 10−3 whereas 400 epochs were required for the FitzHugh-Nagumo system with a ∆t = 1.

heading: B.2 MULTISCALE BEHAVOR OF LEM.
text: LEM (3) is designed to represent multiple scales, with terms ∆tn,∆tn being explicitly designed to learn possible multiple scales. In the following, we will investigate if in practice, LEM learns multiple scales and uses them to yield the observed superior empirical performance with respect to competing models.
To this end, we start by recalling the proposition 4.5 where we showed that in principle, LEM can learn the two underlying timescales of a fast-slow dynamical system (see proposition E.4 for a similar result for the universal approximation of a r-time scale (with r ≥ 2) dynamical system with LEM). Does this hold in practice ? To further investigate this issue, we consider the FitzHughNagumo dynamical system (12) which serves as a prototype for a two-scale dynamical system. We consider this system (12) with the two time-scales being τ = 0.02 and 1 and train LEM for this system. In Fig. 4, we plot the empirical histogram that bins the ranges of learned scales ∆tn,∆tn ≤ ∆t = 2 (for all n and d) and counts the number of occurrences of ∆tn,∆tn in each bin. From this figure, we observe that there is a clear concentration of learned scales around the values 1 and τ = 0.02, which exactly correspond to the underlying fast and slow time scales. Thus, for this model problem, LEM is exactly learning what it is designed to do and is able to learn the underlying time scales for this particular problem.
Nevertheless, one might argue that these learnable mutliple scales ∆tn,∆tn are not necessary and a single scale would suffice to provide good empirical performance. We check this possibility on the FitzHugh-Nagumo data set by simply setting ∆tn,∆tn ≡ ∆t1 (with 1 being the vector with all entries set to 1), for all n and tuning the hyperparameter ∆t. The comparative results are presented in Table 9. We see from this table by not allowing for learnable ∆tn,∆tn and simply setting them to a single scale parameter ∆t and tuning this parameter only leads to results that are comparable to the baseline LSTM model. On the other hand, learning ∆tn,∆tn resulted in an error that is a factor of 6 less than the baseline LSTM test error. Thus, we demonstrate the importance of the ability of the proposed LEM model to learn multiple scales in this example.
Hence, the multiscale resolution of LEM seems essential for the fast-slow dynamical system. Does this multiscale resolution also appear for other datasets and can it explain aspects of the observed empirical performance ? To this end, we consider the Google12 Keyword spotting data set and start by pointing out that given the spatial (with respect to input dimension d) and temporal (with respect to sequence length N ) heterogeneities, a priori, it is unclear if the underlying data has a multiscale structure. We plot the empirical histograms of ∆tn,∆tn in Fig. 5 to observe that even for this problem, the terms ∆tn,∆tn are expressed over a range of scales, amounting to 2 − 3 orders of magnitude. Thus, a range of scales are present in the trained LEM even for this example, but do they affect the empirical performance of LEM ? We investigate this question by performing an ablation study and reporting the results in Fig. 6. In this study, we clip the values of ∆tn,∆tn to lie within the range [2−i, 1], for i = 0, 1, . . . , 7 and plot the statistics of the observed test accuracy of LEM. We observe from Fig. 6 that by clipping ∆tn,∆tn to lie near the default (single scale) value of 1 results in very poor empirical performance of an accuracy of ≈ 65%. Then the accuracy jumps to around 90% when an order of magnitude range for ∆tn,∆tn is considered, before monotonically and slowly increasing to yield the best empirical performance for the largest range of values of ∆tn,∆tn, considered in this study. A closer look at the empirical histograms plotted in Fig. 5 reveal that the proportion of occurrences of ∆tn,∆tn decays as a power law, and not exponentially, with respect to the scale amplitude. This, together with results presented in Fig. 6 suggest that not only do a range of scales occur in learned ∆tn,∆tn, the small scales also contribute proportionately to the dynamics and enable the increase in performance shown in Fig. 6.
Finally, in Fig. 7, we plot the empirical histograms of ∆tn and ∆tn for the learned LEM on the sMNIST data set to observe that again a range of scales are observed and the observed occurrences of ∆tn and ∆tn at each scale decays as a power law with respect to scale amplitude. Hence, we have sufficient empirical evidence to claim that the multi-scale resolution of LEM seems essential to its observed performance. However, further investigation is required to elucidate the precise mechanisms through this multiscale resolution enables superior performance, particularly on problems where the multiscale structure of the underlying data may not be explicit.

heading: B.3 ON GRADIENT-STABLE INITIALIZATION.
text: Specialized weight initialization is a popular tool to increase the performance of RNNs on long-term dependencies tasks. One particular approach is the so-called chrono initialization (Tallec & Ollivier, 2018) for LSTMs, where all biases are set to zero except for the bias of the forget gate as well as the input gate (bf and bi in the LSTM (15)), which are sampled from
bf ∼ log(U [1, Tmax − 1]) bi = −bf ,
where Tmax denotes the maximal temporal dependency of the underlying sequential data. We can see in Table 2 that the chrono initialization significantly improves the performance of LSTM on the EigenWorms task. Hence, we are interested in extending the chrono initialization to LEMs. One possible manner for doing this is as follows: Initialize all biases of LEM to zero except for b1 in (3), which is sampled from
b1 ∼ − log(U [1, Tmax∆t− 1]).
We test the chrono initialization for LEM on the EigenWorms dataset, where we train LEM (without tuning ∆t, i.e., setting ∆t = 1), with and without chrono initialization. We provide the results in Table 10, where we show again the results of LSTM with and without chrono initialization as well as the LEM result with tuned ∆t and without chrono initialization from Table 2 for comparison. We
see from Table 10 that when ∆t is fixed to 1, the chrono initialization significantly improves the result of LEM. However, if we tune ∆t, but do not use the chrono initialization, we significantly improve the performance of LEM again. We further remark that tuning ∆t as well as using chrono initialization for LEM does not improve the results obtained with simply tuning ∆t in LEM. Thus, we conclude that chrono initialization can successfully be adapted to LEM. However, tuning ∆t (which controls the gradients) is still advisable in order to obtain the best possible results.

heading: C RELATION BETWEEN LEM AND THE HODGKIN-HUXLEY EQUATIONS
text: We observe that the multiscale ODEs (2), on which LEM is based, are a special case of the following ODE system,
dz dt = Fz (y, t)−Gz(y, t) z, dy dt = Fy (z, t) H (y, t)−Gy(y, t) y. (13)
As remarked in the main text, it turns out the well-known Hodgkin-Huxley equations Hodgkin & Huxley (1952), modeling the the dynamics of the action potential of a biological neuron can also be written down in the abstract form (13), with dy = 1, dz = 3 and the variables y = y modeling the voltage and z = (z1, z2, z3) modeling the concentration of Potassium activation, Sodium activation and Sodium inactivation channels.
The exact form of the different functions in (13) for the Hodgkin-Huxley equations is given by, Fz(y) = (α1(y), α2(y), α3(y)) ,
Gz(y) = (α1(y) + β1(y), α2(y) + β2(y), α3(y) + β3(y)) ,
α1(y) = 0.01(10 + ŷ − y) e 10+ŷ−y 10 − 1 , α2(y) = 0.1(25 + ŷ − y) e 25+ŷ−y 10 − 1 , α3(y) = 0.07e ŷ−y 20 ,
β1(y) = 0.125e ŷ−y 80 , β2(y) = 4e ŷ−y 18 , β3(y) =
1
1 + e1+ ŷ−y 10
,
Fy(z, t) = u(t) + z 4 1 + z 3 2z3, H(y) = c1(ȳ − y) + c2(ȳ − y), Gy(y) = c3,
(14)
with input current u and constants ŷ, ȳ, c1,2,3, whose exact values can be read from (Hodgkin & Huxley, 1952).
Thus, the multiscale ODEs (2) and the Hodgkin-Huxley equations are special case of the same general family (13) of ODEs. Moreover, the gating functions Gy,z(y), that model voltage-gated ion channels in the Hodgkin-Huxley equations, are similar in form to ∆tn,∆tn in (2).
It is also worth highlighting the differences between our proposed model LEM (and the underlying ODE system (2)) and the Hodgkin-Huxley ODEs modeling the dynamics of the neuronal action potential. Given the complicated form of the nonlinearites Fy,z,Gy,z,H in the Hodgkin-Huxley equations (14), we cannot use them in designing any learning model. Instead, building on the abstract form of (13), we propose bespoke non-linearities in the ODE (2) to yield a tractable learning model, such as LEM (3). Moreover, it should be emphasized that the Hodgkin-Huxley equations only model the dynamics of a single neuron (with a scalar voltage and 3 ion channels), whereas the hidden state dimension d of (2) can be arbitrary.

heading: D RELATION BETWEEN LEM AND LSTM
text: The well-known LSTM (Hochreiter & Schmidhuber, 1997) (in its mainly-used version using a forget gate (Gers et al., 2000)) is given by,
fn = σ̂(Wfhn−1 + Vfun + bf )
in = σ̂(Wihn−1 + Viun + bi)
on = σ̂(Wohn−1 + Voun + bo)
cn = fn cn−1 + in σ(Whn−1 + Vun + b) hn = on σ(cn).
(15)
Here, for any 1 ≤ n ≤ N , hn ∈ Rd is the hidden state and cn ∈ Rd is the so-called cell state. The vectors in, fn,on ∈ Rd are the input, forget and output gates, respectively. un ∈ Rm is the input
signal and the weight matrices and bias vectors are given by W,Wf,i,o ∈ Rd×d,V,Vf,i,o ∈ Rm×d and b,bf,i,o ∈ Rd, respectively. It is straightforward to relate LSTM (15) and LEM (3) by first setting the cell state cn = zn, for all 1 ≤ n ≤ N and the hidden state hn = yn. We further need to assume that the input state in = ∆tn and the forget state has to be fn = 1−∆tn. Finally, the output state of the LSTM (15) has to be
on = ∆tn = 1, ∀1 ≤ n ≤ N. Under these assumptions and by setting ∆t = 1, we can readily observe that the LEM (3) and LSTM (15) are equivalent.
A different interpretation of LEM, in relation to LSTM, is as follows; LEM can be thought of a variant of LSTM but with two cell states yn, zn per unit and no output gate. The input gates are ∆tn and ∆tn and the forget gates are coupled to the input gates. Given that the state zn is fed into the update for the state yn, one can think of one of the cell states sitting above the other, leading to a more sophisticated recursive update for LEM (3), when compared to LSTM (15).
Another key difference between LEM and LSTM is the scaling of the learnable gates in LEM by the small hyperparameter ∆t. It is natural to examine whether scaling LSTM with such a small hyperparameter ∆t will improve its performance on sequential tasks with long-term dependencies. To this end, we propose to scale the input and forget gate of an LSTM with a small hyperparameter in two different ways, where we denote the new forget gate as f̂n and the new input gate as în. The first version is (∆t-LSTM v1)
f̂n = ∆tfn, în = ∆tin,
while the second version is (∆t-LSTM v2)
f̂n = (1−∆t)fn, în = ∆tin.
We can see in Table 11 that both ∆t-scaled versions of the LSTM lead to some improvements over the baseline LSTM for very long-sequence Eigenworms dataset, while still performing very poorly when compared to LEM. Moreover, we can see that standard sub-sampling routines, such as a truncation of the BPTT algorithm or random sub-sampling, applied to LSTMs lead to better improvements than ∆t-scaling the forget and input gate.

heading: E SUPPLEMENT TO THE RIGOROUS ANALYSIS OF LEM
text: In this section, we will provide detailed proofs of the propositions in Section 4 of the main article. We start with the following simplifying notation for various terms in LEM (3),
An−1 = W1yn−1 + V1un + b1,
Bn−1 = W2yn−1 + V2un + b2,
Cn−1 = Wzyn−1 + Vzun + bz,
Dn = Wyzn + Vyun + by.
Note that for all 1 ≤ n ≤ N , An,Bn,Cn,Dn ∈ Rd. With the above notation, LEM (3) can be written componentwise, for each component 1 ≤ i ≤ d as,
zin = z i n−1 + ∆tσ̂(A i n−1)σ(C i n−1)−∆tσ̂(Ain−1)zin−1,
yin = y i n−1 + ∆tσ̂(B i n−1)σ(D i n)−∆tσ̂(Bin−1)yin−1.
(16)
Moreover, we will use the following order-notation,
β = O(α), for α, β ∈ R+ if there exists constants C,C such that Cα ≤ β ≤ Cα. M = O(α), for M ∈ Rd1×d2 , α ∈ R+ if there exists constant C such that ‖M‖ ≤ Cα. (17)
Note that the techniques of proof in this following three sub-sections burrow heavily from those introduced in Rusch & Mishra (2021a).

heading: E.1 PROOF OF PROPOSITION 4.1 OF MAIN TEXT.
text: First, we prove Proposition 4.1, which yields the bound (4) for the hidden states of LEM.
Proof. The proof of the bound (4) is split into 2 parts. We start with the first equation in (16) and rewrite it as,
zin = ( 1−∆tσ̂(Ain−1) ) zin−1 + ∆tσ̂(A i n−1)σ(C i n−1).
Noting that the activation functions are such that 0 ≤ σ̂(x) ≤ 1, for all x and −1 ≤ σ(x) ≤ 1, for all x and using the fact that ∆t ≤ 1, we have from the above expression that,
zin ≤ ( 1−∆tσ̂(Ain−1) ) max ( zin−1, 1 ) + ∆tσ̂(Ain−1) max ( zin−1, 1 ) ,
≤ max ( zin−1, 1 ) .
By a symmetric argument, one can readily show that,
zin ≥ min(−1, zin−1).
Combining the above inequalities yields, min(−1, zin−1) ≤ zin ≤ max ( zin−1, 1 ) . (18)
Iterating (18) over n and using zi0 = 0 for all 1 ≤ i ≤ d leads to,
−1 ≤ zin ≤ 1, ∀n, ∀1 ≤ i ≤ d. (19)
An argument, identical to the derivation of (19), but for the hidden state y yields,
−1 ≤ yin ≤ 1, ∀n, ∀1 ≤ i ≤ d. (20)
Thus, we have shown that the hidden states remain in the interval [−1, 1], irrespective of the sequence length.
Next, we will use the following elementary identities in the proof,
b(a− b) = a 2 2 − b 2 2 − 1 2 (a− b)2 , (21)
for any a, b ∈ R, and also,
ab ≤ a 2 2 + b2 2 , ∀ > 0. (22)
We fix 1 ≤ i ≤ d and multiply the first equation in (16) with zin−1 and apply (21) to obtain,
(zin) 2
2 =
(zin−1) 2
2 + ∆tσ̂(Ain−1)σ(C i n−1)z i n−1 −∆tσ̂(Ain−1)(zin−1)2 + (zin − zin−1)2 2
= (zin−1) 2
2 + ∆tσ̂(Ain−1)σ(C i n−1)z i n−1 −∆tσ̂(Ain−1)(zin−1)2
+ ∆t2
2
( σ̂(Ain−1)σ(C i n−1)− σ̂(Ain−1)zin−1 )2 , (from (16))
≤ (z i n−1) 2
2 + ∆tσ̂(Ain−1)|σ(Cin−1)||zin−1| −∆tσ̂(Ain−1)(zin−1)2
+ ∆t2
2 (σ̂(Ain−1)σ(C i n−1))
2 + ∆t2
2 σ̂(Ain−1) 2(zin−1) 2
+ ∆t2σ̂(Ain−1) 2|σ(Cin−1)||zin−1| (as (a− b)2 ≤ a2 + b2 + 2|a||b|)
We fix = 2−∆t1+∆t in the elementary identity (22) to yield,
|σ(Cin−1)||zin−1| ≤ σ(Cin−1) 2 2 + (zin−1) 2 2
Applying this to the inequality for (zin) 2 leads to,
(zin) 2 2 ≤ (z i n−1) 2 2 + ( ∆tσ̂(Ain−1) + ∆t 2σ̂(Ain−1) 2 ) σ(Cin−1)2 2
−∆tσ̂(Ain−1) [ 1−
2 − ∆tσ̂(A
i n−1)
2 − ∆tσ̂(A
i n−1)
2
] (zin−1) 2.
Using the fact that 0 ≤ σ̂(x) ≤ 1 for all x ∈ R, σ2 ≤ 1 and that = 2−∆t1+∆t , we obtain from the last line of the previous equation that,
(zin) 2 ≤ (zin−1)2 +
∆t+ ∆t2
≤ (zin−1)2 +
∆t(1 + ∆t)2
2−∆t , ∀1 ≤ n.
Iterating the above estimate over n = 1, . . . , n̄, for any 1 ≤ n̄ and setting n̄ = n yields,
(zin) 2 ≤ (zi0)2 + n
∆t(1 + ∆t)2
2−∆t ,
⇒ (zin)2 ≤ tn (1 + ∆t)2
2−∆t as z i 0 = 0, tn = n∆t.
Taking a square root in the above inequality yields,
|zin| ≤ ∆ √ tn, ∀n, ∀1 ≤ i ≤ d. (23)
with ∆ defined in the expression (4).
We can repeat the above argument with the hidden state y to obtain,
|yin| ≤ ∆ √ tn, ∀n, ∀1 ≤ i ≤ d. (24)
Combining (23) and (24) with the pointwise bounds (19) and (20) yields the desired bound (4).

heading: E.2 PROOF OF PROPOSITION 4.2 OF MAIN TEXT.
text: Proof. We can apply the chain rule repeatedly (for instance as in Pascanu et al. (2013)) to obtain,
∂En ∂θ
= ∑
1≤k≤n
∂En ∂Xn ∂Xn ∂Xk ∂+Xk ∂θ︸ ︷︷ ︸
∂E (k) n ∂θ
. (25)
Here, the notation ∂ +Xk ∂θ refers to taking the partial derivative of Xk with respect to the parameter θ, while keeping the other arguments constant.
A straightforward application of the product rule yields,
∂Xn ∂Xk
= ∏
k<`≤n
∂X` ∂X`−1 . (26)
For any k < ` ≤ n, a tedious yet straightforward computation yields the following representation formula,
∂X` ∂X`−1 = I2d×2d + ∆tE `,`−1 + ∆t2F`,`−1. (27)
Here E`,`−1 ∈ R2d×2d is a matrix whose entries are given below. For any 1 ≤ i ≤ d, we have, E`,`−12i−1,2j−1 ≡ 0, j 6= i E`,`−12i−1,2i−1 = −σ̂(Ai`−1),
E`,`−12i−1,2j = (W1)i,j σ̂ ′(Ai`−1) ( σ(Ci`−1)− zi`−1 ) + (Wz)i,j σ̂(A i `−1)σ
′(Ci`−1), ∀1 ≤ j ≤ d E`,`−12i,2j−1 = (Wy)i,j σ̂(B i `−1)σ
′(Di`), ∀1 ≤ j ≤ d E`,`−12i,2j = (W2)i,j σ̂ ′(Bi`−1) ( σ(Di`)− yi`−1 ) , j 6= i
E`,`−12i,2i = −σ̂(Bi`−1) + (W2)i,iσ̂′(Bi`−1) ( σ(Di`)− yi`−1 ) .
(28) Similarly, F`,`−1 ∈ R2d×2d is a matrix whose entries are given below. For any 1 ≤ i ≤ d, we have, F`,`−12i−1,j ≡ 0, ∀ 1 ≤ j ≤ 2d, F`,`−12i,2j−1 = −(Wy)i,j σ̂(Aj`−1)σ̂(Bi`−1)σ′(Di`), 1 ≤ j ≤ d,
F`,`−12i,2j = σ̂(B i `−1)σ ′(Di`) d∑ λ=1 (Wy)i,λ (( σ(Cλ`−1)− zλ`−1 ) σ̂′(Aλ`−1)(W1)λ,j + σ̂(A λ `−1)σ ′(Cλ`−1)(Wz)λ,j ) .
(29) Using the fact that,
sup x∈R
max {|σ(x)|, |σ′(x)|, |σ̂(x)|, |σ̂′(x)|} ≤ 1,
the pointwise bounds (4), the notation tn = n∆t for all n, the definition of η (6) and the definition of matrix norms, we obtain that,
‖E`,`−1‖∞ ≤ max { 1 + ‖Wz‖∞ + (1 + min(1,∆ √ t`))‖W1‖∞, 1 + ‖Wy‖∞ + (1 + min(1,∆ √ t`))‖W2‖∞ } ≤ 1 + (2 + min(1,∆ √ t`))η. (30) By similar calculations, we obtain,
‖F`,`−1‖∞ ≤ ‖Wy‖∞ ( 1 + (1 + min(1,∆ √ t`))‖W1‖∞ + ‖Wz‖∞ ) ≤ η(1 + (2 + min(1,∆ √ t`))η).
(31)
Applying (30) and (31) in the representation formula (27) and observing that ∆t ≤ 1 and ` ≤ n, we obtain.∥∥∥∥ ∂X`∂X`−1 ∥∥∥∥ ∞ ≤ 1 + ( 1 + (2 + min(1,∆ √ t`))η ) ∆t+ η ( 1 + (2 + min(1,∆ √ t`))η ) ∆t2,
≤ 1 + Γ 2 ∆t,
With Γ = 2 (1 + η) (1 + 3η) (32)
Using the expression (26) with the above inequality yields,∥∥∥∥∂Xn∂Xk ∥∥∥∥ ∞ ≤ ( 1 + Γ 2 ∆t )n−k . (33)
Next, we choose ∆t << 1 small enough such that the following holds,( 1 + Γ
2 ∆t
)n−k ≤ 1 + Γ(n− k)∆t, (34)
for any 1 ≤ k < n. Hence applying (34) in (33), we obtain,∥∥∥∥∂Xn∂Xk ∥∥∥∥ ∞ ≤ 1 + Γ(n− k)∆t. (35)
For the sake of definiteness, we fix any 1 ≤ α, β ≤ d and set θ = (Wy)α,β in the following. The following bounds for any other choice of θ ∈ Θ can be derived analogously. Given this, it is straightforward to calculate from the structure of LEM (3) that entries of the vector ∂
+Xk ∂(Wy)α,β
are given by, (
∂+Xk ∂(Wy)α,β ) j
≡ 0, ∀ j 6= 2α,( ∂+Xk
∂(Wy)α,β ) 2α = ∆tσ̂(Bαk−1)σ ′(Dαk )z β k .
(36)
Hence, by the pointwise bounds (4), we obtain from (36) that∥∥∥∥ ∂+Xk∂(Wy)α,β ∥∥∥∥ ∞ ≤ ∆tmin(1,∆ √ tk). (37)
Finally, it is straightforward to calculate from the loss function En = 12‖yn − yn‖2 that ∂En ∂Xn = [ 0,y1n − ȳ1, . . . . . . , 0,ydn − ȳd ] . (38)
Therefore, using the pointwise bounds (4) and the notation Ŷ = ‖ȳ‖∞, we obtain∥∥∥∥ ∂En∂Xn ∥∥∥∥ ∞ ≤ Ŷ + min(1,∆ √ tn). (39)
Applying (35), (37) and (39) in the definition (25) yields,∣∣∣∣∣ ∂E(k)n∂(Wy)α,β ∣∣∣∣∣ ≤ ∆tmin(1,∆√tk)(Ŷ + min(1,∆√tn)) (1 + Γ(n− k)∆t) . (40)
Observing that 1 ≤ k ≤ n, we see that n− k ≤ n and tk ≤ tn. Therefore, (41) can be estimated for any 1 ≤ k ≤ n by,∣∣∣∣∣ ∂E(k)n∂(Wy)α,β ∣∣∣∣∣ ≤ ∆tmin(1,∆√tn)(Ŷ + min(1,∆√tn)) (1 + Γtn) , 1 ≤ k ≤ n. (41) Applying the bound (41) in (25) leads to the following bound on the total gradient,∣∣∣∣ ∂En∂(Wy)α,β ∣∣∣∣ ≤ n∑ k=1 ∣∣∣∣∣ ∂E(k)n∂(Wy)α,β ∣∣∣∣∣
≤ tn min(1,∆ √ tn) ( Ŷ + min(1,∆ √ tn) ) (1 + Γtn)
≤ tn(1 + Ŷ)(1 + Γtn) ≤ (1 + Ŷ)tn + (1 + Ŷ)Γt2n
(42)
which is the desired bound (6) for θ = (Wy)α,β .
Moreover, for long-term dependencies i.e., k << n, we can set tk = k∆t < 1, with k independent of sequence length n, in (40) to obtain the following bound on the partial gradient,∣∣∣∣∣ ∂E(k)n∂(Wy)α,β ∣∣∣∣∣ ≤ ∆t 32 ∆√k (1 + Ŷ) (1 + Γtn) , 1 ≤ k << n. (43) Remark E.1. The bound (6) on the total gradient depends on tn = n∆t, with n being the sequence length and ∆t ≤ 1, a hyperparameter which can either be chosen a priori or determined through a hyperparameter tuning procedure. The proof of the bound (42) relies on ∆t being sufficiently small. It would be natural to choose ∆t ∼ n−s, for some s ≥ 0. Substituting this expression in (6) leads to a bound of the form, ∣∣∣∣∂En∂θ
∣∣∣∣ = O (n2(1−s)) (44) If s = 1, then clearly ∣∣∂En ∂θ ∣∣ = O(1) i.e., the total gradient is bounded. Clearly, the exploding gradient problem is mitigated in this case.
On the other hand, if s takes another value, for instance s = 12 which is empirically observed during the hyperparameter training (see Section B.1, then we can readily observe from (44) that∣∣∂En ∂θ ∣∣ = O(n). Thus in this case, the gradient can grow with sequence length n but only linearly and not exponentially. Thus, the exploding gradient problem is also mitigated in this case.

heading: E.3 PROOF OF PROPOSITION 4.3 OF MAIN TEXT.
text: To mitigate the vanishing gradient problem, we need to obtain a more precise characterization of the gradient ∂E (k) n
∂θ defined in (25). For the sake of definiteness, we fix any 1 ≤ α, β ≤ d and set θ = (Wy)α,β in the following. The following formulas for any other choice of θ ∈ Θ can be derived analogously. Moreover, for simplicity of notation, we set the target function X̄n ≡ 0. Proposition 4.3 is a straightforward corollary of the following, Proposition E.2. Let yn, zn be the hidden states generated by LEM (3), then we have the following representation formula for the hidden state gradient,
∂E (k) n
∂θ = ∆tσ̂(Bαk−1)σ ′(Dαk )z β k (y α n − ȳαn)
+ ∆t2σ̂(Bαk−1)σ ′(Dαk )z β k  d∑ j=1 ( yjn − ȳjn ) n∑ `=k+1 σ̂′(Bj`−1) ( σ(Dj`)− y j `−1 ) (W2)j,2α  + ∆t2σ̂(Bαk−1)σ ′(Dαk )z β k [ n∑
`=k+1
σ̂(Bα`−1) (y α n − ȳαn)
] +O(∆t3).
(45)
Here, the constants in O could depend on η defined in (6) (main text).
Proof. The starting point for deriving an asymptotic formula for the hidden state gradient ∂E (k) n
∂θ is to observe from the representation formula (27), the bound (31) on matrices F`,`−1 and the order notation (17) that,
∂X` ∂X`−1 = I2d×2d + ∆tE `,`−1 +O(∆t2), (46)
as long as η is independent of ∆t.
By using induction and the bounds (30),(31), it is straightforward to calculate the following representation formula for the product,
∂Xn ∂Xk
= ∏
k<`≤n
∂X` ∂X`−1 = I2d×2d + ∆t n∑ `=k+1 E`,`−1 +O(∆t2). (47)
Recall that we have set θ = (Wy)α,β . Hence, by the expressions (38) and (36), a direct but tedious calculation leads to,
∂En ∂Xn I2d×2d ∂+Xk ∂θ = ∆tσ̂(Bαk−1)σ ′(Dαk )z β k (y α n − ȳαn) , (48)
n∑ `=k+1 ∂En ∂Xn E`,`−1 ∂+Xk ∂θ = (49)
∆tσ̂(Bαk−1)σ ′(Dαk )z β k  d∑ j=1 ( yjn − ȳjn ) n∑ `=k+1 σ̂′(Bj`−1) ( σ(Dj`)− y j `−1 ) (W2)j,2α − n∑ `=k+1 σ̂(Bα`−1) (y α n − ȳαn)  . (50)
Therefore, by substituting the above expression into the representation formula (47) yields the desired formula (45).
In order to prove the formula (7) (see Proposition 4.3 of main text), we focus our interest on longterm dependencies i.e., k << n. Then, a closer perusal of the expression in (48), together with the pointwise bounds (4) which implies that yk−1 ≈ O( √ ∆t), results in the following,
∂En ∂Xn I2d×2d ∂+Xk ∂θ
= O ( ∆t 3 2 ) . (51)
Similarly, we also obtain,
∆t n∑ `=k+1 ∂En ∂Xn E`,`−1 ∂+Xk ∂θ = O ( ∆t 3 2 ) . (52)
Combining (51) and (52) results in the desired asymptotic bound (7).
Remark E.3. The upper bound on the gradient (6) and the gradient asymptotic formula (7) impact the choice of the timestep hyperparameter ∆t. For sequence length n, if we choose ∆t ∼ n−s, with s ≥ 0, we see from Remark E.1 that the upper bound on the total gradient scales like O(n2(1−s)). On the other hand, from (7), the gradient contribution from long-term dependencies will scale like O(n−3s2 ). Hence, a small value of s ≈ 0, will ensure that the gradient with respect to long-term dependencies will be O(1). However, the total gradient will behave like O(n2) and possibly blow up fast. Similarly, setting s ≈ 1 leads to a bounded gradient, while the contributions from long-term dependencies decay as fast as n −3 2 . Hence, one has to find a value of s that balances both these requirements. Equilibrating them leads to s = 47 , ensuring that the total gradient grows sub-linearly while long-term dependencies still contribute with a sub-linear decay. This value is very close to the empirically observed value of s = 12 which also ensures that the total gradient grows linearly and the contribution of long-term dependencies decays sub-linearly in the sequence length n.

heading: E.4 PROOF OF PROPOSITION 4.4
text: Proof. To prove this proposition, we have to construct hidden states yn, zn, output state ωn, weight matrices W1,2,y,z,Wy,V1,2,y,z and bias vectors b1,2,y,z such that LEM (3) with output state ωn = Wyyn approximates the dynamical system (8).
Let R∗ > R >> 1 and ∗ < , be parameters to be defined later. By the theorem for universal approximation of continuous functions with neural networks with the tanh activation function σ = tanh (Barron, 1993), given ∗, there exist weight matrices W1 ∈ Rd1×dh , V1 ∈ Rd1×du ,W2 ∈ Rdh×d1 and bias vector b1 ∈ Rd1 such that the tanh neural network defined by,
N1(h, u) = W2σ (W1h+ V1u+ b1) , (53)
approximates the underlying function f in the following manner,
max max(‖h‖,‖u‖)<R∗
‖f(h, u)−N1(h, u)‖ ≤ ∗. (54)
Similarly, one can readily approximate the identity function g(h, u) = h with a tanh neural network of the form, N̄2(h) = W̄2σ ( W̄1h ) , (55)
such that max
‖h‖,‖u‖<R∗ ‖g(h)−N2(h)‖ ≤ ∗. (56)
Next, we define the following dynamical system,
z̄n = W2σ (W1ȳn−1 + V1un + b1) , ȳn = W̄2σ ( W̄1z̄n ) ,
(57)
with initial states z̄0 = ȳ0 = 0.
Using the approximation bound (54), we derive the following bound, ‖φn − ȳn‖ = ∥∥f (φn−1,un)− z̄n + z̄n − ȳn∥∥
≤ ‖f ( φn−1,un ) −W2σ (W1ȳn−1 + V1un + b1) ‖+ ‖g(z̄n)− W̄2σ ( W̄1z̄n ) ‖
≤ ‖f ( φn−1,un ) − f (ȳn−1,un) ‖+ ‖f (ȳn−1,un)−W2σ (W1ȳn−1 + V1un + b1) ‖
+ ‖g(z̄n)− W̄2σ ( W̄1z̄n ) ‖
≤ Lip(f)‖φn−1 − ȳn−1‖+ 2 ∗ (from (54), (56)). Here, Lip(f) is the Lipschitz constant of the function f on the compact set {(h, u) ∈ Rdh×du : ‖h‖, ‖u‖ < R∗}. Note that one can readily prove using the fact that ȳ0 = z̄0 = 0, bounds (54), (56) and the assumption ‖φn‖, ‖un‖ < R, that ‖z̄n‖, ‖ȳn‖ < R∗ = 2R. Iterating the above inequality over n leads to the bound,
‖φn − ȳn‖ ≤ 2 ∗ n−1∑ λ=0 Lip(f)λ. (58)
Hence, using the Lipschitz continuity of the output function o in (8), one obtains,
‖on − o(ȳn)‖ ≤ 2 ∗Lip(o) n−1∑ λ=0 Lip(f)λ, (59)
with Lip(o) being the Lipschitz constant of the function o on the compact set {h ∈ Rdh : ‖h‖ < R∗}. Next, we can use the universal approximation theorem for neural networks again to conclude that given a tolerance ̄, there exist weight matrices W3 ∈ Rd2×dh ,W4 ∈ Rdh×d2 and bias vector b2 ∈ Rd2 such that the tanh neural network defined by,
N3(h) = W4σ (W3h+ b2) , (60)
approximates the underlying output function o in the following manner,
max ‖h‖<R∗
‖o(h)−N3(h)‖ ≤ ̄. (61)
Now defining, ω̄n = W4σ (W3ȳn + b2) , (62)
we obtain from (61) and (59) that,
‖on − ω̄n‖ ≤ ̄+ 2 ∗Lip(o) n−1∑ λ=0 Lip(f)λ. (63)
Next, we introduce the notation, z̃n = σ (W1ȳn−1 + V1un + b1) , ỹn = σ ( W̄1z̄n ) . (64)
From (57), we see that z̄n = W2z̃n, ȳn = W̄2ỹn (65)
Thus from (65) and (63), we have
ω̄n = W4σ (W3W2ỹn + b2) , = W4σ ( W3W2σ ( W̄1W2z̃n ) + b2 ) .
(66)
Define the function R : Rdh × Rdu 7→ Rdo by, R(z) = W4σ ( W3W2σ ( W̄1W2z ) + b2 ) . (67)
The function, defined above, is clearly Lipschitz continuous. We can apply the universal approximation theorem for tanh neural networks to find, for any given tolerance ̃, weight matrices W5 ∈ Rd3×d4 ,W6 ∈ Rdo×d3 , V2 ∈ Rd3×du and bias vector b3 ∈ Rd3 such that the following holds,
max max(‖z‖)<R∗
‖R(z)−W6σ(W5z + b3)‖ ≤ ̃. (68)
Denote ω̃n = W6σ(W5z̃n + b3), then from (68) and (66), we obtain that
‖ω̄n − ω̃n‖ ≤ ̃. Combining this estimate with (63) yields,
‖on − ω̃n‖ ≤ ̃+ ̄+ 2 ∗Lip(o) n−1∑ λ=0 Lip(f)λ. (69)
Now, we collect all ingredients to define the LEM that can approximate the dynamical system (8). To this end, we define hidden states zn,yn ∈ R2dh as
zn = [z̃n, ẑn] , yn = [ỹn, ŷn] ,
with z̃n, ẑn, ỹn, ŷn ∈ Rdh . These hidden states are evolved by the dynamical system,
z⊥n = σ
([ W1W̄2 0
0 0
] y⊥n−1 + [V1un, 0] ⊥ + [b1, 0] ⊥ ) ,
y⊥n = σ ([ W̄1W2 0 W5 0 ] z⊥n + [0, 0] ⊥ + [0, b3] ⊥ ) (70)
and the output state is calculated by,
ω⊥n = [0,W6]y ⊥ n . (71)
Finally, we can recast the dynamical system (70), (71) as a LEM of the form (3) for the hidden states yn, zn, defined in (70), with the following parameters, Now, define the hidden states ȳn, z̄n ∈ Rdy for all 1 ≤ n ≤ N by the LEM (3) with the following parameters,
∆t = 1, dy = 2dh,
W1 = W2 = V1 = V2 = 0
b1 = b2 = b∞,
Wz =
[ W1W̄2 0
0 0
] , Vz = [V1, 0], bz = [b1, 0]
Wy = [ W̄1W2 0 W5 0 ] , Vy = 0, bz = [0, b3].
Wy = [0,W6].
(72)
Here, b∞ ∈ Rdh is defined as b∞ = [b∞, b∞, . . . , . . . , b∞],
with 1 << b∞ is such that |1− σ̂(b∞)| ≤ δ. (73)
The nature of the sigmoid function guarantees the existence of such a b∞ for any δ. As δ decays exponentially fast, we set it to 0 in the following for notational simplicity.
It is straightforward to verify that the output state of the LEM (3) with parameters given in (72) is ωn = ω̃n.
Therefore, from (69) and by setting ̄ < 3 , ̃ < 3 and
∗ < 6Lip(o) N−1∑ λ=0 Lip(f)λ ,
we prove the desired bound (9).

heading: E.5 PROOF OF PROPOSITION 4.5
text: Proof. The proof of this proposition is based heavily on the proof of Proposition 4.4. Hence, we will highlight the main points of difference.
As the steps for approximation of a general Lipschitz continuous output map are identical to the corresponding steps in the proof of proposition 4.4 (see the steps from Eqns. (59) to (69)), we will only consider the following linear output map for convenience herein,
o(ψn) = Wcψn. (74)
Let R∗ > R >> 1 and ∗ < , be parameters to be defined later. By the theorem for universal approximation of continuous functions with neural networks with the tanh activation function σ = tanh, given ∗, there exist weight matrices W f1 ,W f 2 ∈ Rd1×dh , V f1 ∈ Rd1×du ,W f3 ∈ Rdh×d1 and bias vector bf1 ∈ Rd1 such that the tanh neural network defined by,
Nf (h, c, u) = W f 3 σ ( W f1 h+W f 2 c+ V f 1 u+ b f 1 ) , (75)
approximates the underlying function f in the following manner,
max max(‖h‖,‖c‖,‖u‖)<R∗
‖f(h, c, u)−Nf (h, c, u)‖ ≤ ∗. (76)
Next, we define the following map,
G(h, c, u) = g(h, c, u) + ( 1− 1
τ
) c, (77)
for any τ > 0.
By the universal approximation theorem, given ∗, there exist weight matrices W g1 ,W g 2 ∈ Rd2×dh , V g1 ∈ Rd2×du ,W g3 ∈ Rdh×d2 and bias vector bg1 ∈ Rd2 such that the tanh neural network defined by,
Ng(h, c, u) = W g 3 σ (W g 1 h+W g 2 c+ V g 1 u+ b g 1) , (78)
approximates the function G (77) in the following manner,
max max(‖h‖,‖c‖,‖u‖)<R∗
‖G(h, c, u)−Nf (h, c, u)‖ ≤ ∗. (79)
Note that the sizes of the neural network Ng can be made independent of the small parameter τ by simply taking the sum of the neural networks approximating the functions g and ĝ(h, c, u) = c with tanh neural networks. As neither of these functions depend on the small parameter τ , the sizes of the corresponding neural networks are independent of the small parameter too.
Next, as in the proof of proposition 4.4, one can readily approximate the identity function f̂(h, c, u) = h with a tanh neural network of the form,
N̄f (h) = W̄2σ ( W̄1h ) , (80)
such that max
‖h‖,‖c‖,‖u‖<R∗ ‖f̂(h, c, u)−Nf (h)‖ ≤ ∗, (81)
and with the same weights and biases, one can approximate the identity function ĝ(h, c, u) = c with the tanh neural network, N̄g(c) = W̄2σ ( W̄1c ) , (82)
such that max
‖h‖,‖c‖,‖u‖<R∗ ‖ĝ(h, c, u)−Ng(c)‖ ≤ ∗. (83)
Next, we define the following dynamical system, ẑn = W f 3 σ ( W f1 ỹn−1 +W f 2 ŷn−1 + V f 1 un + b f 1 ) ,
z̃n = W̄2σ ( W̄1ŷn−1 ) ,
ŷn = (1− τ)ŷn−1 + τW g3 σ (W g1 ẑn +W g2 z̃n + V g1 un + bg1) , ỹn = W̄2σ ( W̄1ẑn ) ,
(84)
with hidden states ẑn, z̃n, ŷn, ỹn ∈ Rdh and with initial states ẑ0 = z̃0 = ŷ0 = ỹ0 = 0. We derive the following bounds,
‖φn − ẑn‖ = ‖f(φn−1,ψn−1,un)−W f3 σ ( W f1 ỹn−1 +W f 2 ŷn−1 + V f 1 un + b f 1 ) ‖
≤ ‖f(φn−1,ψn−1,un)− f(ỹn−1, ẑn−1,un)‖, + ‖f(ỹn−1, ẑn−1,un)−W f3 σ ( W f1 ỹn−1 +W f 2 ŷn−1 + V f 1 un + b f 1 ) ‖
≤ Lip(f) ( ‖φn−1 − ẑn−1‖+ 2‖ỹn−1 − ẑn−1‖+ ‖ψn−1 − ỹn−1‖ ) + ∗ (by (79))
≤ Lip(f) ( ‖φn−1 − ẑn−1‖+ ‖ψn−1 − ỹn−1‖ ) + (1 + 2Lip(f)) ∗ (by (81), (84)),
and ‖ψn − ŷn‖ = ‖(1− τ)(ψn−1 − ŷn−1) + τ ( G(φn,ψn−1,un)−W g3 σ (W g2 z̃n +W g1 ẑn + V g1 un + bg1) ) ‖
≤ ‖ψn−1 − ŷn−1‖+ τ‖G(φn,ψn−1,un)−G(ẑn, z̃n,un)‖ + τ‖G(ẑn, z̃n,un)−W g3 σ (W g2 z̃n +W g1 ẑn + V g1 un + bg1) ‖ ≤ ‖ψn−1 − ŷn−1)‖+ τLip(G) ( ‖φn − ẑn‖+ ‖z̃n − ŷn−1‖+ ‖ψn−1 − ŷn−1‖ ) + τ ∗, ≤ (1 + τLip(G))(1 + Lip(f))‖ψn−1 − ŷn−1‖+ τLip(G)Lip(f)‖φn−1 − ẑn−1‖ + τ(1 + Lip(G)(2 + 2Lip(f))) ∗,
where the last inequality follows by using the previous inequality together with (84) and (83).
As τ < 1, it is easy to see from (77) that Lip(G) < Lip(g) + 2τ . Therefore, the last inequality reduces to,
‖ψn − ŷn‖ ≤ (3 + τLip(g))(1 + Lip(f))‖ψn−1 − ŷn−1‖+ (2 + τLip(g))Lip(f)‖φn−1 − ẑn−1‖ + (τ + (2 + τLip(g))(2 + 2Lip(f))) ∗.
Adding we obtain, ‖φn − ẑn‖+ ‖ψn − ŷn‖ ≤ C∗ ( ‖φn−1 − ẑn−1‖+ ‖ψn−1 − ŷn−1‖ ) +D∗ ∗, (85)
where, C∗ = max{(3 + Lip(g))Lip(f),Lip(f)(3 + Lip(g))(1 + Lip(f))}, D∗ = 1 + (2 + Lip(g))(2 + 2Lip(f)).
(86)
Iterating over n leads to the bound,
‖φn − ẑn‖+ ‖ψn − ŷn‖ ≤ ∗D∗ n−1∑ λ=0 (C∗)λ. (87)
Here, Lip(f),Lip(g) are the Lipschitz constants of the functions f ,g on the compact set {(h, c, u) ∈ Rdh×dh×du : ‖h‖, ‖c‖, ‖u‖ < R∗}. Note that one can readily prove using the zero values of initial states, the bounds (81), (83) and the assumption ‖φn‖, ‖ψn‖, ‖un‖ < R, that ‖ẑn‖, ‖z̃n‖, ‖ŷn‖, ‖ỹn‖ < R∗ = 2R.
Using the definition of the output function (11) and the bound (87) that,
‖on − o(ŷn)‖ ≤ ‖Wc‖ ∗D∗ n−1∑ λ=0 (C∗)λ. (88)
Defining the dynamical system, z∗n = σ ( W f1 W̄2ȳn−1 +W f 2 W g 3 y ∗ n−1 + V f 1 un + b f 1 ) z̄n = σ ( W̄1W g 3 y ∗ n−1 )
y∗n = (1− τ)y∗n−1 + τσ ( W g1W f 3 z ∗ n +W g 2 W̄2z̄n + V g 1 un + b g 1 ) ,
ȳn = σ ( W̄1W f 3 z ∗ n ) .
(89)
By multiplying suitable matrices to (84), we obtain that,
ẑn = W f 3 z ∗ n, z̃n = W̄2z̄n, ŷn = W g 3 y ∗ n, ỹn = W̄2ȳn. (90)
Finally, in addition to b∞ defined in (58), for any given τ ∈ (0, 1], we introduce bτ ∈ R defined by σ̂(bτ ) = τ. (91)
The existence of a unique bτ follows from the fact that the sigmoid function σ̂ is monotone. Next, we define the two vectors b∞,bτ ∈ R2dh as
bi∞ = b∞, ∀ 1 ≤ i ≤ 2dh, biτ = bτ , ∀ 1 ≤ i ≤ dh, biτ = b∞, ∀ dh + 1 ≤ i ≤ 2dh.
(92)
We are now in a position to define the LEM of form (3), which will approximate the two-scale dynamical system (10). To this end, we define the hidden states zn,yn ∈ R2dh such that zn = [z∗n, z̄n] and yn = [y ∗ n, ȳn]. The parameters for the corresponding LEM of form (3) given by,
∆t = 1, dy = 2dh
W1 = W2 = V1 = V2 ≡ 0, b1 = b∞, b2 = bτ ,
Wz =
[ W f2 W g 3 W f 1 W̄2
W̄1W g 3 0
] , Vz = [V f 1 0], bz = [b f 1 , 0],
Wy =
[ W g1W f 3 W g 2 W̄2
W̄1W f 3 0
] , Vz = [V g 1 0], bz = [b g 1, 0],
(93)
and with following parameters defining the output states, Wy = [WcW g 3 0] , (94) yields an output state ωn = Wyyn.
It is straightforward to observe that ωn ≡ o(ŷn). Hence, the desired bound (11) follows from (87) by choosing,
∗ = D∗ N−1∑ λ=0 (C∗)λ .
The proof of proposition 4.5 can be readily extended to prove the following proposition about a general r-scale dynamical system of the form,
φ1n = τ1f 1(φ1n−1,φ 2 n−1, . . . ,φ r n−1un),
φ2n = τ2f 2(φ1n−1,φ 2 n−1, . . . ,φ r n−1un),
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
φrn = τrf r(φ1n−1,φ 2 n−1, . . . ,φ r n−1un),
on = o(φ 1 n,φ 2 n, . . . ,φ r n).
(95)
Here, τ1 ≤ τ2 . . . . . . ≤ τr ≤ 1, with r > 1, are the r-time scales of the dynamical system (95). We assume that the underlying maps f1,2,...,r are Lipschitz continuous. We can prove the following proposition,
Proposition E.4. For all 1 ≤ n ≤ N , let φ1,2,...,rn ,on be given by the r-scale dynamical system (95) with input signal un. Under the assumption that there exists a R > 0 such that max{‖φ1n‖, |φ2n‖, . . . , |φrn‖, ‖un‖} < R, for all 1 ≤ n ≤ N , then for any given > 0, there exists a LEM of the form (3), with hidden states yn, zn and output state ωn with ωn = Wyn such that the following holds, ‖on − ωn‖ ≤ , ∀1 ≤ n ≤ N. (96) Moreover, the weights, biases and size (number of neurons) of the underlying LEM (3) are independent of the time-scales τ1,2,...,r.

heading: F LEMS EMULATE HETEROGENEOUS MULTISCALE METHODS FOR ODES
text: Following Kuehn (2015), we consider the following prototypical example of a fast-slow system of ordinary differential equations,
φ′(t) = 1
τ (f(ψ)− φ) ,
ψ′(t) = g(φ,ψ). (97)
Here φ,ψ ∈ Rm are the fast and slow variables respectively and 0 < τ << 1 is a small parameter. Note that we have rescaled time and are interested in the dynamics of the slow variable ψ(t) in the time interval [0, T ].
A naive time-stepping numerical scheme for (97) requires a time step size δt ∼ O(τ). Thus, the computation will entail time updates N ∼ O(1/τ). Hence, one needs a multiscale ODE solver to approximate the solutions of the system (97). One such popular ODE solver can be derived by using the Heterogenous multiscale method (HMM); see Kuehn (2015) and references therein. This in turns, requires using two time stepping schemes, a macro solver for the slow variable, with a time step ∆t of the form,
ψn = ψn−1 + ∆̃tg(φn,ψn−1). (98)
Here, the time step ∆̃t < 1 is independent of the small parameter τ .
Moreover, the fast variable is updated using a micro solver of the form,
φ (k) n−1 = φ (k−1) n−1 − δt(f(ψn−1)− φ (k−1) n−1 ), 1 ≤ k ≤ K.
φn = φ K n−1,
φ (0) n−1 = φn−1.
(99)
Note that the micro time step size δt and the number of micro time steps K are assumed to independent of the small parameter τ .
It is shown in Kuehn (2015) (Chapter 10.8) that for any given small tolerance > 0, one can choose a macro time step ∆̃t, a micro time step δt, the number K of micro time steps, the number N of macro time steps, independent of τ , such that the discrete states ψn approximate the slow-variable ψ(tn) (with tn = n∆̃t) of the fast-slow system (97) to the desired accuracy of .
Our aim is to show that we can construct a LEM of the form (3) such that the states φn,ψn, defined in (98), (99) can be approximated to arbitrary accuracy. By combining this with the accuracy of HMM, we will prove that LEMs can approximate the solutions of the fast-slow system (97) to desired accuracy, independent of the small parameter τ in (97).
Proposition F.1. Let φn,ψn ∈ Rm, for 1 ≤ n ≤ N , be the states defined by the HMM dynamical system (98), (99). For any given > 0, there exists a LEM of the form (3) with hidden states [zn,yn], where zn,yn ∈ Rdm and output states ωhn, ωcn such that the following holds,
max { ‖φn − ωhn‖, ‖ψn − ωcn‖ } ≤ , ∀1 ≤ n ≤ N. (100)
Proof. We start by using iteration on the micro solver (99) from k = 1 to k = K to derive the following,
φn = δtφn−1 + (1− δt)f(ψn−1), δt = (1− δt)K .
(101)
As δt < 1, we have that δt < 1.
By the universal approximation theorem for tanh neural networks, for any given tolerance ∗, there exist weight matrices W f1 ∈ Rd1×m,W f2 ∈ Rm×d1 and bias vector bf1 ∈ Rd1 such that the tanh neural network defined by,
Nf (c) = W f 2 σ ( W f1 c+ b f 1 ) , (102)
approximates the underlying function f in the following manner,
max ‖c‖<R∗
‖f(c)−Nf (c)‖ ≤ ∗. (103)
Next, we define the following map,
G(h, c) = g(h, c) + c, (104)
By the universal approximation theorem, given ∗, there exist weight matrices W g1 ,W g 2 ∈ Rd2×m,W g3 ∈ Rm×d2 and bias vector bg1 ∈ Rd2 such that the tanh neural network defined by,
Ng(h, c) = W g 3 σ (W g 1 h+W g 2 c+ b g 1) , (105)
approximates the function G (104) in the following manner,
max max(‖h‖,‖c‖)<R∗
‖G(h, c)−Ng(h, c)‖ ≤ ∗. (106)
Next, as in the proof of propositions 4.4 4.5, one can readily approximate the identity function f̂(h, c) = h with a tanh neural network of the form,
N̄f (h) = W̄2σ ( W̄1h ) , (107)
such that max
‖h‖,‖c‖<R∗ ‖f̂(h, c)−Nf (h)‖ ≤ ∗, (108)
and with the same weights and biases, one can approximate the identity function ĝ(h, c) = c with the Tanh neural network, N̄g(c) = W̄2σ ( W̄1c ) , (109)
such that max
‖h‖,‖c‖<R∗ ‖ĝ(h, c)−Ng(c)‖ ≤ ∗. (110)
Then, we define the following dynamical system, ẑn = δtẑn + (1− δt)W f2 σ ( W f1 ŷn−1 + b f 1 ) ,
z̃n = W̄2σ ( W̄1ŷn−1 ) ,
ŷn = (1− ∆̃t)ŷn−1 + ∆̃tW g3 σ (W g1 ẑn +W g2 z̃n + bg1) , ỹn = W̄2σ ( W̄1ẑn ) ,
(111)
with hidden states ẑn, z̃n, ŷn, ỹn ∈ Rm and with initial states ẑ0 = z̃0 = ŷ0 = ỹ0 = 0. Completely analogously as in the derivation of (87), we can derive the following bound,
‖φn − ẑn‖+ ‖ψn − ŷn‖ ≤ C∗ ∗, (112) with constant C∗ = C∗ (n,Lip(f),Lip(g)).
Defining the dynamical system,
z∗n = δtz ∗ n + (1− δt)σ ( W f1 W g 3 ŷn−1 + b f 1 ) z̄n = σ ( W̄1W g 3 y ∗ n−1 )
y∗n = (1− ∆̃t)y∗n−1 + ∆̃tσ ( W g1W f 3 z ∗ n +W g 2 W̄2z̃n + b g 1 ) ȳn = σ ( W̄1W f 2 z ∗ n ) .
(113)
By multiplying suitable matrices to (113), we obtain that,
ẑn = W f 2 z ∗ n, z̃n = W̄2z̄n, ŷn = W g 3 y ∗ n, ỹn = W̄2ȳn. (114)
In addition to b∞ defined in (58), for δt ∈ (0, 1], we introduce bδ ∈ R defined by
σ̂(bδ) = δt. (115)
Similarly for ∆̃t ∈ (0, 1], we introduce b∆ ∈ R defined by
σ̂(b∆) = ∆̃t. (116)
The existence of unique bδ and b∆ follows from the fact that the sigmoid function σ̂ is monotone.
Next, we define the two vectors b∞,bδ,b∆ ∈ R2m as biδ = bδ, ∀ 1 ≤ i ≤ m, biδ = b∞, ∀m+ 1 ≤ i ≤ 2m, bi∆ = b∆, ∀ 1 ≤ i ≤ m, bi∆ = b∞, ∀m+ 1 ≤ i ≤ 2m.
(117)
We define the LEM of form (3), which will approximate the HMM (98),(99). To this end, we define the hidden states zn,yn ∈ R2m such that zn = [z∗n, z̄∗n] and yn = [y∗n, ȳ∗n]. The parameters for the corresponding LEM of form (3) given by,
∆t = 1, dy = 2m
W1 = W2 = V1 = V2 ≡ 0, b1 = bδ, b2 = b∆,
Wz =
[ W f1 W g 3 0
W̄1W g 3 0
] , Vz = 0, bz = [b f 1 , 0],
Wy =
[ W g1W f 3 W g 2 W̄2
W̄1W f 2 0
] , Vz = 0, bz = [b g 1, 0].
(118)
The output states are defined by,
ωhn = W f 2 z ∗ n, ω h n = W g 3 y ∗ n (119)
It is straightforward to observe that ωhn = ẑn, ω c n = ŷn. Hence, the desired bound (100) follows from (112) by choosing, ∗ =
C∗ .

Review summaries:
This paper introduces a model class of neural networks whose architecture is defined by a circuit that computes a discretized step of a pair of multi-scale ODEs. The key novelty is on the multi-scale aspect of this system. In terms of representation capability, this class can represent LSTMs (and vice-versa) and can also approximate any dynamical system from a very broad class of dynamical systems (as well as multi-scale dynamical systems). The authors also prove upper bounds in l_infinity on both the hidden state values and the gradients, and argue this addresses the exploding gradient problem. They also address the vanishing gradient problem by arguing that the magnitude of the partial gradient corresponding to the contribution that depends on step k of the circuit does not depend on k.  They then present a variety of experiments using this new architectural component for modeling several sequence problems. They include a small study on the multi-scale nature of the datasets they study as well.

The authors propose a new RNN architecture, long expressive memory (LEM), motivated by a system of ODEs with multiple time constants. They prove that it can avoid the vanishing gradient problem while retaining the flexibility to approximate a broad class of dynamical systems. They report comparable or improved prediction performance of LEM-based sequence models across a very wide variety of tasks, as compared to several recent alternatives.

This paper introduces Long Expressive Memory (LEM; Eq. 3) which is a new architecture for recurrent networks derived from discretization of a particular system of multi-scale ODEs (Eq. 2). It presents interesting theoretical results that characterize the properties of the proposed architecture (Sec. 4). The key results are that under certain conditions (e.g. fine discretization) the gradient propagation through time is well behaved, and that it is expressive enough to represent general dynamical systems, and multiscale dynamical systems in particular. Empirical results are presented on the adding problem, a synthetic two-scale dynamical system, sequential MNIST/CIFAR-10 classification, EigenWorms classification of very long time series, heart rate prediction, the Speech Commands dataset, and character-level modeling on Penn Treebank (Sec. 5). Across all tasks, the authors report best performance using the proposed architecture when compared to various methods in the literature. 

The paper tries to propose a new recurrent architecture that could address the well-known issues (of recurrent models) like vanishing gradient and exploding gradient. The architecture is a kind of a realization of numerical discretization of ODEs using an implicit-explicit time-stepping scheme. The authors clearly explain why the designed model architecture can capture multiscale data. The authors also connect the proposed method to vanilla LSTM, Hodgkin-Huxley equations, and heterogeneous multiscale methods for ODEs.  The authors also try to provide theoretical evidence that the proposed methods mitigate the gradient exploding and vanishing issues (under some conditions) while learning informative representations for long/short sequence data.  The authors tried on many tasks across different domains, and the proposed methods consistently outperformed the baselines. The benefits are pretty significant in some tasks.