Recurrent neural networks: vanishing and exploding gradients are not the end of the story

Authors: nicolas  zucchet
antonio  orvieto
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 14:46:39

Abstract: Recurrent neural networks (RNNs) notoriously struggle to learn long-term memories, primarily due to vanishing and exploding gradients. The recent success of state-space models (SSMs), a subclass of RNNs, to overcome such difficulties challenges our theoretical understanding. In this paper, we delve into the optimization challenges of RNNs and discover that, as the memory of a network increases, changes in its parameters result in increasingly large output variations, making gradient-based learning highly sensitive, even without exploding gradients. Our analysis further reveals the importance of the element-wise recurrence design pattern combined with careful parametrizations in mitigating this effect. This feature is present in SSMs, as well as in other architectures, such as LSTMs. Overall, our insights provide a new explanation for some of the difficulties in gradient-based learning of RNNs and why some architectures perform better than others.

Main paper: 
 vanishing and exploding gradients: 
 Let us first introduce the notations we will be using throughout the rest of the paper. We consider a recurrent neural network with hidden state h t , update function f θ parametrized by θ, and input sequence (x t ). The average performance of the network is measured by a loss L. We haveEQUATIONThe gradient of the instantaneous loss L t with respect to the parameters θ is then equal toEQUATIONIn the equation above, we used ∂ to denote partial derivatives and d for total derivatives. Using this notation enables us to distinguish between ∂ ht L t , which corresponds to the error backpropagated from the current loss term to the hidden state through the readout function, and d ht L, which accumulates the errors that are backpropagated through the future hidden state values. In particular, ∂ ht L = ∂ ht L t and d ht L = ∂ ht L t (h t ) + t ′ >t d ht L t ′ (h t ′ ). When stacking several recurrent layers on top of each other, ∂ ht L corresponds to the current error being backpropagated to the hidden state h t through the hierarchy of the network and d ht L to future error signals backpropagated through the recurrence.Early work [5] highlighted the difficulty for gradient descent to make recurrent neural networks remember past inputs that will later be useful to produce a desired behavior. This is due to the fact that error signals flowing backward in time tend to either explode or vanish. The key quantity isEQUATIONOne can remark that this quantity exponentially converges to 0 when the spectral radius of the Jacobian ∂ h f θ is upper bounded by a constant strictly smaller than 1, and can exponentially explode if there exists some component bigger than 1. The error signal at time t backpropagated to time t ′ behaves similarly, as d h t ′ L t = ∂ ht L t d h ′ t h t . Gradient-based learning of long-term memories is thus difficult: the contribution of past hidden states to the current loss becomes either negligible or predominant as the time span considered increases.Since then, the analysis has been refined [6] [7] [8] and the development of recurrent architectures has mostly been driven by the desire to solve this pathological issue. Most famously, the LSTM [3] unit, and later on the GRU [23] , solve this problem by using memory neurons that facilitate direct information storage and retrieval, and thus facilitate error backpropagation. Other approaches to solving this problem, to name a few, involve gradient clipping [24, 8] , activity normalization [25] [26] [27] , careful weight initialization [28, 29] or enforcing architectural constraints such as hierarchical processing [30, 31] , orthogonal weight matrices [32] [33] [34] and oscillations [35] [36] [37] . 
the curse of memory: 
 According to common deep learning wisdom, it is often believed that solving the vanishing and exploding gradients problem enables recurrent neural networks to learn long-term dependencies. We challenge this view and question: is solving those issues really enough to ensure well-behaved loss landscapes? We answer negatively by showing that gradients can explode as the memories of the network are kept for longer, even when the dynamics of the network remains stable.Recurrent neural networks have something special: the very same update function f θ is applied over and over. Therefore, modifying the parameters θ will not only influence one update, as changing the weights of a given layer in a feedforward neural network would, but all of them. As the memory of the network gets longer, the hidden states keep a trace of the effects of more updates. Hidden states thus become increasingly sensitive to parameter changes. This is the curse of memory. We borrow A. Evolution of the second moment of d λ h t as a function of the recurrent parameter λ and of the input x auto-correlation decay rate ρ, when h t+1 = λh t + x t . As the memory of the network increases (λ → 1), h t becomes more sensitive to changes in λ, particularly as the elements in the input sequence are more correlated (ρ → 1). The explosion of d λ h t is faster than the one of h t , as highlighted with the grey line obtained for ρ = 1. See Section 2.2 for more detail. B, C. Illustration of the phenomenon on the toy one-dimensional teacher-student task of Section 4.1, in which the teacher is parametrized by a real number λ * and the student by a complex number λ. In B, λ varies on the real axis, and it varies on the circle of radius λ * parametrized by θ in C. The loss becomes sharper as information is kept longer in memory, making gradient-based optimization nearly impossible.the term from [38, 39] , although we use it in a different context, and note that Martens and Sutskever [40] hypothesized that such a phenomenon could arise in RNNs and hinder their optimization.Let us formalize our intuition by considering the sensitivity of the hidden state h t on the parameters θ:EQUATIONWhen information stays in the network's memory for longer, the number of non-negligible Jacobian d h t ′ h t terms increases. As a result, the magnitude of this sensitivity increases when the network encodes longer-term dependencies, and learning θ becomes trickier. It is crucial to observe that this phenomenon arises even when exploding gradients are removed from the picture by constraining the eigenvalues of the recurrent Jacobian to be smaller than one and ensuring that the network dynamics remains stable. The rest of this section will be dedicated to studying this behavior quantitatively.We study how hidden state and gradient magnitudes evolve as the network encodes longer-term dependencies. Ideally, these quantities do not vanish or explode, as it improves the conditioning of the loss landscape [41] and eases optimization [42, 43] . We operate under the following assumptions: a) Linear diagonal recurrent neural networks. We restrict ourselves to update functions of the form f θ (h t , x t+1 ) = λ ⊙ h t + x t+1 with λ a vector of the size of h t and ⊙ the element-wise product. For ease of exposition, we present results for real-valued λ here; see Appendix B.2 for the complex-valued setting. While this assumption is strong, it allows us to identify some crucial mechanisms and models like S4 [17] , S5 [19] and LRUs [20] satisfy it. We later show that our analysis can model some features of more sophisticated networks. Note that we do not consider the input and readout mappings usually featured in recurrent layers as they are feedforward layers and signal propagation within them is already well understood [e.g., 44, 45] . b) Infinite time horizon. We consider infinite sequences and initialize the network dynamics at t 0 = -∞. It simplifies our calculations while being a reasonable assumption when the sequences considered are longer than the characteristic timescales of the dependencies we want to learn. c) Wide-sense stationarity. We assume the different quantities that the network receives, which include the inputs x t , to be wide-sense stationary (WSS). A random process X t is said to be WSS if its auto-correlation function is independent of time, that is, for all t ∈ Z and ∆ ∈ Z, E X [X t+∆ X t ] =: R X (∆), where E X denotes the expectation over the data. It corresponds to assuming that the statistics of the data are invariant to time shifts. This is a standard assumption when analyzing stochastic processes [46] . It keeps our calculations concise and does not qualitatively affect our conclusions (cf. Section 5). We discuss how to relax it in Appendix B.2.4.We are now equipped to analyze signal propagation in one recurrent layer, both in the forward and backward passes. We show that both hidden states and backpropagated errors explode as |λ| → 1.Forward pass. Here, we are interested in understanding how the hidden state second moment E[h 2 t ] evolves as a function of λ and of the input auto-correlation function R x . After a calculation that we defer to Appendix B.2, we obtainEQUATIONImportantly, this quantity goes to infinity as longer-term dependencies are encoded within the network, that is |λ| → E dh t dλ 2 = 1 + λ 2 (1 -λ 2 ) 3   R x (0) + 2 ∆≥1 λ ∆ R x (∆)   + 2 (1 -λ 2 ) 2   ∆≥1 ∆λ ∆ R x (∆)   . (6)We plot the exact behavior of this quantity when the auto-correlation of x satisfies R x (∆) = ρ |∆| on Figure 1 and refer the interested reader to Appendix B.2 for a derivation of Equation 6. More generally, the hidden state of the network, and thus its final output, becomes increasingly sensitive to changes in recurrent parameters as the network reaches the edge of dynamical stability (|λ| → 1).The last quantity we need to consider is the error that is backpropagated to the inputs x of the recurrent layer. It can be observed that the backward pass is dual to the forward pass in the sense that it is a recurrent process that receives backpropagated errors ∂ ht L and it runs in reverse time:EQUATIONin which we made use of ∂ xt h t = 1. It follows that the analysis we did for the forward pass also holds here. Crucially, this implies that the explosion behavior will be most significant for the recurrent parameters rather than for potential input or readout weights.We now generalize our results to fully connected linear recurrent neural networks of the form h t+1 = Ah t + x t . For the sake of the analysis, we assume that A is complex diagonalizable, that is there exists a complex-valued matrix P and a complex-valued vector λ such that A = P diag(λ)P -1 . Note that this occurs with probability one under random initialization of A [20] . In this case,EQUATIONandEQUATIONFrom the analysis above, we know that the dominating term in the limit |λ| → 1 among ∂ P h t , d λ h t and d -1 P h t is d λ h t , as P and P -1 act as readout and input weights. Given that all other terms do not directly depend on the magnitude of λ, we have thatd A h t ≃ ∂ h diag t h t d λ h diag t ∂ A λ; cf. Appendix B.2.3 for formal statements. This has two consequences: First, the sensitivity of h t on A will explode as longer memories are encoded and this directly comes from the eigenvalues of A. Second, as each entry of A typically impacts all eigenvalues of the matrix, the explosion behavior will be distributed across all entries, whereas it was concentrated on the eigenvalues for the diagonal case. We will later observe that this has significant practical consequences and partly explains why fully connected linear RNNs are difficult to train. As a side note, we remark that enforcing the matrix A to be orthogonal solves vanishing and exploding gradient issues but these weights may remain sensitive to learn because of the curse of memory. 
