Year: 2021

Abstract:
Viewing recurrent neural networks (RNNs) as continuous-time dynamical systems, we propose a recurrent unit that describes the hidden state’s evolution with two parts: a well-understood linear component plus a Lipschitz nonlinearity. This particular functional form facilitates stability analysis of the long-term behavior of the recurrent unit using tools from nonlinear systems theory. In turn, this enables architectural design decisions before experimentation. Sufficient conditions for global stability of the recurrent unit are obtained, motivating a novel scheme for constructing hidden-to-hidden matrices. Our experiments demonstrate that the Lipschitz RNN can outperform existing recurrent units on a range of benchmark tasks, including computer vision, language modeling and speech prediction tasks. Finally, through Hessian-based analysis we demonstrate that our Lipschitz recurrent unit is more robust with respect to input and parameter perturbations as compared to other continuous-time RNNs.

heading: 1 INTRODUCTION
text: Many interesting problems exhibit temporal structures that can be modeled with recurrent neural networks (RNNs), including problems in robotics, system identification, natural language processing, and machine learning control. In contrast to feed-forward neural networks, RNNs consist of one or more recurrent units that are designed to have dynamical (recurrent) properties, thereby enabling them to acquire some form of internal memory. This equips RNNs with the ability to discover and exploit spatiotemporal patterns, such as symmetries and periodic structures (Hinton, 1986). However, RNNs are known to have stability issues and are notoriously difficult to train, most notably due to the vanishing and exploding gradients problem (Bengio et al., 1994; Pascanu et al., 2013).
Several recurrent models deal with the vanishing and exploding gradients issue by restricting the hidden-to-hidden weight matrix to be an element of the orthogonal group (Arjovsky et al., 2016; Wisdom et al., 2016; Mhammedi et al., 2017; Vorontsov et al., 2017; Lezcano-Casado & MartinezRubio, 2019). While such an approach is advantageous in maintaining long-range memory, it limits the expressivity of the model. To address this issue, recent work suggested to construct hidden-tohidden weights which have unit norm eigenvalues and can be nonnormal (Kerg et al., 2019). Another approach for resolving the exploding/vanishing gradient problem has recently been proposed by Kag et al. (2020), who formulate the recurrent units as a differential equation and update the hidden states based on the difference between predicted and previous states.
In this work, we address these challenges by viewing RNNs as dynamical systems whose temporal evolution is governed by an abstract system of differential equations with an external input. The data are formulated in continuous-time where the external input is defined by the function x = x(t) ∈ Rp, and the target signal is defined as y = y(t) ∈ Rd. Based on insights from dynamical systems theory, we propose a continuous-time Lipschitz recurrent neural network with the functional form{
ḣ = AβA,γAh+ tanh(WβW ,γW h+ Ux+ b) ,
y = Dh ,
(1a) (1b)
where the hidden-to-hidden matrices Aβ,γ ∈ RN×N and Wβ,γ ∈ RN×N are of the form{ AβA,γA = (1− βA)(MA +MTA ) + βA(MA −MTA )− γAI WβW ,γW = (1− βW )(MW +MTW ) + βW (MW −MTW )− γW I,
(2a)
(2b)
where βA, βW ∈ [0, 1], γA, γW > 0 are tunable parameters and MA,MW ∈ RN×N are trainable matrices. Here, h = h(t) ∈ RN is a function of time t that represents an internal (hidden) state, and ḣ = ∂h(t)∂t is its time derivative. The hidden state represents the memory that the system has of its past. The function in Eq. (1) is parameterized by the hidden-to-hidden weight matrices A ∈ RN×N and W ∈ RN×N , the input-to-hidden encoder matrix U ∈ RN×p, and an offset b. The function in Eq. (1b) is parameterized by the hidden-to-output decoder matrix D ∈ Rd×N . Nonlinearity is introduced via the 1-Lipschitz tanh activation function. While RNNs that are governed by differential equations with an additive structure have been studied before (Zhang et al., 2014), the specific formulation that we propose in (1) and our theoretical analysis are distinct.
Treating RNNs as dynamical systems enables studying the long-term behavior of the hidden state with tools from stability analysis. From this point of view, an unstable unit presents an exploding gradient problem, while a stable unit has well-behaved gradients over time (Miller & Hardt, 2019). However, a stable recurrent unit can suffer from vanishing gradients, leading to catastrophic forgetting (Hochreiter & Schmidhuber, 1997b). Thus, we opt for a stable model whose dynamics do not (or only slowly do) decay over time. Importantly, stability is also a statement about the robustness of neural units with respect to input perturbations, i.e., stable models are less sensitive to small perturbations compared to unstable models. Recently, Chang et al. (2019) explored the stability of linearized RNNs and provided a local stability guarantee based on the Jacobian. In contrast, the particular structure of our unit (1) allows us to obtain guarantees of global exponential stability using control theoretical arguments. In turn, the sufficient conditions for global stability motivate a novel symmetric-skew decomposition based scheme for constructing hidden-to-hidden matrices. This scheme alleviates exploding and vanishing gradients, while remaining highly expressive.
In summary, the main contributions of this work are as follows:
• First, in Section 3, using control theoretical arguments in a direct Lyapunov approach, we provide sufficient conditions for global exponential stability of the Lipschitz RNN unit (Theorem 1). Global stability is advantageous over local stability results since it guarantees non-exploding gradients regardless of the state. In the special case where A is symmetric, we find that these conditions agree with those in classical theoretical analyses (Lemma 1).
• Next, in Section 4, drawing from our stability analysis, we propose a novel scheme based on the symmetric-skew decomposition for constructing hidden-to-hidden matrices. This scheme mitigates the vanishing and exploding gradients problem, while obtaining highly expressive hidden-to-hidden matrices.
• In Section 6, we show that our Lipschitz RNN has the ability to outperform state-of-theart recurrent units on computer vision, language modeling and speech prediction tasks. Further, our results show that the higher-order explicit midpoint time integrator improves the predictive accuracy as compared to using the simpler one-step forward Euler scheme.
• Finally, in Section 7), we study our Lipschitz RNN via the lens of the Hessian and show that it is robust with respect to parameter perturbations; we also show that our model is more robust with respect to input perturbations, compared to other continuous-time RNNs.

heading: 2 RELATED WORK
text: The problem of vanishing and exploding gradients (and stability) have a storied history in the study of RNNs. Below, we summarize two particular approaches to the problem (constructing unitary/orthogonal RNNs and the dynamical systems viewpoint) that have gained significant attention.
Unitary and orthogonal RNNs. Unitary recurrent units have received attention recently, largely due to Arjovsky et al. (2016) showing that unitary hidden-to-hidden matrices alleviate the vanishing and exploding gradients problem. Several other unitary and orthogonal models have also been proposed (Wisdom et al., 2016; Mhammedi et al., 2017; Jing et al., 2017; Vorontsov et al., 2017; Jose et al., 2018). While these approaches stabilize the training process of RNNs considerably, they also
limit their expressivity and their prediction accuracy. Further, unitary RNNs are expensive to train, as they typically involve the computation of a matrix inverse at each step of training. Recent work by Lezcano-Casado & Martinez-Rubio (2019) overcame some of these limitations. By leveraging concepts from Riemannian geometry and Lie group theory, their recurrent unit exhibits improved expressivity and predictive accuracy on a range of benchmark tasks while also being efficient to train. Another competitive recurrent design was recently proposed by Kerg et al. (2019). Their approach is based on the Schur decomposition, and it enables the construction of general nonnormal hidden-to-hidden matrices with unit-norm eigenvalues.
Dynamical systems inspired RNNs. The continuous time view of RNNs has a long history in the neurodynamics community as it provides higher flexibility and increased interpretability (Pineda, 1988; Pearlmutter, 1995; Zhang et al., 2014). In particular, RNNs that are governed by differential equations with an additive structure have been extensively studied from a theoretical point of view (Funahashi & Nakamura, 1993; Kim et al., 1996; Chow & Li, 2000; Hu & Wang, 2002; Li et al., 2005; Trischler & D’Eleuterio, 2016). See Zhang et al. (2014) for a comprehensive survey of continuous-time RNNs and their stability properties.
Recently, several works have adopted the dynamical systems perspective to alleviate the challenges of training RNNs which are related to the vanishing and exploding gradients problem. For nonsequential data, Ciccone et al. (2018) proposed a negative-definite parameterization for enforcing stability in the RNN during training. Chang et al. (2019) introduced an antisymmetric hidden-tohidden weight matrix and provided guarantees for local stability. Kag et al. (2020) have proposed a differential equation based formulation for resolving the exploding/vanishing gradients problem by updating the hidden states based on the difference between predicted and previous states. Niu et al. (2019) employed numerical methods for differential equations to study the stability of RNNs.
Another line of recent work has focused on continuous-time models that deal with irregular sampled time-series, missing values and multidimensional time series. Rubanova et al. (2019) and De Brouwer et al. (2019) formulated novel recurrent models based on the theory of differential equations and their discrete integration. Lechner & Hasani (2020) extended these ordinary differential equation (ODE) based models and addresses the issue of vanishing and exploding gradients by designing an ODE-model that is based on the idea of long short-term memory (LSTM). This ODE-LSTM outperforms the continuous-time LSTM (Mei & Eisner, 2017) as well as the GRU-D model (Che et al., 2018) that is based on a gated recurrent unit (GRU).
The link between dynamical systems and models for forecasting sequential data also provides the opportunity to incorporate physical knowledge into the learning process which improves the generalization performance, robustness, and ability to learn with limited data (Chen et al., 2019).

heading: 3 STABILITY ANALYSIS OF LIPSCHITZ RECURRENT UNITS
text: One of the key contributions in this work is that we prove that model (1) is globally exponentially stable under some mild conditions on A and W . Namely, for any initial hidden state we can guarantee that our Lipschitz unit converges to an equilibrium if it exists, and therefore, gradients can never explode. We improve upon recent work on stability in recurrent models, which provide only a local analysis, see e.g., (Chang et al., 2019). In fact, global exponential stability is among the strongest notions of stability in nonlinear systems theory, implying all other forms of Lyapunov stability about the equilibrium h∗ (Khalil, 2002, Definitions 4.4 and 4.5). Definition 1. A point h∗ is an equilibrium point of ḣ = f(h, t) if f(h∗, t) = 0 for all t. Such a point is globally exponentially stable if there exists some C > 0 and λ > 0 such that for any choice of initial values h(0) ∈ RN ,
‖h(t)− h∗‖ ≤ Ce−λt‖h(0)− h∗‖, for any t ≥ 0. (3)
The presence of a Lipschitz nonlinearity in (1) plays an important role in our analysis. While we focus on tanh in our experiments, our proof is more general and is applicable to models whose nonlinearity σ(·) is an M -Lipschitz function. Specifically, we consider the general model
ḣ = Ah+ σ(Wh+ Ux+ b) , (4) for which we have the following stability result. In the following, we let σmin and σmax denote the smallest and largest singular values of the hidden-to-hidden matrices, respectively.
Theorem 1. Let h∗ be an equilibrium point of a differential equation of the form (4) for some x ∈ Rp. The point h∗ is globally exponentially stable if the eigenvalues of Asym := 12 (A + A
T ) are strictly negative, W is non-singular, and either (a) σmin(Asym) > Mσmax(W ); or (b) σ is monotone non-decreasing, W +WT is negative definite, and ATW +WTA is positive definite.
The two cases show that global exponential stability is guaranteed if either (a) the matrix A has eigenvalues with real parts sufficiently negative to counteract expanding trajectories in the nonlinearity; or (b) the nonlinearity is monotone, both A and W yield stable linear systems u̇ = Au, v̇ = Wv, and A,W have sufficiently similar eigenvectors. In practice, case (b) occasionally holds, but is challenging to ensure without assuming specific structure onA,W . Because such assumptions could limit the expressiveness of the model, the next section will develop a tunable formulation for A and W with the capacity to ensure that case (a) holds.
In Appendix A.1, we provide a proof of Theorem 1 using a direct Lyapunov approach. One advantage of this approach is that the driving input x is permitted to evolve in time arbitrarily in the analysis. The proof relies on the classical Kalman-Yakubovich-Popov lemma and circle criterion from control theory — to our knowledge, these tools have not been applied in the modern RNN literature, and we hope our proof can illustrate their value to the community.
In the special case whereA is symmetric and x(t) constant, we show that we can also inherit criteria for both local and global stability from a class of well-studied Cohen–Grossberg–Hopfield models. Lemma 1. Suppose that A is symmetric and W is nonsingular. There exists a diagonal matrix D ∈ RN×N , and nonsingular matrices L, V ∈ RN×N such that an equilibrium of (4) is (globally exponentially) stable if and only if there is a corresponding (globally exponentially) stable equilibrium for the system
ż = Dz + Lσ(V z + Ux+ b). (5)
For a thorough review of analyses of (5), see (Zhang et al., 2014). In this special case, the criteria in Theorem 1 coincide with those obtained for the corresponding model (5). However, in practice, we will not choose A to be symmetric.

heading: 4 SYMMETRIC-SKEW HIDDEN-TO-HIDDEN MATRICES
text: In this section we propose a novel scheme for constructing hidden-to-hidden matrices. Specifically, based on the successful application of skew-symmetric hidden-to-hidden weights in several recent recurrent architectures, and our stability criteria in Theorem 1, we propose an effective symmetricskew decomposition for hidden matrices. Our decomposition allows for a simple control of the matrix spectrum while retaining its wide expressive range, enabling us to satisfy the spectral constraints derived in the previous section on both A and W . The proposed scheme also accounts for the issue of vanishing gradients by reducing the magnitude of large negative eigenvalues.
Recently, several methods used skew-symmetric matrices, i.e., S + ST = 0 to parameterize the recurrent weights W ∈ RN×N , see e.g., (Wisdom et al., 2016; Chang et al., 2019). From a stability analysis viewpoint, there are two main advantages for using skew-symmetric weights: these matrices generate the orthogonal group whose elements are isometric maps and thus preserve norms (Lezcano-Casado & Martinez-Rubio, 2019); and the spectrum of skew-symmetric matrices is purely imaginary which simplifies stability analysis (Chang et al., 2019). The main shortcoming of this parametrization is its reduced expressivity, as these matrices have fewer than half of the parameters of a full matrix (Kerg et al., 2019). The latter limiting aspect can be explained from a dynamical systems perspective: skew-symmetric matrices can only describe oscillatory behavior, whereas a matrix whose eigenvalues have nonzero real parts can also encode viable growth and decay information.
To address the expressivity issue, we aim for hidden matrices which on the one hand, allow to control the expansion and shrinkage of their associated trajectories, and on the other hand, will be sampled from a superset of the skew-symmetric matrices. Our analysis in Theorem 1 guarantees that Lipschitz recurrent units maintain non-expanding trajectories under mild conditions on A and W . Unfortunately, this proposition does not provide any information with respect to the shrinkage of paths. Here, we opt for a system whose expansion and shrinkage can be easily controlled. Formally, the latter requirement is equivalent to designing hidden weights S with smallRλi(S), i = 1, 2, . . . , N , where R(z) denotes the real part of z. A system of the form (4) whose matrices A and W exhibit
small spectra and satisfy the conditions of Theorem 1, will exhibit dynamics with moderate decay and growth behavior and alleviate the problem of exploding and vanishing gradients. To this end, we propose the following symmetric-skew decomposition for constructing hidden matrices:
Sβ,γ := (1− β) · (M +MT ) + β · (M −MT )− γI, (6) where M is a weight matrix, and β ∈ [0.5, 1], γ > 0 are tuning parameters. In the case (β, γ) = (1, 0), we recover a skew-symmetric matrix, i.e., S1,0 + ST1,0 = 0. The construction Sβ,γ is useful as we can easily bound its spectrum via the parameters β and γ, as we show in the next proposition. Proposition 1. Let Sβ,γ satisfy (6), and let M sym = 12 (M +M
T ). The real parts <λi(Sβ,γ) of the eigenvalues of Sβ,γ , as well as the eigenvalues of S sym β,γ = Sβ,γ + S T β,γ , lie in the interval
[(1− β)λmin(M sym)− γ, (1− β)λmax(M sym)− γ].
A proof is provided in Appendix A.2. We infer that β controls the width of the spectrum, while increasing γ shifts the spectrum to the left along the real axis, thus enforcing eigenvalues with nonpositive real parts. Choosing our hidden-to-hidden matrices to be AβA,γA and WβW ,γW of the form (6) for different values of βA, βW and γA, γW , we can ensure small spectra and satisfy the conditions of Theorem 1 as desired. Note, that different tuning parameters β and γ affect the stability behavior of the Lipschitz recurrent unit. This is illustrated in Figure 1, where different values for β and γ are used to construct both Aβ,γ and Wβ,γ and applied to learning simple pendulum dynamics.
One cannot guarantee that model parameters will remain in the stability region during training. However, we can show that when β is taken to be close to one, the eigenvalues of Asymβ,γ and W sym β,γ (which dictate the stability of the RNN) change slowly during training. Let ∆δF denote the change in a function F depending on the parameters of the RNN (1) after one step of gradient descent with step size δ with respect to some loss L(y). For a matrix A, we let λk(A) denote the k-th singular value of A. We have the following lemma. Lemma 2. As β → 1−, maxk |∆δλk(Asymβ,γ )|+ maxk |∆δλk(W sym β,γ )| = O(δ(1− β)2).
Therefore, provided both the initial and optimal parameters lie within the stability region, the model parameters will remain in the stability region for longer periods of time with high probability as β → 1. Further empirical evidence of parameters often remaining in the stability region during training are provided alongside the proof of Lemma 2 in the Appendix (see Figure 5).

heading: 5 TRAINING CONTINUOUS-TIME LIPSCHITZ RECURRENT UNITS
text: ODEs such as Eq. (1) can be approximately solved by employing numerical integrators. In scientific computing, numerical integration is a well studied field that provides well understood techniques (LeVeque, 2007). Recent literature has also introduced new approaches which are designed with neural network frameworks in mind (Chen et al., 2018).
To learn the weightsA,W,U and b, we discretize the continuous model using one step of a numerical integrator between sequence entries. In what follows, a subscript t denotes discrete time indices,
∆t represents the time difference between a pair of consecutive data points. Letting f(h, t) = Ah+ tanh(Wh+ Ux(s) + b) so that ḣ(t) = f(h, t), the exact and approximate solutions for ht+1 given ht are given by
ht+1 = ht + ∫ t+∆t t f(h(s), s)ds := ht + ∫ t+∆t t Ah(s) + tanh(Wh(s) + Ux(s) + b) ds (7)
≈ ht + ∆t · scheme [f, ht, ∆t] , (8)
where scheme represents one step of a numerical integration scheme whose application yields an approximate solution for 1∆t ∫ t+∆t t f(h(s), s)ds given ht using one or more evaluations of f .
We consider both the explicit (forward) Euler scheme,
ht+1 = ht + ∆t ·Aht + ∆t · tanh(zt), (9)
as well as the midpoint method which is a two-stage explicit Runge-Kutta scheme (RK2),
ht+1 = ht + ∆t ·Ah̃+ ∆t · tanh(Wh̃+ Uxt + b), (10)
where h̃ = ht + ∆t/2 · Aht + ∆t/2 · tanh(zt) is an intermediate hidden state. The RK2 scheme can potentially improve the performance since the scheme is more accurate, however, this scheme also requires twice as many function evaluations as compared to the forward Euler scheme. Given a β and γ that yields a globally exponentially stable continuous model, ∆t can always be chosen so that the model remains in the stability region of forward Euler and RK2 (LeVeque, 2007).

heading: 6 EMPIRICAL EVALUATION
text: In this section, we evaluate the performance of the Lipschitz RNN and compare it to other state-ofthe-art methods. The model is applied to ordered and permuted pixel-by-pixel MNIST classification, as well as to audio data using the TIMIT dataset. We show the sensitivity with respect to to random initialization in Appendix B. Appendix B also contains additional results for: pixel-by-pixel CIFAR10 and a noise-padded version of CIFAR-10; as well as for character level and word level prediction using the Penn Tree Bank (PTB) dataset. All of these tasks require that the recurrent unit learns long-term dependencies: that is, the hidden-to-hidden matrices need to have sufficient memory to remember information from far in the past.

heading: 6.1 ORDERED AND PERMUTED PIXEL-BY-PIXEL MNIST
text: The pixel-by-pixel MNIST task tests long range dependency by sequentially presenting 784 pixels to the recurrent unit, i.e., the RNN processes one pixel at a time (Le et al., 2015). At the end of the
sequence, the learned hidden state is used to predict the class membership probability of the input image. This task requires that the RNN has a sufficient long-term memory in order to discriminate between different classes. A more challenging variation to this task is to operate on a fixed random permutation of the input sequence.
Table 1 provides a summary of our results. The Lipschitz RNN, with hidden dimension of N = 128 and trained with the forward Euler and RK2 scheme, achieves 99.4% and 99.3% accuracy on the ordered pixel-by-pixel MNIST task. For the permuted task, the model trained with forward Euler achieves 96.3% accuracy, whereas the model trained with RK2 achieves 96.2% accuracy. Hence, our Lipschitz recurrent unit outperforms state-of-the-art RNNs on both tasks and is competitive even when a hidden dimension of N = 64 is used, however, it can be seen that a larger unit with more capacity is advantageous for the permuted task. Our results show that we significantly outperform the Antisymmetric RNN (Chang et al., 2019) on the ordered tasks, while using fewer weights. That shows that the antisymmetric weight paramterization is limiting the expressivity of the recurrent unit. The exponential RNN is the next most competitive model, yet this model requires a larger hidden-to-hidden unit to perform well on the two considered tasks.

heading: 6.2 TIMIT
text: Next, we consider the TIMIT dataset (Garofolo, 1993) to study the capabilities of the Lipschitz RNN for speech prediction using audio data. For our experiments, we used the publicly available implementation of this task by Lezcano-Casado & Martinez-Rubio (2019). This implementation applies the preprocessing steps suggested by Wisdom et al. (2016): (i) downsample each audio sequence to 8kHz; (ii) process the downsampled sequences with a short-time Fourier transform using a Hann window of 256 samples and a window hop of 128 samples; and (iii) normalize the logmagnitude of the Fourier amplitudes. We obtain a set of frames that each have 129 complex-valued Fourier amplitudes and the task is to predict the log-magnitude of future frames. To compare our results with those of other models, we used the common train / validation / test split: 3690 utterances from 462 speakers for training, 192 utterances for validation, and 400 utterances for testing.
Table 2 lists the results for the Lipschitz recurrent unit as well as for several benchmark models. It can be seen that the Lipschitz RNN outperforms other state-of-the-art models for a fixed number of parameters (≈ 200K). In particular, LSTMs do not perform well on this task, however, the recently proposed momentum based LSTMs (Nguyen et al., 2020) have improvemed performance. Interestingly, the RK2 scheme leads to a better performance since this scheme provides more accurate approximations for the intermediate states.

heading: 7 ROBUSTNESS WITH RESPECT TO PERTURBATIONS
text: An important consideration beyond accuracy is robustness with respect to input and parameter perturbations. We consider a Hessian-based analysis and noise-response analysis of different continuous-time recurrent units and train the models on MNIST. Here, we reshape each MNIST thumbnail into sequences of length 98 so that each input has dimension x ∈ R8. We consider this
simpler problem so that all models obtain roughly the same training loss. Here we use stochastic gradient decent (SGD) with momentum to train the models.
Eigenanalysis of the Hessian provides a tool for studying various aspects of neural networks (Hochreiter & Schmidhuber, 1997a; Sagun et al., 2017; Ghorbani et al., 2019). Here, we study the Hessian H spectrum with respect to the model parameters of the recurrent unit using PyHessian (Yao et al., 2019). The Hessian provides us with insights about the curvature of the loss function L. This is because the Hessian is defined as the derivatives of the gradients, and thus the Hessian eigenvalues describe the change in the gradient of L as we take an infinitesimal step into a given direction. The eigenvectors span the (local) surface of the loss function at a given point, and the corresponding eigenvalue determines the curvature in the direction of the eigenvectors. This means that larger eigenvalues indicate a larger curvature, i.e., greater sensitivity, and the sign of the eigenvalues determines whether the curvature will be positive or negative.
To demonstrate the advantage of the additional linear term and our weight parameterization, we compare the Lipschitz RNN to two other continuous-time recurrent units. First, we consider a simple neural ODE RNN (Rubanova et al., 2019) that takes the form
ḣ = tanh(Wh+ Ux+ b), y = Dh, (11)
where W is a simple hidden-to-hidden matrix. As a second model we consider the antisymmetric RNN (Chang et al., 2019), that takes the same form as (11), but uses a skew-symmetric scheme to parameterize the hidden-to-hidden matrix as W := (M −MT )−γI , where M is a trainable weight matrix and γ is a tunable parameter.
Table 3 reports the largest eigenvalue λmax(H) and the trace of the Hessian tr(H).The largest eigenvalue being smaller indicates that our Lipschitz RNN found a flatter minimum, as compared to the simple neural ODE and Antisymmetric RNN. It is known that such flat minima can be perturbed without significantly changing the loss value (Hochreiter & Schmidhuber, 1997a). Table 3 also reports the condition number κ(H) := λmax(H)λmin(H) of the Hessian. The condition number κ(H) provides a measure for the spread of the eigenvalues of the Hessian. It is known that first-order methods can slow down in situations where κ is large (Bottou & Bousquet, 2008). The condition number and trace of our Lipshitz RNN being smaller also indicates improved robustness properties.
Next, we study the sensitivity of the response yT at time T in terms of the test accuracy with respect to a sequence of perturbed inputs {x̃1, . . . , x̃T } ∈ R8. We consider three different perturbations. The results for the artificially constructed perturbations are presented in Table 3, showing that the Lipschitz RNN is more resilient to adversarial perturbation. Here, we have considered the projected gradient decent (PGD) (Goodfellow et al., 2014) method with l∞, and the DeepFool
method (Moosavi-Dezfooli et al., 2016) with l2 and l∞ norm ball perturbations. We construct the adversarial examples with full access to the models, using 7 iterations. The step size for PGD is set to 0.01.
Further, Figure 2 shows the results for white noise and salt and pepper noise. It can be seen that the Lipschitz unit is less sensitive to input perturbations, as compared to the simple neural ODE RNN, and the antisymmetric RNN. In addition, we also show the results for an unitary RNN here.

heading: 7.1 ABLATION STUDY
text: The performance of the Lipschitz recurrent unit is due to two main innovations: (i) the additional linear term; and (ii) the scheme for constructing the hidden-to-hidden matrices A and W in Eq. (6). Thus, we investigate the effect of both innovations, while keeping all other conditions fixed. More concretely, we consider the following ablation recurrent unit
ht+1 = ht + α · ·Aht + · tanh(zt), with zt = Wht + Uxt + b, (12)
where α controls the effect of the linear hidden unit. Both A and W depend on the parameters β, γ.
Figure 3a studies the effect of the linear hidden unit, with β = 0.65 for the ordered task and β = 0.8 for the permuted task. In both cases we use γ = 0.001. It can be seen that the test accuracies of both the ordered and permuted pixel-by-pixel MNIST tasks clearly depend on the linear hidden unit. For α = 0, our models reduces to simple neural ODE recurrent units (Eq. (11)). The recurrent unit degenerates for α > 1.6, since the external input is superimposed by the hidden state. Figure 3b studies the effect of the hidden-to-hidden matrices with respect to β. It can be seen that β = {0.65, 0.70} achieves peak performance for the ordered task, and β = {0.8, 0.85} does so for the permuted task. Note that β = 1.0 recovers an skew-symmetric hidden-to-hidden matrix.

heading: 8 CONCLUSION
text: Viewing RNNs as continuous-time dynamical systems with input, we have proposed a new Lipschitz recurrent unit that excels on a range of benchmark tasks. The special structure of the recurrent unit allows us to obtain guarantees of global exponential stability using control theoretical arguments. In turn, the insights from this analysis motivated the symmetric-skew decomposition scheme for constructing hidden-to-hidden matrices, which mitigates the vanishing and exploding gradients problem. Due to the nice stability properties of the Lipschitz recurrent unit, we also obtain a model that is more robust with respect to input and parameter perturbations as compared to other continuoustime units. This behavior is also reflected by the Hessian analysis of the model. We expect that the improved robustness will make Lipschitz RNNs more reliable for sensitive applications. The theoretical results for our symmetric-skew decomposition of parameterizing hidden-to-hidden matrices also directly extend to the convolutional setting. Future work will explore this extension and study the potential advantages of these more parsimonious hidden-to-hidden matrices in combination with our parameterization in practice. Research code is shared via github.com/erichson/LipschitzRNN.

heading: ACKNOWLEDGMENTS
text: We would like to thank Ed H. Chi for fruitful discussions about physics-informed machine learning and the Antisymmetric RNN. We are grateful to the generous support from Amazon AWS and Google Cloud. NBE and MWM would like to acknowledge IARPA (contract W911NF20C0035), NSF, ONR and CLTC 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.

heading: A PROOFS

heading: A.1 PROOFS OF THEOREM 1 AND LEMMA 1
text: There are numerous ways that one can analyze the global stability of (4) through the related model (5), many of which are discussed in Zhang et al. (2014). Instead, here we shall conduct a direct approach and avoid appealing to diagonalization in order to obtain cleaner conditions, and a more straightforward proof that readily applies in the time-inhomogeneous setting.
Our method of choice relies on Lyapunov arguments summarized in the following theorem, which can be found as (Khalil, 2002, Theorem 4.10). For more details on related Lyapunov theory, see also Hahn (1967); Sastry (2013).
Theorem 2. An equilibrium h∗ for ḣ = f(t, h) is globally exponentially stable if there exists a continuously differentiable function V : [0,∞)×RN → [0,∞) such that for all h ∈ RN and t ≥ 0,
k1‖h− h∗‖α ≤ V (t, h) ≤ k2‖h− h∗‖α, and ∂V ∂t + ∂V ∂h ≤ −k3‖h− h∗‖α,
for some constants k1, k2, k3, α > 0. and V̇ (h) < 0 for h 6= h∗.
To simplify matters, we shall choose a Lyapunov function V : RN → [0,∞) that is independent of time. The most common type of Lyapunov function satisfying the conditions of Theorem 2 is of the form V (h) = (h− h∗)TP (h− h∗), where P is a positive definite matrix. One need only show that V̇ (h) ≤ −(h − h∗)TQ(h − h∗) for some other positive definite matrix Q to guarantee global exponential stability.
The construction of the Lyapunov function V that satisfies the conditions of Theorem 2 is accomplished using the Kalman-Yakubovich-Popov lemma, which is a statement regarding strictly positive real transfer functions. We use the following definition, equivalent to other standard definitions by (Khalil, 2002, Lemma 6.1). Definition 2. A function G : C→ CN×N is strictly positive real if it satisfies the following:
(i) The poles of G(s) have negative real parts.
(ii) G(iω) +G(−iω)T is positive definite for all ω ∈ R, where i = √ −1.
(iii) Either G(∞) + G(∞)T is positive definite or it is positive semidefinite and limω→∞ ω
2MT [G(iω) +G(−iω)T ]M is positive definite for any N × (N − q) full-rank matrix M such that MT [G(∞) +G(∞)T ]M = 0, where q = rank[G(∞) +G(∞)T ].
The following is presented in (Khalil, 2002, Lemma 6.3). Lemma 3 (Kalman-Yakubovich-Popov). Let A,W : RN → RN be full-rank square matrices. There exists a symmetric positive-definite matrix P and matrices L,U and a constant > 0 such that
PA+ATP = −LTL− P P = LTU −WT
UTU = 0,
if and only if the transfer function G(s) = W (sI − A)−1 is strictly positive real. In this case, we may take = 2µ, where µ > 0 is chosen so that G(s− µ) remains strictly positive real.
A shorter proof for case (a) is available to us through the (multivariable) circle criterion — the following theorem is a corollary of (Khalil, 2002, Theorem 7.1) suitable for our purposes. Theorem 3 (Circle Criterion). The system of differential equations
ḣ = Ah+ ψ(t,Wh)
is globally exponentially stable towards an equilibrium at the origin if ‖ψ(t, y)‖ ≤ M‖y‖ for some M > 0 and Z(s) = [I + MG(s)][I −MG(s)]−1 is strictly positive real, where G(s) = W (sI −A)−1.
Both the Kalman-Yakubovich-Popov lemma and the circle criterion are classical results in control theory, and are typically discussed in the setting of feedback systems (Khalil, 2002, Chapter 6, 7). Our presentation here is less general than the complete formulation, but makes clearer the connection to RNNs. With these tools, we state our proof of Theorem 1.
Proof of Theorem 1. To begin, we shall center the differential equation about the equilibrium. By assumption, there exists h∗ such that Ah∗ = −σ(Wh∗ + Ux(t) + b). Letting h̄ = h− h∗, we find that
˙̄h = Ah+ σ(Wh+ Ux(t) + b)
= Ah̄+Ah∗ + σ(Wh̄+Wh∗ + Ux(t) + b)
= Ah̄+ σ(Wh̄+Wh∗ + Ux(t) + b)− σ(Wh∗ + Ux(t) + b). (13) It will suffice to show that (13) is globally exponentially stable at the origin.
Let us begin with case (a). The proof follows arguments analogous to (Khalil, 2002, Example 7.1). Let G(s) = W (A− sI)−1 denote the transfer function for the system (13). Letting
ψ(t, x) = σ(x+Wh∗ + Ux(t) + b)− σ(Wh∗ + Ux(t) + b), since σ is M -Lipschitz, we know that ‖ψ(t, x)‖ ≤ M‖x‖ for any x ∈ RN . Therefore, let Z(s) = [I + MG(s)][I −MG(s)]−1 denote the transfer function in the circle criterion. Our objective is to show that Z(s) is strictly positive real — by Theorem 3, this will guarantee the desired global exponential stability of (4). First, we need to show that the poles of Z(s) have negative real parts. This can only occur when G(s) itself has poles or I −MG(s) is singular. The former case occurs precisely where A − sI is singular, which occurs when s is an eigenvalue of A. Since A + AT is assumed to be negative definite, A must have eigenvalues with negative real part by Lemma 4, and so the poles of G(s) also have negative real parts. The latter case is more difficult to treat. First, since σmax(AB) ≤ σmax(A)σmax(B) and σmax(B−1) = σmin(B)−1,
σmax(G(s)) ≤ σmax(W )
σmin(A− sI) . (14)
Therefore, we observe that
σmin(I −MG(s)) ≥ 1− σmax(MG(s)) ≥ 1−Mσmax(G(s))
≥ 1− Mσmax(W ) σmin(A− sI) .
From the Fan-Hoffman inequality (Bhatia, 2013, Proposition III.5.1), we have that σmin(A− sI) = σmin(sI −A) ≥ λmin ( <(s)I − A+A T
2
) = <(s) + λmin ( −A+A T
2
) ,
and since A+AT is negative definite, for any s with <(s) ≥ 0, σmin(A− sI) ≥ <(s) + σmin ( A+AT
2
) ≥ σmin(Asym). (15)
Since σmin(Asym) > Mσmax(W ), it follows that σmin(I − MG(s)) > 0 whenever s has nonnegative real part, and so the poles of Z(s) must have negative real parts.
Next, we need to show that Z(iω) + Z(−iω)T is positive definite for all ω ∈ R. Observe that
Z(iω) + Z(−iω)T = [I +MG(iω)][I −MG(iω)]−1 + [I −MG(−iω)T ]−1[I +MG(−iω)T ] = 2[I −MG(−iω)T ]−1[I −M2G(−iω)TG(iω)][I −MG(iω)]−1.
From Sylvester’s law of inertia, we may infer that Z(iω) +Z(−iω)T is positive definite if and only if I + Yω is positive definite, where Yω = M2G(−iω)TG(iω). If we can show that the eigenvalues of Yω lie strictly within the unit circle, that is, σmax(Yω) < 1 for all ω ∈ R, then I + Yω will necessarily be positive definite. From (14) and (15), we may verify that
sup ω∈R σmax(G(iω)) ≤ sup ω∈R
σmax(W ) σmin(A− iωI) ≤ σmax(W ) σmin(Asym) .
Therefore, σmax(Yω) ≤M2σmax(G(−iω)T )σmax(G(iω)) ≤ ( Mσmax(W )
σmin(Asym)
)2 < 1,
by assumption. Finally, since Z(∞) + Z(∞)T = 2I is positive definite, Z(s) is strictly positive real and Theorem 3 applies.
Now, consider case (b). The proof proceeds in two steps. First, we verify that the transfer function G(s) = W (A − sI)−1 satisfies the conditions of the Kalman-Yakubovich-Popov lemma. Then, using the matrices P , L, U , and the constant inferred from the lemma, a Lyapunov function is constructed which satisfies the conditions of Theorem 2, guaranteeing global exponential stability. Once again, condition (i) of Lemma 3 is straightforward to verify: G(s) exhibits poles when s is an eigenvalue of A, and so the poles of G(s) also have negative real parts. Furthermore, condition (iii) is easily satisfied with M = I since G(∞) + G(∞)T = 0. To show that condition (ii) holds, observe that for any ω ∈ R, letting A−T = (A−1)T for brevity,
G(iω) +G(−iω)T = W (A− iωI)−1 + (A+ iωI)−TWT
= (A+ iωI)−T [(A+ iωI)TW +WT (A− iωI)](A− iωI)−1.
Since the inner matrix factor is Hermitian, Sylvester’s law of inertia implies that G(iω) +G(−iω)T is positive definite if and only if
Bω := (A+ iωI) TW +WT (A− iωI).
is positive definite. SinceBω is a Hermitian matrix, it has real eigenvalues, with minimal eigenvalue given by the infimum of the Rayleigh quotient:
λmin(Bω) = inf ‖v‖=1
vTBωv
= inf ‖v‖=1
vT (ATW +WTA)v + iωvT (W −WT )v
= inf ‖v‖=1
vT (ATW +WTA)v
= λmin(A TW +WTA).
By assumption, ATW + WTA has strictly positive eigenvalues, and hence Bω and G(iω) + G(−iω)T are positive definite. Therefore, Lemma 3 applies, and we obtain matrices P,L, U and a constant > 0 with the corresponding properties.
Now we may construct our Lyapunov function V . Let v = Wh̄ and
u(t) = σ(v(t) +Wh∗ + Ux(t) + b)− σ(Wh∗ + Ux(t) + b),
so that ˙̄h = Ah̄ + u. Since σ is monotone non-decreasing, σ(x) − σ(y) ≥ 0 for any x ≥ y. This implies that for each i = 1, . . . , N , vi and ui have the same sign. In particular, vTu ≥ 0. Now, let V (h) = hTPh be our Lyapunov function, noting that V is independent of t. Taking the derivative of the Lyapunov function over (13) and using the properties of P,L, U, ,
V̇ (h̄) = h̄TP ˙̄h+ ˙̄hTPh̄
= h̄T (PA+ATP )h̄+ 2h̄TPu
= h̄T (−LTL− P )h̄+ 2h̄T (LTU −WT )u = −(Lh̄)T (Lh̄) + (Lh̄)TUu+ (Uu)T (Lh̄)− uTUTUu− 2vTu = −(Lh̄+ Uu)T (Lh̄+ Uu)− h̄TPh̄− 2vTu.
Since vTu ≥ 0 and (Lh̄ + Uu)T (Lh̄ + Uu) ≥ 0, it follows that V̇ (h̄) ≤ − λmin(P )‖h‖2, and hence global exponential stability follows from Theorem 2 and positive-definiteness of P .
To finish off discussion regarding the results from Sec. 3, we provide a quick proof of Lemma 1 using a simple diagonalization argument.
Proof of Lemma 1. Since A is symmetric and real-valued, by (Horn & Johnson, 2012, Theorem 4.1.5), there exists an orthogonal matrix P and a real diagonal matrix D such that A = PDPT . Letting z = PTh where h satisfies (4), since h = Pz, we see that
ż = PTPDPTh+ PTσ(Wh+ Ux+ b)
= Dz + PTσ(WPz + Ux+ b).
Therefore, z satisfies (5) with L = PT and V = WP , both of which are nonsingular by orthogonality of P . By the same argument, for any equilibrium h∗, taking z∗ = PTh∗,
Dz∗ + PTσ(WPz∗ + Ux+ b) = PT (PDPTh∗ + σ(Wh∗ + Ux+ b))
= PT (Ah∗ + σ(Wh∗ + Ux+ b)) = 0,
implying that z∗ is an equilibrium of (5). Furthermore, since
‖z − z∗‖2 = (PTh− PTh∗)T (PTh− PTh∗) = (h− h∗)TPPT (h− h∗) = ‖h− h∗‖2,
from orthogonality of P . Because every form of Lyapunov stability, both local and global, including global exponential stability, depend only on the norm ‖h − h∗‖ (Khalil, 2002, Definitions 4.4 and 4.5), h∗ is stable under any of these forms if and only if z∗ is also stable.
We remark that the proof of Lemma 1 can extend to matrices A which have real eigenvalues and are diagonalizable. These attributes are implied for symmetric matrices. However, they can be difficult to ensure in practice for nonsymmetric matrices without imposing difficult structural constraints.

heading: A.2 PROOF OF PROPOSITION 1
text: The proof of Proposition 1 relies on the following lemma, which we also have made use of several times throughout this work. Lemma 4. For any matrix A ∈ RN×N , the real parts of the eigenvalues <λi(A) are contained in the interval [λmin(Asym), λmax(Asym)], where Asym = 12 (A+A T ).
Proof. Recall by the min-max theorem, for 〈u, v〉 = u∗v, where u∗ is the conjugate transpose of u, the upper and lower eigenvalues of A+AT satisfy
λmin(A+A T ) = inf v∈CN , ‖v‖=1 〈v, (A+AT )v〉 = inf v∈CN , ‖v‖=1 〈v,Av〉+ 〈Av, v〉,
λmax(A+A T ) = sup v∈CN , ‖v‖=1 〈v, (A+AT )v〉 = sup v∈CN , ‖v‖=1 〈v,Av〉+ 〈Av, v〉.
Let λi(A) = u + iω be an eigenvalue of A with corresponding eigenvector v satisfying ‖v‖ = 1. Since Av = (u+ iω)v,
〈v,Av〉+ 〈Av, v〉 = 〈v,Av〉+ 〈v,Av〉 = 2<〈v,Av〉 = 2u‖v‖2 = 2u.

heading: Hence, λmin(A+AT ) ≤ u ≤ λmax(A+AT ).
text: Proof of Proposition 1. By construction, Ssymβ,γ = Sβ,γ + S T β,γ = (1− β)M sym − γI, and so from Lemma 4, both the real parts <λi(Sβ,γ) of the eigenvalues of Sβ,γ as well as the eigenvalues of Ssymβ,γ lie in the interval
[λmin(S sym β,γ ), λmax(S sym β,γ )] = [λmin((1− β)M sym − γI), λmax((1− β)M sym − γI)]. If β < 1, for any eigenvalue λ of Ssymβ,γ with corresponding eigenvector v,
(1− β)M symv − γv = λv, and so M symv = λ+ γ 1− β v
implying that λ+γ1−β is an eigenvalue of M sym, and therefore contained in [λmin(M sym), λmax(M
sym)]. In particular, we find that [λmin(S sym β,γ ), λmax(S sym β,γ )] ⊆ [(1− β)λmin(M
sym)− γ, (1− β)λmax(M sym)], (16) as required. Finally, if β = 1, then (16) still holds, since both intervals collapse to the single point {−γ}.
Figure 4 illustrates the effect of β onto the eigenvalues of Aβ,γ with the largest and smallest real parts. It can be seen, both empirically and theoretically, that the real part of the eigenvalues converges towards zero as β tends towards one, i.e., we yield a skew-symmetric matrix with purely imaginary eigenvalues in the limit. Thus, for a sufficiently large parameter β we yield a system that approximately preserves an “energy” for a limited time-horizon
Rλi(Aβ,γ) ≈ 0, for i = 1, 2, . . . , N. (17)

heading: A.3 PROOF OF LEMMA 2
text: First, it follows from Gronwall’s inequality that the norm of the final hidden state ‖h(T )‖ is bounded uniformly in β. From Weyl’s inequalities and the definition of Aβ,γ ,
max k |∆δλk(Asymβ,γ )| ≤ ‖∆δA sym β,γ ‖ = (1− β)‖∆δM sym A ‖.
By the chain rule, for each element M ijA of the matrix MA,
∂L
∂M ijA =
∂L
∂y(T )
∂y(T )
∂h(T )
∂h(T ) ∂M ijA = ∂L ∂y(T ) D ∂h(T ) ∂M ijA .
Now, for any collection of parameters θi,
d
dt ∑ i ∂h ∂θi = A ∑ i ∂h ∂θi + ∑ i ∂A ∂θi h+ sech2 (Wh+ Ux+ b) ( W ∑ i ∂h ∂θi + ∑ i ∂W ∂θi h ) ,
and from Gronwall’s inequality,∥∥∥∥∥∑ i ∂h(T ) ∂θi ∥∥∥∥∥ ≤ (∥∥∥∥∥∑ i ∂Aβ,γ ∂θi ∥∥∥∥∥+ ∥∥∥∥∥∑ i ∂Wβ,γ ∂θi ∥∥∥∥∥ ) ‖h‖ e(‖Aβ,γ‖+‖Wβ,γ‖)T .
Since ∆δM sym A = δ ∂L ∂MA
+ δ (
∂L ∂MA
)T ,
‖∆δM symA ‖ ≤ ‖∆δM sym A ‖F
≤ δ √√√√∑ i,j ( ∂L ∂M ijA + ∂L ∂M ijA )2
≤ δ ∥∥∥∥∂L∂y ∥∥∥∥ ‖D‖ ‖h‖ e(‖Aβ,γ‖+‖Wβ,γ‖)T √√√√∑
i,j ∥∥∥∥∥∂Aβ,γ∂M ijA + ∂Aβ,γ∂M jiA ∥∥∥∥∥ 2 .
Since ∂(MAh) ∂MijA
= ∂(MTAh)
∂MjiA , it follows that
∂Aβ,γ ∂M ijA + ∂Aβ,γ ∂M jiA = 2(1− β)
( ∂(MAh)
∂M ijA + ∂(MTAh) ∂M jiA
) ,
and so ‖∆δM symA ‖ = O(δ(1−β)), and therefore maxk |∆δσk(A sym β,γ )| = O(δ(1−β)2). Similarly, for the matrix MW ,
max k ∣∣∣∆δλk(W symβ,γ )∣∣∣ ≤ (1− β) ‖∆δM symW ‖ ≤ δ(1− β) ∥∥∥∥∂L∂y ∥∥∥∥ ‖D‖ ‖h‖ e(‖Aβ,γ‖+‖Wβ,γ‖)T √√√√∑ i,j ∥∥∥∥∥∂Wβ,γ∂M ijW + ∂Wβ,γ∂M jiW ∥∥∥∥∥ 2
= 2δ(1− β)2 ∥∥∥∥∂L∂y ∥∥∥∥ ‖D‖ ‖h‖ e(‖Aβ,γ‖+‖Wβ,γ‖)T √√√√∑
i,j
( ∂(MWh)
∂M ijW + ∂(MTWh) ∂M jiW
)2 ,
and hence maxk |∆δλk(W symβ,γ )| = O(δ(1− β)2).
In Figure 5, we plot the most positive real part of the eigenvalues of Aβ,γ and Wβ,γ during training for the ordered MNIST task. As β increases, the eigenvalues change less during training, remaining in the stability region provided by case (b) of Theorem 1 for more of the training time.

heading: B ADDITIONAL EXPERIMENTS

heading: B.1 SENSITIVITY TO RANDOM INITIALIZATION FOR MNIST AND TIMIT
text: The hidden matrices are initialized by sampling weights from the normal distributionN (0, σ), where σ is the variance, which can be treated as a tuning parameter. In our experiments we typically chose a small σ; see the Table 8 for details. To show that the Lipschitz RNN is insensitive to random initialization, we have trained each model with 10 different seeds. Table 4 shows the maximum, average and minimum values obtained for each task. Note that higher values indicate better performance on the ordered and permuted MNIST tasks, while lower values indicate better performance on the TIMIT task.

heading: B.2 ORDERED PIXEL-BY-PIXEL AND NOISE-PADDED CIFAR-10
text: The pixel-by-pixel CIFAR-10 benchmark problem that has recently been proposed by (Chang et al., 2019). This task is similar to the pixel-by-pixel MNIST task, yet more challenging due to the increased sequence length and the more difficult classification problem. Similar to MNIST, we flatten the CIFAR-10 images to construct a sequence of length 1024 in scanline order, where each element of the sequence consists of three pixels (one from each channel).
A variation of this problem is the noise-padded CIFAR-10 problem (Chang et al., 2019), where we consider each row of an image as input at time step t. The rows from each channel are stacked so that we obtain an input of dimension x ∈ R96. Then, after the 32 time step which process the 32
row, we start to feed the recurrent unit with independent standard Gaussian noise for 968 time steps. At the final point in T = 1000, we use the learned hidden state for classification. This problem is challenging because only the first 32 time steps contain signals. Thus, the recurrent unit needs to recall information from the beginning of the process.
Table 5 provides a summary of our results. Our Lipschitz recurrent unit outperforms both the incremental RNN (Kag et al., 2020) and the antisymmetric RNN (Chang et al., 2019) by a significant margin. This impressively demonstrates that the Lipschitz unit enables the stable propagation of signals over long time horizons.

heading: B.3 PENN TREE BANK (PTB)

heading: B.3.1 CHARACTER LEVEL PREDICTION
text: Next, we consider a character level language modeling task using the Penn Treebank Corpus (PTB) (Marcus et al., 1993). Specifically, this task studies how well a model can predict the next character in a sequence of text. The dataset is composed of a train / validation / test set, where 5017K characters are used for training, 393K characters are used for validation and 442K characters are used for testing. For our experiments, we used the publicly available implementation of this task by Kerg et al. (2019), which computes the performance in terms of mean bits per character (BPC).
Table 6 shows the results for back-propagation through time (BPTT) over 150 and 300 time steps, respectively. The Lipschitz RNN performs slightly better then the exponential RNN and the nonnormal RNN on this task. (Kerg et al., 2019) notes that orthogonal hidden-to-hidden matrices are not particular well-suited for this task. Thus, it is not surprising that the Lipschitz unit has a small advantage here.
For comparison, we have also tested the Antisymmetric RNN (Chang et al., 2019) on this task. The performance of this unit is considerably weaker as compared to our Lipschitz unit. This suggests that the Lipschitz RNN is more expressive and improves the propagation of meaningful signals over longer time scales.

heading: B.3.2 WORD-LEVEL PREDICTION
text: In addition to character-level prediction, we also consider word-level prediction using the PTB corpus. For comparison with other state-of-the-art units, we consider the setup by Kusupati et al. (2018), who use a sequence length of 300. Table 7 shows results for back-propagation through time (BPTT) over 300 time steps. The Lipschitz RNN performs slightly better than the other RNNs on this task and the baseline LSTM for the test perplexity metric reported by Kusupati et al. (2018).

heading: C TUNING PARAMETERS
text: For tuning we utilized a standard training procedure using a non-exhaustive random search within the following plausible ranges for the our weight parameterization β = 0.65, 0.7, 0.75, 0.8, γ = [0.001, 1.0]. For Adam we explored learning rates between 0.001 and 0.005, and for SGD we considered 0.1. For the step size we explored values in the range 0.001 to 1.0. We did not perform an automated grid search and thus expect that the models can be further fine-tuned.
The tuning parameters for the different tasks that we have considered are summarized in Table 8.
For pixel-by-pixel MNIST and CIFAR-10, we use Adam for minimizing the objective. We train all our models for 100 epochs, with scheduled learning rate decays at epochs {90}. We do not use gradient clipping during training. Figure 6 shows the test accuracy curves for our Lipschitz RNN for the ordered and permuted MNIST classification tasks.
For TIMIT we use Adam with default parameters for minimizing the objective. We also tried Adam using betas (0.0, 0.9) as well as RMSprop with α = 0.9, however, Adam with default values worked best in our experiments. We train the model for 1200 epochs without learning-rate decay. Similar to Kerg et al. (2019) we train our model with gradient clipping, however, we observed that the performance of our model is relatively insensitive to the clipping value.
For the character level prediction task, we use Adam with default parameters for minimizing the objective, while we use RMSprop with α = 0.9 for the word level prediction task. We train the model for 200 epochs for the character-level task, and for 500 epochs for the word-level task.

Review summaries:
Considering a continuous time RNN with Lipschitz-continuous nonlinearity, the authors formulate sufficient conditions on the parameter matrices for the network to be globally stable, in the sense of a globally attracting fixed point. They provide a specific parameterization for the hidden-to-hidden weight matrices to control global stability and error gradients, consisting of a weighted combination of a symmetric and a skew-symmetric matrix (and some diagonal offset). The authors discuss numerical integration by forward-Euler and RK2, and thoroughly benchmark their approach against a large set of other state-of-the-art RNNs on various tasks including versions of MNIST and TIMIT.  Finally, they highlight improved stability of their RNN against parameter and input perturbations.