Year: 2021

Abstract:
Recurrent neural networks (RNNs) are wide-spread machine learning tools for modeling sequential and time series data. They are notoriously hard to train because their loss gradients backpropagated in time tend to saturate or diverge during training. This is known as the exploding and vanishing gradient problem. Previous solutions to this issue either built on rather complicated, purpose-engineered architectures with gated memory buffers, or more recently imposed constraints that ensure convergence to a fixed point or restrict (the eigenspectrum of) the recurrence matrix. Such constraints, however, convey severe limitations on the expressivity of the RNN. Essential intrinsic dynamics such as multistability or chaos are disabled. This is inherently at disaccord with the chaotic nature of many, if not most, time series encountered in nature and society. Here we offer a comprehensive theoretical treatment of this problem by relating the loss gradients during RNN training to the Lyapunov spectrum of RNN-generated orbits. We mathematically prove that RNNs producing stable equilibrium or cyclic behavior have bounded gradients, whereas the gradients of RNNs with chaotic dynamics always diverge. Based on these analyses and insights, we adapt the old idea of teacher forcing to yield an effective yet simple training technique for chaotic data, and offer guidance on how to choose relevant hyperparameters according to the Lyapunov spectrum.

text: Recurrent neural networks (RNNs) are wide-spread machine learning tools for modeling sequential and time series data. They are notoriously hard to train because their loss gradients backpropagated in time tend to saturate or diverge during training. This is known as the exploding and vanishing gradient problem. Previous solutions to this issue either built on rather complicated, purpose-engineered architectures with gated memory buffers, or - more recently - imposed constraints that ensure convergence to a fixed point or restrict (the eigenspectrum of) the recurrence matrix. Such constraints, however, convey severe limitations on the expressivity of the RNN. Essential intrinsic dynamics such as multistability or chaos are disabled. This is inherently at disaccord with the chaotic nature of many, if not most, time series encountered in nature and society. Here we offer a comprehensive theoretical treatment of this problem by relating the loss gradients during RNN training to the Lyapunov spectrum of RNN-generated orbits. We mathematically prove that RNNs producing stable equilibrium or cyclic behavior have bounded gradients, whereas the gradients of RNNs with chaotic dynamics always diverge. Based on these analyses and insights, we adapt the old idea of teacher forcing to yield an effective yet simple training technique for chaotic data, and offer guidance on how to choose relevant hyperparameters according to the Lyapunov spectrum.

heading: 1 INTRODUCTION
text: Recurrent neural networks (RNNs) are widely used across various fields in engineering and science for learning sequential tasks or modeling and predicting time series (Lipton et al., 2015). Yet, they struggle when long-term temporal dependencies, very slow, or hugely varying time scales are involved (Hochreiter, 1991; Bengio et al., 1994; Schmidt et al., 2021; Li et al., 2018; Rusch & Mishra, 2021a). Time series or sequential data with such properties are, however, very common in fields like climate physics (Thomson, 1990), neuroscience (Fusi et al., 2007; Russo & Durstewitz, 2017), ecology (Turchin & Taylor, 1992), or language processing (Cho et al., 2014b). Training RNNs on such data is hard because the loss gradients backpropagated in time easily saturate or diverge in this process. This is commonly referred to as the exploding and vanishing gradient problem (EVGP) (Hochreiter, 1991; Bengio et al., 1994; Pascanu et al., 2013).
One solution to the EVGP is based on specifically designed RNN architectures with gating mechanisms, such as long short-term memory (LSTM) (Hochreiter & Schmidhuber, 1997) or gated recurrent units (GRU) (Cho et al., 2014a). These architectures allow states at earlier time steps to more easily influence activity much later through a kind of protected memory buffer, thus alleviating the EVGP by structural design. In practice, such models need to be backed up by further techniques like gradient clipping to keep the gradients in check (Pascanu et al., 2013). The relatively complex architectural design of these networks impedes their mathematical analysis and requires reverse engineering after training (Maheswaranathan et al., 2019; Monfared & Durstewitz, 2020a;b; Schmidt et al., 2021). Partly to forego these complications, a variety of other solutions has been proposed recently, imposing restrictions on the recurrence matrix to bound the gradients (Arjovsky et al., 2016; Chang et al., 2019), or enforcing global stability by design or regularization (Erichson et al., 2021; Kolter & Manek, 2019). Often these procedures dramatically curtail the expressivity of the RNN (Kerg et al., 2019; Orhan & Pitkow, 2020; Schmidt et al., 2021); in particular, they rule out chaotic dynamics (see further discussion below).
This is at odds with the plethora of chaotic phenomena in nature, engineering, and society. Chaotic dynamics are commonplace, almost default in any complex physical or biological system. This includes scientific areas as diverse as neuroscience (Durstewitz & Gabriel, 2007; van Vreeswijk & Sompolinsky, 1996), physiology (Kesmia et al., 2020), geophysics (Sivakumar, 2004), climate systems (Tziperman et al., 1997), astrophysics (Laskar & Robutel, 1993), ecology (Duarte et al., 2010), chemical reactions (Field & Györgyi, 1993), cell (Olsen & Degn, 1977) or population (May, 1987) biology. Chaotic phenomena are also crucial for the understanding of societal and epidemiological processes, such as the spread of diseases (Mangiarotti et al., 2020), or in economics (Faggini, 2014). They are further relevant in purely technical contexts such as electrical engineering (Tchitnga et al., 2019; Kamdjeu Kengne et al., 2021) or laser optics (Kantz et al., 1993). They have even been suggested to play an up to now largely neglected, but potentially very significant role in speech recognition (Sabanal & Nakagawa, 1996) and natural language processing (Inoue et al., 2021). Hence, in almost any practical setting, chaotic phenomena abound. They cannot, in general, be ignored when devising RNN training algorithms.
Here we offer a comprehensive theoretical treatment of the relation between RNN dynamics and the behavior of the loss gradients during training. We find a close connection between an RNN’s loss gradients and the largest Lyapunov exponent of its freely generated orbits. We mathematically prove that RNNs producing stable fixed point or cyclic behavior have bounded gradients. Crucially, however, the loss gradients of RNNs producing chaotic dynamics always diverge. Hence, the chaotic nature of many time series data induces a principle problem, and, despite significant efforts in the past to solve the EVGP, training RNNs on such data remains an open issue. We illustrate the implications of our theory for RNN training on several simulated and empirical chaotic time series, and adapt the old idea of sparsely forced BPTT as a simple yet effective remedy that enables to learn the underlying dynamics despite exploding gradients.

heading: 2 RELATED WORKS
text: Exploding and vanishing gradients. While ‘classical’ remedies of the EVGP (Hochreiter, 1991; Bengio et al., 1994; Pascanu et al., 2013) rest on purpose-tailored architectures with gating mechanisms, which safeguard information flow across longer temporal distances (Hochreiter & Schmidhuber, 1997; Cho et al., 2014a), the focus has recently shifted to simpler RNNs that address the EVGP by restricting the recurrence matrix to be orthogonal (Henaff et al., 2016; Helfrich et al., 2018; Jing et al., 2019), unitary (Arjovsky et al., 2016), or antisymmetric (Chang et al., 2019), or by ensuring globally stable fixed point solutions (Kag et al., 2020; Kag & Saligrama, 2021a), for example through co-trained Lyapunov functions (Kolter & Manek, 2019). However, all these approaches impose strong limitations on the dynamical repertoire of the RNN, enforcing global convergence to fixed points or simple cycles.1 In doing so, they drastically reduce the expressiveness of these models (Kerg et al., 2019; Orhan & Pitkow, 2020). To address this problem, Erichson et al. (2021) somewhat relaxed the constraints on the recurrence matrix by introducing a skew-symmetric decomposition combined with a Lipschitz condition on the activation function. Another recent approach discretizes oscillator ODEs to arrive at a stable system of coupled (Rusch & Mishra, 2021a) or independent (Rusch & Mishra, 2021b) oscillators which increase the RNN’s expressiveness while bounding its gradients. By design (and as acknowledged by the authors), neither of these architectures is capable of producing chaotic dynamics, however, as the underlying ODEs do not allow for exponential divergence of close-by trajectories (a prerequisite for chaos). Given these often principle limitations of parametrically or dynamically strongly constrained models, a fruitful direction may be to modify the training process itself, e.g. through modified or auxiliary loss functions (Trinh et al., 2018; Schmidt et al., 2021), or special procedures for parameter updating (Kag & Saligrama, 2021b) or loss truncation (Williams & Zipser, 1989; Menick et al., 2021). Our empirical evaluation will follow up on such ideas, but also highlight that simple loss truncation, windowing, or architectural solutions like LSTMs are not sufficient.
Learning dynamical systems. Surprisingly disconnected from the work on the EVGP and learning long-term dependencies, a huge and long-standing literature deals with training RNNs on nonlinear dynamical systems (DS) (Pearlmutter, 1990; Trischler & D’Eleuterio, 2016; Vlachas et al., 2020), including chaotic systems like the famous Lorenz equations (Lorenz, 1963) or chaotic turbulence in fluid dynamics (Li et al., 2021). Teacher forcing (TF; Williams & Zipser (1989), Pearlmutter (1990),
1We make this point more formal in Appx. A.1.6.
Doya (1992), see also Goodfellow et al. (2016)) is one of the earliest techniques introduced to keep RNN trajectories on track while training. The idea behind TF is to simply replace RNN states by observations when available, thereby also effectively cutting off the gradients. TF essentially derives from ideas in dynamical control theory, and adaptive schemes that increasingly hand over control to the RNN throughout training have been devised (Abarbanel, 2013; Abarbanel et al., 2018; Bengio et al., 2015). A related technique from the control theory literature is “multiple shooting” (Voss et al., 2004): Here the whole observed time series is chopped into chunks, and for each chunk of trajectory a new initial condition is estimated. Explicit constraints ensure continuity between the separate trajectory bits during optimization. State space models and the Expectation-Maximization algorithm became popular particularly in the 90es for uncovering the latent dynamics underlying a set of time series observations (Ghahramani & Roweis, 1999), and remain an important tool until today (Durstewitz, 2017; Koppe et al., 2019). Most recently, approaches based on variational inference and the reparameterization trick (Kingma & Welling, 2014), like sequential variational autoencoders (SVAE), gained in popularity for DS approximation (Hernandez et al., 2020). “Deterministic” RNNs (i.e., with latent states not treated as random variables), like conventional LSTMs (Vlachas et al., 2018), remain top choices for DS reconstruction, however.
Although connections between DS ideas and loss gradients have been drawn early on (Bengio et al., 1994), so far only particular scenarios (like fixed point attractors) have been considered. Closest to our work is recent work by Schmidt et al. (2021), where non-divergence of loss gradients is established when RNNs converge to fixed points or cycles. However, this was done only for the particular class of piecewise-linear RNNs (PLRNNs), more restrictive conditions for cycles were imposed than assumed here, and - most importantly - the chaotic case on which we focus here was not considered. Another recent study (Engelken et al., 2020) also points out the general connections between Lyapunov exponents and loss gradients that we develop in sect. 3.1, but does not provide any in-depth theoretical treatment, proofs, or empirical evaluation, as we do here. Thus, a systematic theoretical framework that relates RNN dynamics more generally, and across a range of different RNN architectures, to the behavior of its training gradients, is still lacking so far.