mitigating the curse of memory: 
 We have discussed the sensitivity of recurrent networks to parameter updates. Given this problem, how can it be mitigated? We show that recurrent networks with diagonal connectivity are particularly well suited for this purpose. Besides enabling control over the Jacobian and avoiding exploding gradients, they facilitate the mitigation of the curse of memory. We additionally highlight that deep state-space models and gated RNNs inherently incorporate such mechanisms.Both forward and backward passes explode as the network encodes longer memories. When h t+1 = λh t + x t+1 , we argue that it is relatively straightforward to mitigate this effect. We aim to keepE[h 2 t ], E[(d λ h t ) 2 ] and E[(d xt h t ) 2 ]independent of λ, similar to initialization schemes that maintain the magnitude of neural activity constant in deep networks [44, 45] , regardless of the layer width [42, 47, 43] . ] over all λ values when the input auto-correlation satisfies R x (∆) = ρ |∆| with ρ = 0, but does not manage do to so for other type of distributions (ρ ̸ = 0). Here, we use γ(λ) = √ 1 -λ 2 , decouple it from λ when differentiating, and take λ = exp(-exp(ν)), as in [20] . The grey line indicates the value the two quantities take without any normalization and reparametrization, when ρ = 1.t ] to stay constant is to introduce a scaling factor γ(λ) applied to the inputs a neuron receives, that satisfiesγ(λ) 2 E[h 2 t ] = Θ(1). Given that the backward propagation of output errors to inputs is dual to the forward pass, the role of γ in the backward pass will be similar. The value γ needs to take therefore both depends on the input distribution to normalize the forward pass, as well as on the output error distribution to normalize the backward pass. Perfect normalization is likely unrealistic, but some normalization can help, as shown in Figure 2 .A.Eigenvalue reparametrization. We are now left with keeping the gradient of the loss with respect to λ under control. Input normalization partly reduces the memoryinduced exploding effect, but not entirely as the variance of d λ h t is much larger than the one of h t (cf. Fig. 1 .A). Reparametrization can close that gap. Indeed, if λ is parametrized by ω, we have that d ω h t = d λ h t d ω λ. Choosing a parameterization that is more and more granular as λ goes to 1 thus helps in keeping the magnitude ofd ω h t constant. Assuming γ is independent of λ for simplicity, achieving E[(d ω h t ) 2 ] = Θ(1) requires solving the differ- ential equation γ(λ) 2 λ ′ (ω) 2 E[(d λ h t ) 2 ] = 1.While deriving a universal optimal parametrization is again unrealistic due to dependency on the input distribution, reparametrization definitely helps, as shown in Figure 2 .B. Figure 6 illustrates how it affects the loss landscape.The case of complex numbers. We have not yet discussed the case λ ∈ C, relevant for SSMs such as S4 [17] . We extend our analysis to complex-valued λ in Appendix B.3.2. Briefly, it reveals that changes in the magnitude of λ have a similar impact as in the real case, but this similarity does not extend to the angle. To keep the sensitivity on the angle constant, its parametrization must depend on the magnitude of λ. However, doing so hurts learning, particularly far from optimality, as we exemplify in Appendix B.3.2. A key implication of this analysis is that the sensitivity of the hidden state on the angle of λ explodes as its magnitude approaches 1.Deep state-space models, as well as gated RNNs, feature some form of normalization and reparametrization which help keeping signal propagation under control. We discuss how below.Deep state-space models. Deep SSMs were originally motivated as discretizations of the differential equation ḣ = Ah + Bx [16] . Naïve discretization of the differential equation yields h t+1 = (Id + dtA)h t + dtBx t+1 which already provides some input normalization. More elaborate discretization schemes, such as the zero-order hold, effectively reparametrize the A matrix, e.g. with exp(dtA). Here, diagonalization arises from computational efficiency and simplicity reasons [18] . While such models can approximate any smooth mappings [48, 49] , their expressivity remains limited [50] . The next generation of these models, including Mamba [21] , incorporates input-dependent gates which modulate dt depending on the input x t . The theory we developed above does not strictly apply to this setting as dt is no longer constant. However, since the recurrence Jacobian remains diagonal, we expect the qualitative behaviors we analyzed to remain. Gated RNNs. While the original motivation behind gated RNNs such as LSTMs [3] or GRUs [23] largely differs from the one of SSMs, they share similar mechanisms. In these networks, the memories stored in hidden neurons can be erased through a forget gate, and incoming inputs can selectively be written in memory through an input gate. Mathematically, this corresponds to hidden state updates of the form h t+1 = f t+1 ⊙h t +i t+1 ⊙x t+1 , with the forget f t+1 and input i t+1 gates being independent non-linear functions of x t+1 and h t . The forget gate is akin to λ and usually involves a sigmoid non-linearity, which has a similar effect as reparametrizing λ in the backward pass. The input gate can act as an input normalization depending on the initialization of the network or if is coupled to the forget gate as in the GRU (f t = 1 -i t ) [29] . Importantly, the gates here depend on the hidden states and thus make the Jacobian ∂ ht h t+1 non diagonal. Yet, we argue that these architectures still have a bias towards diagonality. Indeed, the contributions of the hidden state through the forget and input gates are indirect, and they can be ignored when the weights connecting the hidden states to the gates are small. Empirically, we find that GRUs lie in this regime at initialization, cf. Section D.2, so that our theory accurately captures signal propagation in GRUs. We additionally confirm that signal propagation is well behaved in gated RNNs in Section 5. In regimes in which this approximation does not hold, studying signal propagation requires a much more sophisticated analysis than the one we have done here [51] . 
a linear teacher-student analysis: 
 We start our empirical analysis with a teacher-student task using linear recurrent networks [52] . This is arguably the simplest setting in which one can train recurrent networks. Yet, as we shall see, it is already remarkably complex and captures some of the differences between different architectures observed in more realistic settings [20] . It additionally makes it possible to control the characteristic time constants of the data, which is only possible with synthetic data.Our investigation starts with the one-dimensional setting to provide an intuitive illustration of the consequences of the curse of memory on the loss landscape. We then address the general setting and observe that linear networks indeed suffer from the curse of memory, and that the remedies we studied in the last section are effective. We additionally find that diagonality greatly modifies the structure of the loss landscape and helps optimizers with adaptive learning rates to compensate for eventual increased sensitivities.We first consider a student and a teacher following the one-dimensional dynamics h t+1 = λh t + x t+1 , with complex-valued parameter λ for the student and λ * for the teacher. For simplicity, we independently draw each x t+1 from a unit normal distribution (its autocorrelation function is R x (∆) = δ ∆=0 ) and note that other input distributions do not qualitatively change the results. The performance of the student is measured by a loss L that averages the per time-step losses L t := This simple model already captures two key difficulties of gradient-based learning of recurrent neural networks. In Figure 1 , we plot the resulting loss landscape for different λ * values, when λ evolves on the positive part of the real axis (Fig. 1 .B) and when it evolves on the circle of radius |λ * | in the complex plane (Fig. 1 .C). We restrict λs to have absolute values smaller than one: exploding gradients are out of the picture. Still, two difficulties for gradient-based learning appear here. On one side, vanishing gradients lead to flat loss regions that are hard to escape. On the other side, the loss sharpens as the student encodes longer memories because of the curse of memory. As a consequence, gradient-based optimization is extremely tedious, already in this simple example.We now move to the general case in which the teacher evolves according toEQUATION)with h t ∈ R n , x t ∈ R drawn i.i.d. from N (0, 1), A ∈ R n×n , B ∈ R n×1 , C ∈ R 1×n and D ∈ R 1×1 .Here both inputs and outputs are scalars.Given the intuition we have developed so far, we expect fully connected linear recurrent neural networks to struggle solving the task when the teacher encodes longer memories, not only because of exploding gradients but also due to the curse of memory. Conversely, diagonality facilitates the eigenvalue reparametrization needed to avoid exploding gradients and make them better behaved. We run the following experiment to verify this intuition. We draw random teachers with hidden dimension n = 10 and transform the complex eigenvalues of the recurrent matrix A to have magnitudes close to a value ν 0 that we control1 . The larger ν 0 is, the longer the memories encoded by the teacher are. We train a linear RNN, as well as an LRU [20] , with hidden dimension 64 on this task. The students are therefore largely overparametrized. We chose the LRU architecture to represent deep SSMs due to its simplicity. This architecture uses input normalization and an exponential reparametrization of the eigenvalues, similar to what we analyze in Section 3. Both networks are trained using the Adam optimizer [53] and cosine annealing schedule for 10k steps, on batches of size 128. To ensure that we are in the infinite sequence length regime, we take the sequences to be of length 300, that is three times longer than the characteristic time scale of the teacher. Learning rates are tuned separately for each method and training distribution. The results, which we plot in Figure 3 .A, confirm our intuition: LRUs significantly outperform linear RNNs when long memories have to be learned, despite having 10 times fewer parameters.Next, we wonder which design choices behind the LRU architecture are crucial to this performance improvement. To this end, we interpolate between a linear RNN and an LRU in the following way: First, we restrict the weight matrix of the linear RNN to a block diagonal with blocks of size 2. Each block can represent a complex number, so the network can represent 32 complex numbers in total.We additionally double the number of hidden neurons. Second, we change those 2 × 2 blocks (and their input and output weights) to be complex numbers. Finally, we add the γ input normalization and the exponential parametrization to obtain the final LRU architecture. We report the results of this experiment in Figure 3 .B. Surprisingly, we find the gap comes from making the weight matrix block diagonal (4 × 4 blocks are here enough, cf. Figure 11 in the Appendix). Interestingly, this change reduces the number of parameters the model has and slightly reduces the model expressivity. An explanation of this behavior is therefore likely to be related to the optimization properties of those models. We confirm this hypothesis in the next section. Hessian of the loss at optimality, its 10 eigenvectors with greatest eigenvalues and its eigenspectra for a fully connected RNN (A) and a complex diagonal one (B). The spectra are almost the same. However, the top eigenvectors are concentrated on few coordinates for the complex diagonal one but not for the fully connected one. C, D. This structure makes it possible for Adam to efficiently deal with the extra sensitivity, as shown with the effective learning rates that it uses at the end of learning. For the fully connected one (C), Adam uses small learning rates to compensate for the sensitivity, whereas it can use larger ones for the complex diagonal one without hindering training stability. The horizontal grey line shows the learning rate used, which is here 10 -3 .So far, our results highlight the importance of having a close to diagonal recurrent connectivity matrix.In this section, we show that this parametrization alone does not mitigate any exploding behavior but modifies the structure of the loss landscape, making it possible for optimizers with adaptive learning rates to compensate for these behaviors.To demonstrate this, we consider the Hessian of the loss:EQUATION)If the network can perfectly fit the target data, which is the case in the experiments above, the second term vanishes at optimality. We plot the Hessian at optimality in Figure 4 .A and B for a standard linear recurrent network and one with complex diagonal parametrization, both with 4 hidden neurons (ν 0 = 0.99). We observe that the eigenvalue spectra are similar for the two architectures, both exhibiting large terms that are characteristic of the curse of memory, which makes learning with stochastic gradient descent almost impossible2 . However, their structures differ. For the fully connected linear RNN, the top eigenvectors are distributed over many coordinates, whereas they are concentrated on a few coordinates for the complex diagonal one. This feature aids adaptive optimization [e.g., 57]: adapting to large curvature is much easier for Adam when the pathological directions are aligned to the canonical basis. This is what we observe in practice.In Figure 4 .C and D, we compare the effective learning rate used by Adam, which we compute by providing a vector of ones to the optimizer. For the dense linear RNN, the adaptive learning rates cannot compensate for the intricate coupling between components, resulting in small learning rates overall. Conversely, the sensitive directions of complex diagonal RNNs are concentrated on few ] for the different recurrent layers and specific groups of parameters. For the complex diagonal RNN, the gradients of all parameters explode and in particular the ones of the recurrent parameters, whereas only the ones of the angle of λ explode for the LRU, consistently with the theory. Error signal propagation in GRUs is under control: the magnitude of the gradients is independent of ν 0 . The GRU-specific parameters exhibit smaller gradients than the feedforward (ff) ones. C. Layer normalization keeps the overall gradient magnitude under control in cRNNs. Batch normalization yields similar results. parameters, which adaptive learning rates can compensate for. This leads to more targeted and overall larger learning rates, significantly speeding up learning. As a side note, the complex eigenvalues of the teacher come in conjugate pairs. However, during training, the complex values of the complex RNN are not conjugates of each other, thereby increasing Hessian diagonality. Finally, performing this analysis for the LRU, we find that the Hessian spectrum is similar to the diagonal setting and that the exploding dimensions of the Hessian are almost exclusively due to the angle parameter, consistently with our theoretical analysis; see Figure 9 .A and C. The loss landscape for S4 model [17] can be qualitatively similar to the complex diagonal RNN or to the LRU, depending on which regime it is in; see Figure 9 .B and D.Before concluding this section, we investigate whether certain eigenvalue distributions can break the diagonal structure of the Hessian, thereby complicating optimization and increasing the need for eigenvalue reparametrization. In Appendix C.2, we provide a theoretical quantification of the intuitive result that more concentrated eigenvalues lead to less diagonal Hessian. Consequently, the performance gap between complex-valued diagonal networks and LRUs widens, although the former still greatly outperform their fully-connected counterpart (see Figure 10 ). An important corollary is that increasing the number of hidden neurons breaks the diagonal structure of the loss landscape, thus reducing the effectiveness of optimizers with adaptive learning rates in mitigating the curse of memory. This observation may explain why Orvieto et al. [20] reported a more substantial performance improvement from eigenvalue reparametrization than what we observe in our study (cf. block diagonal vs. LRU in Figure 3 .B). 
