On Lyapunov Exponents for RNNs:
Understanding Information Propagation Using Dynamical Systems Tools

The propagation of error gradients in deep learning leads to the study of recursive compositions and their stability [1]. Vanishing and exploding gradients arise from long products of Jacobians of the hidden state dynamics whose norm exponentially grows or decays, hindering training [2]. To mitigate this sensitivity, much effort has been made to mathematically understand the link between model parameters and the the eigen- and singular-value spectra of these long products [3, 4, 5]. For architectures used in practice, this approach appears to have limited a scope [6, 7] likely due to spectra having non-trivial properties reflecting complicated long-time dependencies within the trajectory. Fortunately, the theory of dynamical systems has formulated this problem for general dynamical systems, i.e., for arbitrary architectures. The Lyapunov spectrum (LS) formed by the Lyapunov exponents (LEs) is precisely the measurement that characterizes rates of expansion and contraction of nonlinear system dynamics.

While there are results using LEs in machine learning [8, 9, 10], most only look at the largest exponent, the sign of which indicates the presence or absence of chaos. Additionally, to our knowledge, current machine learning uses of LEs are for the autonomous case, i.e. , analyzing RNN dynamics in the absence of inputs, a regime at odds with the most common use of RNNs as input-processing systems. Thus, there is likely much to be gained from both the study of other features of the LS and from including input-driven regimes into LEs analysis. Indeed, a concurrent work coming from theoretical neuroscience takes up this approach [11] and outlines the ongoing work using LEs to study models of neural networks in the brain. However, the machinery to compute the LS and the theory surrounding its interpretation is still relatively unknown in the machine learning community.

In this exposition-style paper, we aim to fill this gap by presenting an overview of LEs in the context of RNNs, and discuss their usefulness for studying training dynamics. We present encouraging results which suggest that LEs could serve in efficient hyperparameter tuning by allowing early detection of performance. To compute LEs, we present a novel algorithm for various RNN architectures, taking advantage of modern deep learning environments. Finally, we highlight key directions of research that have the potential to leverage LEs for improvement of RNNs robustness and performance.

Here we develop the problem of spectral constraints for robust gradient propagation. For transparent exposition, in this section we will consider the "vanilla" RNN defined as

where $\mathbf{V}$ is the recurrent weight matrix that couples elements of the hidden state vector $\mathbf{h}_{t}\in\mathbb{R}^{N}$, $\mathbf{U}$ projects the input, $\mathbf{x}_{t}$, into the network, $\mathbf{b}$ is a constant bias vector, and $\phi$ applies a nonlinear scalar transformation element-wise. The readout weights, $\mathbf{W}$, output the activity into the output variable, $\mathbf{o}_{t}$. The loss over $T$ iterates is $L=\sum_{t=1}^{T}L_{t}$, with $L_{t}=f(\mathbf{y}_{t},\hat{\mathbf{y}}_{t})$, with $f$ some scalar loss function (e.g., cross-entropy loss, $f(\mathbf{y},\hat{\mathbf{y}})={\mathbf{y}}\cdot\log\hat{\mathbf{y}}$), $\hat{\mathbf{y}}_{t}$ the prediction (e.g., $\hat{\mathbf{y}_{t}}=\mathrm{softmax}(\mathbf{o}_{t})$), and $\mathbf{y}_{t}$ is a one-hot binary target vector. The parameters, $\Theta=(\mathbf{V},\mathbf{U},\mathbf{b},\mathbf{W})$, are learned by following the gradient of the loss, e.g. in the space of recurrent weights $\mathbf{V}$,

where $\textrm{diag}(\mathbf{x})$ is the diagonal matrix formed by the vector $\mathbf{x}$ and $\phi^{\prime}$ is the derivative of $\phi$ and ^⊤ denotes transpose. Here,

where $\nabla_{\mathbf{o}_{s}}L$ is some simple expression (e.g. $\hat{\mathbf{y}}-\mathbf{y}_{t}$ for cross-entropy loss) and $\mathbf{J}_{t}=\frac{\partial\mathbf{h}_{t}}{\partial\mathbf{h}_{t-1}}$ is the Jacobian of the hidden state dynamics,