heading: 3 THEORETICAL ANALYSIS: RELATION BETWEEN RNN DYNAMICS AND LOSS GRADIENTS
text: In our analysis, we will cover all major types of system dynamics (fixed points, cycles, chaos, and quasi-periodicity), and mathematically investigate their implications for the loss gradients. We will do this for all major classes of RNNs, including standard RNNs with largely arbitrary activation function, LSTMs, GRUs, and PLRNNs. The next section will first develop and illustrate the basic intuition behind the relations between RNN dynamics and loss gradients.

heading: 3.1 PRELIMINARIES: RNN DYNAMICS AND LOSS GRADIENTS
text: Formally, all popular RNN architectures, including LSTMs, GRUs, or PLRNNs, are discrete time DS, defined by a (first-order-Markovian) recursive prescription for the temporal evolution of the latent states zt ∈ RM of the form
zt = Fθ(zt−1, st), (1)
where st ∈ RN is the input at time t and θ are RNN parameters. For instance, for standard RNNs we have Fθ(zt−1, st) = f ( Wzt−1 + Bst + h ) , where f is an element-wise activation function like tanh or a rectified linear unit (ReLU).
Assuming we start at some initial value z1 ∈ RM , and given a sequence of external inputs S = {st}, we can recursively rewrite eq. (1) as
zT = Fθ(Fθ(Fθ(...Fθ(z1, s2)...))) =: F T−1 θ (z1, s2). (2)
In DS theory, we characterize the long-term behavior of such sequences by its spectrum of Lyapunov exponents. The Lyapunov exponents estimate the exponential growth rates in different local directions of the system’s state space, and the largest Lyapunov exponent gives the dominant exponential behavior. Let us denote the system’s Jacobian at time t by
Jt := ∂Fθ(zt−1, st)
∂zt−1 = ∂zt ∂zt−1 . (3)
For instance, for standard RNNs we would have Jt = W diag ( f ′ ( Wzt−1 + Bst + h )) , where diag denotes a diagonal matrix for which the i-th diagonal entry is the derivative of f w.r.t. zi,t−1. Then, the maximal Lyapunov exponent along an RNN trajectory {z1, z2, · · · , zT , · · · } is defined as
λmax := lim T→∞
1 T log ∥∥∥∥∥ T−2∏ r=0 JT−r ∥∥∥∥∥ , (4) where ∥ · ∥ denotes the spectral norm (or any subordinate norm) of a matrix. If λmax < 0 nearby trajectories will ultimately converge to a fixed point or cycle, while for λmax > 0 (a necessary condition for chaos) initially nearby trajectories will exponentially separate, i.e. we will have divergence along one (or more) directions in state space. This accounts for the sensitive dependence on initial conditions in chaotic systems. Now let L(W ,B,h) be some loss function employed for RNN training that decomposes in time as L = ∑Tt=1 Lt. Suppose we fancy BPTT as our training algorithm (similar derivations could be performed for RTRL), we recursively develop the loss gradients w.r.t. some RNN parameter θ in time (i.e., across layers of the RNN unrolled in time) as
∂L ∂θ
= T∑
t=1
∂Lt ∂θ with ∂Lt ∂θ
= t∑
r=1
∂Lt ∂zt ∂zt ∂zr ∂+zr ∂θ , (5)
and ∂zt ∂zr = ∂zt ∂zt−1 ∂zt−1 ∂zt−2 · · · ∂zr+1 ∂zr = t−r−1∏ k=0 ∂zt−k ∂zt−k−1 = t−r−1∏ k=0 Jt−k, (6)
where ∂+ denotes the immediate derivative. Now observe that the behavior of the loss gradients crucially depends on the product series of Jacobians in eqn. (6) : If the maximum absolute eigenval-
ues of the Jacobians Jt will, in the geometric mean, be larger than 1 (i.e., ∥∥∥∏T−2r=0 JT−r∥∥∥1/T > 1),
gradients will explode as T → ∞, while they will saturate if ∥∥∥∏T−2r=0 JT−r∥∥∥1/T < 1. Thus, the key point to note is that the same terms that occur in the definition of the Lyapunov spectrum, eqn. (4), resurface in the loss gradients, eqn. (5) & (6). This accounts for the tight links between system dynamics and gradients.

heading: 3.2 FIXED POINTS AND CYCLIC DYNAMICS
text: Let us start by considering the simplest types of dynamics that can occur in RNNs (or any discretetime DS): fixed points and cycles. In fact, by far most of the literature on global stability in RNNs and on loss gradients focused on just fixed points (Chang et al., 2019; Kolter & Manek, 2019; Erichson et al., 2021), with only few authors who recently started to also connect cyclic behavior to loss gradients (Schmidt et al., 2021; Rusch & Mishra, 2021a). Recall that a fixed point of a recursive map zt = F (zt−1) is defined as a point z∗ for which we have z∗ = F (z∗).2 Likewise, a k-cycle (k > 1) is a set of temporally consecutive periodic points Pk := {zt1 , zt2 , . . . ,ztk} = {zt1 , F (zt1), . . . , F k−1(zt1)} that we obtain from recursive application of the map such that each of the cyclic points ztr ∈ Pk is a fixed point of the k times iterated map F k (with k being the smallest positive integer for which this holds). To simplify the subsequent treatment, we will collectively refer to fixed points and cycles as k-cycles (k ≥ 1). Further recall that a fixed point or k-cycle is called stable if the maximum absolute eigenvalue of the Jacobian evaluated at that point is smaller than 1, neutrally stable if exactly 1, and unstable otherwise. Although the results we develop in this and the following sections will hold more widely, we will restrict our attention to recursive maps Fθ from the class of RNNs R = {standardRNN, LSTM, GRU, PLRNN} (see Appx. A.1 for details). Based on the observations made in the previous sections we can state the following theorem that links RNN dynamics and loss gradients: Theorem 1. Consider an RNN Fθ ∈ R parameterized by θ, and assume that it converges to a stable fixed point or k-cycle Γk (k ≥ 1) with BΓk as its basin of attraction3. Then for every z1 ∈ BΓk (i)
2From here on we will suppress the explicit dependence on external inputs st notation-wise, see Remark 2. 3The basin of attraction is defined as the set of all points from which orbits converge to the resp. attractor.
the Jacobian ∂zT∂z1 exponentially vanishes as T → ∞; (ii) for Γk the tangent vectors ∂zT ∂θ and thus the gradient of the loss function, ∂LT∂θ , will be bounded from above, i.e. will not diverge for T → ∞; and (iii) for the PLRNN (27) both ∥∥∂zT ∂θ
∥∥ and ∥∥∂LT∂θ ∥∥ will remain bounded for every z1 ∈ BΓk as T → ∞. Proof. (i) Assume that Γk is a stable k-cycle (k ≥ 1) denoted by Γk = {z1, z2, · · · , zT , · · · } = {zt∗k , zt∗k−1, · · · , zt∗k−(k−1), zt∗k , zt∗k−1, · · · , zt∗k−(k−1), · · · }.
(7)
Then, the largest Lyapunov exponent of Γk is given by
λΓk = limt→∞
1 t ln ∥∥J∗t J∗t−1 · · · J∗2∥∥ = lim j→∞ 1 jk ln ∥∥∥∥∥ ( k−1∏
s=0
Jt∗k−s )j∥∥∥∥∥ . (8) By assumption of stability of Γk we have λΓk < 0 and also ρ (∏k−1 s=0 Jt∗k−s ) < 1, which implies
lim t→∞
J∗t J ∗ t−1 · · · J∗2 = lim
j→∞ ( k−1∏ s=0 Jt∗k−s )j = 0. (9)
Now suppose that Oz1 is an orbit of (1) converging to Γk, i.e. z1 ∈ BΓk . Since Oz1 and Γk have the same largest Lyapunov exponent, we have
λOz1 = limT→∞
1 T ln ∥JT JT−1 · · · J2∥ = λΓk < 0, (10)
and hence for z1 ∈ BΓk
lim T→∞ ∥∥∥∥∂zT∂z1 ∥∥∥∥ = limT→∞ ∥JT JT−1 · · · J2∥ = 0. (11)
(ii)& (iii) See Appx. A.2.1.
Remark 1. The result of Theorem 1 part (i) will be generally true for any first-order-Markovian recursive map (1), but the conclusions in part (ii) may hinge on its specific definition. Remark 2. None of the results above and throughout sect. 3 require the dynamics to be autonomous, the theory applies whether there is external input or not. In fact, mathematically, non-autonomous (externally forced) systems can always be rewritten as autonomous dynamical systems (Alligood et al., 1996; Perko, 2001; Zhang et al., 2009), see Appx. A.1.1 for details.
The results above ensure that loss gradients will not diverge (explode) as T → ∞ in RNNs that are “well-behaved” in the sense that they converge to a fixed point or cycle. This is a generalization of the results given in Theorem 1 in Schmidt et al. (2021), where this was shown only a) for the specific class of PLRNNs and b) for specific constraints imposed on the eigenvalue spectrum of the RNN’s Jacobians which were relaxed in our theorem above.
While our treatment above is centered on the “exploding-gradients” case, various architectural modifications or regularization techniques can ensure that gradients do not vanish either, i.e. remain bounded from below as well. This was established, for instance, in Schmidt et al. (2021) for PLRNNs using ’manifold attractor regularization’. In Appx. A.2.1 we show that the results from Theorem 2 from Schmidt et al. (2021) on doubly bounded (from below and above) loss gradients can indeed be extended to the more general case covered by Theorem 1 above.