signal propagation in deep recurrent networks at initialization: 
 The ultimate goal of our theoretical quest is to gain insights into the training of practical recurrent network architectures. Specifically, we aim to verify whether the trends established theoretically and in controlled experiments hold in practice, by studying signal propagation at initialization on a realistic next-token prediction natural language processing task.To that matter, we provide sentences as input to deep recurrent networks that contain four blocks and use a next-token prediction loss to measure their performance. Each block consists of a recurrent layer followed by a feedforward gated linear unit [58] . By default, there are no normalization layers in this architecture. More details can be found in Appendix D.1. We empirically study how E[h 2 t ] and E[(d θ L) 2 ] evolve when the characteristic time scale of the recurrent layers, controlled through ν 0 , increases. We compare three different recurrent layers: a complex-valued diagonal RNN (cRNN), a LRU and a GRU initialized with the Chrono initialization [29] .The results are consistent with our theory. Complex-valued RNNs suffer from the curse of memory and recurrent parameters grow faster to infinity than the rest as ν 0 goes to 1. Perhaps more surprisingly, a finer grain analysis reveals that the gradient magnitude is independent of the layer; see Figure 13 . LRUs almost perfectly mitigate exploding behaviors in the forward pass (Figure 5 .A) as well as in the backward pass (Figure 5 .B), except for the angle parameter, consistently with our previous analysis. We also wonder whether layer normalization can replace the input normalization and reparametrization of the LRU. We find that it mitigates the memory-induced gradient explosion at the macroscopic level (Figure 5 .C), but it likely kills any learning signal for the smallest eigenvalues [20] . Finally, the GRU manages to keep the gradient magnitude constant over different characteristic time constants, consistently with the intuition we developed in Section 3.2. Preliminary experiments revealed that same trends also hold for LSTMs.The results presented above for the GRU align qualitatively with the intuition developed throughout the paper. We now consider how well our theory can quantitatively explain this behavior. The primary difference between our simple model and GRUs is that the λ values, referred to as forget gates in the GRU terminology, depend on both inputs and hidden states, and are therefore not constant. Interestingly, we find that GRUs almost behave like the diagonal linear RNNs we have focused on in this paper, particularly for slowly decaying recurrent neurons with high ν 0 values (see Appendix D.2). Consequently, applying our theory as if this context-dependency does not exist only introduces minor approximation errors, which we confirm empirically in Appendix D.3. Given the similarity of the Chrono initialization to those used in modern architectures like Mamba [21] and Hawk [22] , we expect our theory to also serve as a good proxy for studying signal propagation in these models at initialization. 