$\mathbf{J}_{t}$ varies in time with $\mathbf{x}_{t}$ and $\mathbf{h}_{t-1}$ via $\mathbf{a}_{t}$ and so is treated as a random matrix with ensemble properties arising from the specified input statistics and the emergent hidden state statistics.

Gradients can vanish or explode with $T$ according to the singular value spectra of the products of Jacobians appearing in Eq. (3) [1]. This fact motivated the study of hyperparameter constraints that control the average of the latter over the input distribution. Specifically, the scale parameter of Gaussian or orthogonal initializations of $\mathbf{V}$ is chosen such that the condition number (ratio of largest to smallest singular value) of $\mathbf{J}_{t}$ is bounded around 1 for all $t$ [2], or such that the first and second moments of the distribution of squared singular values (averaged over inputs) of the Jacobian products are chosen near 1 and 0, respectively according to dynamical isometry [4]. Parameter conditions for the latter have been derived for i.d.d. Gaussian input. The approach has been extended to RNNs under the assumption of untied weights with good empirical correspondence in vanilla and minimalRNN [3], but larger discrepancy in LSTMs [12], suggesting there are other contributions to stability in more complex, state-of-the-art architectures [6, 7].

Better understanding how robust gradient propagation emerges in these less idealized settings that include high-dimensional and temporally correlated input sequences demands a more general theory. These Jacobian products appearing in Eq. (3) can each be expressed as the transpose of a forward sequence, $\prod_{r=t+1}^{s}\mathbf{J}^{\top}_{r}=\left(\mathbf{J}_{s}\cdots\mathbf{J}_{t+1}\right)^{\top}$. Note that the backward time gradient dynamics shares properties with the forward time hidden state dynamics. This is an important insight that suggests the use of advanced tools from dynamical systems theory may help better understand gradient propagation.

Here we describe the stochastic Lyapunov exponents from the ergodic theory of non-autonomous dynamical systems and outline their connection to the conditions that support gradient-based learning in RNNs. Jacobians are linear maps that evolve state perturbations forward along a trajectory, $\mathbf{u}_{t}=\mathbf{J}_{t}\mathbf{u}_{t-1}$, e.g. for some initial perturbation $\mathbf{h}_{0}(\epsilon)=\mathbf{h}_{0}+\epsilon\mathbf{u}_{0}$ with $\epsilon\ll 1$ and a given vector $\mathbf{u}_{0}$. Thus, perturbations can vanish or explode under the hidden state dynamics according to the singular value spectrum of $\mathbf{T}_{t}=\mathbf{J}_{t}\cdots\mathbf{J}_{1}$. The linear stability of the dynamics is obtained taking $\epsilon\to 0$ and then $t\to\infty$. In particular, the Lyapunov exponents, $\{\lambda_{i}\}_{i=1}^{N}$, are the exponential growth rates associated with the singular values of $\mathbf{T}_{t}^{1/t}$ for $t\to\infty$, i.e. the logarithms of the eigenvalues of $\lim_{t\to\infty}\left(\mathbf{T}_{t}^{\top}\mathbf{T}_{t}\right)^{\frac{1}{2t}}$. As a property of the stationary dynamics, the Lyapunov exponents are independent of initial condition for ergodic systems, i.e. those with only one or the same type of attractor. If the maximum Lyapunov exponent $\lambda_{\textrm{max}}$ is positive, the stationary dynamics is chaotic and small perturbations explode, otherwise it is stable and small perturbations vanish. The theory was initially developed for autonomous dynamical systems, where stable dynamics implies limit-cycle or fixed-point attractors. How the shape of the Lyapunov spectrum varies with network model parameters has provided unprecedented insight into autonomous ( evolving in the absence of inputs) neural network models in theoretical neuroscience [13, 11, 14].