heading: 3.3 CHAOTIC DYNAMICS
text: We will now consider the all-important chaotic case. Let F be a recursive map and Oz1 = {z1, z2, z3, · · · } be an orbit of F . The orbit is chaotic if (i) it is not asymptotically periodic and (ii) has at least one positive Lyapunov exponent (Glendinning & Simpson, 2021; Meiss, 2007). If the system’s invariant set is bounded, condition (ii) is considered a standard signature of chaos, as in this case two nearby orbits separate exponentially fast, but at the same time their mutual separation cannot go to infinity so that there are also folds. The following theorem states the sufficient condition for exploding gradients:
Theorem 2. Suppose that an RNN Fθ ∈ R (parameterized by θ) has a chaotic attractor Γ∗ with BΓ∗ as its basin of attraction. Then, for every orbit with z1 ∈ BΓ∗ , (i) the Jacobians connecting temporally distal states zT and zt ( T ≫ t), ∂zT∂zt , will exponentially explode for T → ∞, and (ii) the tangent vector ∂zT∂θ and so the gradients of the loss function, ∂LT ∂θ , will diverge as T → ∞.
Proof. Let the RNN Fθ ∈ R have a chaotic orbit denoted by Γ∗ = {z∗1 , z∗2 , · · · , z∗T , · · · }. Then, denoting by J∗T the Jacobian of (1) at z ∗ T ∈ Γ∗, the largest Lyapunov exponent of Γ∗ is given by
λ = lim T→∞
1 T ln ∥∥J∗T J∗T−1 · · · J∗2∥∥ . (12)
Since Γ∗ is chaotic, so λ > 0. Hence, from (12), it is concluded that
lim T→∞ ∥∥J∗T J∗T−1 · · · J∗2∥∥ = lim T→∞ ∥∥∥∥∂z∗T∂z∗t ∥∥∥∥ = ∞ , T ≫ t. (13)
Now, according to Oseledec’s multiplicative ergodic Theorem, nearly all the points in the basin of attraction of Γ∗ have the same largest Lyapunov exponent λ. Thus, (13) holds for every z1 ∈ BΓ∗ . (ii) See Appx. A.2.2.
Remark 3. The first part of Theorem 2 holds for all first-order-Markovian recursive maps (1). Note that for LSTMs, ∂zT∂zt (z := (h, c) T) denotes the full Jacobian of both hidden and cell states.
We collect some further mathematical results and remarks related to Theorem 2 in Appx. A.3.1.
Hence, the essential result is that for all popular RNNs R and activation functions, loss gradients will inevitably diverge if the RNN latent states converge to a chaotic attractor.

heading: 3.4 QUASI-PERIODICITY
text: Quasi-periodicity is a long-term behavior which occurs on a torus and, superficially, bears some similarity to chaos in the sense that, strictly speaking, orbits are also aperiodic. That is, as T → ∞, trajectories will never close up with themselves. Moreover, every trajectory becomes arbitrarily close to any point on the torus, that is, it is dense. One important difference between quasi-periodic and chaotic systems is, however, that in a quasi-periodic system, as time passes, two close initial conditions are linearly diverging, while in a chaotic system the divergence is exponential. Theorem 3. Assume that an RNN Fθ ∈ R (parameterized by θ) has a quasi-periodic attractor Γ with BΓ as its basin of attraction. Then, for every z1 ∈ BΓ
∀ 0 < ϵ < 1 ∃T0 > 1 s.t. ∀T ≥ T0 =⇒ (1− ϵ)T−1 < ∥∥∥∥∂zT∂z1 ∥∥∥∥ < (1 + ϵ)T−1. (14) Proof. See Appx. A.2.3.
According to Theorem 3, for every orbit converging to a quasi-periodic attractor, the Jacobians ∂zT∂zt may diverge or vanish as T → ∞, but this will not occur exponentially fast as T → ∞. Thus, even for bounded non-chaotic RNNs we may sometimes stumble into the problem of diverging gradients. Although this may be a less common scenario, we point out it may occur if we train RNNs on real data from oscillatory systems with incommensurate frequencies, as for instance encountered in electronic engineering.
In Appx. A.3.2 we have collected further mathematical results on the connection between RNN dynamics and loss gradients that hold regardless of the RNN’s limiting behavior.

heading: 4 EMPIRICAL EVALUATION
text: Our theoretical results imply that chaotic time series pose a principle challenge for RNN training that cannot easily be circumvented through specifically designed architectures, constraints, or regularization criteria. If the underlying DS we aim to capture is chaotic, loss gradients propagated back in time will inevitably explode. Here we will work out some implications for RNN training and potential remedies empirically.