conclusion: 
 Vanishing and exploding gradients complicate the learning of recurrent networks, but solving these problems is not enough. We uncovered yet another difficulty of training such networks, which is rooted in their iterative nature and arises at the edge of dynamical stability. Reparametrizations and adaptive learning rates can effectively mitigate this behavior in practice, and diagonalizing the recurrence simplifies both. Our analysis additionally reveals the complexity of learning the angle of complex eigenvalues, which may explain why complex numbers were not found to be useful in most recent state-space model architectures [21, 22] .A side finding of our study is the symbiosis between independent modules, which are here neurons and can be more generally small heads, with adaptive learning rate optimizers in linear recurrent networks. Such a design pattern has promising properties: it facilitates online learning [59] and compositional generalization [60] , allows for high level of parallelization [22] , and matches, at a high level, the modular organization of the cortex in cortical columns [61] . Understanding how to increase the expressivity of small linear modules while keeping their great optimization properties constitutes a promising avenue for future research. 
 
 Appendix: The theory we introduced, with its focus on signal propagation, only addresses the training dynamics of recurrent neural networks. Consequently, it does not provide insights into other important questions such as generalization abilities or memory capacities of these networks. The main assumption underlying this analysis is that the recurrence is both diagonal and linear. While this approach offers valuable insights, it can only approximate signal propagation in more sophisticated architectures. Given the simplicity of the analytical tools employed, there is little hope that this framework can be extended to more general settings without significant modifications.In this section, we rigorously define all the architectures we use in the main text and in the appendix, as well as precisely describe how we initialize them.Let us start by introducing the linear recurrent neural network we are using. It satisfiesh 0 = 0 (12) h t+1 = Ah t + Bx t+1 (13) y t = Ch t + Dx t (14) with x t ∈ R din , h t ∈ R n , y t ∈ R dout , A ∈ R n×n , B ∈ R n×din , C ∈ R dout×n and D ∈ R dout×din .We draw the entries of B, C and D from independent truncated normal distributions with fan_in scaling. We draw each entry of A from N (0, 1/ √ n) independently and then apply the following postprocessing to it: First we complex diagonalize A, which we can do almost surely. Note λ its eigenvalues. We then transform them according toEQUATION)with ν 0 and θ 0 two scalars that we control. This transformation has several benefits: we are guaranteed that the magnitude of λ is within [ν 0 , 1] (and in [ν 0 , ν 0 + (1 -ν 0 ) tanh(1)] in the limit n → ∞ as the eigenvalues of A stay within the unit circle in that limit), and conjugate pairs of eigenvalues remain conjugate. This last point ensures that the resulting matrix remains real without having to change the eigenvectors. When we split the connectivity matrix A into several independent blocks, we initialize each block separately with the scheme described above.Both architectures have recurrence of the form h 0 = 0 (16)EQUATIONwith ⊙ the element-wise product,x t ∈ R din , h t ∈ C n , y t ∈ R dout , λ ∈ C n , B ∈ C n×din , γ ∈ R n , C ∈ C dout×n and D ∈ R dout×din .For the complex-valued linear RNN, we takeEQUATIONso that it can be considered as parametrizing the diagonalized version of the linear RNN.For the LRU [20] , we takeEQUATION)For both architectures, we use the LRU initialization so that λ is uniformly distributed on the ring between the circles of radii ν 0 and 1, with absolute angle restricted to be below θ 0 . For the LRU, we initialize γ to 1 -|λ| 2 .We consider a slightly modified version of S4 here:EQUATIONEQUATIONwithx t ∈ R din , h t ∈ C n , y t ∈ R dout , ω re A + iω im A ∈ C n , B ∈ C n×din , ω ∆ ∈ R n , C ∈ C dout×n and D ∈ R dout×din .The main difference with the standard architecture is that we do not consider any state expansion so that there is no parameter sharing. The makes this architecture closer to the linear RNN architecture we mostly focus on in this paper, while keeping the kind of parametrization it uses. In the same spirit, we do not couple the input and the readout matrices in any way.The initialization we use also differs from existing ones as we initialize ∆ to 1 and exp ∆(ω re A + iω im A ) in the same way as λ in the LRU.The GRU version we use is the following:EQUATIONwith σ the sigmoid function, h t ∈ R n , x t ∈ R din and all the other matrices appropriately sized.We initialize the parameters with the Chrono initialization [29] : all parameters are initialized using standard practice (orthogonal initialization for the weights taking h as input, Lecun initialization for the rest, biases initialized at 0) except for b f , which is initialized asEQUATIONHere, T min and T max denote the minimum and maximal characteristic time scale. In the experiments of Section 5 we takeEQUATIONT max = 2 1 -ν 0 (33) to enable the comparison with other architectures. Indeed, when ignoring the dependence of the forget gate f on the input x and on the hidden state h, this corresponds to having f ∈ ν 0 , 1+ν0 2 .This section introduces all the theoretical results we directly or indirectly mention in the main text, as well as provides a proof for them.Most, if not all the calculations, that we will be doing in this section involves infinite sums. We state and prove two useful lemmas to simplify later calculations. Lemma 1. For α, β ∈ C satisfying |α| < 1 and |β| < 1, and (u n ) n∈Z a bounded sequence satisfying u -n = u n , we have n,m≥0EQUATIONProof. The proof naturally comes from separating the indices n and m in three sets: one in which the two are equals, one in which n is larger and one in which m is larger. This gives n,m≥0EQUATIONEQUATIONEQUATIONEQUATION)Lemma 2. In the same conditions as Lemma 1, we haveEQUATIONProof. This follows from remarking thatEQUATIONand using Lemma 1 to get the final result.We recall the assumptions that we stated in Section 2.2: a) Linear diagonal recurrent neural networks. We restrict ourselves to networks satisfying h t+1 = λ ⊙ h t + x t+1 with λ, h t and x t complex numbers. Without loss of generality, we focus on the one dimensional setting. We additionally consider λs with absolute values smaller than 1.b) Infinite time horizon. We consider infinite sequences and initialize the network dynamics at t 0 = -∞.c) Wide-sense stationarity. We assume the different quantities that the network receives, which includes the inputs x t , to be wise-sense stationary (WSS). A random process X t is said to be WSS if its auto-correlation function is independent of time, that is, for all t ∈ R and ∆ ∈ R, E X t+∆ Xt := R X (∆).Without loss of generality, we can take t = 0 given the wide-sense stationarity and infinite time horizon assumptions. Let us first remark that we haveEQUATIONso thatEQUATION= n,m≥0EQUATION)We used Lemma 1 to obtain the last equality. In Section 2.2, we focused on the real case λ = λ, so this formula becomes Equation 5. If we further assume that the auto-correlation of x decreases exponentially with decay rate ρ, that is R x (∆) = ρ |∆| , we can further simplify the last expression:EQUATION)It follows that if the inputs are i.i.d. (ρ = 0), we haveE[|h 0 | 2 ] = (1 -|λ| 2 ) -1, and if the inputs are constant equal to 1 (ρ = 1), we haveE[|h 0 | 2 ] = |1 -λ| -2 .Differentiating the update h t+1 = λh t + x t+1 with respect to λ givesEQUATIONso thatEQUATIONEQUATION= n,m≥0EQUATIONEQUATIONNote that some extra technicalities are needed to justify these equations as λ and h t are complex valued: these formulas hold as they would in the real-valued case as h t is an holomorphic function of λ.We can now compute the variance of the sensitivity of the hidden state with respect to the parameters.EQUATIONUsing Lemma 2 givesEQUATIONDifferentiating this quantity as a product gives (56) which then simplifies asE dh t dλ 2 =   1 + αβ (1 -αβ) 3   R x (0) + ∆≥1 (α ∆ + β ∆ )R x (∆)   + 0 + α (1 -αβ) 2   ∆≥1 ∆α ∆-1 R x (∆)   + β (1 -αβ) 2   ∆≥1 ∆β ∆-1 R x (∆)     α=λ,β= λ,EQUATION)Note that Equation 6 in the main text is the real-valued version of that formula.Let us now further simplify this equation when R x (∆) = ρ |∆| . If we use this in the differentiated quantity before differentiating it, we getEQUATION)Calculating this quantity manually is painful. Instead, we use the following trick. Its denominator is rather easy to compute, it is equal to (1 -αβ) 3 (1 -ρα) 2 (1 -ρβ) 2 . We thus multiply it to the derivative of the function we want to compute in order to obtain a polynomial with unknown factors, and use polynomial regression tools to derive the resulting coefficients. Massaging the obtained expression to make it easier to compute the closed-form value of this quantity when ρ = 0 and ρ = 1, we getE dh t dλ 2 = (1 -ρ)(1 + |λ| 2 ) + ρ 2 (1 -|λ| 2 ) 3 + ρ(1 -ρ)|λ| 2 (ρ|λ| 2 (1 + |λ| 2 ) -2λ -2 λ) (1 -|λ| 2 ) 3 |1 -ρλ| 4 .(59) This is the quantity we plot on Figure 1 .A, when λ is real-valued. When ρ = 0, this quantity becomesEQUATION)and it is equal toEQUATIONwhen ρ = 1. Additionally, it will diverge whenever |λ| → 1 when ρ < 1, and when λ → 1 when ρ = 1.Regarding the backpropagation of errors to the inputs, the analysis we did in the main text also holds for complex number given that h t is an holomorphic function of x t and it thus behaves as the forward pass once replacing the input distribution with the one of output errors ∂ ht L t .We now turn to the non-diagonal case. For the sake of simplicity, we assume that recurrent matrix is complex diagonalizable and that its eigenvalues are all different. This will enable us to differentiate the eigenvalues and the eigenvectors. We consider dynamics of the form h t+1 = Ah t + x t+1 (62) As A is complex diagonalizable, there exists a complex-valued matrix P and a complex-valued vector λ such that A = P diag(λ)P -1 (63)EQUATION)The linear recurrent neural network considered above is equivalent to its diagonal versionEQUATIONEQUATIONWe now differentiate h t w.r.t. to A using the diagonal parametrization and obtainEQUATION)d A P , d A P -1 and d A λ can be considered constant. Intuitively, the eigenvalues and eigenvectors move smoothly as we restricted ourselves to the case in which eigenvalues are singular. If this is not the case, math becomes trickier as the eigenvectors are not uniquely defined. We can study the behavior of those quantities in more detail, following Boeddeker et al. [62] :EQUATIONEQUATIONThe F introduced in the last equation is equal toEQUATIONImportantly, those two quantities do not grow to infinity as the absolute value of the eigenvalues goes to 1, which means that we can consider those derivatives to be independent of |λ| for the sake of our analysis. Note that the previous argument assumes that eigenvalues do not collapse.d λ h diag tis the dominating term in d A h t . We wonder which of the three different terms that appear in d A h t (Equation 67) will be the dominating one as |λ| (or λ) goes to 1. In the previous paragraph, we have shown that the derivative of P -1 , P and λ can be considered constant for the sake of our analysis. We thus focus on the other terms.First, we have∂h t,l ∂P ij = h diag t,i 1 j=l (71)so the magnitude of this quantity is roughly the one of h diag t , which corresponds to the low pass filtering of the inputs with different λ values.Second, we know that ∂ h diag t h t does not change in magnitude as λ changes, as P remains bounded.So, for the third term of the sum, we are left to study the behavior of d P -1 h diag t . We can show that it evolves according toEQUATION)It follows that the third term in the sum also corresponds to a low pass filtering of the inputs.Finally, we know that the second term, the one in d λ h diag t will grow faster to infinity as it corresponds to two consecutive low pass filters with the same λ values (c.f. calculation above). It will thus be the dominating term in the infinite memory limit.In our analysis, we make the assumption that the different quantities given to the network are widesense stationary, that is their statistics are invariant to a temporal shift. In practice, this assumption will likely never be satisfied, as sequences are finite and as parts of the sequence (e.g., the beginning of a text) can have different statistics than other parts (e.g., the end of a text).It should be noted that the analysis we have done in the wide-sense stationary case can provide an upper bound on the different quantity we study. Indeed, if there exists a function U such thatEQUATIONminor modifications to our analysis enable to upper bound the different quantities we study, replacing R x by U .In our experiments of Section 5, we plug the empirical covariance defined asEQUATIONinto our analytical expressions. It comes with an approximation as it averages the autocorrelation for all positions, which are not equal when wide-sense stationarity is not met.In this section, we consider a diagonal linear recurrent neural network of the formEQUATIONwith γ(λ) the input normalization factor and λ parametrized by a vector ω. Next, we study the effect of input normalization and reparametrization, first in the real-valued setting and then in the complex-valued one.Let us start with the forward pass: asEQUATIONγ rescales the value the hidden state takes. To avoid any explosion behavior, we thus ideally want γ to be the inverse of value of E[(h t ) 2 ] without normalization, which we have computed in Equation 45.The same behavior holds for the backpropagation of errors to the inputs asEQUATION)We now move to the impact of the parametrization. To simplify the calculation, we will ignore the dependency of γ on λ when differentiating it. This can easily be done in automatic differentiation software by removing this dependency from the computational graph with γ(stop_gradient(λ)). We then haveEQUATIONand d λ h t which is rescaled by γ compared to the calculation we did above. As a consequence, both the input normalization and the parametrization can help to mitigate the curse of memory.We now extend the previous analysis to the complex case, and take a polar parametrization of λ: λ(ω) = ν(ω) exp(iθ(ω)). The effect of the input normalization does not change when moving to complex numbers. The role of the reparametrization is however a bit more subtle. To simplify the analysis, we will further assume that E[|d λ h t | 2 ] is only a function of ν. This asumption holds in the case of ρ = 0 and ρ = 1, c.f. Section B.2.2, but not necessarily otherwise. To ensure that that this quantity does not depend on λ, we thus wantd ω ν 2 E(ν) = Θ(1) and ν 2 d ω θ 2 E(ν) = Θ(1).The second means that the angle parametrization must depend on the value ν takes. Let us take the ρ = 0 example to get an idea of what the ideal parametrization should be. First, we haveEQUATION)We are left with the differential equation ν ′ = Θ(1 -ν 2 ), which is for example solved with ν = tanh(ω ν ). Now let us look at the parametrization of θ. If we ignore the ν 2 term for simplicity, the approximate differential equation it needs to solve is d ω θ = Θ(1 -ν 2 ), which can be solved by θ = stop_gradient(1 -ν 2 )ω θ . The exact detail of this calculation do not really matter as this is heavily input distribution dependent. However, the interesting part here is that the angle parameter must be rescaled by a function of ν. This makes intuitive sense when looking looking at the sharpness of the loss around optimality in Figure 1 .C, but this also makes the loss even flatter further away from optimality. We will come back to this point in Section C. 1.4 , showing that in practice, such a parametrization complicates the learning of the θ. Learning complex numbers is thus difficulty, because of the angle.This section is dedicated to detail the theoretical results behind our analysis of the teacher-student task, present all the details necessary to reproduce our empirical experiments, and provide additional analysis.C.1 1D settingIn this toy example, we are interested in learning a simple 1-dimensional linear recurrent neural network which follows the dynamicsEQUATIONto reproduce the hidden state h * t of a teacher with recurrent parameter λ * . Note that we here allow all variables to be complex-valued. We take the loss to beEQUATIONWe assume x to be drawn from a wide-sense stationary distribution so that we can focus on studying the behavior of oneL t (λ, λ * ) := 1 2 E x |h t -h * t |2 to understand the behavior of the full loss L, in the limit of infinitely long sequences (T → ∞). Moreover, to further simplify the calculations, we assume that x is real-valued and that R x (∆) = ρ |∆| .Let us now proceed with the calculation:EQUATIONWe have shown in Section B.2 that in the limit of t → ∞,EQUATIONSimilar derivations hold for the other three terms in the loss. Grouping them gives the exact value of the loss. We omit the formula as it is not particularly insightful. In the case of constant inputs (ρ = 1), we haveL t (λ, λ * ) = 1 2 1 1 -λ - 1 1 -λ * 2 .(90)In the case of i.i.d. inputs (ρ = 0), we haveEQUATION)This is the loss we plot on Figure 1 .B and C.Having a simple closed-form solution for the value the loss takes gives us the possibility to investigate in more detail what an optimal normalization and parametrization should be. We focus on the case ρ = 0.For ρ = 0, the optimal normalization is γ(λ) = 1 -|λ| 2 . Given that we now add an input normalization to the student, we must also add it to the teacher for the student to be able to fit it. The loss becomesEQUATION)Next, we parametrize λ as λ = ν(ω ν ) exp(iθ(ω θ )) and seek to find a parametrization such that, at optimality, E[(d ων h t ) 2 ] = 1 and E[(d ω θ h t ) 2 ] = 1. Given that the student perfectly fit the teachergenerated data at optimality and that the loss we use is the mean-squared error, this corresponds to having d 2 ων L t = 1 and d 2 ω θ L t = 1.Deriving the optimal ν parametrization. We now compute the Hessian of the loss w.r.t. ω ν . First, we can simplify our calculations by restricting ourselves to the case θ = θ * . The loss becomesEQUATION)Differentiating this function a first time, we obtainEQUATION)Differentiating it a second time givesEQUATION)Leveraging the fact thatEQUATION)we finally get, whenEQUATION)Given that we are at optimality, we haveEQUATION)To keep that quantity constant, we thus have to solve the differential equationEQUATIONwhich gives ν(ω ν ) = tanh(ω ν ).Deriving the optimal θ parametrization. We now move to the parametrization of θ. We haveEQUATION)For notational convenience, we denoteEQUATION)We haveEQUATIONandEQUATIONAt optimality (θ = θ * and ν = ν * ), we have α(0) = (1 -ν 2 ) 2 , α ′ (0) = 0 and α ′′ (0) = 2ν 2 , so thatEQUATION)The optimal parametrization thus has to satisfyEQUATION)that isθ(ω θ ) = ω θ 1 -ν 2 ν √ 1 + ν 2 (107)There are two things we can remark:-First, the parametrization that we derived for the general case in Section B.3.2, which additionally ignored the dependence of γ on λ, is relatively accurate. The only difference is the apparition of the extra ν term, which becomes insignificant in the long memory limit ν → 1.-Second, the optimal θ parametrization has to be a function of ν, and thus ω ν , so the differential equation ν needs to satisfy changes. Yet, this considerably simplifies the calculation and there is no simple solution to that problem. One could still argue that the initial choice we made, that is to use a polar parametrization, is the issue. It could be, but most practical models end up using that choice so highlighting the limitations of this choice has important practical consequences.In the rest of this section, we ignore the dependency of θ on ν, and consider the optimal parametrization in this setting to beEQUATIONWe now visualize the effect of input normalization and reparametrization on the loss landscape. We focus on two such reparametrizations:-the one used in the LRU [20, 22] with γ(λ) = √ 1 -λ 2 , ν = exp(-exp(ω exp ν )) and θ = exp(ω exp θ ).-the optimal one we derived in the previous Section (c.f. Equations 108 and 109), which is tailored to this specific setting.We use this one-dimensional teacher-student setting to test whether having a parametrization that avoids exploding behaviors at optimality, such as the one we derived in Section C.1.2, facilitates learning. Figure 6 already hints towards the fact the basin of attraction of the global minima is either extremely narrow or that their number decreases as longer memories are considered, making learning more tedious. Figure 7 confirms it. In this figure, we plot the learning dynamics obtained using the Adam optimizer with a learning rate of 10 -3 for 50k steps, starting from λ 0 = 0.99 exp(iπ/4). We consider three different parametrizations of the angle:EQUATIONThe first one does not reparametrize the angle, the second one is the one used in the LRU and the third one is the optimal one we derived above. We use ν = tanh(ω ν ) to parametrize the magnitude in the three cases. We set λ * to λ * = 0.99 exp(iπ/100). The θ landscape when ν is correct therefore corresponds to the ones plotted in the last two columns of Figure 6 . This example shows that efforts to reduce the sharpness of the loss at optimality, as done in the last parametrization, inevitably make the loss flatter elsewhere and optimization impossible. Each line corresponds to a different parametrization: the first line uses a polar parametrization (λ = ν exp(iθ)), the second line uses the double exponential parametrization used in the LRU (exp) and the third one is the optimal parametrization for that task (tanh). Overall, both reparametrizations enable to control the explosion of the Hessian. However, the size of basins of attraction around optimality, or their number, shrinks as |λ * | goes to 1 for the angle, but not for the absolute value, highlighting how difficult learning the angle can be.In Section 4, we argue that the Hessian at optimality is an important object to understand the learning dynamics in the linear teacher-student task we consider. We here provide some theoretical analysis of its structure in the complex diagonal setting, that is we consider a recurrent network of the formEQUATION)EQUATIONwith λ, b and c complex vectors of size n, with n the number of hidden neurons, and d a scalar. We additionally take the loss to be the mean-square error, which is also the one we use in our numerical experiments. Note that, as in our theoretical analysis of Section 2, we consider infinitely long sequences and wide-sense stationary inputs.Recall that the Hessian of the loss is equal toEQUATION)At optimality, only the first term remains, as ∂ ht L t is 0 for all data points. Given that we have shown earlier, e.g. in Section B.2, that the most sensitive parameters to learn are the recurrent ones λ, we focus on the Hessian with respect to these parameters in the following. Figure 7 : Learning the angle is difficult, even in a simple one-dimensional task. The target λ value is equal to λ * = 0.99 exp(iπ/100) and is plotted in yellow. The black lines correspond to the Adam learning dynamics. A. When the angle is not reparametrized (θ = ω θ ), the loss landscape is extremely sharp in the ω θ direction, but Adam compensates for it. B. When the angle is parametrized exponentially (θ = exp(ω θ )), the loss landscape becomes smoother. However, this only hold when the considered angles are small enough, as the exponential parametrization does not bring extra granularity elsewhere. C. When reparametrizing the angle to reduce the gradient explosion as |λ| → 1, the loss becomes extremely tricky to navigate. The parameters are first attracted to a nearby valley, which is flat on the ω θ direction and only indirectly connected to the global minimum. Such a reparametrization thus hinders optimization far away from optimality. See Section C.1.4 for more detail.Before delving into more specific calculations, we make a few remarks on how to deal the Hessian when having complex-valued parameters. We will mostly leverage the fact that the loss L is realvalued.Before that, we recall a few facts about Wirtinger derivatives:-For f (z) a complex-valued function of z, the Wirtinger derivatives are defined as:EQUATIONEQUATION-We haveEQUATION-Leveraging the fact that L is real-valued so that L = L, we haveEQUATIONand, similarly, d λ dλL = dλd λ L. Second derivatives are symmetric, so we additionally have d λ dλL = dλd λ L ⊤ , which means that the complex Hessian is a Hermitian matrix.Taken all together, this shows that the full complex Hessian, which contains all cross derivatives, has a similar structure to the real case.In this section, we compute the full complex Hessian with respect to the recurrent eigenvalue λ and defer the analysis of reparametrization to the next section.First, let us remark thatEQUATIONso thatEQUATIONWe assumed that we are at optimality so that the network perfectly fits the target trajectories and ∂ yt L t = 0. Additionally, L t is the mean-squared error loss so that ∂ 2 yt L t = Id. It follows thatEQUATIONIn the last equation, we made use of the fact that the parameter λ i only affects the hidden state h t,i and not the others, sod λj h t,i = 0 if i ̸ = j.The previous calculation applied to one sequence, we now take the expectation over the data:EQUATIONNote that we introduced a slight abuse of notation in the previous equation as d λ h t is in general a matrix. However, given that the hidden neurons are independent here due to the diagonal connectivity, it is effectively a vector, and we treat it that way. Let us now compute the expectation, using similar calculation techniques to the one we used in Section B.2:EQUATIONWe can now remark that this quantity is very similar to the one we have encountered in Section B.2, up to the presence of b i b j , and can be simplified using Lemma 2. For conciseness, we note S(λ i , λ j ) the right-hand side of the last equation without the b i b j factor. Putting this result back in the Hessian, we getEQUATIONTo gain further intuition of the behavior of this quantity, we take R x (∆) = ρ |∆| , ρ being a real number. A similar calculation to the one we did in Section B.2 givesEQUATION)Given the complexity of this formula, we visualize the magnitude of S(λ i , λ j ) on Figure 8 . Interestingly, we observe this quantity is large when λ i and λ j are conjugate to each other and inputs are uncorrelated. However, as elements in the input sequence get more correlated (ρ → 1), this effect disappears and |S| increases as one of the two eigenvalue gets closer to 1 in the complex plane. In both cases, the effect gets amplified as the magnitude of the eigenvalue increases. ). This term measures how "similar" eigenvalues are in the Hessian. When ρ = 0, eigenvalues are mostly "similar" when they are conjugate to each other. As ρ increases, this effect decreases and eigenvalues become more "similar" if one of them gets close to 1.We also need to compute dλd λ L to get the full complex Hessian. Similarly to the previous calculation, we first haveEQUATIONIt follows thatEQUATION)Using the symmetry with the complex Hessian matrix, we now have all its components.So far, we have computed the complex Hessian, which is not of direct use as we end up optimizing real numbers in practice. Here, we study the impact of different parametrizations of λ on the Hessian. Given that this parametrization only affects λ and not the other parameters in the network and that we only consider the Hessian at optimality here, computing the Hessian of those parameters reduces to left and right multiplying the Hessian by derivatives of λ and λ with respect to these parameters. For future reference, we introduceEQUATION)withA ij := b i b j c i c j S(λ i , λ j ) and B ij = b i bj c i cj S(λ i , λj ).Real-imaginary parametrization: λ = ω re + ω im . We aim at computing the matrixEQUATIONwhich is the building block to compute the full Hessian. First, let us remark that d ωre,i λ i = 1/2, d ωre,i λi = 1/2, d ωim,i λ i = i/2 and d ωim,i λi = -i/2. It follows thatEQUATIONOnce again we emphasize that the first line only holds as we are at optimality. Similar calculations give the rest of the elements of H RI ij :EQUATION)Given the intuition we gained on the structure of S previously, and the fact that A ij ∝ S(λ i , λ j ) and B ij ∝ S(λ i , λj ), we know that this block will have large components if the two corresponding eigenvalues are conjugate of each other or aligned to each other, or if one of them is close to 1.One other quantity that we can calculate is the trace of the Hessian H RI , which is equal to the sum of its eigenvalues. Note that this does not correspond to the eigenvalues of the full Hessian matrix, as it additionally contains entries for other parameters. Yet it already provides some idea of how large the Hessian will be, as the value of this trace appears in the value of the full trace. We haveEQUATIONEQUATIONwhere we used that S(λ i , λi ) is real-valued in the last line. As a side note, this formula partly justifies why studying the expected squared magnitude of d λ h t in Section 2 makes general sense, asEQUATIONMagnitude-angle parametrization: λ = ν(ω ν ) exp(iθ(ω θ )). The calculations for this parametrization are similar to the previous one, with the following differences:EQUATIONEQUATIONEQUATION)After some calculations we obtaind 2 L dω ν,i dω ν,j = ν ′ (ω ν,i )ν ′ (ω ν,j ) 2Re[e i(θ(ω θ,i )+θ(ω θ,j ) A ij + e i(θ(ω θ,i )-θ(ω θ,j )) B ij ] (152)d 2 L dω ν,i dω θ,j = ν ′ (ω ν,i )ν(ω ν,j )θ ′ (ω θ,j ) 2Im[-e i(θ(ω θ,i )+θ(ω θ,j ) A ij + e i(θ(ω θ,i )-θ(ω θ,j )) B ij ] (153)d 2 L dω θ,i dω ν,j = ν(ω ν,i )θ ′ (ω θ,i )ν ′ (ω ν,j ) 2Im[-e i(θ(ω θ,i )-θ(ω θ,j ) A ij + e i(θ(ω θ,i )-θ(ω θ,j )) B ij ] (154)d 2 L dω θ,i dω θ,j = ν(ω θ,i )θ ′ (ω θ,i )ν(ω θ,j )θ ′ (ω θ,j ) 2Re[-e i(θ(ω θ,i )+θ(ω θ,j ) A ij + e i(θ(ω θ,i )-θ(ω θ,j ))EQUATIONWe use the linear RNN architecture defined in Appendix A.1 as teacher and implement our experiments in JAX [63] , using the default Flax [64] implementation of RNNs and the LRU implementation of Zucchet et al. [65] . Code is available here.We initialize RNNs in the same way we initialized the teacher, and initialize the eigenvalues of the LRU and other complex-valued networks with magnitude in [ν 0 , 1] and angle within [-θ 0 , θ 0 ].Given that we are interested in the optimization properties of the different architectures, we only report training losses and do not perform any cross validation.Here are additional details related to the different figures:-Figure 3 : see Tables 1 and 2 .-Figure 4 and 9: for panels A and B, we use ν 0 = 0.99 and draw A in a slightly different manner to the one described above (we directly draw the eigenvalues and eigenvectors so that we have two pairs of complex eigenvalues). We use automatic differentiation to compute the Hessian. For panels C and D, we use the same setup as described in Table 2 , but keep the learning rate constant over the course of learning. We report the effective learning rate at the end of learning.-Figure 10 : for panels A, B and C, we draw the magnitude and angle of 10 λ values independently, uniformly in [ν 0 , 1+ν0 2 ] and [-θ 0 , θ 0 ]. Importantly, this means that there are no conjugate pairs, which leads to more diagonal Hessian matrices at optimality than in Figure 4 . For panel D, see Table 3 .-Figure 11 : same setup as for Figure 3 .As a rule of thumb, each LRU (or complex-valued diagonal network) experiment takes 3 minutes on a consumer-scale GPU (NVIDIA GeForce RTX 3070) and each RNN experiment takes 10 minutes on a CPU. The scans behind the results reported in the different figures require on the order of few hundreds run each. Including our preliminary explorations, the results we report in this section required 30 days of compute, one third of it on GPUs and two thirds on CPUs.In the main text, we only provide an analysis of the loss landscape for the fully connected linear recurrent neural network and its complex-valued diagonal counterpart. We here complete this result by performing the same analysis for the LRU and S4. Given that S4 involves some form of parameter sharing between the magnitude and the angle of the recurrence complex eigenvalues through ∆, which is required to avoid the explosion of the angle gradients, we are particularly interested in observing whether it fully mitigates or not the gradient explosion effect. The results of Figure 9 .B and D do not reveal such benefits: the loss landscape does not have high curvature on the ω im A direction, but it is moved in the ω ∆ direction. We haven't investigated whether qualitative changes arise when changing the input distribution. We chose the learning rates for the two architectures on preliminary scans and verified that non of the extreme learning rates were optimal in the final scan. For the RNN, we found that initializing with ν 0 = 0 gave better results than initializing with the same distribution the teacher has, so we included this choice in the scan. In the LRU, the exponential parametrization of the magnitude ν = exp(-exp(ω ν )) efficiently mitigates the Hessian explosion but not the one of the angle θ = exp(ω θ ), consistently with the theoretical and empirical evidence we have accumulated so far. In B, ∆ is set to 0.01 in the S4 architecture. Setting it to 1 (not plotted here) leads to an Hessian at optimality that has large entries on all recurrent parameters ω re A , ω im A and ω ∆ , similarly to the behavior we observed for the complex diagonal RNN studied in the main text. For panel D, we initialized ∆ at 1 given that this is the initialization that yielded best performance. We note that we didn't find qualitative changes in this plot when changing the initialization scheme of S4.The goal of this experiment is to better understand how the concentration of eigenvalues λ affect the learning dynamics. For fully connected RNNs, there is no reason to expect a major change in behavior. However, it is different for diagonal RNNs. The theoretical analysis we have done in Section C.2 provides us with the following insights. When the elements in the input sequence are uncorrelated, as it is the case here, the entries in the Hessian corresponding to two different eigenvalues increase if they are aligned or conjugate to each other, and if their magnitude is large. We therefore expect that, as the interval on which the angle of the teacher's eigenvalues shrinks (θ 0 → 0), those eigenvalues will be more likely to be "similar" to each other. This results in large non-diagonal terms, as we confirm in Figure 10 .A, B and C. The LRU suffers less from this problem thanks to its reparametrization, which reduces the overall magnitude of Hessian entries related to the magnitude, and partly the one of angle parameters (when it is a small positive number). As a consequence, the performance between these two architectures increases as θ 0 → 0, as seen on Figure 10 .D. In Figure 3 , we have shown that constraining the connectivity matrix to be block diagonal with blocks of size 2 × 2 lead to a critical boost in performance. Further analysis revealed that this arises as the Hessian becomes more diagonal and Adam can thus better compensate for gradients explosion. Here, we study this behavior in more detail by interpolating between the fully connected case and the block diagonal one. This can be achieved by increasing the number of independent heads from 1 (fully connected case) to 32 (2 × 2 block-diagonal connectivity matrix, as we have 64 hidden neurons). In particular, we are interested in understanding how big heads can be while keeping this performance boost. We plot the final performance, as well as the evolution of the effective learning rate for the A matrix over the course of learning on Figure 11 . We find that slightly bigger heads, until 4 × 4 (which corresponds to 16 heads), yield similar benefits. Additionally, the learning rate analysis reveals that the adaptive learning rates of Adam can more selectively compensate for potential gradient explosion cases as the number of heads increases, allowing for bigger learning rates overall and better performance. × . We detail the experimental setup used in Section 5. We select the first 512 tokens from 1024 random sequences in the Wikipedia dataset [66] and pass them through the BERT [67] tokenizer and embedding layer. This yields a dataset of 1024 examples with length 512 and feature dimension 724. Figure 12 shows the autocorrelation function of these inputs, revealing that the i.i.d. assumption serves (ρ = 0) as a good first order approximation. This validates the relevance of our toy experiments for studying signal propagation in more realistic settings. To refine this approximation, we can include a high correlation term (ρ close to 1).We examine realistic networks comprising 4 blocks with the following structure: a recurrent layer, a non-linearity, a gated linear unit [58, GLU] and a skip connection. By default, we omit normalization layers, but when included, as in Figure 5 .C, we place one normalization layer before the recurrent layer and another one before the GLU. All the layers involved contain 256 neurons. We also incorporate a linear encoder at the beginning of the network and a linear decoder at its end.The loss that we use is a next-token mean-squared error, defined asEQUATIONwhere xt (x 1:t-1 ) represents the prediction of the network. Figure 5 reports the average squared value of the hidden state or the gradient. This average is computed over sequences, but also over all neurons / parameters and over all time steps. We compute gradients using batches of size 8.In Figure 5 we vary ν 0 , which controls the magnitude of the eigenvalues of the recurrent Jacobian. Specifically, we sample those magnitudes in the interval [ν 0 , (1 + ν 0 )/2]. For the complex-valued diagonal RNN and the LRU, we apply the LRU initialization. For the LSTM, we use the Chrono Varying strength of hidden state to gates connections C Figure 14 : GRUs behave like diagonal linear networks. This figure illustrates the evolution of the recurrent Jacobian dh∆ dh0 of a GRU, when provided with the BERT embeddings of a sentence extracted from the Wikipedia dataset. A. In the first row, we take the forget gates to be independent of x and h and we sample them with the Chrono initialization [29] for different intervals. The resulting network is linear and diagonal, similar to what we have studied in the theory. These plots therefore serve as visual reference for comparison with the realistic case. B. The second row shows the same plots as the previous row, except that the gates are now dependent on the inputs x and on the hidden states h. As mentioned in A.4, we initialize all the linear layers taking x as input with LeCun initialization and the layers taking h as input with orthogonal initialization. The recurrent Jacobian evolves similarly than in the constant case, particularly on slowly decaying dimensions (large T values). C. In the row, we aim to break the diagonality of the model by increasing the strength σ h of the hidden state to gate connections, for T ∈ [1, 16] . The case σ h = 0 corresponds to gates that only depend on x, similar to architectures like Mamba or Hawk. The plot with σ h = 1 is the same as the middle one in B. For σ h = 3 and higher, the diagonality progressively breaks and the recurrent Jacobian eventually explodes. We note that these plots were obtained from a single example. Yet, we have found the behavior we report here to be typical of the general behavior.Having established that gated RNNs behave similarly to the linear diagonal RNNs considered in our theoretical investigation, a question arises: How well can our theory describe signal propagation in gated RNNs on realistic data? To address this, we study a simplified version of the GRU network (similar to the one studied by Chen et al. [51] ), which incorporates a realistic gating mechanism:EQUATION)As in the rest of this section, we provide BERT embeddings of sentences from the Wikipedia dataset as inputs.To apply our theory to this architecture, we must address two main challenges:1. The gate f t+1 depends on both h t+1 and x t , making it non-constant. Based on our empirical results from the previous section, we reasonably approximate f t+1 ≈ σ(b f ) =: λ, ignoring this dependency. 2. Our theoretical derivations have overlooked cases where the recurrence strength λ normalizes the input x t . When considered, we detached the normalization factor from the computational graph (as in Section B.3.2). However, we can extend our calculations from Section B.2 to accommodate this setting: With these adjustments in place, we can now empirically test our theoretical predictions. For simplicity, we approximate the auto-correlation function R x (blue line in Figure 12 ) as R x (∆) ≈ 0.332δ ∆=0 + 0.044 × 0.997 ∆ . Figure 15 presents our results, which reveal:EQUATIONEQUATION• An almost perfect match between theory and practice for constant gates, confirming that our sample size is large enough. • A very precise, though not perfect, match for context-dependent gates.• The variance of h t and ∂ht ∂b f shows minimal dependence on λ, indicating that the magnitude of error signals received by b f , W f h , and W f x are largely independent of the time constants encoded in the network. Input and hidden-dependent gates For each of these models / sequences, we measure h 2 512,i and d b f,i h 2 512,i (i being the index of one of the hidden neurons). We report this measurement as a function of the time constant T encoded by the neuron (λ = T /1 + T ). The empirical mean is obtained with a kernel regression with the Gaussian kernel K(a, b) := exp(-(a -b) 2 /100). The theoretical prediction comes from the approach described in Section D.3.We draw each entry of A from N (0, 1/ √ n), complex diagonalize it, and apply the transformation x → ν0 + (1 -ν0)tanh(x) to the absolute values of the eigenvalues.The gradient Lipschitz constant L of the loss equals the maximum Hessian eigenvalue[54]. This quantity sets a bound 2/L for the maximum globally stable learning rate. While convergence might happen in a subspace, it is generally aligned with the top Hessian eigenspace near the solution[55,56]. 


Decision: The paper analyses what they call the "curse of memory" in recurrent neural networks (a term borrowed from prior research), i.e., the fact that the variance of the activations can exponentially diverge with respect to the size of the memory for very long sequences. In a sense, the analysis is dual to the one over gradients, as both exploit long chain of Jacobians over the sequence time steps.

The paper had a very live rebuttal period, at the end of which it moved from a potential reject (3/4/4/5/6), to a scenario where one reviewer is recommending acceptance and all the other reviewers remain non-committal with borderline scores. Some parts of the discussion have been correctly summarized by the authors here on OpenReview. Before proceeding to the decision, I would like to weigh on some remaining concerns:

* Reviewer `zt6Z` is concerned that the analysis only involves the transition matrix of SSMs ($\mathbf{A}$ in common notation), and not a full-fledged model. In addition, they were concerned about the use of an artificial scenario to validate the results. From what I understand, both points were addressed but the reviewer is not convinced they can be integrated easily in a camera-ready version. After going through the paper carefully, I agree on this with the authors that both points are relatively minor modifications to the paper, and I respectfully disagree with the reviewer (who, in any case, voted for a borderline accept).

* Reviewer `uPdc` is concerned about the relation of the paper with the LRU paper by Orvieto et al. In their opinion, the paper should be made more "LRU-centric", and they are also concerned that "*LRUs are not mentioned at all until the experriment*". The authors have argued that this paper complements the LRU paper, since most of the techniques of LRU (e.g., input normalization) are introduced empirically, and they are motivated from a signal propagation point of view here. I note that the authors have been very careful with prior literature (e.g., page 3, "*Martens and Sutskever [40] hypothesized that such a phenomenon could arise in RNNs*"), and LRUs are indeed cited in Section 2.2a as an example of application of the paper. However, the analysis does hold for other types of architectures. Hence, after careful analysis, I side with the authors also on this point.

* Reviewer `Bq7Y` has a few concerns (e.g., typos, discussion of the limitations) which were solved in the rebuttal period. They are now in a position similar to `zt6Z` concerning a potential inability to fit the material for the camera ready. On this, I have the same opinion as for `zt6Z`, since most modifications appear to be small.

* Reviewer `xzKN` is concerned about the assumption of wide sense stationarity made by the authors. This does not invalidate their claim, but they propose a different condition that should be less restrictive. This does not appear to be especially critical with respect to the results, as the authors also argue.

Importantly, all reviewers agree on the novelty of the paper and the potential for further investigation for RNNs, both in the original reviews and during the rebuttal period. I believe this is a strong point in favor of the paper. In addition, my own evaluation is that the paper is very well written, easy to follow (despite being fundamentally a theoretical paper), well organized, and with clear and easy to read figures and experiments that complement the analysis. None of the modifications discussed above would modify this substantially, and I do not see a strong reason for waiting another review cycle for incorporating them. Hence, I suggest acceptance for NeurIPS due to the importance of the topic and the results.

Reviews: review: paper_summary: paper_summary: The paper provides a theoretical and practical analysis on the difficulties of training recurrent neural networks. In particular, given the current rise of interest in leveraging recurrent mechanisms for long sequence processing -- due to novel architectural components and solutions (deep state-space models, linear/diagonal recurrences, exponential parametrizations), identifying the crucial components allowing to preserve long range information is a very important research avenue.
strength_weakness: strength_weakness: I have some concerns with some of the paper initial framing and setting (see the Questions section) and the experimental analysis, which has been carried on very synthetic tasks which I am not sure are capable to represent real learning scenarios, thus possibly hindering the analysis validity.
questions: questions: 1) (Minor) While I understand that recent literature has completely overriden this notion, the term "state-space model" refers [1,2,3] to a very general broad concept which is common in control theory and dynamical systems -- which simply refers to a model equipped with state variables controlled by first-order differential equations or difference equations. Such definition  covers the typical feedback connection of Recurrent Models. Thus, stating in the abstract  and in the main text "state-space models, a subclass of RNNs" seem a bit too strong to me -- even if several recent works highlighted the connections. I would be completely fine with specifying "Deep" SSMs (just to highlight that the authors are referring to modern models) or "Structured" SSMs and to mention the suggested references, to better frame the paper to the readers. 

2) In Eq. 1, the authors consider a model not equipped with an output mapping (that was instead considered in Eq. 10, i.e. y_t = Ch_t + Dx_t) -- the hidden state h_t is directly provided to the loss function.  Why is that the case, and do the considered analysis still hold when considering such mapping (that is common in RNNs and SSMs)?

3) The authors refers to an Infinite time horizon/infinite sequences (lines 98/100) as the basis for their analysis. Practically, how long the learning process necessitates to be to the considered issue to emerge/is this something that also the experimental analysis considers? In such case, does this have relation with recent literature highlighting the difficulties of sequential models when dealing with Infinite length sequences [3, 4]? I would like a comment by the authors on this.
 
4) Line (40): the notation (x_t)_t is not very clear to me.

5) Could the authors give an intuition on lines (113/116)?  Why an higher correlation in the input data should imply an higher variance in the model state (in my understanding, similar inputs imply similar loss function, thus smaller gradients)