RNN dynamics are non-autonomous because the hidden state dynamics $\mathbf{h}_{t}$ are driven by inputs $\mathbf{x}_{t}$. The theory of random dynamical systems generalizes stability analyses to non-autonomous dynamics driven by an input sequence sampled from a stationary distribution [15]. Analytical results typically employ uncorrelated Gaussian inputs, but the framework is expected to apply to a wider range of well-behaved input statistics. This includes those with finite, low-order moments and finite correlation times like character streams from written language and sensor data from motion capture systems. The time course of the driven dynamics depends on the specific input realization but, critically, the theory guarantees that the stationary dynamics for all input realizations share the same stability properties, which will in general depend on the input distribution (e.g., its variance). Stability here is quantified using the same definition of Lyapunov exponents in the autonomous case (now called stochastic Lyapunov exponents; we hereon will omit stochastic for brevity). In Fig. 1, we show vanilla RNN dynamics in stable, marginally stable, and chaotic regimes, as measured by the sign of the largest LE. For stable driven dynamics, the stationary activity on the time-dependent attractor (called a random sink) is independent of initial condition. This holds true also in the marginally stable case, $\lambda_{max}\approx 0$, where in addition there are directions along which the sizes of perturbations (error gradients) are maintained over many iterates. For chaotic driven dynamics, the activity variance over initial conditions fluctuates in time in a complex way. Understanding stability properties in this non-autonomous setting is at the frontier of analysis of neural network dynamics in both theoretical neuroscience [11, 16, 17] and machine learning [18, 19].

The practical calculation of the Lyapunov spectra (see Sec. 4) is fast and offers more robust numerical behaviour and generality over singular value decomposition. Together, this suggests the Lyapunov spectra as a useful diagnostic for RNN sequence learning.

There are a handful of features of the Lyapunov spectrum that dictate gradient propagation and thus, influence training speed and performance. The maximum LE determines the linear stability. The mean exponent determines the rate of contraction of full volume elements. The LE variance measures heterogeneity in stability across different directions and can reflect the conditioning of the product of many Jacobians. We use the first two in the analysis of our experiments later in the paper.

We remind readers that much of the theory regarding high-dimensional dynamics is derived in the so-called thermodynamic limit $N\to\infty$. Here, the number of exponents diverges and the spectrum over $i/N$ becomes stationary (so called extensive dynamics), i.e., it is insensitive to $N$ so long as $N$ is large enough (e.g. [13, 11]). Here, self-averaging effects can take hold and enable accurate analytical results based on Gaussian assumptions justified by central limit theorem arguments. ‘Large enough’ $N$ is typically in the hundreds (though this should be checked for any given application) and so these results can be useful for studying RNNs.

Finally, we note that obtaining the Lyapunov exponents from their definition relates to the singular values of $\mathbf{T}_{t}^{1/t}$ for large $t$, and thus is numerically impractical. The standard approach [20] is to exploit the fact that $m$-dimensional volume elements grow at a rate $\lambda^{(m)}=\sum_{i=1}^{m}\lambda_{i}$ and so the desired rates, $\lambda_{1}=\lambda^{(1)}$, $\lambda_{2}=\lambda^{(1)}-\lambda_{1}$, … arise from volumes obtained by projecting orthogonal to subspaces of increasing dimension. These rates are a direct output of computationally efficient orthonormalization procedures such as QR-decomposition. Next, we discuss the practical aspects of computing the Lyapunov exponents in data-driven RNNs.

We extend this framework of Lyapunov exponents to recurrent neural networks by calculating the asymptotic trajectories of the hidden states of the networks when driven by the same signal.

To illustrate the computation and application of Lyapunov exponents, we applied the approach to two tasks with distinct task constraints: character prediction from sentences, where performance depends on capturing low-dimensional and long temporal correlations, and signal prediction from motion capture data, where performance relies also on signal correlations across the many dimensions of the input.