heading: 4.1 TRAINING ON SYSTEMS WITH EXPLODING GRADIENTS BY SPARSE TEACHER FORCING
text: To illustrate the connections between theory and RNN training, we revive the old idea of TF (Williams & Zipser, 1989) as a mechanism for truncating error gradients while training. However, we would like to do this such that important information about the system dynamics does not get lost, for which Lyapunov theory offers some guidance. Specifically, we should not force the system back onto the true trajectory all or most of the time (as in “classical TF”), but should effectively “re-calibrate” it only at certain time points chosen wisely according to the system’s local divergence rates. This procedure will be referred to as sparsely forced BPTT in the following.
Assume we want to train an RNN with hidden states zt ∈ RM and linear (or affine) output layer on a time-series {x1,x2, · · · ,xT } generated by a chaotic system.4 The linear output layer x̂t = Bzt, B ∈ RN×M , maps the RNN hidden states into the observation space. This allows us to modify the original TF procedure by constructing a control series {z̃1, z̃2, · · · , z̃T } from the observations by “inverting” the linear output mapping 5
z̃t = (B TB)−1BTxt. (15)
The idea is to supply this control signal only sparsely, separated by the learning interval τ between consecutive forcings. Hence, defining T = {nτ + 1}n∈N0 as the set of all time points at which we force the RNN onto the ‘true’ values, the RNN updates can be written as
zt+1 = { RNN(z̃t) if t ∈ T RNN(zt) else . (16)
This forcing is applied after calculation of the loss, such that Lt = ∥xt −Bzt∥22 irrespective of whether t is in T or not (and of course it is applied only during training, not at test time!). Replacing hidden states zt with their teacher-forced signals z̃t simply breaks divergence between true and predicted trajectories at time points t ∈ T , and also cuts off the Jacobians by breaking the temporal contingency (for details see Appx. A.7). The learning interval τ hence controls how many time steps are included in the gradient calculation and has to be chosen with care such as to balance the effects of exploding gradients vs. those of loosing relevant time scales and long-term dependencies. While it is general wisdom that an optimal batch size will facilitate training, the point here is thus much more specific: Ideally τ should be chosen in accordance with the system’s Lyapunov spectrum, for instance based on the predictability time (Bezruchko & Smirnov, 2010)
τpred = ln 2
λmax . (17)
We emphasize that such a simple recipe for addressing the exploding gradient problem is based on modifying the training routine, and is thus in principle applicable to any model architecture.6

heading: 4.2 EXAMPLE 1: LORENZ SYSTEM AND EXTERNALLY FORCED DUFFING OSCILLATOR IN CHAOTIC REGIME
text: Let us illustrate these ideas on two classical textbook examples of chaotic DS, the chaotic Lorenz attractor as an autonomous system, and the chaotically forced Duffing oscillator as an example with explicit external input (see Appx. A.4 for details). Trajectories were repeatedly drawn from these systems, on which we trained a PLRNN, a vanilla RNN with tanh activation function, and a LSTM by stochastic gradient descent (SGD) to minimize the MSE loss between predicted and actual observations. As optimizer we used Adam (Kingma & Ba, 2015) from PyTorch (Paszke et al., 2017) with a learning rate of 0.001. For all models, training proceeded solely by sparsely forced BPTT and did not employ gradient clipping or any other technique that may interfere with optimal loss truncation.
In nonlinear DS reconstruction, we are mainly interested in reproducing invariant properties of the underlying system such as the attractor geometry (or topology; Takens (1981); Sauer et al. (1991)) or
4Note that in DS reconstruction one usually considers the data as observations (unsupervised problem); according to Remark 2, however, it is mainly a matter of the scientific question addressed whether we include certain variables as explicit inputs or as observations.
5To ensure invertibility, one could add a regularizer λI to BTB in eqn. (15), as in ridge regression, but we did not find this necessary in any of our examples.
6All code produced here is available at [placeholder].
the frequency composition (i.e., time-averaged properties), while measures like ahead-prediction errors are less meaningful especially on chaotic time series (Wood, 2010; Koppe et al., 2019). Thus, in evaluating training performance, here we follow Koppe et al. (2019) in using a Kullback-Leibler divergence Dstsp to quantify the agreement between observed and generated probability distributions across state-space to asses the overlap in attractor geometry (Appx. A.5). Moreover, we employ a dimension-wise frequency correlation measure (PSC) to quantify the agreement of power-spectra of the observed and generated time-series (Appx. A.5).
Fig. 1 shows the dependence of the reconstruction quality on the learning interval τ for all RNN architectures on (a) the Lorenz and (b) the externally forced Duffing system. Fig. 2 provides particular examples of reconstructions for τ chosen too small, too large, or about right. For all models we find a system-dependent range for the optimal learning interval that agrees well with the predictability time defined in eqn. (17), where estimates for the maximal Lyapunov exponent were taken from the literature (Rosenstein et al., 1993; Gilpin, 2021). As a reference, dashed lines represent the reconstruction performance for all architectures when trained with classical BPTT and gradient clipping. The training procedure was the same as for sparsely forced BPTT, except that instead of supplying a control-signal, gradients were normalized to 1 prior to each parameter update. As evidenced by the much worse performance, gradient clipping does not effectively address the EVGP, even for LSTMs. As further shown in Fig. 9 in Appx. A.6.4, using the optimal window length τ but resetting the initial condition to zero (instead of its control value z̃t) for each chunk equally destroys performance. This suggests that neither mere gradient normalization nor simple windowing are sufficient, but will wipe out essential information about the dynamics.
In Appx. A.6 we collect further results on the chaotic Rössler attractor (Fig. 5), high-dimensional Mackey-Glass equations (Fig. 7), and the Lorenz attractor with partial observations (Fig. 8).

heading: 4.3 EXAMPLE 2: CHAOTIC WEATHER DATA
text: As for an empirical example, we trained all RNNs (vanilla RNN, PLRNN, LSTM) on a temperature time series recorded at the Weather Station at the Max Planck Institute for Biogeochemistry in Jena, Germany. To expose the chaotic behavior and obtain a robust estimate of the maximal Lyapunov exponent, trends and yearly cycles were removed, and nonlinear noise-reduction was performed (Kantz et al. (1993); Appx. A.4). The maximal Lyapunov exponent was determined with the TISEAN package (Hegger et al., 1999), as shown in Figure 3 (a). The value obtained is in close agreement with the literature (Millán et al., 2010).
Figure 3 shows that also for these empirical data the optimal training interval τ agrees well with the predictability time, eqn. (17), for all trained RNNs. Furthermore, as was the case for the DS benchmarks, gradient clipping was not able to satisfactorily tackle the EVGP, even when paired with architectures like LSTMs explicitly designed for alleviating this problem. Similar results are reported for another real-world dataset, electroencephalogram (EEG) recordings, in Appx. A.11.

heading: 5 DISCUSSION AND CONCLUSIONS
text: In this paper we proved that RNN dynamics and loss gradients are intimately related for all major types of RNNs and activation functions. If the RNN is “well behaved” in the sense that its dynamics converges to a fixed point or cycle, loss gradients will remain bounded, and established remedies (Hochreiter & Schmidhuber, 1997; Schmidt et al., 2021) can be used to refrain them from vanishing. However, if the dynamics are chaotic, gradients will always explode. This constitutes a principle problem in RNN training that cannot easily be mastered through architectural design or gradient clipping. It is furthermore a practically highly relevant one, as most time series we encounter in nature, and many from man-made systems as well, are inherently chaotic. While we do not offer a full solution to this problem here, we suggest it might be tackled in training by taking a system’s local divergence rates as measured through the Lyapunov spectrum into account. Hence, rather than conquering the EVGP by structural design or specific constraints or regularization terms, we recommend to put the focus more on the training process itself. We illustrated this point empirically using sparsely forced BPTT, a training technique that pulls trajectories back on track at times determined by the maximal Lyapunov exponent. Doing so leads to optimal reconstruction results for a variety of simulated and real-world benchmarks, regardless of the specific RNN architecture employed in training. We stress that our goal above all was to provide a mathematically grounded perspective on the problem, with the empirical section focused on elucidating the practical implications of the theoretical results. Empirically, for instance, precise Lyapunov exponents may sometimes be hard to obtain, although our empirical examples confirm that estimates based on log-divergence plots may work sufficiently well.

heading: ACKNOWLEDGMENTS

heading: ETHICS STATEMENT
text: The current work deals with theoretical-mathematical aspects of RNN training and performs basic research on limitations of training algorithms. As such it has no direct ethical implications. Since many real world applications depend on accurate forward predictions of time series, however, this
paper may raise awareness for an important issue that will also be relevant in practical settings, including sensitive domains like medical time series, weather forecasts, or traffic control.

heading: REPRODUCIBILITY STATEMENT
text: All theoretical results in this paper were carefully and thoroughly proven, with all proofs and detailed derivations available in the Appendix. Likewise, we will make available all code used in the empirical section in a way that will allow others to easily reproduce the results from this paper. This means we will include everything, starting with the code for benchmark simulations and simulated time series data used for evaluation, code for our model and training algorithms, up to the metafiles that produce the actual figures in this work, on our lab github site. All of this will be clearly documented.

heading: A APPENDIX

heading: A.1 THEOREMS: PRELIMINARIES

heading: A.1.1 TRANSFORMING NON-AUTONOMOUS INTO AUTONOMOUS DISCRETE-TIME DS
text: Following (Zhang et al., 2009), and based on similar reasoning as for continuous time (ODE-based) DS (Alligood et al., 1996; Perko, 2001), let us consider the non-autonomous discrete-time DS
xt+1 = F (xt, t), x ∈ Rn. (18) Defining zt = (xt, t)T and G(zt) = (F (xt, t), t + 1)T , system (18) can be rewritten as the autonomous system
zt+1 = G(zt), z ∈ Rn+1. (19) Hence, in all our theoretical treatment we can confine our attention to systems of the form eqn. 18.

heading: A.1.2 RNN DERIVATIVES
text: Considering the loss function L = ∑Tt=1 Lt of an RNN Fθ ∈ R parameterized by θ, we have ∂L ∂θ = T∑ t=1 ∂Lt ∂θ , (20)
where ∂Lt ∂θ = ∂Lt ∂zt ∂zt ∂θ . (21)
The tangent vector ∂zT∂θ has the form
∂zT ∂θ = ∂+zT ∂θ + T−2∑ t=1 ( t−1∏ r=0 JT−r ) ∂+zT−t ∂θ , (22)
where ∂+ denotes the immediate partial derivative. Since for an RNN Fθ ∈ R the activation function is element-wise, with θ the m-th element of a parameter vector θ (or belonging to the m-th row of a parameter matrix θ), we have
∂+zT ∂θ
= ( 0 · · · 0 ∂
+zm,T ∂θ 0 · · · 0
)T . (23)
For instance, let θ = W be a weight matrix, then
∂L ∂W =  ∂L ∂w11 ∂L ∂w12 · · · ∂L∂w1M ∂L ∂w21 ∂L ∂w22 · · · ∂L∂w2M ...
∂L ∂wM1 ∂L ∂wM2 · · · ∂L∂wMM
 . (24)
In this case, for the standard RNN we have ∂+zT ∂wmk = (0 · · · 0 zk,T−1 ξmk(zT−1) 0 · · · 0)T = 1(m,k) ξmk(zT−1) zT−1, (25)
where ξmk(zT−1) = f ′wm,k (∑M j=1 wmj zj,T−1 + ∑M j=1 bmj sj,T + hm )
, and f ′wm,k stands for the derivative of f with respect to wm,k.
Therefore, for standard RNNs, (22) becomes
∂zT ∂wmk = 1(m,k) ξmk(zT−1) zT−1 + T−2∑ t=1 ( t−1∏ r=0 JT−r ) 1(m,k) ξmk(zT−t−1) zT−t−1. (26)

heading: A.1.3 PIECEWISE-LINEAR RNN (PLRNN)
text: The PLRNN has the generic form (Koppe et al., 2019; Schmidt et al., 2021)
zt = F (zt−1) = Azt−1 + Wϕ(zt−1) + Cst + h+ εt, (27)
where ϕ(zt−1) = max(zt−1, 0) is the element-wise rectified linear unit (ReLU) function, zt ∈ RM is the neural state vector, A ∈ RM×M is a diagonal matrix of auto-regression weights, W ∈ RM×M is a matrix of connection weights, h ∈ RM is the bias vector, st ∈ RK the external input weighted by C ∈ RM×K , and εt ∼ N (0,Σ) a Gaussian noise term with diagonal covariance matrix Σ. Equation (27) can be rewritten as
zt = (A+WDΩ(t−1))zt−1 + Cst + h+ εt =: WΩ(t−1) zt−1 + Cst + h+ εt, (28)
where DΩ(t) := diag(dΩ(t)) is a diagonal matrix with dΩ(t) := (d1, d2, · · · , dM ) an indicator vector such that dm(zm,t) =: dm = 1 whenever zm,t > 0, and zeros otherwise.
For the PLRNN (28) we have
Jt = ∂zt
∂zt−1 = WΩ(t−1), (29)
and ∥∥WΩ(t−1)∥∥ ≤ ∥A∥ + ∥W ∥.
Furthermore, the derivatives (22) for the PLRNN (28) are
∂zT ∂wmk = 1(m,k)DΩ(T−1) zT−1 + T−1∑ j=2 ( j−1∏ i=1 WΩ(T−i) ) 1(m,k)DΩ(T−j) zT−j . (30)

heading: A.1.4 LONG SHORT-TERM MEMORY (LSTM)
text: The LSTM is defined by the equations it = σ ( Wiist + Whiht−1 + bi ) ft = σ ( Wifst + Whfht−1 + bf
) gt = tanh ( Wigst + Whght−1 + bg
) ot = σ ( Wiost + Whoht−1 + bo
) ct = ft ⊙ ct−1 + it ⊙ gt ht = ot ⊙ tanh (ct) (31)
where {st} is the input sequence, W denotes weight matrices, b bias terms, it,ft, gt,ot demonstrate the input, forget, cell, and output gates, ht and ct are the hidden and cell states at time t respectively, σ is the sigmoid activation function, and ⊙ represents the element-wise (Hadamard) product (see (Hochreiter & Schmidhuber, 1997; Graves et al., 2016; Vlachas et al., 2018) for further information on LSTMs).
Defining zt := (ht, ct)T, the LSTM (31) can be represented as the first-order recursive map
zt = Fθ(zt−1) =
( ot ⊙ tanh (ft ⊙ ct−1 + it ⊙ gt)
ft ⊙ ct−1 + it ⊙ gt
) . (32)
The term ∂Lt∂θ in (20) for some LSTM parameter θ can be written as
∂Lt ∂θ = t∑ r=1 ∂Lt ∂ht ∂ht ∂zt ∂Lt ∂zt ∂zt ∂zr ∂zr ∂θ . (33)
A necessary condition for LSTMs to have a chaotic orbit is given by:
Proposition 1. Let the LSTM given by (31) have a chaotic attractor Γ∗ with BΓ∗ as its basin of attraction. Then for every z1 = (h1, c1)T ∈ BΓ∗
γ := lim T→∞
T √√√√√ ∥∥∥∥∥∥ ∂hT∂h1 ∂hT∂c1
∂cT ∂h1 ∂cT ∂c1 ∥∥∥∥∥∥ > 1. (34) Proof. The Jacobian matrix of (32) for t > 1 can be written in the block form
∂zt ∂zt−1 = Jt =  ∂ht∂ht−1 ∂ht∂ct−1 ∂ct
∂ht−1 ∂ct ∂ct−1  . (35) Further, due to the chain rule, we have
Jt Jt−1 =  ∂ht∂ht−1 ∂ht−1∂ht−2 + ∂ht∂ct−1 ∂ct−1∂ht−2 ∂ht∂ht−1 ∂ht−1∂ct−2 + ∂ht∂ct−1 ∂ct−1∂ct−2 ∂ct
∂ht−1 ∂ht−1 ∂ht−2 + ∂ct∂ct−1 ∂ct−1 ∂ht−2
∂ct ∂ht−1 ∂ht−1 ∂ct−2 + ∂ct∂ct−1 ∂ct−1 ∂ct−2

=  ∂ht∂ht−2 ∂ht∂ct−2 ∂ct
∂ht−2 ∂ct ∂ct−2  , (36) and by induction we obtain
∂zt ∂z1 = Jt Jt−1 Jt−2 · · · J2 =  ∂ht∂h1 ∂ht∂c1 ∂ct ∂h1 ∂ct ∂c1  . (37) Now assume that (32) has a chaotic orbit given by
Γ∗ = {z∗1 , z∗2 , · · · , z∗T , · · · }. (38)
According to (37), the largest Lyapunov exponent of Γ∗ is given by
λΓ∗ = lim T→∞
1 T ln ∥∥J∗T J∗T−1 · · · J∗2∥∥ = lim T→∞ 1 T ln ∥∥∥∥∥∥ ∂h∗T∂h∗1 ∂h∗T∂c∗1
∂c∗T ∂h∗1 ∂c∗T ∂c∗1 , ∥∥∥∥∥∥ . Since Γ∗ is chaotic, so λΓ∗ > 0, which gives
lim T→∞
T √√√√√ ∥∥∥∥∥∥ ∂h∗T∂h∗1 ∂h∗T∂c∗1
∂c∗T ∂h∗1 ∂c∗T ∂c∗1 ∥∥∥∥∥∥ > 1. (39) Based on Oseledec’s multiplicative ergodic Theorem, (39) holds for every z1 ∈ BΓ∗ . This completes the proof.

heading: A.1.5 GATED RECURRENT UNIT (GRU)
text: A GRU network is defined by the equations
zt = σ ( Wz st + Uzht−1 + bz ) rt = σ ( Wr st + Urht−1 + br
) ht = (1− zt) ⊙ tanh ( Wh st + Uh(rt ⊙ ht−1) + bh ) + zt ⊙ ht−1, (40)
where rt represents the reset gate, zt the update gate, st and ht denote the inputs and the hidden state respectively, Wz,Wr,Wh ∈ RM×N and Uz,Ur,Uh ∈ RM×M are weight matrices, bz, br, bh ∈ RM are bias vectors, and σ is the element-wise logistic sigmoid function (for more details about GRUs see Cho et al. (2014a)).

heading: A.1.6 UNITARY EVOLUTION RNN (URNN)
text: The uRNN, proposed in (Arjovsky et al., 2016), is defined as the nonlinear DS zt = σb ( Wzt−1 + V st ) , (41)
for which W ∈ U(M) is an unitary matrix, V ∈ CM×N , b ∈ RM is the bias parameter, st is the real- or complex-valued input of dimension N , and
[σb(z)]i = [σmodReLU(z)]i =
{( |zi|+ bi ) zi |zi| if |zi|+ bi ≥ 0
0 if |zi|+ bi < 0 . (42)
Proposition 2. The uRNN given by (41) cannot have any chaotic orbit.
Proof. For any arbitrary orbit Oz1 of (41) we have
∥JT JT−1 · · · J2∥ = ∥∥∥∥∥ T−2∏ k=0 DT−k W T ∥∥∥∥∥ , (43) where Dt = diag ( σ′b ( Wzt−1 + V st )) . Since W is unitary and so a norm preserving matrix, it is concluded that∥∥∥∥∥ T−2∏ k=0 DT−k W T ∥∥∥∥∥ ≤ T−2∏ k=0 ∥∥DT−k WT∥∥ = T−2∏ k=0 ∥DT−k∥ = 1, (44)
which implies
λmax = lim T→∞
1 T ln ∥JT JT−1 · · · J2∥ ≤ 0. (45)
This rules out the existence of chaos (since λmax > 0 is a necessary condition for Oz1 to be chaotic).
Note that, more generally, any RNN which is constrained such as to exhibit global convergence to a fixed point or cycle, by definition must have a maximum Lyapunov exponent λmax ≤ 0 (in accordance with Theorem 1), hence cannot exhibit chaotic behavior by definition.

heading: A.2 THEOREMS: PROOFS

heading: A.2.1 PROOF OF THEOREM 1, PARTS (II) & (III)
text: Proof. (ii) If J is the Jordan normal form of ∏k−1 s=0 Jt∗k−s , then ∏k−1 s=0 Jt∗k−s = P J P −1, where
J =  Jm1(λ1) 0 0 · · · 0 0 Jm2(λ2) 0 · · · 0 ... . . . . . . . . . ...
0 . . . 0 Jmp−1(λp−1) 0 0 . . . . . . 0 Jmp(λp)  , (46) and mi is the algebraic multiplicity of each eigenvalue λi. Since ρ (∏k−1 s=0 Jt∗k−s ) < 1, so the eigenvalue λi associated with each Jordan block satisfies |λi| < 1 (i = 1, · · · , p). Moreover, every mi ×mi Jordan block has the form
Jmi(λi) =  λi 1 0 · · · 0 0 λi 1 · · · 0 ... ... . . . . . .
... 0 0 . . . λi 1 0 0 . . . 0 λi  . (47) Accordingly ∥∥∥∥∥ ( k−1∏ s=0 Jt∗k−s
)j∥∥∥∥∥ = ∥∥P J j P−1∥∥ ≤ p ∥∥J j∥∥ , (48) in which p = ∥P ∥
∥∥P−1∥∥. Furthermore, for j ∈ N, J j is a block diagonal matrix of the form J j =  J jm1(λ1) 0 0 · · · 0 0 J jm2(λ2) 0 · · · 0 ... . . . . . . . . . ...
0 . . . 0 J jmp−1(λp−1) 0 0 . . . . . . 0 J jmp(λp)
 , (49)
in which every mi ×mi Jordan block has the form
J jmi(λi) =  λji ( j 1 ) λj−1i
( j 2 ) λj−2i · · ·
( j mi−1 ) λj−mi+1i
0 λji ( j 1 ) λj−1i · · ·
( j mi−2 ) λj−mi+2i
... ...
. . . . . . ...
0 0 . . . λji ( j 1 ) λj−1i
0 0 . . . 0 λji
 . (50)
In addition, for every block J jmi(λi), we have∥∥J jmi(λi)∥∥ ≤ √mi ∥∥J jmi(λi)∥∥∞ = √mi mi∑ q=1 ∣∣∣(J jmi(λi))1q∣∣∣ = √ mi
mi∑ q=1 ( j q − 1 ) |λi|j−q+1 = |λi|j √ mi ( |λi|1−mi mi∑ q=1 ( j q − 1 ) |λi|mi−q )
≤ |λi|j jmi √ mi ( |λi|1−mi mi∑ q=1 |λi|mi−q ) =: |λi|j jmi Nλi . (51)
Moreover, for any 1 < r̃i < 1|λi| , there exists some li such that j mi < r̃ ji for j ≥ li. This means for j ≥ li ∥∥J jmi(λi)∥∥ ≤ Nλi |r̃i λi|j , (52) such that |r̃i λi| = r̃i|λi| < 1.
Besides, for J j = J jm1(λ1)⊕ J jm2(λ2)⊕ · · · ⊕ J jmp(λp)∥∥J j∥∥ = max 1≤i≤p ∥∥J jmi(λi)∥∥ =: ∥∥J jm(λ)∥∥ . (53) Hence, from (48), (52) and (53), it is deduced that for j ≥ l∥∥∥∥∥( k−1∏ s=0 Jt∗k−s
)j∥∥∥∥∥ ≤ pNλ |r̃ λ|j =: p̄ rj , (54) in which r = |r̃ λ| < 1.
Furthermore, let for Γk
max T≥1
{ ∥J∗T ∥ } = max
0≤s≤k−1
{ ∥Jt∗k−s∥ } = m̄,
max T≥1 {∥∥∥∥∂+zT∂θ ∥∥∥∥} = max0≤s≤k−1{ ∥∥∥∥∂+zt∗k−s∂θ ∥∥∥∥} = ξ,
max T≥1
{ ∥zT ∥ } = max
0≤s≤k−1
{ ∥zt∗k−s∥ } = q̄. (55)
Hence, defining z0 = 0, for this k-cycle∥∥∥∥∂zT∂θ ∥∥∥∥ = ∥∥∥∥∥∂+zT∂θ + T−2∑ t=1 ( t−1∏ r=0 J∗T−r ) ∂+zT−t ∂θ ∥∥∥∥∥ =
∥∥∥∥∥∂+zT∂θ + T−1∑ t=1 ( t−1∏ r=0 J∗T−r ) ∂+zT−t ∂θ ∥∥∥∥∥ ≤ q̄ ξ ( 1 +
T−1∑ t=1 ∥∥∥∥∥ t−1∏ r=0 J∗T−r ∥∥∥∥∥ ) . (56)
On the other hand, for T = kj , from (54) and (55) we have
T−1∑ t=1 ∥∥∥∥∥ t−1∏ r=0 J∗T−r ∥∥∥∥∥ = kj−1∑ t=1 ∥∥∥∥∥ t−1∏ r=0 J∗kj−r ∥∥∥∥∥ = k−1∑ t=1 ∥∥∥∥∥ t−1∏ r=0 J∗kj−r ∥∥∥∥∥ + 2k−1∑ t=k ∥∥∥∥∥ t−1∏ r=0 J∗kj−r ∥∥∥∥∥ +
3k−1∑ t=2k ∥∥∥∥∥ t−1∏ r=0 J∗kj−r ∥∥∥∥∥+ · · ·+ (j−1)k−1∑ t=(j−2)k ∥∥∥∥∥ t−1∏ r=0 J∗kj−r ∥∥∥∥∥+ kj−1∑ t=(j−1)k ∥∥∥∥∥ t−1∏ r=0 J∗kj−r ∥∥∥∥∥ =
k−1∑ t=1 ∥∥∥∥∥ t−1∏ r=0 J∗kj−r ∥∥∥∥∥ + j∑ i=2 ik−1∑ t=(i−1)k ∥∥∥∥∥ t−1∏ r=0 J∗kj−r ∥∥∥∥∥ ≤ ( m̄+ m̄2 + · · ·+ m̄k−1 ) +
j∑ i=2 p̄ ( 1 + m̄+ m̄2 + · · ·+ m̄k−1 ) ri−1.
(57)
Thus, considering ( m̄+ m̄2 + · · ·+ m̄k−1 ) = M, it is deduced that
lim T→∞ ∥∥∥∥∂zT∂θ ∥∥∥∥ = limj→∞ ∥∥∥∥∂zkj∂θ ∥∥∥∥ ≤ q̄ ξ(1 + M+ p̄ r(1 +M)1− r ) = M̄ < ∞, (58)
which, by (21), implies ∂LT∂θ will be bounded for T → ∞.
(iii) Consider the PLRNN given by (27), where for simplicity we ignore the external inputs and noise terms. Let {zt1 , zt2 , zt3 , . . .} be an orbit which converges to Γk. Hence
lim n→∞ d(ztn ,Γk) = 0, (59)
which implies there exists a neighborhood U of Γk and k sub-sequences {ztkm}∞m=1, {ztkm+1}∞m=1, · · · , {ztkm+(k−1)}∞m=1 of the sequence {ztn}∞n=1 such that all these sub-sequences belong to U and
a) ztkm+s = F k(ztk(m−1)+s), s = 0, 1, 2, · · · , k − 1,
b) lim m→∞ ztkm+s = zt∗k−s, s = 0, 1, 2, · · · , k − 1,
c) for every ztn ∈ U there is some s ∈ {0, 1, 2, · · · , k − 1} such that ztn ∈ {ztkm+s}∞m=1.
In this case, for every ztn ∈ U with ztn ∈ {ztkm+s}∞m=1 , there exists some ñ ∈ N such that ztn = ztkñ+s and lim
ñ→∞ ztkñ+s = zt∗k−s . Therefore, continuity of F results in
lim ñ→∞ F (ztkñ+s) = F (zt∗k−s), (60)
and so by (28)
lim ñ→∞
( WΩ(tkñ+s) ztkñ+s + h ) = WΩ(t∗k−s) zt∗k−s + h, (61)
which implies
lim ñ→∞ WΩ(tkñ+s) ztkñ+s = WΩ(t∗k−s) zt∗k−s. (62)
Assuming lim ñ→∞ WΩ(tkñ+s) = L, since (62) holds for every zt∗k−s, substituting zt∗k−s = e T 1 = (1, 0, · · · , 0)T in (62), we can prove that the first column of L equals the first column of WΩ(t∗k−s). Performing the same procedure for zt∗k−s = eTi , i = 2, 3, · · · ,M , yields
lim ñ→∞ WΩ(tkñ+s) = WΩ(t∗k−s). (63)
According to (59), U contains an infinite number of terms of the sequence {ztn}∞n=1, i.e.
∃N ∈ N s.t. n ≥ N =⇒ ztn ∈ U. (64)
Suppose that ztn ∈ U for some n ≥ N . Thus, there exists some s ∈ {0, 1, 2, · · · , k − 1} such that ztn ∈ {ztkm+s}∞m=1. Without loss of generality let s = 0. Hence, there is some ñ ∈ N such that
ztn = ztkñ and lim ñ→∞ ztkñ = zt∗k . In this case, moving forward in time gives ztn = ztkñ ( ztn ∈ {ztkm}∞m=1 ) , lim
ñ→∞ ztkñ = zt∗k ,
ztn+1 = ztkñ+1 ( ztn+1 ∈ {ztkm+1}∞m=1 ) , lim
ñ→∞ ztkñ+1 = zt∗k−1,
ztn+2 = ztkñ+2 ( ztn+2 ∈ {ztkm+2}∞m=1 ) , lim
ñ→∞ ztkñ+2 = zt∗k−2,
...
ztn+k−1 = ztkñ+k−1 ( ztn+(k−1) ∈ {ztkm+k−1}∞m=1 ) , lim
ñ→∞ ztkñ+k−1 = zt∗k−(k−1),
ztn+k = ztk(ñ+1) ( ztn+k ∈ {ztkm}∞m=1 ) , lim
ñ→∞ ztk(ñ+1) = zt∗k ,
ztn+k+1 = ztk(ñ+1)+1 ( ztn+k+1 ∈ {ztkm+1}∞m=1 ) , lim
ñ→∞ ztk(ñ+1)+1 = zt∗k−1,
...
ztn+2k−1 = ztk(ñ+1)+k−1 ( ztn+2k−1 ∈ {ztkm+k−1}∞m=1 ) , lim
ñ→∞ ztk(ñ+1)+k−1 = zt∗k−(k−1),
ztn+2k = ztk(ñ+2) ( ztn+2k ∈ {ztkm}∞m=1 ) , lim
ñ→∞ ztk(ñ+2) = zt∗k ,
... (65)
Consequently, for n ≥ N and j ∈ N, we can write kj−1∏ i=0 WΩ(tn+kj−1−i)
= ( k∏
i=1
WΩ(tk(ñ+j)+k−i) )( k∏ i=1 WΩ(tk(ñ+j−1)+k−i) ) · · · ( k∏ i=1 WΩ(tk(ñ)+k−i) )
= j∏ l=0 k∏ i=1 WΩ(tk(ñ+j−l)+k−i). (66)
On the other hand, in equation (28), there are different configurations for matrix DΩ(t−1) and hence different forms for matrix WΩ(tkñ+s) . In this case, the phase space of the system is divided into different sub-regions by some borders; see (Monfared & Durstewitz, 2020a;b) for more details. Also, since the system (28) is a linear map in each sub-region, the k periodic points of Γk must belong to different sub-regions (at least two different sub-regions). Accordingly, based on (63) and (65), there exists some Ñ ∈ N such that for every ñ ≥ Ñ both ztkñ+s and zt∗k−s belong to the same sub-region and so the matrices WΩ(tkñ+s) and WΩ(t∗k−s) (s ∈ {0, 1, 2, · · · , k − 1}) are identical. Hence, for n ≥ N , ñ ≥ Ñ and j ∈ N, equation (66) becomes
kj−1∏ i=0 WΩ(tn+kj−1−i) = j∏ l=0 k∏ i=1 WΩ(tk(ñ+j−l)+k−i) = ( k−1∏ s=0 WΩ(t∗k−s) )j . (67)
Therefore, similar to the part (ii), we can prove for every z1 ∈ BΓk , ∂zT∂θ and ∂LT∂θ will also remain bounded.