6) In line 135, the authors use yet another form for the state update that do not consider the B linear mapping. Why is that the case?

7) The teacher-student task should be better described and justified. Why did the authors chose to tackle such task instead of other synthetic tasks (such as selective copy, induction heads etc. see recent surveys [3] for references on synthetic datasets ). And do the authors believe that the teacher-student task is general enough (i.e., points sampled from a normal distribution) to represent a setting where the state has to learn meaningful long range dependencies? Do the conclusions still hold when dealing with tasks  requiring to store more informative data? I believe that  experiments could benefit from experimentations on more "difficult" settings (selective copy, Long Range Arena) to confirm the theoretical considerations with real learning scenarios. 
 
[1]  Genshiro Kitagawa. A self-organizing state-space model. Journal of the American Statistical Association, pages 1203–1215, 1998. 4

[2]  Anton Schafer and Hans Georg Zimmermann. Recurrent neural networks are universal approximators. In Proceedings of ICANN,  pages 632–640. Springer, 2006.

[3] Tiezzi, Matteo, et al. "State-Space Modeling in Long Sequence Processing: A Survey on Recurrence in the Transformer Era." arXiv preprint arXiv:2406.09062 (2024).

[4] Tiezzi, Matteo, et al. "On the resurgence of recurrent models for long sequences: Survey and research opportunities in the transformer era." arXiv preprint arXiv:2402.08132 (2024).
limitations: limitations: The paper does not discuss limitations.
paper_summary: paper_summary: The paper discusses "the curse of memory" which is the hypersensitivity of hidden states to parameters as the memory increases. This can lead to optimization issues, even if the exploding/vanishing gradient issue is addressed. The authors discuss solutions to this problem: complex diagonalization, normalization, and reparametrization. To this end, the authors derive closed-form solutions for the variance of the hidden state and its sensitivity to the parameters and discuss how the solutions can prevent the variances from diverging. This can give insights into why certain families of RNNs are performant. In particular, they show that both the variance of the hidden state and its derivative wrt to the parameters diverges as the memory goes to infinity (although at different rates). Therefore, normalizing the hidden state and reparametrizing the variable representing the memory's temporal scale can avoid the divergence issue.
strength_weakness: strength_weakness: - The main weakness of the paper is the presentation. It is unclear what the contribution of papers is. At many points in the manuscript, it is unclear if the authors are discussing previous results on LRUs or putting forward new results. It would be good if the authors explicitly spell out their contributions.
- Although the paper is presented as discussing the curse of memory and the solutions to it, it reads more like "why LRUs are performant" The solution discussed: exponential reparametrization, normalization, and diagonalization, are all introduced and employed in LRU architecture [20]. However, the authors only discuss SSMs and gated RNNs as networks that implicitly address "the curse of memory". Meanwhile, LRUs employ the same solution discussed in the paper explicitly. Although LRUs can be thought of as a subclass of SSMs, vanilla SSMs [17] do not use the normalization, diagonalization, and exponential reparametrization discussed in the paper. Therefore, discussing LRUs only in experiments "to represent SSMs due to its simplicity" is not justified. 
- If the solution to the curse of memory is not the contribution of the paper, then the novelty and contribution become limited. As I mentioned earlier, LRUs employ the same solution and sometimes with the same justification, for example, exponential parametrization is used for higher granularity around 1 [20]. 
- To motivate their results, the authors mention that "We answer negatively by showing that gradients can explode as the memory of the network increases, even when the dynamics of the network remains stable." but it is unclear where this is shown. 