In this paper, we presented an exposition and example application of Lyapunov exponents for understanding training stability in RNNs. We motivated them as a natural quantity related to stability of dynamics and useful in a natural complementary approach to existing mathematical approaches for understanding training stability focused on the singular value spectrum. We adapted the standard algorithm for their computation to RNNs, and showed how it could be made more efficient in standard machine learning development environments. We demonstrated that even the basic implementation runs almost 2 orders of magnitude faster than typical training time (see [11] for further precision and efficiency considerations) and that it converges relatively quickly with learning time. The latter implies it could be useful as a readout of performance early in training. To test this application, we studied it in two tasks with distinct constraints. In both we find interpretable variation with performance reflecting these distinct constraints.

There are a few points to be made about LEs in RNNs. First, one should consider carefully the degrees of freedom in the architecture under consideration. We have considered the hidden state dynamics of an LSTM such that there are as many Lyapunov exponents as there are units. LSTMs have gating variables that can also carry their own dynamics and could be considered as degrees of freedom, adding another $N$ exponents to the spectrum. While we were able to obtain evidence from the spectra of hidden states alone, recent work [25, 11] suggests that studying the stability in the subspace of gating variables can also be informative. A second point is that one should also consider the role of the input weights $\mathbf{U}$ (c.f. Eq. 1). The strength of inputs (e.g. input signal variance) that drive high-dimensional dynamics has long been known to have a stabilizing effect [26, 16, 27]. Thus, the input statistics of the task’s data can play a role that is modulated by $\mathbf{U}$, making the latter a target for gradient-based algorithms aiming to decrease loss. Due to limited space, we have not analyzed these weights here, but believe it is important to do so in applications. Finally, we raise perhaps the most glaring adaptation: that of approximating this asymptotic quantity to settings of finite-length sequences. We proposed to use averaging over inputs, in addition to averaging over initial conditions, as a way to take advantage of batch tensor mechanics to achieve faster LE estimates. This is in contrast to long-time averages used in the LE definition, although our method is justified under assumptions of ergodicity and short transient dynamics. Analyzing the influence of input batch size on this precision is a topic of future research.

The task analysis we have performed suggests a use case for the Lyapunov spectrum as having features that serve as an early readout of performance useful to speed up hyperparameter search. We looked the maximum and mean Lyapunov exponent as two features highlighted by hypotheses we made based on our knowledge of the task constraints. Presumably, there are others relevant here and in other tasks. The search for generic features that do not require knowledge of task constraints (e.g. dynamical isometry) is important, as is demonstrating how useful these features are in practise with complex sequence data. One limitation is how quickly they converge with learning. For example, as a composite quantity, the mean LE converges rather quickly, while the maximum LE is significantly slower to converge (see supplemental material). In [11], it is shown that loss function/learning rule combinations can also sculpt the spectra and alter its statistics, demonstrating that these are more sources of variation to understand. Interestingly, Full Force [28], an alternative to Back-propagation-through-time, has a qualitatively similar spectra (same max, min, and mean LE) but with much lower variance, presumably a desirable feature for generalization. This method also demands more information about performance than just a gradient, suggesting the interesting hypothesis that the precision with which the features of the Lyapunov spectra can be sculpted can scale with the amount of information provided in training. Finally, recent work has given theoretical grounding to why LSTMs avoid vanishing gradients [25]. Of course this was the reason for the design of LSTM so it only recapitulates the original design intuition. Looking forward however as more complex architectures are developed with less interpretable design, this approach based on LEs can still provide novel insight into why they work.

We close with a short discussion of open problems as this is a prominent area in understanding network dynamics at the intersection of machine learning and computational neuroscience. Having supporting analytical results is essential for robust control of complex problems. To this end, understanding how and where the assumption of untied weights breaks down and how forward propagation differs from backward propagation will be important in extending analytical results on spectral constraints. The Lyapunov spectrum can guide this work. Last, we believe insightful parallels in theoretical neuroscience work. For example, Lyapunov spectrum have been derived for additional degrees of freedom in the unit dynamics, different connectivity ensembles, and more. Also, a discrepancy between the loss of linear stability and the onset of chaotic dynamics in driven systems must be understood [27]. Making these connections explicit will serve both fields.