heading: A.2.2 PROOF OF THEOREM 2, PART (II)
text: Proof. (ii) Let for every T > 2
LT := J ∗ T J ∗ T−1 · · · J∗2 . (68)
{LT }T∈N, T>2 is a sequence of matrices LT = [ l (T ) ij ] 1≤i,j≤M and, due to (13), limT→∞ ∥LT ∥ = ∞. Hence, there is at least one sub-sequence {l(Tn)mk }Tn∈N, Tn>2 (for some m, k ∈ {1, 2, · · · ,M}) such that limTn→∞ l(Tn)mk = ∞.
On the other hand
∂z∗T ∂θ = ∂+z∗T ∂θ + T−2∑ t=1 ( t−1∏ r=0 J∗T−r ) ∂+z∗T−t ∂θ . (69)
Moreover, there exists some N > 2 such that (for t = T −N + 1) ∂+z∗N−1
∂θ ̸= 0. (70)
For θ as the k-th element of a parameter vector θ (or belonging to the k-th row of a parameter matrix θ), the term ( T−N∏
r=0
J∗T−r
) ∂+z∗N−1
∂θ (71)
is a vector in which the i-th element is l(T )ik ∂+z∗k,N−1 ∂θ .
Since limTn→∞ l (Tn) mk = ∞, due to (70) limTn→∞ l (Tn) mk ∂+z∗k,N−1 ∂θ = ∞, which implies ∂z∗T ∂θ will diverge as T → ∞. Similarly, by (21), we can prove ∂L ∗ T
∂θ is divergent for T → ∞. By Oseledec’s multiplicative ergodic Theorem, the results also hold for every z1 ∈ BΓ∗ .