- minor comments:
L117: the paragraph on backward pass is confusing and hard to understand.
L139: check wording. Could not understand it.
questions: questions: - What are the main contributions of the paper? Which results and which sections?
- What are other types of normalization and reparametrization, other than the ones employed by LRUs, that can address the curse of memory?
limitations: limitations: NA
paper_summary: paper_summary: The paper explores challenges RNNs face in learning long-term dependencies. While generally attributed to the exploding and vanishing gradients problem (EVGP), the authors reveal that as a network's memory increases, parameter changes cause an explosion of the second moment of hidden states and their gradients, making gradient-based learning highly ill-posed. They study different design choices in recurrent architectures such as element-wise recurrence, gating mechanisms and careful parameterization with respect to the curse of memory. They apply their theoretical findings to successful architectures such as state-space models (SSMs) and LSTMs, offering new insights into why some RNN architectures perform better in gradient-based learning.
strength_weakness: strength_weakness: - To add to EVGP mitigation techniques (l. 61 - 64): There also exist regularization techniques based on dynamical systems theory, see especially [1]
- There are quite some typos across the document, see the list below. I recommend running the paper through some sort of language checker to improve the manuscript from a readability and presentation side.
- The authors repeatedly call $\mathbb{E}[X^2]$ a variance (where $X$ here is a placeholder for any RV, i.e.  $h_t$ and $d_\theta h_t$ in the paper), which is technically not correct. It is the second moment. I’d consider fixing this or explaining why the missing term for the variance, $\mathbb{E}[X]^2$ (first moment squared), can be neglected.