heading: A.2.3 PROOF OF THEOREM 3
text: Proof. Let Γ = {z1, z2, . . . zT , · · · } be a quasi-periodic attractor. Then, the largest Lyapunov exponent of Γ is
λ = lim T→∞
1 T ln ∥JT JT−1 · · · J2∥ = lim T→∞ 1 T ln ∥∥∥∥∂zT∂z1 ∥∥∥∥ = 0. (72)
We prove for every 0 < ϵ < 1
lim T→∞ (1− ϵ)T−1 < lim T→∞ ∥∥∥∥∂zT∂z1 ∥∥∥∥ < limT→∞(1 + ϵ)T−1. (73)
For this purpose, we show ∀ 0 < ϵ < 1
(I) limT→∞(1− ϵ)T−1 < limT→∞ ∥∥∥∂zT∂z1 ∥∥∥, and
(II) limT→∞ ∥∥∥∂zT∂z1 ∥∥∥ < limT→∞(1 + ϵ)T−1.
Assume for the sake of contradiction that (I) does not hold. Then there exists some 0 < ϵ < 1 such that
lim T→∞ (1− ϵ)T−1 ≥ lim T→∞ ∥∥∥∥∂zT∂z1 ∥∥∥∥ . (74)
Therefore
∃T0 > 1 s.t. ∀T ≥ T0 =⇒ (1− ϵ)T−1 ≥ ∥∥∥∥∂zT∂z1 ∥∥∥∥ , (75)
and so
∃T0 > 1 s.t. ∀T ≥ T0 =⇒ ln(1− ϵ)T−1 T − 1 ≥ ln ∥∥∥∂zT∂z1 ∥∥∥ T − 1 . (76)
Consequently, due to (72), for T → ∞ we have ln(1 − ϵ) ≥ 0. This implies ϵ ≤ 0, which is a contradiction.
Similarly if we assume (II) is not true, then there exists some 0 < ϵ < 1 such that
lim T→∞ ∥∥∥∥∂zT∂z1 ∥∥∥∥ ≥ limT→∞(1 + ϵ)T−1. (77)
Thereby
∃T0 > 1 s.t. ∀T ≥ T0 =⇒ ∥∥∥∥∂zT∂z1 ∥∥∥∥ ≥ (1 + ϵ)T−1, (78) and thus
∃T0 > 1 s.t. ∀T ≥ T0 =⇒ ln ∥∥∥∂zT∂z1 ∥∥∥ T − 1 ≥ ln(1 + ϵ)T−1 T − 1 . (79)
This means ln(1 + ϵ) ≤ 0 as T → ∞, i.e. ϵ ≤ 0, which is a contradiction. Therefore (14) holds for Γ and also, according to Oseledec’s multiplicative ergodic Theorem, for every z1 in the basin of attraction of Γ.