Typos / text-bugs:
- l. 119: “[...] a low a pass filtered [...]
- l. 191: “[...] does not hurts [...]”
- l. 192 “Several RNN architectures implicitly alleviates [...]”
- l. 200 “While such models can can approximate any smooth mappings [...]” 
- l. 219 “anaylsis”
- l. 257 “Both network are trained [...]”
- word missing in l. 283 “[...] that characteristic of the curse of memory, [...]”?
- l. 298 consistently -> consistent

[1] Schmidt et al., Identifying nonlinear dynamical systems with multiple time scales and long-range dependencies (ICML 2021)
questions: questions: - in l. 180 you assume $\gamma$ is independent of $\lambda$ yet in the Diff. eq. in l. 182 you introduce the dependency again ($\gamma(\lambda)$)? Is this connected to what you write in the caption of Fig. 2 where you decouple $\gamma$ from $\lambda$ when differentiating? Can you explain the reason for this in more detail?
- in l. 289 you mention that you probe ADAM’s effective learning rate by providing a vector of ones to the optimizer - does that mean you simply read out ADAM’s (corrected) second-moment estimates and report $\frac{\eta}{\sqrt{v_n} + \epsilon}$, where $\eta$ is the global learning rate and $v_n$ the second moment estimates at iteration $n$? If so, at which iteration/epoch do you query the effective learning rate?
limitations: limitations: While an explicit section on limitations is missing, the authors mention their limitations and specific assumptions throughout the manuscript.
paper_summary: paper_summary: In this work, the authors present an analysis of Recurrent Neural Networks (RNNs), focusing on problems which impede optimization as the length of input sequences and model memory increase. The authors first focus on the well-studied problem of vanishing/exploding gradients, arguing that this problem by itself does not explain the difficulties in optimization. The authors turn to the less studied "curse of memory", and argue that this issue persists for RNNs even when gradients are controlled to not explode or vanish. Finally, the authors study how some of the alternative recurrent architectures such as LSTMs and SSMs mitigate the curse of memory, and provide some suggestions for how optimization can be improved.
strength_weakness: strength_weakness: I have three major issues with the paper in its current form. If issue 1 is addressed, I am willing to increase the score of the presentation from fair to excellent, which would bring the paper to a borderline reject. Issues 2 and 3 are more relevant to the soundness of the work, and I've lowered the score to a 2. If the authors address these, I will consider raising the score; however, I can't go above a 3 without a more substantial discussion of limitations, which I think might not be feasible within the rebuttal period; however, I am still open to adjusting this score depending on the discussion.

1. One of the more glaring issues with this work is with the writing. The authors need to take some time to carefully check for typos, incomplete sentences, and correct grammar. I have listed a few particular issues here:
  * line 30: remains -> remain
  * line 69: You say here that you are showing the gradients can explode, but in the next section say you are controlling the gradients so that they do not explode (lines 83-86). Am I missing something conceptually here that makes this not an inconsistency?
   * line 70: remains -> remain
   * Figure 1: the last sentence of the caption is unclear and seems to be missing some words.
   * line 86: should be "studying this behavior quantitatively", not the other way around
   * line 192: alleviates -> alleviate
   * line 283: characteristic -> are characteristic

if these issues are all fixed and the authors will take some time to carefully edit, I will update my score to reflect this effort. Overall the typos only affect clarity in a few spots.

2. In section 2.2 - you justify assumptions A) and B), but I couldn't find any justification of assumption C). Unless I am missing something, this needs to be provided, as assuming Wide-sense stationarity is not necessarily standard. At the very least, this assumption could be acknowledged as a limitation and room for future work. 

3. There is no real discussion section at all, and there is particularly no discussion of limitations.
questions: questions: See above - why is the assumption of Wide sense stationarity needed, and how is your analysis affected by not adopting the assumption?
limitations: limitations: No limitations are presented or discussed in the paper as far as I could find. 