heading: A.3 ADDITIONAL RESULTS ON RELATION BETWEEN DYNAMICS AND GRADIENTS

heading: A.3.1 FURTHER RESULTS AND REMARKS RELATED TO THEOREM 2
text: Remark 4. The result of Theorem 2 also holds for unstable orbits {z1, z2, z3, · · · } with positive largest Lyapunov exponent. Trivially, for such orbits that diverge to infinity (unbounded latent states) gradients of the loss function will explode as T → ∞. Remark 5. For RNNs with ReLU activation functions there are finite compartments in the phase space each with a different functional form. In such a case, to define the largest Lyapunov exponent of Γ∗, in the proof of Theorem 2 we assume that Γ∗ never maps to the points of the borders.
Based on Theorem 2, we can also formulate the necessary conditions for chaos and diverging gradients in standard RNNs with particular activation functions by considering the norms of their recurrence matrix, for which the following Corollary provides the basis: Corollary 1. Let for a standard RNN∥∥diag(f ′(Wzt−1 + Bst + h))∥∥ ≤ γ < ∞. (80)

heading: If the RNN is chaotic, then ∥W ∥ γ > 1 .
text: Proof. Assume for the sake of contradiction that ∥W ∥ γ ≤ 1 . From∥∥∥∥∥∥ ∏
2<t≤T
W diag ( f ′ ( Wzt−1 + Bst + h ))∥∥∥∥∥∥ ≤ ∏
2<t≤T ∥∥W diag(f ′(Wzt−1 + Bst + h))∥∥ ≤ (∥W ∥ γ)T−2, (81)
it is concluded that limT→∞ ∥∥∥∏2<t≤T W diag(f ′(Wzt−1 + Bst + h))∥∥∥ < ∞ , which contradicts (13). This means ∥W ∥ γ > 1 is a necessary condition for the standard RNN to be chaotic.
Remark 6. For RNN with the tanh and sigmoid activation functions γ = 1 and γ = 14 , respectively. Thus, by Corollary 1, the necessary conditions for chaos in these two cases are ∥W ∥ > 1 and ∥W ∥ > 4 , respectively.

heading: A.3.2 OTHER CONNECTIONS BETWEEN DYNAMICS AND GRADIENTS
text: As sect. 3 elucidated, there is a direct link between the norms of the Jacobians of the RNN along trajectories and the EVGP. By observing this link, we can formulate some general conditions that will have implications for the behavior of the gradients regardless of the limiting behavior of the RNN, as collected in the following theorem: Theorem 4. Let Oz1 = {z1, z2, . . . zT , · · · } be a sequence (orbit) generated by an RNN Fθ ∈ R parameterized by θ, and PT := JT − I, T = 2, 3, · · · .
(i) Assume that Oz1 is an orbit for which ∥∥∥∂+zT∂θ ∥∥∥ ≤ ξ ∀t. If ∑∞T=2 ∥JT ∥ < ∞, then the
Jacobian ∂zT∂z1 , the tangent vector ∂zT ∂θ and thus the gradient of the loss function, ∂LT ∂θ , will be bounded for T → ∞.
(ii) If ∑∞
T=2 ∥PT ∥ < ∞, then the Jacobian ∂zT∂z1 will neither vanish nor explode as T → ∞.
Proof. Let ∥.∥ be any matrix norm satisfying ∥A1A2∥ ≤ ∥A1∥ ∥A2∥.
(i) By boundedness of ∂ +zT ∂θ we have∥∥∥∥∂zT∂θ ∥∥∥∥ = ∥∥∥∥∥∂+zT∂θ + T−2∑ t=1 ( t−1∏ r=0 JT−r ) ∂+zT−t ∂θ ∥∥∥∥∥ ≤ ξ ( 1 +
T−2∑ t=1 ∥∥∥∥∥ t−1∏ r=0 JT−r ∥∥∥∥∥ ) ≤ ξ ( 1 + T−2∑ t=1 t−1∏ r=0 ∥JT−r∥ ) . (82)
Moreover,
1 + T−2∑ t=1 t−1∏ r=0 ∥JT−r∥ ≤ 1 + ∑ p ∥Jp∥ + ∑ p<q ∥Jp∥ ∥Jq∥ + ∑ p<q<r ∥Jp∥ ∥Jq∥ ∥Jr∥ + · · ·
= ( 1 + ∥JT ∥ )( 1 + ∥JT−1∥ ) · · · ( 1 + ∥J2∥ ) =: T∏ t=2 ( 1 + ∥Jt∥ ) . (83)
Since ∑∞ T=2 ∥JT ∥ converges, according to (Wedderburn, 1964), the infinite products ∏∞ T=2 ( 1 +
∥JT ∥ ) in (83) converge to a finite number K̃ ̸= 0. Consequently, by (82) and (83)
lim T→∞ ∥∥∥∥∂zT∂θ ∥∥∥∥ ≤ K̃ < ∞, (84)
which implies ∂LT∂θ will be bounded for T → ∞.
Furthermore
lim T→∞ ∥∥∥∥∂zT∂z1 ∥∥∥∥ ≤ ∞∏
T=2
∥JT ∥ := lim T→∞
( ∥JT ∥ ∥JT−1∥ · · · ∥J2∥ ) ≤ ∞∏ T=2 ( 1 + ∥JT ∥ ) ≤ K̃,
(85) which completes the proof. (ii) Since ∑∞
T=1 ∥PT ∥ < ∞, due to (Wedderburn, 1964) the infinite product ∞∏
T=2
( I + PT ) = ∞∏ T=2 JT := lim T→∞ JT JT−1 · · · J2, (86)
converges to a matrix K ̸= O, which implies
0 < lim T→∞ ∥∥∥∥∂zT∂z1 ∥∥∥∥ = ∥K∥ < ∞. (87)
Part (i) of Theorem 4 relaxes some of the conditions required in Theorem 1 for bounded gradients by imposing a Lipschitz condition on the immediate derivatives. Part (ii) generalizes conditions satisfied, for instance, in orthogonal (unitary) RNNs (Arjovsky et al., 2016; Henaff et al., 2016) or fully regularized PLRNNs (Schmidt et al., 2021).
Proposition 3. Let Oz1 = {z1, z2, . . . zT , · · · } be an orbit generated by an RNN Fθ ∈ R (parameterized by θ), and ∥JT ∥ ≠ 0, T ≥ 2. If ∑∞ T=2 ln ∥JT ∥ diverges to −∞, then the Jacobian ∂zT∂z1 vanishes as T tends to infinity.
Proof. For ∥JT ∥ ̸= 0, T ≥ 2, we have 0 ≤ ∥∥∥∥∂zT∂z1 ∥∥∥∥ ≤ ∥JT ∥ ∥JT−1∥ · · · ∥J2∥ = eln ∥JT ∥ eln ∥JT−1∥ · · · eln ∥J2∥ = e∑Tt=2 ln ∥Jt∥. (88)
Hence if ∑∞
T=2 ln ∥JT ∥ → −∞, then
lim T→∞ ∂zT ∂z1 = O. (89)

heading: A.4 EMPIRICAL EVALUATION: DATASETS
text: Lorenz attractor The Lorenz system (Lorenz, 1963) is a simplified model for atmospheric convection, given by
dx
dt = σ(y − x),
dy
dt = x(ρ− z)− y, (90)
dz dt = xy − βz.
The system is of particular interest for its chaotic regime and was studied here for σ = 16, ρ = 45.92 and β = 4. For these parameters the Lorenz system is known to have a maximal Lyapunov exponent λmax = 1.5 (Rosenstein et al., 1993). To generate a time series, the ODEs were integrated with a step size ∆t = 0.01 using scipy.integrate. Accordingly, the prediction time is τpred =
ln(2) ∆t λmax = 46.2.
Duffing oscillator The Duffing oscillator (Duffing, 1918) is an example of a periodically forced oscillator with nonlinear elasticity
ẍ+ δẋ+ βx+ αx3 = γ cos(ωt). (91)
Note that this system is non-autonomous, that is externally forced due to the r.h.s. of eqn. 91. The following parameters were chosen to arrive at a chaotically forced oscillator: α = 1.0, β = −1.0, δ = 0.1, γ = 0.35, and ω = 1.4. For these parameters the Duffing oscillator has a maximum Lyapunov exponent of λmax = 0.0995. The dataset used here was created with the code from (Gilpin, 2021) as a three dimensional embedding with step size ∆t = 0.17. The prediction time is τpred = 39.28.
Rössler system Another prime textbook example for a chaotic system is the Rössler system (Rössler, 1976) given by:
dx dt = −y − z,
dy dt = x+ ay, (92)
dz dt = b+ z(x− c).
For the parameters a = 0.15, b = 0.2 and c = 10, the maximal Lyapunov exponent is λmax = 0.09 (Rosenstein et al., 1993). To arrive at a time series, a step size of ∆t = 0.1 was chosen for integration. This gives us a prediction time of τpred = 77.0 for this system.
Mackey-Glass equation The Mackey-Glass equation (Glass & Mackey, 1979) is a nonlinear time delay differential equation
ẋ = β xρ
1 + xnρ − γx with β, γ, ρ > 0. (93)
Here xρ represents the value of the variable x at time t − ρ (note that strictly, mathematically, this makes the system infinite-dimensional). Choosing the parameters to be β = 2, γ = 1.0, n = 9.65, and ρ = 2.0, leads to chaotic behavior with a maximum Lyapunov exponent of λmax = 0.21. The dataset was created as a 10-dimensional embedding with the code from (Gilpin, 2021) using ∆t = 0.04. This yields a prediction time of τpred = 82.2.
Empirical temperature time series This time series was recorded at the Weather Station at the Max Planck Institute for Biogeochemistry in Jena, Germany, spanning the time period between 2009 and 2016, and reassembled by François Chollet for the book Deep Learning with Python. The data set can be accessed at https://www.kaggle.com/pankrzysiu/weather-archive-jena. To expose the underlying chaotic dynamics of the time series, trends and yearly cycles were removed, and nonlinear noise-reduction was performed (using ghkss from TISEAN, see also Kantz et al. (1993)). Fig. 4 (a) shows a snippet of the temperature data in comparison with the denoised time-series. High-frequency noise was further reduced through Gaussian kernel smoothing (σ = 200), and the resulting time series was sub-sampled (every 5th data point was retained). Fig. 4 (b) clearly reveals a fractional dimension of Deff = 2.8 for the de-noised and smoothed time-series. This strongly suggests that the dynamics governing the time series are chaotic. We created a time delay embedding (Kantz & Schreiber, 2003) with m = 5 (estimated by the false nearest neighbor technique, see Kennel et al. (1992)) and delay ∆t = 500 (obtained as the first minimum of the mutual information). The first three embedding dimensions are shown in Fig. 4(c). The maximal Lyapunov exponent of this time series was determined with lyap r from TISEAN (Hegger et al., 1999) to be λmax = 0.016, see Fig. 3(a). This value is in close agreement with the literature (Millán et al., 2010). The predictability time of this system is estimated to be τpred = 43.3.
All datasets used were standardized (i.e., centered with unit variance) prior to training.

heading: A.5 EMPIRICAL EVALUATION: MEASURES OF RECONSTRUCTION QUALITY
text: Attractor overlap To asses the geometrical similarity of the chaotic attractor produced by the RNN to the one underlying the observations, we calculate the Kullback-Leibler divergence of the ground truth distribution ptrue(x) and the distribution pgen(x|z) generated by RNN simulation. To do so in practice, we employ a binning approximation (see (Koppe et al., 2019))
Dstsp (ptrue(x), pgen(x | z)) ≈ K∑
k=1
p̂ (k) true(x) log
( p̂ (k) true(x)
p̂ (k) gen(x | z)
) ,
where K is the total number of bins, and p̂(k)true (x) and p̂ (k) gen (x | z) are estimates obtained as relative frequencies through sampling trajectories from the observed time-series and the trained RNN, respectively.
Power-spectral correlations Since in DS reconstruction we aim to capture invariant (timeindependent) properties of the underlying system, besides the geometrical agreement, we compare the similarity in true and RNN-reconstructed power spectra. To do so, we generate a time series of length 100, 000 from the RNN and calculate its power spectrum using the fast Fourier transform (scipy.fft). To reduce the influence of noise we apply Gaussian kernel smoothing and cut off the long high-frequency tails of the spectra. The dimension-wise correlation between observed and generated spectra are then averaged to give the PSC.

heading: A.6 FURTHER EMPIRICAL EVALUATIONS
text: A.6.1 RECONSTRUCTION: RÖSSLER SYSTEM

heading: A.6.2 RECONSTRUCTION: HIGH-DIMENSIONAL MACKEY-GLASS SYSTEM

heading: A.6.3 RECONSTRUCTION: PARTIALLY OBSERVED LORENZ SYSTEM
text: For this evaluation we trained models only on the variables {y, z} of the Lorenz system, eqn. 90. In order to compute the attractor overlap (Dstsp) in the true state space, however, after training the observation matrix B was recomputed by linearly regressing the first 10 latent states onto the first 10 observations from all three Lorenz variables in eqn. 90.

heading: A.6.4 OTHER INITIALIZATION PROCEDURES: STANDARD BATCHING WITH ZERO RESETTING
text: A common procedure in training RNNs is partitioning the time series into chunks (as we did based on the Lyapunov spectrum), but then simply resetting the hidden states at the beginning of each chunk (window) to 0. Formally this would mean that we do not force the trajectory back on track as in our approach, but instead may kick it off the track. To illustrate this, here we trained an LSTM on chunks (windows) with a length given by the optimal τ (τopt = 30 for the Lorenz system), but then initialized the hidden states to 0 at the beginning of each window. The performance obtained this way is indicated by the green dashed line in Fig. 9 below.

heading: A.6.5 ELECTROENCEPHALOGRAM (EEG) DATA
text: We used EEG data recorded by Schalk et al. (2004) and provided on PhysioNet (Goldberger et al., 2000), from which we took the baseline recording of the first patient for our analysis. Preprocessing was performed as outlined above for the temperature time-series, i.e. we applied nonlinear noise-reduction (see Fig.10 (a)) and Gaussian kernel smoothing (σ = 30). Fig. 10 (b) indicates a fractional dimension Deff = 1.8 for the de-noised and smoothed times series. We created a time delay embedding with an embedding dimension of m = 10 and a delay time of ∆t = 500. The maximal Lyapunov exponent for this time series was determined to be λmax = 0.007, see Fig. 11 (a). To ease training, the time-series was sub-sampled and only every third data point was retained. With this, we obtain a predictability time τpred = 33.01.
Full reconstruction of the dynamics from complex and noisy EEG signals is a very ambitious. Here we therefore focused on the system’s short-term behavior and combined multiple shorter trajectory bits (of length T = 70) in the calculation of our measure for geometrical agreement in state space.

heading: A.7 SPARSELY FORCED BPTT
text: Loss truncation One implicit consequence of the teacher forcing, eqn. (16), is the interruption of the hidden-to-hidden connections at these time points. More specifically, if the system is forced at time t ∈ T , then there is no connection between zt and zt+1, that is
Jt+1 = ∂zt+1 ∂zt = ∂RNN(z̃t) ∂zt = 0. (94)
To see how these vanishing Jacobians truncate the loss gradients w.r.t to some parameter θ, let us focus on the loss gradients immediately after the forcing,
∂Lt+1 ∂θ = ∂Lt+1 ∂zt+1 t+1∑ k=1 ∂zt+1 ∂zk ∂+zk ∂θ
= ∂Lt+1 ∂zt+1 ( ∂+zt+1 ∂θ + t∑ k=1
∂zt+1 ∂zk︸ ︷︷ ︸
=0 , because of (94)
∂+zk ∂θ )
= ∂Lt+1 ∂zt+1 ∂+zt+1 ∂θ . (95)
Eqn. (95) shows that sparsely forced BPTT implicitly truncates the loss gradients because it interrupts the hidden-to-hidden connection from zt to zt+1 for t ∈ T . More generally, defining t̃ := max{t′ ∈ T : t′ ≤ t}, the overall loss gradients are truncated to
∂L ∂θ = T∑ t=1 ∂Lt ∂zt t∑ k=1 ∂zt ∂zk ∂+zk ∂θ
tr. = T∑ t=1 ∂Lt ∂zt t∑ k=t̃ ∂zt ∂zk ∂+zk ∂θ . (96)

Review summaries:
The paper connects the gradient behavior of recurrent neural networks with its dynamics through the maximum Lyapunov exponent of a trajectory, which is generated by the respective RNN. In particular, several connections between RNNs which produce specific dynamics (e.g. with stable fixed points or limit-cycles) and the mitigation of the exploding gradient problem are proved. Moreover, the paper shows that RNNs generating a chaotic trajectory always suffer from the exploding gradient problem.  That being said, the paper emphasizes that in order to approximate chaotic systems, the RNN has to be able to produce chaotic trajectories, which is ruled out by design by several state-of-the-art RNNs that are constructed to solve the exploding and vanishing gradient problem. However, given that these RNNs - at least the moment the optimizer finds parameters such that the resulting RNN produces chaotic trajectories - will always run into the exploding gradient problem, the paper suggests an algorithmic fix for that, by forcing the hidden-states of the RNN to be the pseudo-inverse of the underlying ground-truth trajectory at equidistantly distributed points in time during training. This cuts off the gradients at these points in time and forces the output states to be close to the ground-truth trajectories. The interval for that is chosen based on the maximum Lyapunov exponent. In the paper, three experiments are provided, where two of them are based on approximating a chaotic dynamical system (Lorenz and Rössler), while one is based on real-world data from temperature forecasting.

This paper studies the exploding and vanishing gradient problem (EVGP) in the sequential modelling with chaotic dynamics. Theoretical analysis of the EVGP, based  on the relationship between loss gradients during RNN training and Lyapunov spectrum of RNN-generated orbits, is provided. Inspired by this analysis, an alternative to BPTT training algorithm is proposed, named sparsely forced BPTT, that forces the diverging dynamics to conform to the true trajectory, at regular time interval provided by the  Lyapunov spectrum.

The paper looks at the asymptotic behavior of the Jacobian of various RNN variants (standard, LSTM, GRU, PWL-RNN) when realizing stable vs chaotic dynamics. In particular, the paper shows that in the chaotic setting, the gradients asymptotically explode, i.e., when learning chaotic dynamics, the gradients have to explode (asymptotically). The paper proposes to overcome this limitation by truncating the BP length and applying a teacher forcing method that periodically projects the observation onto the hidden state during training.

The authors provide a study of gradients of the recurrent neural networks using Lyapunov exponents. They claim to offer a simple and effective training for the chaotic data. Some simulation results are provided.