I see no ethical issues with this work.
paper_summary: paper_summary: This paper studies the vanishing and exploding gradient phenomenon in RNNs. In particular, it answer the question: "is solving those issues really enough to ensure well-behaved loss landscapes?" This paper reveals the importance of the element-wise recurrence design pattern combined with careful parametrizations in mitigating the scale issue of hidden states.
strength_weakness: strength_weakness: 1. The assumption of inputs being wide-sense stationarity seems strong. Does it include the setting such as the Cifar10 dataset? 
2. The argument of input normalization seems hand-waving. 
	1. When $\lambda \to 1$, the forward pass and backward pass are exploding with a different speed: Forward is asymptotically $1/(1-\lambda)$ while backward is (at least) asymptotically $1/(1-\lambda)^2$. 
	2. It should be proper to say normalization can relax one scale issue but it cannot achieve the scaling both at the same time. 
3. The paragraph "On the Importance of Adaptive Learning Rates" addresses an essential topic but lacks direct training curves.
	1. While the local Hessian analysis provided is valuable, it doesn’t assure improved global performance.
	2. It would be more effective and indeed necessary to include training curves comparing results using SGD and adaptive learning rates (such as AdamW). This would clearly illustrate the advantages and need for adaptive optimizers.
4. Related works on parameterizations and orthogonal matrix is limited:
	1. https://proceedings.mlr.press/v80/zhang18g.html
	2. https://arxiv.org/abs/1909.09501: This work proposes a parameterised approach to train orthogonal matrix. 
	3. https://arxiv.org/abs/1901.08428
	4. https://openreview.net/pdf?id=ryxepo0cFX
	5. https://arxiv.org/abs/2311.14495: This work proposes the relationship between parameterization and gradient in RNNs/SSMs. 
	6. https://arxiv.org/pdf/2006.12070
5. I am willing to increase the score if the above issues are addressed.
questions: questions: 1. On Page 4, Equation 5 focuses on the second-order expectation of $h_t$​. It would be helpful to understand the specific reason for choosing a second-order expectation over a first-order one. Generally, implementations do not involve the second moment of hidden states. Why is the first-order expectation not considered in this case? Does first-order suffer from the same problem? 
	1. A similar forward-pass blow-up has been proved in the LRU(Orvieto, 2023) paper, proposition 3.3. What is the difference in terms of assumption and technique between your result and the LRU result? 
2. On Page 5, in the paragraph titled "What About Complex Numbers?", you discuss the concept of polar parametrization. This raises an important question: do all parametrizations encounter the same issues as the polar parametrization? If they do, demonstrating this could effectively highlight the limitations of using complex numbers in a polar setting. If not, it may be worthwhile to find the good parameterization.  Currently, the exploration of polar parametrization seems quite restricted.
3. Your title is "vanishing and exploding gradients are not the end of the story". A question based on the saying in Page 9 Figure 5 Panel C is: 
	1. Does layer normalization and batch normalization resolves the all problems?
	2. If not, I think it's hard to claim "Layer normalization keeps the overall gradient magnitude under control".
limitations: limitations: 1. This paper investigates how the hidden states of recurrent models contribute to learning long-term memory. There is very little research available to explain why the mamba model, which uses both gating and state-space models, is successful.
	1. The discussion on Gated RNNs in lines 205-220 lacks rigorous substantiation, relying primarily on speculative assertions to illustrate the impact of gating mechanisms.
	2. The paper does not clearly establish whether the proposed methodology is effective within a gated framework.

Comments: title: Re: Rebuttal
comment: I want to thank the authors and other reviewers for their active discussion. To clarify my first suggestion, I think the paper is LRU-specific both in theoretical section and emperical section. In the theoretical section, the form of the normalization and reparameteraztion discussed are the exact form used in LRUs. But LRUs are not mentioned at all until the experriment. Therefore, my suggestion was to EITHER
- make the paper *more LRU centric* from the beginning and make the paper about LRUs OR
- make the paper *less LRU centric* by discussing more forms of normalization and reparameterizations and including more architectures in the experiments. 

Regarding the contribution list, I suggest to include it in the paper. In the current form, it is not clear what the contribution of the paper is. Including the list provided in the global rebuttal in the manuscript will help the reader understand the paper contributions.
We thank the reviewer for updating their score.

The bounded input assumption the reviewer is suggesting is an interesting one that we did not consider. Such an approach gives an upper bound on the quantities we study that would correspond to the case $\rho=1$ in our analysis (up to a multiplicative constant of course). Our approach has two main benefits compared to this one:
- It also provides a lower bound on the second-moments considered given that it is an equality
- Our upper bound is more precise when the auto correlation function decreases quickly (i.e. $\rho < 1$, c.f. Fig 1.A), which is the more realistic case (see Figure 1 in our rebuttal PDF for an example)

To get those results, we have to assume stationarity (which is arguably stronger than boundedness), which indeed does not hold in practice. However, given that we empirically found that our results almost perfectly hold when stationarity is not met, we believe that the benefits of this assumption largely outweigh the fact that it will never be exactly achieved in practice. We will include this discussion to the appendix in the next version of the paper. We hope that this argument convinces the reviewer of the interest of the stationarity assumption.
We thank the reviewer for their confidence and for increasing their score.

We join the reviewer in thanking all the reviewers for the active discussion and in helping us sharpen the contributions of our work.
We thank the reviewer for their positive feedback on our work and for increasing their score!
Thank you for the response.
Given the lack of substantial discussion in the original submission and the extent of the rebuttal provided, I still have some concerns regarding the ability to fit this additional material within the revision period. Given the feedback provided to other reviewers and since Neurips does not support revision uploads, I will have to take the authors at their word that their revisions and the inclusion of the material from the rebuttal will not affect the clarity and organization of the revised manuscript.

I will raise my score to a borderline accept in light of this consideration, and wish you luck for the final decision. 

To the Neurips area chairs and board - I would encourage in the future allowing authors to submit revised manuscripts during the rebuttal period so that reviewers can guarantee that issues of clarity and promised changes to the manuscript are actually delivered. 

To the authors and other reviewers, I have rarely seen an active discussion during this period between authors and reviewers and I was happy to see the substantial time and effort put into revising this work on both sides.
I have carefully reviewed both the reviews and the rebuttals.

However, I believe that the section on the wide-sense stationary assumption may not be essential for deriving the key results. Estimating hidden states could potentially be approached deterministically with bounded inputs, which might more effectively highlight the advantages of specific parameterizations, initializations, or normalizations.

While the inclusion of the stationary property facilitates a discussion of the second-order moment, this approach seems somewhat forced. The bound on the second-order moment within a random process context does not provide a significantly more valuable insight. Although generalizing to more diverse input distributions is a noteworthy contribution, it is crucial to clearly define the boundaries of this generalization. Without specifying the limits of input expansion, the contribution may not be fully strengthened. It remains unclear how this paper’s contribution stands apart from the LRU paper and what kinds of stronger understanding can be derived for complex model such as gated.

Therefore, I will increase my score from 4 to 5 and encourage the authors to address the issues related to the stationary assumption and to clarify the improvements over the LRU paper in the next version.
I thank the authors for the clarifications, elaborate rebuttal and additional experiments in the general response. I think this work is of great value for the community. Hence, I increased my score to 7.
We here provide additional details on the B, C and D mappings, in the hope that they will clarify the remaining concerns of the reviewer. These mappings operate token-wise whereas the recurrence mixes information between tokens. This fundamental difference justifies why B, C and D should be considered separately from the recurrence and similarly to the rest of the feedforward layers of the network. On the more technical side, integrating B and C in our theory only requires changing the quantities we give as input to our theory, e.g. we should consider the statistics of Bx instead of x ($C\delta$ instead of $\delta$ for the backward pass through the recurrence). Our theory does not aim to capture the role of $D$ given that it side steps the recurrence. To summarize, integrating these mapping does not affect our conclusions (Sections 4 and 5 confirm it empirically).

We believe that the concerns previously raised by the reviewer (as well as the ones mentioned by other reviewers) only require minor modifications to the corresponding sections of the manuscript (see also our answers to reviewers uPdc and Bq7Y). We would appreciate additional input from the reviewer to help us better identify which parts require substantial modifications before the next version of the paper.
I acknowledge that I have read the reviews and rebuttals and I thank the authors for taking their time to answer my questions, clarifiying some aspects of the work. I still have some concerns regarding questions 2 and 6 (the assumptions by the authors on analyzing simplified versions of the models -- do the results still hold when considering models equipped with B,C,D?), which I believe are important and require a more in depth justification.  I agree with other reviewers that the work is promising but requires some additional work to better organize the contribution (i.e., regarding my points: the aforementioned setting simplification on B,C,D mapping and the teacher-student experimental setting should be better framed in the context of the paper). Regardless, I will keep my positive score.
We thank the reviewer for their time and valuable input. Based on our understanding, the reviewer is suggesting two major revisions: first, make the paper less LRU-centric, and second, use the contribution statement from our rebuttal to reorganize the paper. 

Regarding the first point, we agree that this criticism applies to the empirical part of the submitted paper, but we believe that Figure 3 from the PDF already addresses this concern. We additionally believe that focusing on one architecture in most of the paper, while discussing how our findings apply to other architecture (both theoretically and empirically), helps us in keeping the exposition simple. 

As for the second point, our rebuttal’s contribution list follows the exact same structure the paper. 

As a consequence, we are uncertain how to structure the paper differently. Could the reviewer provide more specific feedback on what we could improve or clarify what additional experiments/results they would like to see in our paper?
We appreciate that the reviewer took the time to read our rebuttal and thank them for updating their score.

However, we have to admit we are not quite sure why the reviewer still believes the paper requires substantial revision. From their review, we had understood that clarifying our contributions (as we did in the rebuttal) at the end of the introduction, adding a discussion of the wide-sense stationarity in Section 2 and the discussion of the limitations of our work in a new paragraph of the conclusion would sufficiently address the concerns raised by the reviewer. Can the reviewer help us understand their viewpoint by clarifying what is missing in our rebuttal? We appreciate further feedback so that we can continue improving our paper.
I thank the authors for addressing my comments. After reviewing the authors' response and the remaining concerns of other reviewers I believe this work would benefit from substantial revision and resubmission at a later date. I have changed by score to a borderline reject per my initial review as the authors have adequately acknowledge the particular typos mentioned here. 

The discussion of limitations provided in the rebuttal is interesting, and I agree with the feedback of other reviewers that this work is valuable; however, I think a more substantial revision which works a discussion of limitation and highlights the more particular contributions of this work would serve as a stronger submission in the future.
I thank the authors for taking the time to respond to my comments and other reviewers comments. I will comment on the issues raised by myself here.

 I appreciate the importance of the topic and the fact that the authors study the RNNs and LRUs from a more theoretical standpoint. My general concern is presentation and I believe that the paper will greatly benefit from major reorganization. 

While the general solution to the curse of the memory is discussed, the specific solutions, i.e., the specific form of reparameterrization and specific form of normalization are the forms employed by LRUs. But the LRUs are not really discussed, or mentioned, in the preceding sections. I would suggest to either make the paper “a theoretical study of LRUs” or discuss more specific solutions, e.g., S4, and make it more general. The list of the contributions stated in the global rebuttal is a good starting point. I encourage the authors to resubmit the manuscript after reorganization and revision.  Therefore, I will keep my rating as “Technically solid paper where reasons to reject outweigh reasons to accept.”
Dear Reviewers,

The authors have provided a comprehensive rebuttal. I encourage all of you to respond before the deadline, to have sufficient time to discuss outstanding issues. Please change your scores accordingly, with a small motivation based on the points that were addressed or not (especially concerning the relation with the LRU paper, the novelty of the analysis, and the technical assumptions underneath the methods).

Thanks, the AC