Year: 2021

Abstract:
We study the approximation properties and optimization dynamics of recurrent neural networks (RNNs) when applied to learn input-output relationships in temporal data. We consider the simple but representative setting of using continuoustime linear RNNs to learn from data generated by linear relationships. Mathematically, the latter can be understood as a sequence of linear functionals. We prove a universal approximation theorem of such linear functionals and characterize the approximation rate. Moreover, we perform a fine-grained dynamical analysis of training linear RNNs by gradient methods. A unifying theme uncovered is the non-trivial effect of memory, a notion that can be made precise in our framework, on both approximation and optimization: when there is long-term memory in the target, it takes a large number of neurons to approximate it. Moreover, the training process will suffer from slow downs. In particular, both of these effects become exponentially more pronounced with increasing memory a phenomenon we call the “curse of memory”. These analyses represent a basic step towards a concrete mathematical understanding of new phenomenons that may arise in learning temporal relationships using recurrent architectures.

heading: 1 INTRODUCTION
text: Recurrent neural networks (RNNs) (Rumelhart et al., 1986) are among the most frequently employed methods to build machine learning models on temporal data. Despite its ubiquitous application (Baldi et al., 1999; Graves & Schmidhuber, 2009; Graves, 2013; Graves et al., 2013; Graves & Jaitly, 2014; Gregor et al., 2015), some fundamental theoretical questions remain to be answered. These come in several flavors. First, one may pose the approximation problem, which asks what kind of temporal input-output relationships can RNNs model to an arbitrary precision. Second, one may also consider the optimization problem, which concerns the dynamics of training (say, by gradient descent) the RNN. While such questions can be posed for any machine learning model, the crux of the problem for RNNs is how the recurrent structure of the model and the dynamical nature of the data shape the answers to these problems. For example, it is often observed that when there are long-term dependencies in the data (Bengio et al., 1994; Hochreiter et al., 2001), RNNs may encounter problems in learning, but such statements have rarely been put on precise footing.
In this paper, we make a step in this direction by studying the approximation and optimization properties of RNNs. Compared with the static feed-forward setting, the key distinguishing feature †Equal contribution ‡Corresponding author
here is the presence of temporal dynamics in terms of both the recurrent architectures in the model and the dynamical structures in the data. Hence, to understand the influence of dynamics on learning is of fundamental importance. As is often the case, the key effects of dynamics can already be revealed in the simplest linear setting. For this reason, we will focus our analysis on linear RNNs, i.e. those with linear activations. Further, we will employ a continuous-time analysis initially studied in the context of feed-forward architectures (E, 2017; Haber & Ruthotto, 2017; Li et al., 2017) and recently in recurrent settings (Ceni et al., 2019; Chang et al., 2019; Lim, 2020; Sherstinsky, 2018; Niu et al., 2019; Herrera et al., 2020; Rubanova et al., 2019) and idealize the RNN as a continuous-time dynamical system. This allows us to phrase the problems under investigation in convenient analytical settings that accentuates the effect of dynamics. In this case, the RNNs serve to approximate relationships represented by sequences of linear functionals. On first look the setting appears to be simple, but we show that it yields representative results that underlie key differences in the dynamical setting as opposed to static supervised learning problems. In fact, we show that memory, which can be made precise by the decay rates of the target linear functionals, can affect both approximation rates and optimization dynamics in a non-trivial way.
Our main results are:
1. We give a systematic analysis of the approximation of linear functionals by continuoustime linear RNNs, including a precise characterization of the approximation rates in terms of regularity and memory of the target functional.
2. We give a fine-grained analysis of the optimization dynamics when training linear RNNs, and show that the training efficiency is adversely affected by the presence of long-term memory.
These results together paint a comprehensive picture of the interaction of learning and dynamics, and makes concrete the heuristic observations that the presence of long-term memory affects RNN learning in a negative manner (Bengio et al., 1994; Hochreiter et al., 2001). In particular, mirroring the classical curse of dimensionality (Bellman, 1957), we introduce the concept of the curse of memory that captures the new phenomena that arises from learning temporal relationships: when there is long-term memory in the data, one requires an exponentially large number of neurons for approximation, and the learning dynamics suffers from exponential slow downs. These results form a basic step towards a mathematical understanding of the recurrent structure and its effects on learning from temporal data.

heading: 2 RELATED WORK
text: We will discuss related work on RNNs on three fronts concerning the central results in this paper, namely approximation theory, optimization analysis and the role of memory in learning. A number of universal approximation results for RNNs have been obtained in discrete (Matthews, 1993; Doya, 1993; Schäfer & Zimmermann, 2006; 2007) and continuous time (Funahashi & Nakamura, 1993; Chow & Xiao-Dong Li, 2000; Li et al., 2005; Maass et al., 2007; Nakamura & Nakagawa, 2009). Most of these focus on the case where the target relationship is generated from a hidden dynamical system in the form of difference or differential equations. The formulation of functional approximation here is more general, albeit our results are currently limited to the linear setting. Nevertheless, this is already sufficient to reveal new phenomena involving the interaction of learning and dynamics. This will be especially apparent when we discuss approximation rates and optimization dynamics. We also note that the functional/operator approximation using neural networks has been explored in Chen & Chen (1993); Tianping Chen & Hong Chen (1995); Lu et al. (2019) for nonrecurrent structures and reservoir systems for which approximation results similar to random feature models are derived (Gonon et al., 2020). The main difference here is that we explicitly study the effect of memory in target functionals on learning using recurrent structures.
On the optimization side, there are a number of recent results concerning the training of RNNs using gradient methods, and they are mostly positive in the sense that trainability is proved under specific settings. These include recovering linear dynamics (Hardt et al., 2018) or training in overparameterized settings (Allen-Zhu et al., 2019). Here, our result concerns the general setting of learning linear functionals that need not come from some underlying differential/difference equations, and is also away from the over-parameterized regime. In our case, we discover on the contrary
that training can become very difficult even in the linear case, and this can be understood in a quantitative way, in relation to long-term memory in the target functionals.
This points to the practical literature in relation to memory and learning. The dynamical analysis here puts the ubiquitous but heuristic observations - that long-term memory negatively impacts training efficiency (Bengio et al., 1994; Hochreiter et al., 2001) - on concrete theoretical footing, at least in idealized settings. This may serve to justify or improve current heuristic methods (Tseng et al., 2016; Dieng et al., 2017; Trinh et al., 2018) developed in applications to deal with the difficulty in training with long-term memory. At the same time, we also complement general results on “vanishing and explosion of gradients” (Pascanu et al., 2013; Hanin & Rolnick, 2018; Hanin, 2018) that are typically restricted to initialization settings with more precise characterizations in the dynamical regime during training.
The long range dependency within temporal data has been studied for a long time in the time series literature, although its effect on learning input-output relationships is rarely covered. For example, the Hurst exponent (Hurst, 1951) is often used as a measure of long-term memory in time series, e.g. fractional Brownian motion (Mandelbrot & Ness, 1968). In contrast with the setting in this paper where memory involves the dependence of the output time series on the input, the Hurst exponent measures temporal variations and dependence within the input time series itself. Much of the time series literature investigates statistical properties and estimation methods of data with long range dependence (Samorodnitsky, 2006; Taqqu et al., 1995; Beran, 1992; Doukhan et al., 2003). One can also combine these classic statistical methodologies with the RNN-like architectures to design hybrid models with various applications (Loukas & Oke, 2007; Diaconescu, 2008; Mohan & Gaitonde, 2018; Bukhari et al., 2020).

heading: 3 PROBLEM FORMULATION
text: The basic problem of supervised learning on time series data is to learn a mapping from an input temporal sequence to an output sequence. Formally, one can think of the output at each time as being produced from the input via an unknown function that depends on the entire input sequence, or at least up to the time at which the prediction is made. In the discrete-time case, one can write the data generation process
yk = Hk(x0, . . . , xk−1), k = 0, 1, . . . (1)
where xk, yk denote respectively the input data and output response, and {Hk : k = 0, 1, . . . } is a sequence of ground truth functions of increasing input dimension accounting for temporal evolution. The goal of supervised learning is to learn an approximation of the sequence of functions {Hk} given observation data.
Recurrent neural networks (RNN) (Rumelhart et al., 1986) gives a natural way to parameterize such a sequence of functions. In the simplest case, the one-layer RNN is given by
hk+1 = σ(Whk + Uxk), ŷk = c >hk. (2)
Here, {hk} are the hidden/latent states and its evolution is governed by a recursive application of a feed-forward layer with activation σ, and ŷk is called the observation or readout. We omit the bias term here and only consider a linear readout or output layer. For each time step k, the mapping {x0, . . . , xk−1} 7→ ŷk parameterizes a function Ĥk(·) through adjustable parameters (c,W,U). Hence, for a particular choice of these parameters, a sequence of functions {Ĥk} is constructed at the same time. To understand the working principles of RNNs, we need to characterize how {Ĥk} approximates {Hk}. The model (2) is not easy to analyze due to its discrete iterative nature. Hence, here we employ a continuous-time idealization that replaces the time-step index k by a continuous time parameter t. This allows us to employ a large variety of continuum analysis tools to gain insights to the learning problem. Let us now introduce this framework.
Continuous-time formulation. Consider a sequence of inputs indexed by a real-valued variable t ∈ R instead of a discrete variable k considered previously. Concretely, we consider the input space
X = C0(R,Rd), (3)
which is the linear space of continuous functions from R (time) to Rd that vanishes at infinity. Here d is the dimension of each point in the time series. We denote an element in X by x := {xt ∈ Rd : t ∈ R} and equip X with the supremum norm ‖x‖X := supt∈R ‖xt‖∞. For the space of outputs we will take a scalar time series, i.e. the space of bounded continuous functions from R to R:
Y = Cb(R,R). (4) This is due to the fact that vector-valued outputs can be handled by considering each output separately. In continuous time, the target relationship (ground truth) to be learned is
yt = Ht(x), t ∈ R (5) where for each t ∈ R, Ht is a functional Ht : X → R. Correspondingly, we define a continuous version of (2) as a hypothesis space to model continuous-time functionals
d dt ht = σ(Wht + Uxt), ŷt = c >ht, (6)
whose Euler discretization corresponds to a discrete-time residual RNN. The dynamics then naturally defines a sequences of functionals {Ĥt(x) = ŷt : t ∈ R}, which can be used to approximate the target functionals {Ht} via adjusting (c,W,U).
Linear RNNs in continuous time. In this paper we mainly investigate the approximation and optimization property of linear RNNs, which already reveals the essential effect of dynamics. The linear RNN obeys (6) with σ being the identity map. Notice that in the theoretical setup, the initial time of the system goes back to −∞ with limt→−∞ xt = 0,∀x ∈ X , thus by linearity (Ht(0) = 0) we specify the initial condition of the hidden state h−∞ = 0 for consistency. In this case, (6) has the following solution
ŷt = ∫ ∞ 0 c>eWsUxt−sds. (7)
Since we will investigate uniform approximations over large time intervals, we will consider stable RNNs, where W ∈ Wm with Wm = {W ∈ Rm×m : eigenvalues of W have negative real parts}. (8) Owing to the representation of solutions in (7), the linear RNN defines a family of functionals
Ĥ := ∪m≥1Ĥm, Ĥm := { {Ĥt(x), t ∈ R} : Ĥt(x) = ∫ ∞ 0 c>eWsUxt−sds,W ∈ Wm, U ∈ Rm×d, c ∈ Rm } .
(9) Here, m is the width of the network and controls the complexity of the hypothesis space. Clearly, the family of functionals the RNN can represent is not arbitrary, and must possess some structure. Let us now introduce some definitions of functionals that makes these structures precise. Definition 3.1. Let {Ht : t ∈ R} be a sequence of functionals.
1. Ht is causal if it does not depend on future values of the input: for every pair of x,x′ ∈ X such that xs = x′s for all s ≤ t, we have Ht(x) = Ht(x′).
2. Ht is linear and continuous if Ht(λx + λ′x′) = λHt(x) + λ′Ht(x′) for any x,x′ ∈ X and λ, λ′ ∈ R, and supx∈X ,‖x‖X≤1 |Ht(x)| <∞, in which case the induced norm can be defined as ‖Ht‖ := supx∈X ,‖x‖X≤1 |Ht(x)|.
3. Ht is regular if for any sequence {x(n) ∈ X : n ∈ N} such that x(n)t → 0 for Lebesgue almost every t ∈ R, limn→∞Ht(x(n)) = 0.
4. {Ht : t ∈ R} is time-homogeneous if Ht(x) = Ht+τ (x(τ)) for any t, τ ∈ R, where x(τ)s = xs−τ for all s ∈ R, i.e. x(τ) is x whose time index is shifted to the right by τ .
Linear, continuous and causal functionals are common definitions. One can think of regular functionals as those that are not determined by values of the inputs on an arbitrarily small time interval, e.g. an infinitely thin spike input. Time-homogeneous functionals, on the other hand, are those where there is no special reference point in time: if the time index of both the input sequence and the functional are shifted in coordination, the output value remains the same. Given these definitions, the following observation can be verified directly and its proof is immediate and hence omitted.
Proposition 3.1. Let {Ĥt : t ∈ R} be a sequence of functionals in the RNN hypothesis space Ĥ (see (9)). Then for each t ∈ R, Ĥt is a causal, continuous, linear and regular functional. Moreover, the sequence of functionals {Ĥt : t ∈ R} is time-homogeneous.

heading: 4 APPROXIMATION THEORY
text: The most basic approximation problem for RNN is as follows: given some sequence of target functionals {Ht : t ∈ R} satisfying appropriate conditions, does there always exist a sequence of RNN functionals {Ĥt : t ∈ R} in Ĥ such that Ht ≈ Ĥt for all t ∈ R? We now make an important remark with respect to the current problem formulation that differs from previous investigations in the RNN approximation: we are not assuming that the target functionals {Ht : t ∈ R} are themselves generated from an underlying dynamical system of the form
Ht(x) = yt where ḣt = f(ht, xt), yt = g(ht) (10)
for any linear or nonlinear functions f, g. This differs from previous work where it is assumed that the sequence of target functionals are indeed generated from such a system. In that case, the approximation problem reduces to that of the functions f, g, and the obtained results often resemble those in feed-forward networks.
In our case, however, we consider general input-output relationships related by temporal sequences of functionals, with no necessary recourse to the mechanism from which these relationships are generated. This is more flexible and natural, since in applications it is often not clear how one can describe the data generation process. Moreover, notice that in the linear case, if the target functionals {Ht} are generated from a linear ODE system, then the approximation question is trivial: as long as the dimension of ht in the approximating RNN is greater than or equal to that which generates the data, we must have perfect approximation. However, we will see that in the more general case here, this question becomes much more interesting, even in the linear regime. In fact, we now prove precise approximation theories and characterize approximation rates that reveal intricate connections with memory effects, which may be otherwise obscured if one considers more limited settings.
Our first main result is a converse of Proposition 3.1 in the form of an universal approximation theorem for certain classes of linear functionals. The proof is found in Appendix A. Theorem 4.1 (Universal approximation of linear functionals). Let {Ht : t ∈ R} be a family of continuous, linear, causal, regular and time-homogeneous functionals on X . Then, for any > 0 there exists {Ĥt : t ∈ R} ∈ Ĥ such that
sup t∈R ‖Ht − Ĥt‖ ≡ sup t∈R sup ‖x‖X≤1 |Ht(x)− Ĥt(x)| ≤ . (11)
The proof relies on the classical Riesz-Markov-Kakutani representation theorem, which says that each linear functional Ht can be uniquely associated with a signed measure µt such that Ht(x) =∫ R x > s dµt(s). Owing to the assumptions of Theorem 4.1, we can further show that the sequence of representations {µt} are related to an integrable function ρ : [0,∞) → Rd such that {Ht} admits the common representation
Ht(x) = ∫ ∞ 0 x>t−sρ(s)ds, t ∈ R, x ∈ X . (12)
Comparing this representation with the solution (7) of the continuous RNN, we find that the approximation property of the linear RNNs is closely related to how well ρ(t) can be approximated by the exponential sums of the form (c>eWtU)>. Intuitively, (12) says that each output yt = Ht(x) is simply a convolution between the input signal and the kernel ρ. Thus, the smoothness and decay of the input-output relationship is characterized by the convolution kernel ρ. Due to this observation, we will hereafter refer to {Ht} and ρ interchangeably.
Approximation rates. While the previous result establishes the universal approximation property of linear RNNs for suitable classes of linear functionals, it does not reveal to us which functionals can be efficiently approximated. In the practical literature, it is often observed that when there is
some long-term memory in the inputs and outputs, the RNN becomes quite ill-behaved (Bengio et al., 1994; Hochreiter et al., 2001). It is the purpose of this section to establish results which make these heuristics statements precise. In particular, we will show that the rate at which linear functionals can be approximated by linear RNNs depends on the former’s smoothness and memory properties. We note that this is a much less explored area in the approximation theory of RNNs.
To characterize smoothness and memory of linear functionals, we may pass to investigating the properties of their actions on constant input signals. Concretely, let us denote by ei (i = 1, . . . , d) the standard basis vector in Rd, and ei denotes a constant signal with ei,t = ei1{t≥0}. Then, smoothness and memory is characterized by the regularity and decay rate of the maps t 7→ Ht(ei), i = 1, . . . , d, respectively. Our second main result shows that these two properties are intimately tied with the approximation rate. The proof is found in Appendix B. Theorem 4.2 (Approximation rates of linear RNN). Assume the same conditions as in Theorem 4.1. Consider the output of constant signals yi(t) := Ht(ei), i = 1, . . . , d. Suppose there exist constants α ∈ N+, β, γ > 0 such that yi(t) ∈ C(α+1)(R), i = 1, . . . , d, and
eβty (k) i (t) = o(1) as t→ +∞, and sup t≥0 β−k|eβty(k)i (t)| ≤ γ, k = 1, . . . , α+ 1, (13)
where y(k)i (t) denotes the k th derivative of yi(t). Then, there exists a universal constant C(α) only depending on α, such that for any m ∈ N+, there exists a sequence of width-m RNN functionals {Ĥt : t ∈ R} ∈ Ĥm such that
sup t∈R ‖Ht − Ĥt‖ ≡ sup t∈R sup ‖x‖X≤1
|Ht(x)− Ĥt(x)| ≤ C(α)γd
βmα . (14)
The curse of memory in approximation. For approximation of non-linear functions using linear combinations of basis functions, one often suffers from the “curse of dimensionality” (Bellman, 1957), in that the number of basis functions required to achieve a certain approximation accuracy increases exponentially when the dimension of input space d increases. In the case of Theorem 4.2, the bound scales linearly with d. This is because the target functional possesses a linear structure, and hence each dimension can be approximated independently of others, resulting in an additive error estimate. Nevertheless, due to the presence of the temporal dimension, there enters another type of challenge, which we coin the curse of memory. Let us now discuss this point in detail.
The key observation is that the rate result requires exponential decay of derivatives ofHt(ei), but the density result (Theorem 4.1) makes no such assumption. The natural question is thus, what happens when no exponential decay is present? We assume d = 1 and consider an example in which the target functional’s representation satisfies ρ(t) ∈ C(1)(R) and ρ(t) ∼ t−(1+ω) as t → +∞. Here ω > 0 indicates the decay rate of the memory effects in our target functional family. The smaller its value, the slower the decay and the longer the system memory. For any ω > 0, the system’s memory vanishes more slowly than any exponential decay. Notice that y(1)(t) = ρ(t) and in this case there exists no β > 0 making (13) true, and no rate estimate can be deduced from it.
A natural way to circumvent this obstacle is to introduce a truncation in time. With T ( 1) we can define ρ̃(t) ∈ C(1)(R) such that ρ̃(t) ≡ ρ(t) for t ≤ T , ρ̃(t) ≡ 0 for t ≥ T + 1, and ρ̃(t) is monotonically decreasing for T ≤ t ≤ T + 1. With the auxiliary linear functional H̃t(x) :=∫ t
0 xt−sρ̃(s)ds, we can have an error estimate (with technical details provided in Appendix B)
sup t∈R ‖Ht − Ĥt‖ ≤ sup t∈R ‖Ht − H̃t‖+ sup t∈R ‖H̃t − Ĥt‖ ≤ C
( T−ω + ω
m T 1−ω
) . (15)
In order to achieve an error tolerance , according to the first term above we require T ∼ − 1ω , and then according to the second term we have
m = O ( ωT 1−ω/ ) = O ( ω −1/ω ) . (16)
This estimate gives us a quantitative relationship between the degree of freedom needed and the decay speed. With ρ(t) ∼ t−(1+ω), the system has long memory when ω is small. Denote the minimum number of terms needed to achieve an L1 error as m(ω, ). The above estimate shows
an upper-bound of m(ω, ) goes to infinity exponentially fast as ω → 0+ with fixed . This is akin to the curse of dimensionality, but this time on memory, which manifests itself even in the simplest linear settings. A stronger result would be that the lower bound of m(ω, )→∞ exponentially fast as ω → 0+ with fixed , and this is a point of future work. Note that this kind of estimates differ from the previous results in the literature (Kammler, 1979b; Braess & Hackbusch, 2005) regarding the order of m(ω, ) ∼ log(1/ ) as → 0 with fixed ω = 1 or 2 in the L∞ or L1 sense. Note that the L1 result has not been proved.

heading: 5 FINE-GRAINED ANALYSIS OF OPTIMIZATION DYNAMICS
text: According to Section 4, memory plays an important role in determining the approximation rates. The result therein only depends on the model architecture, and does not concern the actual training dynamics. In this section, we perform a fine-grained analysis on the latter, which again reveals an interesting interaction between memory and learning.
The loss function (for training) is defined as
ExJT (x; c,W,U) := Ex|ĤT (x)−HT (x)|2 = Ex ∣∣∣∣∣ ∫ T
0
[c>eWtU − ρ(t)>]xT−tdt ∣∣∣∣∣ 2 . (17)
Without loss of generality, here the input time series x is assumed to be finitely cut off at zero, i.e. xt = 0 for any t ≤ 0 almost surely. Training the RNN amounts to optimizing ExJT with respect to the parameters (c,W,U). The most commonly applied method is gradient descent (GD) or its stochastic variants (say SGD), which updates the parameters in the steepest descent direction.
We first show that the training dynamics of ExJT exhibits very different behaviors depending on the form of target functionals. Take d = 1 and consider learning different target functionals with white noise x. We first investigate two choices for ρ: a simple decaying exponential sum and a scaled Airy function. The Airy function target is defined as ρ(t) = Ai(s0[t− t0]), where Ai(t) is the Airy function of the first kind, given by the improper integral Ai(t) = 1π limξ→∞ ∫ ξ 0 cos ( u3 3 + tu ) du. Note that the effective rate of decay is controlled by the parameter t0: for t ≤ t0, the Airy function is oscillatory. Hence for large t0, a large amount of memory is present in the target.
Observe from Figure 1 that the training proceeds efficiently for the exponential sum case. However, in the second Airy function case, there are interesting “plateauing” behaviors in the training loss, where the loss decrement slows down significantly after some time in training. The plateau is sustained for a long time before an eventual reduction is observed.
As a further demonstration of that this behavior may be generic, we also consider a nonlinear forced dynamical system, the Lorenz 96 system (Lorenz, 1996), where the similar plateauing behavior is observed even for a non-linear RNN model trained with the Adam optimizer (Kingma & Ba, 2015). All experimental details are found in Appendix C.3.1.
The results in Figure 1 hint at the fact that there are certainly some functionals that are much harder to learn than others, and it is the purpose of the remaining analyses to understand precisely when and why such difficulties occur. In particular, we will again relate this to the memory effects in the target functional, which shows yet another facet of the curse of memory when it comes to optimization.
Dynamical analysis. To make analysis amenable, we will make a series of simplifications to the loss function (17), by assuming that x is white noise, d = 1, T →∞, and the recurrent kernel W is diagonal. This allows us to write (see Appendix C.1 for details) the optimization problem as
min a∈Rm,w∈Rm+ J(a,w) := ∫ ∞ 0 ( m∑ i=1 aie −wit − ρ(t) )2 dt. (18)
We will subsequently see that these simplifications do not lose the key features of the training dynamics, such as the plateauing behavior. We start with some informal discussion on a probable reason behind the plateauing. A straightforward computation shows that, for k = 1, 2, . . . ,m,
∂J
∂wk (a,w) = 2ak ∫ ∞ 0 (−t)e−wkt ( m∑ i=1 aie −wit − ρ(t) ) dt. (19)
A similar expression holds for ∂J∂ak . Write the (simplified) linear functional representation for linear RNNs as ρ̂(t; a,w) := ∑m i=1 aie
−wit, which serves to learn the target ρ. Observe that plateauing under the GD dynamics occurs if the gradient ∇J is small but the loss J is large. A sufficient condition is that the residual ρ̂(t; a,w) − ρ(t) is large only for large t (meaning the exponential multiplier to the residual is small). That is, the learned functional differs from the target functional only at large times. This again relates to the long-term memory in the target.
Based on this observation, we build this memory effect explicitly into the target functional by considering ρ of the parametrized form
ρω(t) := ρ̄(t) + ρ0,ω(t), (20)
where ρ̄ is the function which can be well-approximated by the model ρ̂, e.g. the exponential sum ρ̄(t) = ∑m∗ j=1 a ∗ je −w∗j t (with w∗j > 0, j = 1, · · · ,m∗). On the other hand, ρ0,ω(t) := ρ0(t − 1/ω) controls the target memory, with ρ0 as any bounded template function in L2(R) ∩ C2(R) with subGaussian tails. As ω → 0+, the support of ρ0,ω shifts towards large times, modelling the dominance of long-term memories. In this case, if the initialization satisfies ρ̂ ≈ ρ̄, the sufficient condition informally discussed above is satisfied as ω → 0+, which heuristically leads to the plateauing.
A simple example of (20) can be ρω(t) = a∗e−w ∗t+ c0e
− (t−1/ω) 2
2σ2 , where a∗, c0, σ 6= 0 and w∗ > 0 are fixed constants. This corresponds to the simple case that m∗ = 1 and ρ0 is the Gaussian density. Observe that as ω → 0+, the memory of this sequence of functionals represented by ρω increases. It can be numerically verified that this simple target functional gives rise to the plateauing behavior, which gets worse significantly as ω → 0+ (see Figure 2 in Appendix C.3.2). Our main result on training dynamics quantifies the plateauing behavior theoretically for general functionals possessing the decomposition (20). For rigorous statements and detailed proofs, see Theorem C.1 in Appendix C.2. Theorem 5.1. Define the loss function Jω as in (18) with the target ρ = ρω as defined in (20). Consider the gradient flow training dynamics
d
dτ θω(τ) = −∇Jω(θω(τ)), θω(0) = θ0, (21)
where θω(τ) := (aω(τ), wω(τ)) ∈ R2m for any τ ≥ 0, and θ0 := (a0, w0). For any ω > 0, m ∈ N+ and θ0 ∈ Rm × Rm+ , 0 < δ 1, define the hitting time
τ0 = τ0(δ;ω,m, θ0) := inf {τ ≥ 0 : |Jω(θω(τ))− Jω(θ0)| > δ} . (22) Assume that m > m∗, and the initialization is bounded and satisfies ρ̂(t; θ0) ≈ ρ̄(t). Then
τ0(δ;ω,m, θ0) & ω 2ec0/ω min { δ√ m , ln(1 + δ) } (23)
for any ω > 0 sufficiently small, where c0 and & hide universal positive constants independent of ω, m and θ0.
Let us sketch the intuition behind Theorem 5.1. Suppose that we currently have a good approximation ρ̂ of the short-term memory part ρ̄, then we can show that the loss is large (J = O(1)) since the long-term memory part ρ0,ω is not accounted for. However, the gradient now is small (∇J = o(1)), since the gradient corresponding to the long-term memory part is concentrated at large t, and thus modulated by exponentially decayed multipliers (see (19)). This implies slowdowns in the training dynamics in the region ρ̂ ≈ ρ̄. It remains to estimate a lower bound on the timescale of escaping from this region, which depends on the curvature of the loss function. In particular, we show that ∇2J is positive semi-definite when ω = 0, but has O(1) positive eigenvalues and multiple o(1) (can be exponentially small) eigenvalues for any 0 < ω 1. Hence, a local linearization analysis implies an exponentially increasing escape timescale, as indicated in (23).
While the target form (20) may appear restrictive, we emphasize that some restrictions on the type of functionals is necessary, since plateauing does not always occur (see Figure 1). In fact, a goal of the preceding analysis is to establish a family of functionals for which exponential slowdowns in training provably occurs, and this can be related to memories of target functionals in a precise way.
The curse of memory in optimization. The timescale proved in Theorem 5.1 is verified numerically in Figure 3 in Appendix C.3.3, where we also show that the analytical setting here is representative of general cases, where plateauing occurs even for non-linear RNNs trained with accelerated optimizers, as long as the target functional has the memory structure imposed in (20).
Theorem 5.1 reveals another aspect of the curse of memory, this time in optimization. When ω → 0+, the influence of target functionalHt does not decay, much like the case considered in the curse of memory in approximation. However, different from the approximation case where an exponentially large number of hidden states is required to achieve approximation tolerance, here in optimization the adverse effect of memory comes with the exponentially pronounced slowdowns of the gradient flow training dynamics. While this is theoretically proved under sensible but restrictive settings, we show numerically in Appendix C.3.3 (Figure 4) that this is representative of general cases.
In the literature, a number of results have been obtained pertaining to the analysis of training dynamics of RNNs. A positive result for training by GD is established in Hardt et al. (2018), but this is in the setting of identifying hidden systems, i.e. the target functional comes from a linear dynamical system, hence it must possess good decay properties provided stablity. On the other hand, convergence can also be ensured if the RNN is sufficiently over-parameterized (large m; Allen-Zhu et al. (2019)). However, both of these settings may not be sufficient in reality. Here we provide an alternative analysis of a setting that is representative of the difficulties one may encounter in practice. In particular, the curse of memory that we established here is consistent with the difficulty in RNN training often observed in applications, where heuristic attributions to memory are often alluded to Hu et al. (2018); Campos et al. (2017); Talathi & Vartak (2015); Li et al. (2018). The analysis here makes the connection between memories and optimization difficulties precise, and may form a basis for future developments to overcome such difficulties in applications.

heading: 6 CONCLUSION
text: In this paper, we analyzed the basic approximation and optimization aspects of using RNNs to learn input-output relationships involving temporal sequences in the linear, continuous-time setting. In particular, we coined the concept curse of memory, and revealed two of its facets. That is, when the target relationship has the long-term memory, both approximation and optimization become exceedingly difficult. These analyses make concrete heuristic observations of the adverse effect of memory on learning RNNs. Moreover, it quantifies the interaction between the structure of the model (RNN functionals) and the structure of the data (target functionals). The latter is a much less studied topic. Here, we adopt a continuous-time approach in order gain access to more quantitative tools, including classical results in approximation theory and stochastic analysis, which help us derive precise results in approximation rates and optimization dynamics. The extension of these results to discrete time may be performed via numerical analysis in subsequent work. More broadly, this approach may be a basic starting point for understanding learning from partially observed time series data in general, including gated RNN variants (Hochreiter & Schmidhuber, 1997; Cho et al., 2014) and other methods such as transformers and convolution-based approaches (Vaswani et al., 2017; Oord et al., 2016). These are certainly worthy of future exploration.

heading: A UNIVERSAL APPROXIMATION THEOREM OF LINEAR FUNCTIONALS BY LINEAR RNNS
text: A key simplification of considering linear functionals is due to the classical representation result below, which allows us to pass from the approximation of functionals to the approximation of functions. In short, this theorem says that while for each measure µ, x 7→ ∫ R x > s dµ(s) defines a linear functional, this is in fact the only way to define a linear functional: for any linear functional H there exists a unique measure µ such that H(x) = ∫ R x > s dµ(s). Theorem A.1 (Riesz-Markov-Kakutani representation theorem). Let H : X → R be a continuous linear functional. Then, there exists a unique, vector-valued, regular, countably additive signed measure µ on R such that
H(x) = ∫ R x>s dµ(s) = d∑ i=1 ∫ R xs,idµi(s). (24)
Moreover, we have
‖H‖ := sup ‖x‖X≤1 |H(x)| = ‖µ‖1(R) := ∑ i |µi|(R). (25)
Proof. Well-known, see e.g. Bogachev (2007), CH 7.10.4.
We will use the representation theorem to prove Theorem 4.1. First, we prove some lemmas. The first shows that there is in fact a common representation of a sequence of linear functionals satisfying the assumptions of Theorem 4.1. Lemma A.1. Let {Ht} be a family of continuous, linear, regular, causal and time homogeneous functionals on X . Then, there exists a measurable function ρ : [0,∞)→ Rd that is integrable, i.e.
‖ρ‖L1([0,∞)) := d∑ i=1 ∫ ∞ 0 |ρi(s)|ds <∞ (26)
and
Ht(x) = ∫ ∞ 0 x>t−sρ(s)ds, t ∈ R. (27)
In particular, {Ht} is uniformly bounded with supt ‖Ht‖ = ‖ρ‖L1([0,∞)) and t 7→ Ht(x) is continuous for all x ∈ X .
Proof. By the Riesz-Markov-Kakutani representation theorem (Theorem A.1), for each t there is a unique regular signed Borel measure µt such that
Ht(x) = ∫ R x>s dµt(s), (28)
and ∑ i |µt,i|(R) = ‖Ht‖. Since {Ht} is causal, we must have ∫∞ t x>s dµt(s) = 0 for any x and thus
Ht(x) = ∫ t −∞ x>s dµt(s). (29)
Now, by time homogeneity we have∫ t −∞ x>s dµt(s) = Ht(x) = Ht+τ (x (τ)) = ∫ t+τ −∞ x>s−τdµt+τ (s). (30)
Take τ = −t and set µ = −µ0 to get
Ht(x) = ∫ ∞ 0 x>t−sdµ(s). (31)
Note that we have ‖µ‖1([0,∞)) = ‖µ0‖1([0,∞)) = ‖H0‖ = ‖Ht‖, and continuity follows from the fact that |Ht+δ(x)−Ht(x)| = ∣∣∣∣∫ ∞
0
(xt+δ−s − xt−s)>dµ(s) ∣∣∣∣
≤ ∑ i ∫ ∞ 0 ‖xt+δ−s − xt−s‖∞d|µi|(s), (32)
which converges to 0 as δ → 0 by dominated convergence theorem. Finally, we will show that each µi is absolutely continuous with respect to λ (Lebesgue measure). Take a measurable E ⊂ [0,∞) such that λ(E) = 0 and set E′ = [0,∞)\E. For each n ≥ 0 set Kn ⊂ E,K ′n ⊂ E′ where Kn,K ′n are closed and µi(E \Kn) ≤ 1/n, µi(E′ \K ′n) ≤ 1/n. For a fixed i ∈ {1, . . . , d}, define x(n) to be such that x(n)t−s,j = 0 for all j 6= i and all s. For j = i, we set x (n) t−s,i = 1 if s ∈ Kn and 0 if s ∈ K ′n, which can then be continuously extended to [0,∞). Observe that by construction, x(n)t−s → 0 for λ-a.e. s, thus by dominated convergence theorem
0 = lim n→∞
Ht(x (n)) = µi(E). (33)
This shows that µi is absolutely continuous with respect to λ, and by the Radon-Nikodym theorem there exists a measurable function ρi : [0,∞)→ R such that for any measurable A ⊂ R we have∫
A
dµi(s) = ∫ A ρi(s)ds, (34)
for i = 1, . . . , d. Hence, we have
Ht(x) = ∫ ∞ 0 x>t−sρ(s)ds (35)
with ‖ρ‖L1([0,∞)) = ∑ i ∫∞ 0 |ρi(s)|ds = ‖µ‖1([0,∞)) <∞.
Lemma A.2. Let ρ : [0,∞)→ R a Lebesgue integrable function, i.e. ‖ρ‖L1([0,∞)) <∞. Then, for any > 0, there exists a polynomial p with p(0) = 0 such that∥∥ρ− p(e−·)∥∥
L1([0,∞)) = ∫ ∞ 0 |ρ(t)− p(e−t)|dt ≤ . (36)
Proof. The approach here is similar to that of the approximation of functions using exponential sums (Kammler, 1976; Braess, 1986). Alternatively, one may also appeal to the density of phase type distributions (He & Zhang, 2007; O’Cinneide, 1990) in the space of positive distributions, and generalizing them to signed measures.
Fix > 0. Define
R(u) = { 1 uρ(− log u), u ∈ (0, 1], 0, u = 0.
(37)
Then, we can check that ‖R‖L1([0,1]) = ‖ρ‖L1([0,∞)) <∞. (38)
By density of continuous functions in L1 there exists a continuous function R̃ on [0, 1] with R̃(0) = 0 such that
‖R− R̃‖L1([0,1]) ≤ /2. (39) By Müntz-Szász theorem (Müntz, 1914; Szász, 1916), there exists a polynomial p with p(0) = 0 such that
‖q − R̃‖L1([0,1]) ≤ /2, (40) and q(u) := p(u)/u is also a polynomial. Therefore, we have∥∥ρ− p(e−·)∥∥
L1([0,∞)) = ∫ 1 0 |R(u)− p(u)/u|du
≤ ∫ 1
0
|R(u)− R̃(u)|du+ ∫ 1
0
|R̃(u)− p(u)/u|du ≤ . (41)
We are now ready to present the proof of Theorem 4.1.
Proof of Theorem 4.1. By (7), for each {Ĥt} ∈ Ĥ we can write
Ĥt(x) = ∫ ∞ 0 x>t−s(U >[eWs]>c)ds. (42)
By Lemma A.1, we can write
Ht(x) = ∫ ∞ 0 x>t−sρ(s)ds, (43)
where ρ is integrable. Thus, we can apply Lemma A.2 to conclude that there exists polynomials pi, i = 1, . . . , d with pi(0) = 0 such that∑
i
‖ρi − pi(e−·)‖L1([0,∞)) ≤ . (44)
Notice that we can write each pi(u) = ∑m j=1 αiju
j for somem equaling the maximal order of {pi}. Taking W = diag(−1, . . . ,−m), c = (1, . . . , 1) and Uij = αji, we have
(U>[eWs]>c)i = pi(e −s), i = 1, . . . , d. (45)
Consequently, we have for any x with ‖x‖∞ ≤ 1, |Ht(x)− Ĥt(x)| = ∣∣∣∣∫ ∞
0
x>t−sρ(s)ds− ∫ ∞
0
x>t−sp(e −s)ds ∣∣∣∣ ≤ ∑ i ∫ ∞ 0 |xt−s,i| ∣∣ρi(s)− pi(e−s)∣∣ ds ≤∑ i ‖ρi − pi(e−·)‖L1([0,∞))
≤ .
(46)

heading: B APPROXIMATION RATES AND THE CURSE OF MEMORY IN APPROXIMATION
text: We first give the proof of Theorem 4.2.
Proof of Theorem 4.2. We fix i ∈ {1, . . . , d} below until the last part of the proof. By Lemma A.1, there exists ρi(t) ∈ Cα[0,∞) such that
yi(t) = Ht(ei) = ∫ t 0 ρi(r)dr, t ≥ 0. (47)
By the assumption,
ρ (k) i (t) = o(e −βt) as t→∞, k = 0, . . . , α. (48)
Consider the transform
qi(s) = { 0 s = 0, ρi(−(α+1) log sβ )
s s ∈ (0, 1]. (49)
For k = 0, . . . , α, one can prove by induction that
q (k) i (s) = (−1) k k∑ j=0 c(j, k) ( α+ 1 β )j ρ(j)i (−(α+1) log sβ ) sk+1 , (50)
where c(j, k) are some integer constants. Together with the assumption, we have∣∣∣q(k)i (e− βα+1 t)∣∣∣ = ∣∣∣∣∣∣ k∑ j=0 c(j, k) ( α+ 1 β )j ρ (j) i (t) e− (k+1)β α+1 t ∣∣∣∣∣∣ ≤ k∑ j=0
∣∣c(j, k)|(α+ 1)j∣∣ γ ≤ C(α)γ, (51) where C(α) is a universal constant only depending on α. Note that for j = 0, . . . , α,
lim s→0+
ρ (j) i
( −(α+1) log s
β ) sk+1 = lim t→∞ ρ (j) i (t)
e− (k+1)β α+1 t
= lim t→∞
ρ (j) i (t)
e−βt e− (α−k)β α+1 t = 0, (52)
hence qi(s) ∈ Cα[0, 1] with qi(0) = q(1)i (0) = · · · = q (α) i (0) = 0. By Jackson’s theorem (Jackson, 1930), for m = 1, 2, . . . , there exists a polynomial Qi,m of degree m− 1 such that
‖qi −Qi,m‖L∞([0,1]) ≤ C(α)γ
mα . (53)
Denote the polynomial Qi,m as
Qi,m(s) = m−1∑ j=0 αi,js j , (54)
and define
φi,m(t) = e − βα+1 tQi,m(e − βα+1 t). (55) Then we have
φi,m(t) = c >eWtui, (56)
where c = (1, 1, . . . , 1), (57)
W =  − βα+1 − 2βα+1 . . .
− mβα+1  , (58) ui = (αi,0, αi,1, . . . , αi,m−1). (59)
With a change of variable s = e− β α+1 t, we have the estimate ‖ρi − φi,m‖L1([0,∞)) = ∫ ∞
0
|ρi(t)− φi,m(t)|dt
= ∫ 1 0 ∣∣∣∣ρi(−(α+ 1) log sβ ) − sQi,m(s) ∣∣∣∣ α+ 1βs ds = α+ 1
β ∫ 1 0 |qi(s)−Qi,m(s)|ds
≤ C(α)γ βmα .
(60)
Finally we define U = [u1, . . . , ud] ∈ Rm×d and have c>eWtU = (φ1,m(t), . . . , φd,m(t)). (61)
Recall the dynamical system (6) d
dt ht = σ(Wht + Uxt), ŷt = c
>ht, (62)
which is together determined by the parameters c,W,U . Similar to the argument in the proof of Theorem 4.1, for any x with ‖x‖∞ ≤ 1 and t, we have
|Ht(x)− Ĥt(x)| ≤ ∑ i ‖ρi − φi,m‖L1([0,∞)) ≤ C(α)γd βmα . (63)
The curse of memory in approximation. Now we explain more technical details of why Theorem 4.2 implies the curse of memory in approximation, as pointed out in the main text. We assume d = 1 and consider an example
Ht(x) := ∫ t 0 xt−sρ(s)ds, t ≥ 0, (64)
in which the density satisfies ρ(t) ∈ C(1)(R) and
ρ(t) ∼ t−(1+ω) as t→ +∞. (65) Here ω > 0 indicates the decay rate of the memory effects in our target functional family {Ht}. The smaller its value, the slower the decay and the longer the system memory. Notice that y(1)(t) = ρ(t), and in this case there exists no β > 0 making the following condition ((13) in the main text) true:
eβty (k) i (t) = o(1) as t→ +∞, and sup t≥0 β−k|eβty(k)i (t)| ≤ γ, (66)
and no rate estimate can be deduced from it.
A natural way to circumvent this obstacle is to introduce a truncation in time. With T ( 1) we can define ρ̃(t) ∈ C(1)(R) such that ρ̃(t) ≡ ρ(t) for t ≤ T , ρ̃(t) ≡ 0 for t ≥ T + 1, and ρ̃(t) is monotonically decreasing for T ≤ t ≤ T + 1. Considering the linear functional H̃t(x) :=∫ t
0 xt−sρ̃(s)ds, we have the truncation error estimate |Ht(x)− H̃t(x)| ≤ ‖x‖X (∫ ∞
T
|ρ(s)|ds ) ∼ ‖x‖XT−ω. (67)
Now the conclusion of Theorem 4.2 (i.e. (63)) is applicable to the truncated {H̃t} (with α = 1), and we have for ∀β > 0, there is a linear RNN (U,W, c) such that the associated functionals {Ĥt} ∈ Ĥm satisfy
sup t∈R ‖H̃t − Ĥt‖ ≤
Cγ βm := C βm sup t≥0 |eβty(1)(t)| β = Cω m eβT β2Tω+1 . (68)
It is straightforward to verify that when β = 2/T , the right-hand side of (68) achieves the minimum, which gives us
sup t∈R ‖H̃t − Ĥt‖ ≤
Cω
m T 1−ω. (69)
Combining (67) and (69) gives us
sup t∈R ‖Ht − Ĥt‖ ≤ C
( T−ω + ω
m T 1−ω
) . (70)
In order to achieve an error tolerance , according to the first term above we require T ∼ − 1ω , and then according to second term we have
m = O ( ωT 1−ω ) = O ( ω − 1 ω ) . (71)
This estimate gives us a quantitative relationship between the degree of freedom needed and the decay speed. When ω is small, i.e., the system has long memory, the size of the RNN required grows exponentially.

heading: C FINE-GRAINED ANALYSIS OF OPTIMIZATION DYNAMICS

heading: C.1 SIMPLIFICATIONS ON THE LOSS FUNCTION
text: Recall the loss function
ExJT (x; c,W,U) = Ex ∣∣∣∣∣ ∫ T
0
[c>eWtU − ρ(t)>]xT−tdt ∣∣∣∣∣ 2 . (72)
The simplifications on ExJT are listed as follows.
1. Take the input data x to be the white noise, so that
xT−tdt in distribution = dBt, (73)
where Bt is the standard d-dimensional Wiener process. As a consequence, simplifying (72) via Itô’s isometry gives
JT (c,W,U) := ExJT (x; c,W,U) = ∫ T
0
∥∥c>eWtU − ρ(t)>∥∥2 2 dt. (74)
2. We focus on the temporal dimension and take the spatial dimension d = 1 in (74).1 Moreover, to investigate the effect of long-term memory, it is necessary to consider the training on large time horizons. Hence, we take T →∞ to get
J∞(c,W, b) := ∫ ∞ 0 ( c>eWtb− ρ(t) )2 dt, (75)
where b is the sole column of U in (74) and ρ(t) becomes a scalar-valued target. This corresponds to the so-called single-input-single-output (SISO) system.
3. Due to the difficulty of directly analyzing ∇W eWt and ∇2W eWt, we consider a further simplified ansatz. Assume that W is a diagonal matrix with negative entries (to guarantee the stability of the model). That is, W = −diag(w) with w ∈ Rm+ . Then we can combine a = b ◦ c (where ◦ denotes the Hadamard product) and rewrite the model as
ρ̂(t; c,W, b) := c>eWtb = m∑ i=1 aie −wit , a>e−wt , ρ̂(t; a,w). (76)
The optimization problem (75) becomes
min a∈Rm,w∈Rm+ J(a,w) := ∫ ∞ 0 ( m∑ i=1 aie −wit − ρ(t) )2 dt. (77)
Here we omit the subscript∞. 2
4. We apply a continuous-time idealization of the gradient descent dynamics by considering the gradient flow with respect to J(a,w). That is,{
a′(τ) = −∇aJ(a(τ), w(τ)), w′(τ) = −∇wJ(a(τ), w(τ)), (78)
with some initial value a(0) = a0 ∈ Rm, w(0) = w0 ∈ Rm+ .
As we will show later, applying the training dynamics (78) to optimize (77) is able to serve as a starting point in the fine-grained dynamical analysis, since it still preserves the plateauing behavior observed in the optimization process (see Figure 1), provided additional structures related to memories as discussed next.

heading: C.2 CONCRETE DYNAMICAL ANALYSIS AND THE CURSE OF MEMORY IN OPTIMIZATION
text: We prove Theorem 5.1 in this section. The basic insight is, by adding long-term memories in targets, one can increase the loss with little effect on the gradient and Hessian, which leads to a significant slow down of the training dynamics near the short-term memory parts of the targets. Therefore, Theorem 5.1 is proved subsequently in the following procedure.
1. We prove that Jω has a large value but small gradient when ρ̂(t; a,w) ≡ ρ̄(t). 1One can see that the spatial dimension plays little role in the previous approximation analysis (see the proof of Theorem 4.2), since each spatial dimension can be handled separately. 2The time horizon is always taken as ∞ in the whole analysis. Note that here we also omit an index m (width of the network, relating to the model capacity), since it remains unchanged in the following content if not specified.
2. We prove that when ρ̂(t; a,w) ≡ ρ̄(t), the Hessian ∇2Jω is positive semi-definite for ω = 0, but for finite, small ω > 0, ∇2Jω has O(1) positive eigenvalues and multiple o(1) eigenvalues.
3. Based on these results, we perform a local linearization analysis on the gradient flow (21) initialized by ρ̂(t; a0, w0) ≡ ρ̄(t), from which and by continuity the timescale of plateauing is derived.
(1) Preliminary Results
As discussed around (20), we consider the target functional with a parametrized representation
ρω(t) = ρ̄(t) + ρ0,ω(t) = m∗∑ j=1 a∗je −w∗j t + ρ0(t− 1/ω).

heading: Here, a∗j := a
text: ∗(w∗j ) 6= 0, w∗j > 0 and w∗i 6= w∗j for any i 6= j, i, j ∈ [m∗], 3 and m∗ < m. The former requirements are just non-degenerate conditions, and the last requirement ensures that the model can perfectly represent the well-approximated part of the target, ρ̄(t). The memory in the target is controlled by ρ0,ω(t) = ρ0(t − 1/ω), with ρ0 as a fixed template function which satisfies the following assumptions.
Assumptions on ρ0. (i) ρ0(t) 6≡ 0; (ii) ρ0 ∈ L2(R) ∩ C2(R); (iii) ρ0 is bounded on R, i.e. ‖ρ0‖L∞(R) <∞; (iv) lim
t→−∞ ρ0(t) = 0.
Remark C.1. The above assumptions (i)(ii)(iii) are rather natural, and (iv) only restricts the single side tail of ρ0 to be zero. In the following analysis, we further focus on ρ0 with light tails, e.g. the sub-Gaussian tail
|ρ0(t)| ≤ c0e−c1t 2
, ∀t : |t| ≥ t0 (79) for some fixed positive constants c0, c1, t0. Obviously, Gaussian densities and continuous functions with compact supports are sub-Gaussian functions.
We begin by the following preliminary estimate that is used throughout the subsequent analysis. Lemma C.1. For any n ∈ N, ω > 0 and w > 0, let
∆n,ω(w) := ∫ ∞ 0 tne−wt|ρ0,ω(t)|dt. (80)
Then
• ∆n,ω(w) is monotonically decreasing on (0,∞);
• lim ω→0+ ∆n,ω(w) = 0;
• In particular, if ρ0 is sub-Gaussian, we further have
∆n,ω(w) . ω −ne−w/ω ( cw 2 2 + c w 3 ) , ω ∈ (0,min{1/2, 1/t0, 2c1/w}). (81)
Here c2 = e 1
4c1 > 1, c3 = et0 > 1, and . hides universal constants only depending on n and ρ0, t0, c0, c1.
Proof. (i) Obviously ∆n,ω(w1) ≤ ∆n,ω(w2) for any w1 > w2 > 0. (ii) By the assumptions on ρ0, we get
lim ω→0+ |ρ0,ω(t)| = lim ω→0+ |ρ0(t− 1/ω)| = lim s→−∞ |ρ0(s)| = 0, ∀t ≥ 0,
and M0 := ‖ρ0‖L∞(R) < ∞, which gives tne−wt|ρ0,ω(t)| ≤ M0tne−wt ∈ L1([0,∞)) for any n ∈ N, ω > 0 and w > 0. By Lebesgue’s dominant convergence theorem, we have
lim ω→0+ ∆n,ω(w) = ∫ ∞ 0 tne−wt · lim ω→0+ |ρ0,ω(t)|dt = 0, ∀n ∈ N, ∀w > 0. (82)
3For any n ∈ N+, [n] := {1, 2, · · · , n}.
(iii) Now we estimate ∆n,ω(w) under the sub-Gaussian condition (79). Suppose 0 < ω < 1/t0, we have |∆n,ω(w)| = ∫ ∞
0
tne−wt|ρ0(t− 1/ω)|dt
= ∫ 1/ω+t0 1/ω−t0 tne−wt|ρ0(t− 1/ω)|dt+ ∫ 1/ω−t0 0 tne−wt|ρ0(t− 1/ω)|dt
+ ∫ ∞ 1/ω+t0 tne−wt|ρ0(t− 1/ω)|dt , I1 + I2 + I3.
Then we bound I1, I2 and I3 respectively: I1 ≤M0 ∫ 1/ω+t0
1/ω−t0 tne−wtdt ≤M0e−w(1/ω−t0) ∫ 1/ω+t0 1/ω−t0 tndt
= M0e wt0 · e−w/ω · (1 + ωt0)
n+1 − (1− ωt0)n+1
(n+ 1)ωn+1
.M0e wt0ω−ne−w/ω(t0 + ω),
where ω ∈ (0, 1/2), and . hides universal constants only related to n and t0. Let 1/c1 := 2σ2, we have
I2 = e −w/ω ∫ −t0 −1/ω (s+ 1/ω)ne−ws|ρ0(s)|ds ≤ e−w/ω ∫ −t0 −1/ω (s+ 1/ω)ne−ws · c0e−c1s 2 ds,
where∫ −t0 −1/ω (s+ 1/ω)ne−wse−c1s 2 ds = e w2 4c1 ∫ −t0 −1/ω (s+ 1/ω)ne−c1(s+ w 2c1 )2ds
≤ eσ 2w2/2 ∫ R (|t|+ |1/ω − σ2w|)n · e− t2 2σ2 dt
= eσ 2w2/2 n∑ k=0 Ckn(1/ω − σ2w)n−k · 2 ∫ ∞ 0 tke− t2 2σ2 dt
≤ eσ 2w2/2 n∑ k=0 Ckn(1/ω) n−k(√2σ)k+1Γ(k + 1 2 ) holds for any ω ∈ (0, 2c1/w). Here the last inequality is due to the Mellin Transform of absolute moments of the Gaussian density (see Proposition 1 in Li et al. (2005)). The argument is similar for I3, which gives the same bound as I2.
Combining all the estimates gives the desired conclusion. The proof is completed.
The main idea to analyze plateauing behaviors is to investigate the local dynamics of the gradient flow (21) when ρ̂ = ρ̄, then extend the results to the setting ρ̂ ≈ ρ̄ by continuity. Recall that both of them are exponential sums, we can obtain the relation of parameters between ρ̂ and ρ̄, according to the following lemma. Lemma C.2. For any m ∈ N+, let λ = (λ1, · · · , λm) with λi 6= λj for any i 6= j, i, j ∈ [m]. Then the series of functions { eλit }m i=1 is linear independent on any interval I ⊂ R.
Proof. The aim is to show m∑ i=1 cie λit = 0, t ∈ I ⇒ ci = 0, ∀i ∈ [m]. (83) (83) holds trivially for m = 1. Assume that (83) holds for m− 1, then m∑ i=1 cie λit = 0, t ∈ I ⇒ m−1∑ i=1 cie (λi−λm)t + cm = 0, t ∈ I (84)
⇒ m−1∑ i=1 ci(λi − λm)e(λi−λm)t = 0, t ∈ I.
By induction, we get ci(λi − λm) = 0 for any i = 1, · · · ,m − 1. Since λ1, · · · , λm are distinct, we have ci = 0 for any i = 1, · · · ,m− 1. Together with (84), we get cm = 0, which completes the proof.
Definition C.1. Let m ≥ m∗. For any partition P: [m] = ∪m∗j=0Ij with Ij1 ∩ Ij2 = ∅ for any j1 6= j2, j1, j2 ∈ {0} ∪ [m∗], and I0 = ∪i0r=1I0,r with I0,r1 ∩ I0,r2 = ∅ for any r1 6= r2, r1, r2 ∈ [i0], where Ij 6= ∅ for any j ∈ [m∗] and I0,r 6= ∅ for any r ∈ [i0] (if I0 6= ∅), define the affine space (with respect to P):
M∗P := { (a,w) ∈ Rm × Rm+ : ∑ i∈Ij ai = a ∗ j , wi = w ∗ j for any i ∈ Ij , j ∈ [m∗];
∑ i∈I0,r ai = 0, wi = vr 6= w∗j for any i ∈ I0,r, r ∈ [i0] and j ∈ [m∗]
} .
Denote the collection of all such affine spaces byM∗ := ⋃ PM∗P .
The following lemma characterizes the relation of parameters, by showing thatM∗ is exactly the set of equivalent points to (a∗, w∗) for the purpose of representation via exponential sums.
Lemma C.3. For any (a,w) ∈ Rm × Rm+ , ρ̂(t; a,w) ≡ ρ̄(t)⇔ (a,w) ∈M∗.
Proof. (i) (⇐) Since (a,w) ∈M∗, there exists P such that (a,w) ∈M∗P . Then for any t ≥ 0,
ρ̂(t; a,w) = m∑ i=1 aie −wit = m∗∑ j=0 ∑ i∈Ij aie −wit = i0∑ r=1 ∑ i∈I0,r aie −wit + m∗∑ j=1 ∑ i∈Ij aie −wit
= i0∑ r=1 ( ∑ i∈I0,r ai ) e−vrt + m∗∑ j=1 (∑ i∈Ij ai ) e−w ∗ j t = m∗∑ j=1 a∗je −w∗j t = ρ̄(t).
(ii) (⇒) Let Ij = { i ∈ [m] : wi = w∗j } for any j ∈ [m∗], and I0 = { i ∈ [m] : wi 6= w∗j for any j ∈
[m∗] } . Recall that ρ̄(t) = ∑m∗ j=1 a ∗ je −w∗j t is non-degenerate: a∗j 6= 0, w∗j > 0 and w∗i 6= w∗j for any i 6= j, i, j ∈ [m∗], we get [m] = ∪m∗j=0Ij , Ij1 ∩ Ij2 = ∅ for any j1 6= j2, j1, j2 ∈ {0} ∪ [m∗]. Combining Lemma C.2 and the non-degeneracy of ρ̄ , Ij 6= ∅ for any j ∈ [m∗]. Assume that there are i0 different components in (wi)i∈I0 , say v1, · · · , vi0 , then vr 6= w∗j for any r ∈ [i0] and j ∈ [m∗]. Let I0,r = { i ∈ I0 : wi = vr } for any r ∈ [i0], we get I0,r 6= ∅ for any r ∈ [i0] (if I0 6= ∅), and I0 = ∪i0r=1I0,r, and I0,r1 ∩ I0,r2 = ∅ for any r1 6= r2, r1, r2 ∈ [i0]. Hence [m] = ∪m∗j=0Ij with I0 = ∪ i0 r=1I0,r forms a P defined in Definition C.1, and
0 ≡ ρ̂(t; a,w)− ρ̄(t) = m∑ i=1 aie −wit − m∗∑ j=1 a∗je −w∗j t
= m∗∑ j=0 ∑ i∈Ij aie −wit − m∗∑ j=1 a∗je −w∗j t
= i0∑ r=1 ∑ i∈I0,r aie −wit + m∗∑ j=1 ∑ i∈Ij aie −wit − m∗∑ j=1 a∗je −w∗j t  =
i0∑ r=1 ( ∑ i∈I0,r ai ) e−vrt + m∗∑ j=1 (∑ i∈Ij ai − a∗j ) e−w ∗ j t.
Again by Lemma C.2, we have ∑ i∈Ij ai = a ∗ j for any j ∈ [m∗] and ∑ i∈I0,r ai = 0 for any r ∈ [i0], which gives (a,w) ∈M∗P . The proof is completed.
Remark C.2. Let I0 = ∅. 4 It is straightforward to check that for any partition P , the dimension ofM∗P is ∑m∗ j=1(|Ij | − 1) = m −m∗. In addition, it can be verified that the cardinality ofM∗ is
m∗! { m m∗ } , where { m m∗ } is the Stirling number of the second kind. 5
(2) Loss and Gradient
Now we show that as ω → 0+, the loss Jω remains lower bounded, while ‖∇Jω‖2 converges to 0. This implies that the loss will saturate at a high value when learning long-term memories.
Proposition C.1. For any (a,w) ∈ Rm × Rm+ satisfying ρ̂(t; a,w) ≡ ρ̄(t), we have
Jω(a,w) ≥ C(ρ0) > 0, ∀ω ∈ (0, C ′(ρ0)), (85)
where C(ρ0), C ′(ρ0) > 0 are constants only depending on ρ0. That is, the loss is lower bounded uniformly in ω.
Proof. Recall the assumptions on ρ0, we have ρ0(t) 6≡ 0 and ρ0 ∈ C(R). Let t1 ∈ R such that ρ0(t1) 6= 0. By continuity, there exists δ0 > 0 such that |ρ0(t)| ≥ |ρ0(t1)|/2 for any t ∈ [t1 − δ0, t1 + δ0]. Hence, for any ω > 0 sufficiently small such that −1/ω < t1 − δ0, we have
‖ρ0,ω‖2L2[0,∞) = ∫ ∞
0 ρ20(t− 1/ω)dt = ∫ ∞ − 1ω ρ20(t)dt ≥ ∫ t1+δ0 t1−δ0 ρ20(t)dt ≥ 1 2 δ0|ρ0(t1)|2 > 0.
Fix any (â, ŵ) ∈ Rm × Rm+ such that ρ̂(t; â, ŵ) ≡ ρ̄(t). Then
Jω(â, ŵ) = ‖ρ̂(t; â, ŵ)− ρ̄(t)− ρ0,ω(t)‖2L2[0,∞) = ‖ρ0,ω‖ 2 L2[0,∞) ≥
1 2 δ0|ρ0(t1)|2 > 0,
which completes the proof.
Proposition C.2. For any (a,w) ∈ Rm × Rm+ satisfying ρ̂(t; a,w) ≡ ρ̄(t), we have
lim ω→0+
‖∇Jω(a,w)‖2 = 0. (86)
In particular, if ρ0 has the sub-Gaussian tail (79), the estimate
‖∇Jω(a,w)‖2 . √ mω−1e−wmin/ω ( c w2min 2 + c wmin 3 ) (1 + ‖a‖∞) (87)
holds for any ω ∈ (0,min{1/2, 1/t0, 2c1/wmin}). Here wmin := mini∈[m] wi > 0, c2, c3 > 1 are constants only related to c1, t0, and . hides universal constants only depending on ρ0, t0, c0, c1.
Proof. A straightforward computation shows, for k = 1, 2, · · · ,m,
∂Jω ∂ak (a,w) = 2  m∑ i=1 ai wk + wi − m∗∑ j=1 a∗j wk + w∗j − 2∆0,ω(wk), (88) ∂Jω ∂wk (a,w) = −2ak  m∑ i=1 ai (wk + wi)2 − m∗∑ j=1 a∗j (wk + w∗j ) 2
+ 2ak∆1,ω(wk). (89) Fix any (â, ŵ) ∈ Rm × Rm+ satisfying ρ̂(t; â, ŵ) ≡ ρ̄(t). By Lemma C.3, we have (â, ŵ) ∈ M∗. Recall Definition C.1, there exists a partition P: [m] = ∪m∗j=0Ij with I0 = ∪ i0 r=1I0,r, where Ij 6= ∅ for any j ∈ [m∗] and I0,r 6= ∅ for any r ∈ [i0] (if I0 6= ∅), such that (â, ŵ) ∈ M∗P , which gives
4That is, the non-degenerate case. Obviously, I0 6= ∅ implies an uncountable M∗P , but they are all degenerate.
5The result follows from basic knowledge of combinatorics. See details in the proof of Theorem D.1.
that ∑ i∈Ij âi = a ∗ j , ŵi = w ∗ j for any i ∈ Ij and j ∈ [m∗]; ∑ i∈I0,r âi = 0, ŵi = vr 6= w ∗ j for any i ∈ I0,r, r ∈ [i0] and j ∈ [m∗]. Therefore, for any n ∈ N+, we have
m∑ i=1 âi (ŵk + ŵi)n − m∗∑ j=1 a∗j (ŵk + w∗j ) n
= i0∑ r=1 ∑ i∈I0,r âi (ŵk + ŵi)n + m∗∑ j=1 ∑ i∈Ij âi (ŵk + ŵi)n − m∗∑ j=1 a∗j (ŵk + w∗j ) n
= i0∑ r=1
∑ i∈I0,r âi (ŵk + vr)n + m∗∑ j=1 ∑ i∈Ij âi (ŵk + w∗j ) n − m∗∑ j=1 a∗j (ŵk + w∗j ) n = 0. (90)
This yields
∂Jω ∂ak (â, ŵ) = −2∆0,ω(ŵk), ∂Jω ∂wk (â, ŵ) = 2âk∆1,ω(ŵk),
and hence
‖∇Jω(â, ŵ)‖22 = 4 m∑ k=1 [ ∆20,ω(ŵk) + â 2 k∆ 2 1,ω(ŵk) ] .
By Lemma C.1, we get limω→0+ ‖∇Jω(â, ŵ)‖2 = 0. If ρ0 has the sub-Gaussian tail (79), again by Lemma C.1, the estimate
∆n,ω(ŵk) ≤ ∆n,ω(ŵmin) . ω−ne−ŵmin/ω ( c ŵ2min 2 + c ŵmin 3 ) (91)
holds for any n ∈ N, ω ∈ (0,min{1/2, 1/t0, 2c1/ŵmin}) and k ∈ [m]. Here ŵmin := mini∈[m] ŵi > 0, c2, c3 > 1 are constants only related to c1, t0, and . hides universal constants only depending on n and ρ0, t0, c0, c1. Therefore
‖∇Jω(â, ŵ)‖2 . √ mω−1e−ŵmin/ω ( c ŵ2min 2 + c ŵmin 3 ) (1 + ‖â‖∞) , ω ∈ (0, 1].
The proof is completed.
(3) Eigenvalues of Hessian
Now we show that minimal eigenvalues of∇2Jω also converges to 0 as ω → 0+.
Proposition C.3. For any (a,w) ∈ Rm × Rm+ satisfying ρ̂(t; a,w) ≡ ρ̄(t), denote the eigenvalues of∇2Jω(a,w) by λ1(ω) ≥ λ2(ω) ≥ · · · ≥ λ2m(ω). If m > m∗, we have
λk(ω) > 0, k = 1, 2, · · · ,m′, (92) lim ω→0+ λk(ω) = 0, k = m ′ + 1,m′ + 2, · · · , 2m (93)
for ω > 0 sufficiently small, where m′ ≤ 2m∗ + |I0| ≤ m + m∗. In particular, if ρ0 has the sub-Gaussian tail (79), the estimate
|λk(ω)| . ω−2e−wmin/ω ( c w2min 2 + c wmin 3 ) (1 + ‖a‖∞) k = m′ + 1,m′ + 2, · · · , 2m (94)
holds for any ω ∈ (0,min{1/2, 1/t0, 2c1/wmin}). Here wmin := mini∈[m] wi > 0, c2, c3 > 1 are constants only related to c1, t0, and . hides universal constants only depending on ρ0, t0, c0, c1.
Proof. A straightforward computation shows, for k, j = 1, 2, · · · ,m, ∂2Jω ∂ak∂aj (a,w) = 2 wk + wj , (95)
∂2Jω ∂ak∂wj (a,w) = −2aj (wk + wj)2 , k 6= j, (96)
∂2Jω ∂ak∂wk (a,w) = −2  m∑ i=1 ai (wk + wi)2 − m∗∑ j′=1
a∗j′
(wk + w∗j′) 2 − ak 2w2k + 2∆1,ω(wk), (97)
∂2Jω ∂wk∂wj (a,w) = 4akaj (wk + wj)3 , k 6= j, (98)
∂2Jω ∂wk∂wk (a,w) = 4ak  m∑ i=1 ai (wk + wi)3 − m∗∑ j′=1
a∗j′
(wk + w∗j′) 3 + a2k 2w3k − 2ak∆2,ω(wk). (99)
Fix any (â, ŵ) ∈ Rm × Rm+ satisfying ρ̂(t; â, ŵ) ≡ ρ̄(t). By (90), we have
∂2Jω ∂ak∂aj (â, ŵ) = 2 ŵk + ŵj ,
∂2Jω ∂ak∂wj (â, ŵ) = −2âj (ŵk + ŵj)2 (k 6= j), ∂ 2Jω ∂ak∂wk (â, ŵ) = − âk 2ŵ2k + 2∆1,ω(ŵk),
∂2Jω ∂wk∂wj (â, ŵ) = 4âkâj (ŵk + ŵj)3 (k 6= j), ∂ 2Jω ∂wk∂wk (â, ŵ) = â2k 2ŵ3k − 2âk∆2,ω(ŵk).
Let
J̄(a,w) := ‖ρ̂(t; a,w)− ρ̄(t)‖2L2[0,∞) (100)
= ∥∥∥∥∥∥ m∑ i=1 aie −wit − m∗∑ j=1 a∗je −w∗j t ∥∥∥∥∥∥ 2
L2[0,∞)
,
and
Eω(a,w) := [
Om×m Diag(∆1,ω(w)) Diag(∆1,ω(w)) −Diag(a ◦∆2,ω(w))
] ,
where ∆n,ω(w) (n = 1, 2) is performed element-wisely. One can verify that
∇2Jω(a,w) = ∇2J̄(a,w) + 2Eω(a,w). (101)
Then we analyze∇2J̄(â, ŵ) and Eω(â, ŵ) respectively. (i) ∇2J̄(â, ŵ). Obviously (â, ŵ) is a global minimizer of J̄(a,w) due to J̄(â, ŵ) = 0. Hence ∇J̄(â, ŵ) = 0 and∇2J̄(â, ŵ) is positive semi-definite. We further show that∇2J̄(â, ŵ) has multiple zero eigenvalues when m > m∗. In fact, since
∂2J̄
∂ak∂aj (â, ŵ) =
2
ŵk + ŵj ,
∂2J̄
∂ak∂wj (â, ŵ) = −2âj (ŵk + ŵj)2 , ∂2J̄ ∂wk∂wj (â, ŵ) = 4âkâj (ŵk + ŵj)3 ,
it is straightforward to verify that for any i, j ∈ Ip, p ∈ [m∗] and any i, j ∈ I0,r, r ∈ [i0],
∇2J̄(â, ŵ)i,: = ∇2J̄(â, ŵ)j,:, âj · ∇2J̄(â, ŵ)m+i,: = âi · ∇2J̄(â, ŵ)m+j,:, where Ai,: denotes the i-th row of matrix A. Notice that ∑ i∈I0,r ∇
2J̄(â, ŵ)m+i,: = 0 for any r ∈ [i0], we conclude that the Hessian ∇2J̄(â, ŵ) has at most 2m∗ + i0 + i2 ≤ 2m∗ + |I0| ≤ m + m∗ different rows, 6 which yields rank(∇2J̄(â, ŵ)) ≤ 2m∗ + |I0| ≤ m + m∗. Therefore,
6Here i2 := |{r ∈ [i0] : |I0,r| ≥ 2}|. When I0 = ∅, the upper bound is 2m∗; when I0 6= ∅, since I0,r 6= ∅ for any r ∈ [i0], let i1 := |{r ∈ [i0] : |I0,r| = 1}| and i2 defined as before. Then i0 = i1 + i2, |I0| = ∑i0 r=1 |I0,r| ≥ i1+2i2 = i0+ i2. The last inequality follows from |I0| = m− ∑m∗ j=1 |Ij | ≤ m−m ∗.
the number of zero eigenvalues of ∇2J̄(â, ŵ) ≥ dim { x ∈ R2m : ∇2J̄(â, ŵ) · x = 0 } = 2m − rank(∇2J̄(â, ŵ)) ≥ 2(m−m∗)− |I0| ≥ m−m∗. Since ∇2J̄(â, ŵ) is positive semi-definite, all the non-zero eigenvalues must be positive.
(ii) Eω(â, ŵ). Let
G (1) k := {y ∈ R : |y| ≤ |∆1,ω(ŵk)|} ,
G (2) k := {y ∈ R : |y + âk∆2,ω(ŵk)| ≤ |∆1,ω(ŵk)|} .
By Gershgorin’s circle theorem, for any eigenvalue of Eω(â, ŵ), say λ(ω), we have λ(ω) ∈⋃m k=1(G (1) k ∪G (2) k ). Combining with Lemma C.1, we get
|λ(ω)| ≤ max k∈[m]
( |âk||∆2,ω(ŵk)|+ |∆1,ω(ŵk)| ) → 0, ω → 0+. (102)
If ρ0 has the sub-Gaussian tail (79), by (91), we further have
|λ(ω)| ≤ max k∈[m]
( |âk||∆2,ω(ŵk)|+ |∆1,ω(ŵk)| ) . ω−2e−ŵmin/ω ( c ŵ2min 2 + c ŵmin 3 ) (1 + ‖â‖∞), ω ∈ (0, 1], (103)
where ω ∈ (0,min{1/2, 1/t0, 2c1/ŵmin}), ŵmin := mini∈[m] ŵi > 0, c2, c3 > 1 are constants only related to c1, t0, and . hides universal constants only depending on ρ0, t0, c0, c1.
Combining (i), (ii) and applying Weyl’s theorem gives the desired result.
(4) Local Linearization Analysis
The previous analysis can now be tied directly to a quantitative dynamics via linearization arguments. It is shown that under mild assumptions, the gradient flow (21) can become trapped in plateaus with an exponentially large timescale. That is, the curse of memory occurs, this time in optimization dynamics instead of approximation rates.
Theorem C.1 (Restatement of Theorem 5.1). For any ω > 0, m ∈ N+ and θ0 = (a0, w0) ∈ Rm × Rm+ , 0 < δ 1, define the hitting time
τ0 = τ0(δ;ω,m, θ0) := inf {τ ≥ 0 : ‖θω(τ)− θ0‖2 > δ} , (104) τ ′0 = τ ′ 0(δ;ω,m, θ0) := inf {τ ≥ 0 : |Jω(θω(τ))− Jω(θ0)| > δ} . (105)
Assume that m > m∗, and the initialization satisfies ρ̂(t; θ0) ≈ ρ̄(t). Then we have
lim ω→0+ τ0(δ;ω,m, θ0) = lim ω→0+
τ ′0(δ;ω,m, θ0) = +∞. (106)
In particular, if ρ0 has the sub-Gaussian tail (79), and the initialization is bounded as (a0, w0) ∈ [a0l , a 0 r] m × [w0l , w0r ]m with constants a0l < a0r , 0 < w0l < w0r , we further have
τ ′0(δ;ω,m, θ0) ≥ τ0(δ;ω,m, θ0) & ω2ew 0 l /ω min { δ√ m , ln(1 + δ) } (107)
for any ω ∈ (0,min{1/2, 1/t0, 2c1/w0r}) sufficiently small, where & hides universal constants only depending on ρ0, t0, c0, c1, w0r , a 0 l and a 0 r .
Proof. Consider the asymptotic expansion with the form
θω(τ) = θ 0 ω(τ) + ∞∑ i=1 δiθiω(τ) = θ 0 ω(τ) + δθ 1 ω(τ) + δ 2θ2ω(τ) + o(δ 2), (108)
for some δ ∈ (0, 1) (with δ 1) and θiω(τ) = O(1) (τ ≥ 0, i = 0, 1, · · · ). 7 For consistency, we have θ0ω(0) = θ0 and θ i ω(0) = 0 for i = 1, 2, · · · . By continuity, τ0 > 0 and ‖θω(τ)− θ0‖2 ≤ δ for any τ ∈ [0, τ0]. The aim is to quantify the scale of τ0. 7Here θiω(τ) denotes the i-th term in the asymptotic expansion of θω(τ), not the i-th power.
Let g0 := ∇Jω(θ0) and H0 := ∇2Jω(θ0). The local linearization on (21) shows
d
dτ θω(τ) = −g0 −H0(θω(τ)− θ0) +O(δ2), τ ∈ [0, τ0].
Combining with (108), we have
d
dτ θ0ω(τ) = −g0 −H0(θ0ω(τ)− θ0), θ0ω(0) = θ0, at O(1) scale,
d
dτ θ1ω(τ) = −H0θ1ω(τ), θ1ω(0) = 0, at O(δ) scale,
d
dτ θ2ω(τ) = −H0θ2ω(τ) +O(1), θ2ω(0) = 0, at O(δ2) scale.
Therefore
θ0ω(τ) = θ0 − (∫ τ
0
e−H0sds ) g0,
θ1ω(τ) = e −H0τθ1ω(0) = 0,
which gives
θω(τ) = θ0 − (∫ τ
0
e−H0sds ) g0 +O(δ2), τ ∈ [0, τ0]. (109)
To achieve a parameter separation gap δ0, i.e. ‖θω(τ)− θ0‖2 = δ0 with δ0 = cδ, c ∈ (0, 1], we need to take τ such that ∥∥∥∥(∫ τ
0
e−H0sds ) g0 ∥∥∥∥ 2 ≥ δ0 2 . (110)
Let H0 = P>ΛP be the eigenvalue decomposition with P orthogonal and Λ = diag(λ1, · · ·λ2m) consisting of the eigenvalues of H0 with λ1 ≥ · · · ≥ λ2m. Then∥∥∥∥(∫ τ
0
e−H0sds ) g0 ∥∥∥∥ 2 = ∥∥∥∥P>(∫ τ 0 e−Λsds ) Pg0 ∥∥∥∥ 2 ≤ ∥∥∥∥∫ τ 0 e−Λsds ∥∥∥∥ 2 ‖g0‖2
≤ ‖g0‖2 ·max {
max i∈[2m],λi 6=0
1
|λi| |e−λiτ − 1|, τ
} .
It is straightforward to verify that h(τ ;λ) := 1|λ| |e −λτ − 1|, τ ≥ 0 monotonically decreases on λ ∈ R for any τ ≥ 0. 8 Hence∥∥∥∥(∫ τ 0 e−H0sds ) g0 ∥∥∥∥ 2 ≤ ‖g0‖2 · { 1 −λ2m (e −λ2mτ − 1), λ2m < 0, τ, λ2m ≥ 0, (111)
and the right hand side monotonically increases on τ ≥ 0. Combining (110), (111) gives
δ0 2 ≤ ‖g0‖2 ·
{ 1
−λ2m (e −λ2mτ − 1), λ2m < 0,
τ, λ2m ≥ 0. (112)
We discuss for different cases:
(i) ‖g0‖2 = 0. Obviously the inequality (112) fails since (110) fails for any τ ≥ 0, which gives τ0 = +∞; (ii) ‖g0‖2 6= 0 and λ2m ≥ 0. By (112), we get τ ≥ δ02‖g0‖2 ;
(iii) ‖g0‖2 6= 0 and λ2m < 0. By (112), we get
τ ≥ 1 −λ2m ln
( 1 + δ0
−λ2m 2‖g0‖2
) .
8With the convention that h(τ ; 0) = τ for any τ > 0, and h(0;λ) ≡ 0 for any λ ∈ R.
If −λ2m ≤ 2‖g0‖2, we have
τ ≥ 1 −λ2m · δ0 −λ2m 2‖g0‖2
· ln ( 1 + δ0 −λ2m 2‖g0‖2 ) δ0 −λ2m 2‖g0‖2
= δ0
2‖g0‖2
( 1 +O ( δ0 −λ2m 2‖g0‖2 )) = δ0 2‖g0‖2 (1 +O(δ0));
if −λ2m > 2‖g0‖2, we have τ ≥ ln(1+δ0)−λ2m .
Combining (i), (ii), (iii) gives
τ0 = τ0(δ;ω,m, θ0) & min
{ δ
‖g0‖2 ,
ln(1 + δ)
|λ2m|
} . (113)
Let the initialization satisfy ρ̂(t; θ0) ≡ ρ̄(t), and assume m > m∗. According to Propositions C.2 and C.3, we have
lim ω→0+ ‖g0‖2 = 0, lim ω→0+ λ2m = 0⇒ lim ω→0+ τ0(δ;ω,m, θ0) = +∞. (114)
If ρ0 has the sub-Gaussian tail (79), again by Propositions C.2 and C.3, we further have
τ0(δ;ω,m, θ0) & ω 2ew0,min/ω min { δ√ m , ln(1 + δ) } 1(
c w20,min 2 + c w0,min 3 ) (1 + ‖a0‖∞)
(115)
for any ω ∈ (0,min{1/2, 1/t0, 2c1/w0,min}), where w0,min := mini∈[m] w0,i > 0, c2, c3 > 1 are constants only related to c1, t0, and & hides universal constants only depending on ρ0, t0, c0, c1. Since the initialization is bounded as (a0, w0) ∈ [a0l , a0r]m×[w0l , w0r ]m with a0l < a0r , 0 < w0l < w0r , let c0a = max{|a0l |, |a0r|}, we get
τ0(δ;ω,m, θ0) & ω 2ew 0 l /ω min { δ√ m , ln(1 + δ) } 1(
c (w0r) 2 2 + c w0r 3 ) (1 + c0a)
& ω2ew 0 l /ω min { δ√ m , ln(1 + δ) } , (116)
where & hides universal constants only related to w0r , a 0 l and a 0 r .
The last task is to show the dynamics of the loss is much slower than the parameter separation when there is plateauing. The argument is trivial since for any τ ∈ [0, τ0],
Jω(θω(τ))− Jω(θ0) = g>0 (θω(τ)− θ0) + (θω(τ)− θ0)>H0(θω(τ)− θ0) + o(δ2) ≥ −‖g0‖2‖θω(τ)− θ0‖2 + λ2m‖θω(τ)− θ0‖22 + o(δ2) = o(1)O(δ) + o(1)O(δ2) + o(δ2) = o(δ2), ω → 0+.
By continuity, the proof is completed.
Remark C.3. The estimate in Theorem C.1 shows a lower bound on the escape time, hence it does not appear to preclude the situation that the plateauing lasts forever. However, in the proof above, if one supposes τ0 = +∞ in (104), i.e. the hypothetical situation where the parameters are trapped forever, and write g̃0 := Pg0 = (g̃0,1, · · · , g̃0,2m), we have∥∥∥∥(∫ τ
0
e−H0sds ) g0 ∥∥∥∥2 2 = g̃>0 (∫ τ 0 e−Λsds )2 g̃0 = 2m∑ i=1 (g̃0,i) 2(h(τ ;λi)) 2 ≥ (g̃0,j)2(h(τ ;λj))2
for any j such that λj < 0. If g̃0,j 6= 0, (109) gives ‖θω(τ)− θ0‖2 ≥ ∥∥∥∥(∫ τ
0
e−H0sds ) g0 ∥∥∥∥ 2 +O(δ2)
≥ |g̃0,j | −λj (e−λjτ − 1) +O(δ2)→ +∞, τ →∞,
which is a contradiction. That is to say, the parameter separation has to achieve the gap δ within a finite time, even if it is exponentially large.
Remark C.4. Recall Lemma C.3, ρ̂(t; a0, w0) ≡ ρ̄(t) if and only if (a0, w0) ∈ M∗ = ⋃ PM∗P , where P is a partition over [m] as defined in Definition C.1. That is, as a union of affine spaces, M∗ is in fact an equivalent set for qualified initializations. As discussed in Remark C.2, when there
is no degeneracy, the cardinality of M∗ is m∗! { m m∗ } (i.e. the number of P), with each M∗P an (m − m∗)-dimensional affine space; when there is degeneracy in some M∗P , it then becoms an uncountable set. Certainly, initializations sufficiently nearM∗ are also qualified by continuity. Remark C.5. Motivated by the idea of weights degeneracy (see Definition C.1), we can further apply similar methods to a (global) landscape analysis on the loss function Jω . The results there show that the plateaus are all over the landscape, even provided general targets (without memory structures). See details in Appendix D.

heading: C.3 NUMERICAL EXPERIMENTS

heading: C.3.1 MOTIVATING TESTS
text: In this section we give details of Figure 1.
(1) Learning linear functionals using linear RNNs (with GD optimizer) The target functional is HT (x) = ∫ T
0 ρ(t)xT−tdt with white noise x, while the representation ρ is
selected as the exponential sum or the scaled Airy function:
1. Exponential sum: ρ(t) = [c∗]>eW ∗tb∗, where c∗, b∗ are standard normal random vectors
of m∗ dimensions, and W = −I − Z>Z with Z ∈ Rm∗×m∗ is a Gaussian random matrix with i.i.d. entries having variance 1/m∗.
2. Airy function: ρ(t) = Ai(s0[t − t0]), where Ai(t) is the Airy function of the first kind, given by the improper integral
Ai(t) = 1
π lim ξ→∞ ∫ ξ 0 cos ( u3 3 + tu ) du. (117)
Note that in the first example, the memory of target functional decays quickly. However, for the second example, the effective rate of decay is controlled by the parameter t0. For t ≤ t0, the Airy function is oscillatory. Hence for large t0, a large amount of memory is present in the target. In Figure 1 we set m∗ = 8 for exponential sums, and t0 = 3, s0 = 2.25 for Airy functions.
In Figure 1 (a) and (b), we plot the gradient descent dynamics on training the linear RNNs (discretized using Euler method, hence equivalent to residual RNNs). We observe an efficient training process for the exponential sum case, while “plateauing” behaviors in the Airy function case. This causes a severe slow down of training. In addition, we also find that the plateauing effect gets worse as t0 (or s0) is increased, which corresponds to more complex Airy functions in the sense of more memory effects. That is, the long-term memory adversely affects the optimization process via gradient descent.
(2) Learning nonlinear functionals using nonlinear RNNs (with Adam optimizer)
To show that the plateauing behavior may be generic, we also consider a nonlinear forced dynamical system target, the Lorenz 96 system (Lorenz, 1996):
ẏ = −y + x+ K∑ k=1 zk/K,
żk = 2[zk+1(zk−1 − zk+2)− zk + y], k = 1, · · · ,K,
(118)
with cyclic indices zk+K = zk, and x is an external stochastic noise. When the unresolved variables zk are unknown, the dynamics of the resolved variable y driven by x is a nonlinear dynamical system with memory effects (but not a linear functional). We use a standard nonlinear RNN (with the tanh activation) to learn the sequence-to-sequence mapping x0:T 7→ y0:T with the Adam optimizer. Figure 1 (c) shows that the training of the Lorenz 96 system with the presence of memory also exhibits the plateauing phenomenon.
In all cases, the trained RNN model has a hidden dimension of 16 and the total length of the path is T = 6.4. The continuous-time RNN is discretized using the Euler method with step size 0.1.

heading: C.3.2 LONG-TERM MEMORY SIGNIFICANTLY CONTRIBUTES TO PLATEAUS
text: Recall the simple example of target with memory
ρ(t) = a∗e−w ∗t + c0e
− (t−1/ω) 2
2σ2 , w∗ > 0. (119)
We aim to learn ρ with the exponential sum
ρ̂m(t; a,w) = m∑ i=1 aie −wit, a ∈ Rm, w ∈ Rm+ . (120)
That is, the simplified linear RNN model with a diagonal recurrent kernel (see (76) and Appendix C.1). The optimization is performed by the gradient flow (21), i.e. a continuous-time idealization of the gradient descent dynamics. The ODE (21) is numerically solved by the AdamsBashforth-Moulton method. The results are illustrated in Figure 2. It is shown that the plateauing time increases rapidly as the memory 1/ω becomes longer.

heading: C.3.3 NUMERICAL VERIFICATIONS
text: (1) The timescale estimate
We first numerically verify the timescale proved in Theorem 5.1 (or Theorem C.1). That is, the time of plateauing (and also parameter separation ‖θ(τ) − θ(0)‖2) is exponentially large as the memory 1/ω → +∞. The results are shown in Figure 3, where we observe good agreement with the predicted scaling.
(2) General cases
To facilitate mathematical tractability, the analysis so far is done on the restrictive cases of the diagonal recurrent kernel W with negative entries, linear activations and the gradient flow training dynamics. However, we show here that the plateauing behavior - which we now understood as a generic feature of long-term memory of the target functional and its interaction with the optimization dynamics - is present even for general cases, and hence our simplified analytical setting is representative of the general situation.
In Figure 4, we still take the target functional as defined in (119), but apply more general models to learn it, including using full (non-diagonal) recurrent kernels of the RNN with no restrictions on
entries, using non-linear (tanh) activations and using the Adam optimizer. Furthermore, we do not take the Itô isometry simplification and instead use actual input sample paths of finite time horizons, just as one would do in training RNNs in practice. We observe that the plateauing behavior is present in all cases. Moreover, in the last case of Adam (which can be thought of as a momentumbased optimization method), the plateauing behavior is somewhat alleviated, although the separation of timescales is still present. This is consistent with our supplemental analysis in Appendix E, where we show that the momentum-based methods will speed up training based on our dynamical (Appendix C.2) and landscape analysis (Appendix D) of plateauing.

heading: D LANDSCAPE ANALYSIS
text: As mentioned in Remark C.5, we can perform a global landscape analysis on the loss function based on the idea of weights degeneracy, which arises from Definition C.1. Recall that the loss function reads
min a∈Rm,w∈Rm+ Jm(a,w) := ∫ ∞ 0 ( m∑ i=1 aie −wit − ρ(t) )2 dt. (121)
The main results of the appendix are summarized as follows.
• In Theorem D.1, we prove that the loss function has infinitely many critical points, which form a factorial number of affine spaces;
• In Theorem D.2, we prove that such (critical) affine spaces are much more than global minimizers provided the target being an exponential sum; 9
• In Theorem D.3, we prove that on such (critical) affine spaces, the Hessian is singular in the sense of processing multiple zero eigenvalues;
• In Proposition D.1, we prove that the (critical) affine spaces contain both saddles and degenerate stable points which are not global optimal.
Instead of a local dynamical analysis in Appendix C.2, we generalize similar methods to a global landscape analysis here, and the results hold for the loss function associated with general targets. More specifically, these results complement our main results (see Theorem 5.1 or Theorem C.1) in the following aspects.
• It is shown that the weights degeneracy is quite common in the whole landscape of the loss function. Unfortunately, the weights degeneracy often worsens the landscape to a large extent;
• It is shown that the weights degeneracy leads to a large number of stable areas (i.e. critical affine spaces), but most of them contribute to non-global minimizers;
• It is shown that these stable areas can also be quite flat, which often connect with local plateaus;
• For the structure of these stable areas, there are both saddles and degenerate critical points (not global optimal). In certain regimes, even saddles can be rather difficult to escape from (see Theorem 5.1 or Theorem C.1).
As a consequence, the optimization problem of linear RNNs is globally and essentially difficult to solve.

heading: D.1 SYMMETRY ANALYSIS ON THE LANDSCAPE
text: This subsection consists of two parts: in Appendix D.1.1, we give main results provided the existence of weights degeneracy; in Appendix D.1.2, we give sufficient conditions to guarantee the existence. Since the key observation to utilize weights degeneracy is to notice the permutation symmetry of coordinates of gradients, we also called it “symmetry analysis”.

heading: D.1.1 GENERIC THEORY
text: We begin with the following definition, which describes the concept of weights degeneracy in a natural and rigorous way. Definition D.1. (coincided critical solutions and affine spaces) Let d ∈ N+ and 1 ≤ d ≤ m. We say (a,w) is a d-coincided critical solution of Jm, if ∇Jm(a,w) = 0, and w = (wi) ∈ Rm+ has d different components. The coincided critical affine spaces are defined as coincided critical solutions that form affine spaces.
9The global minimizers are distinct when the target is an exponential sum. Here we compare the number of (critical) affine spaces with the number of global minimizers (both of them are finite). When the target is not an exponential sum, the same conclusion holds if there are still finite number of global minimizers. See Remark D.3 in Appendix D.1.1 for details.
To guarantee the existence of such solutions, it is necessary to have the following definition. Definition D.2. (â, ŵ) ∈ Rm ⊗ Rm+ is called the non-degenerate global minimizer of Jm, if and only if
Jm(â, ŵ) = inf a∈Rm,w∈Rm+ Jm(a,w), (122)
and (â, ŵ) takes a non-degenerate form
âi 6= 0, ŵi 6= ŵj for i 6= j, i, j = 1, 2, · · · ,m. (123)
For convenience, we also define an index set
N := {n ∈ N+ : Jm has non− degenerate global minimizers for any m ≤ n} , (124)
which is used frequently in the following analysis. For any f ∈ L2[0,∞), let L[f ] be the Laplace transform of f , i.e. L[f ](s) = ∫∞ 0 e−stf(t)dt, s > 0. We begin with the following lemma. Lemma D.1. Assume that ρ is smooth and √ w |L[ρ](w)| → 0 as w → 0+ and w → ∞. Then we have 1 ∈ N and thus N 6= ∅.
Proof. We aim to show that there exists â 6= 0 and ŵ > 0, such that J1(â, ŵ) = inf
a∈R,w∈R+ J1(a,w). (125)
The basic idea is to limit the unbounded domain a ∈ R, w ∈ R+ to a compact set without effecting the minimization of J1(a,w). We have
min a,w>0 J1(a,w) = min w>0 min a
{ 1
2w · a2 − 2L[ρ](w) · a+ ‖ρ‖2L2[0,∞) } = min
w>0 min a
{ 1
2w
( a− 2wL[ρ](w) )2 + [ ‖ρ‖2L2[0,∞) − 2w(L[ρ](w)) 2 ]}
= min w>0
{ ‖ρ‖2L2[0,∞) − 2w(L[ρ](w)) 2 } = J1(a(w), w),
where a(w) := 2wL[ρ](w). Write h(w) := J1(a(w), w), then h(0+) = h(∞) = ‖ρ‖2L2[0,∞). Obviously h(w) < ‖ρ‖2L2[0,∞) for any w > 0, hence
min w>0 h(w) = min w∈[wlb,wub]
h(w), 0 < wlb < wub <∞,
which implies
min a,w>0 J1(a,w) = min w>0 J1(a(w), w) = min w∈[wlb,wub] J1(a(w), w).
That is to say, the minimization of J1(a,w) can be equivalently performed on a 2-dimensional smooth curve (w, a(w))w∈[wlb,wub], which is certainly a compact set. By continuity, J1(a,w) has global minimizers, say (â, ŵ). Obviously ŵ > 0 and â = a(ŵ) 6= 0 (since â = 0 implies J1(â, w) = ‖ρ‖2L2[0,∞), certainly not a minimum), which completes the proof.
Remark D.1. If the target is an exponential sum, i.e. ρ(t) = ∑m∗ j=1 a ∗ je −w∗j t, we know ρ is smooth and √ w |L[ρ](w)| → 0 as w → 0+ and w → +∞; hence 1 ∈ N by Lemma D.1, and thus
N 6= ∅. In fact, L[ρ](w) = ∑m∗ j=1 a∗j w+w∗j
implies that L[ρ](w) = O(1) when w → 0+, and L[ρ](w) = O(1/w) when w → +∞. Theorem D.1. Assume that N 6= ∅ with N defined as (124). Let M := sup N . Then for any
m, d ∈ N+, 1 ≤ d ≤ min{m,M}, there exists at least d! { m d } d-coincided critical affine spaces
of Jm, 10 where { m d } ∈ N+ is called the Stirling number of the second kind.
10The affine spaces degenerate to distinct points when d = m. For sufficient conditions to guarantee M > 1 (to avoid vacuous results), see Theorem D.6 and Remark D.7 in Appendix D.1.2.
Proof. (i) Existence. The key observation is the permutation symmetry of∇Jm: by (19), if ai = aj and wi = wj for some i 6= j, then ∂Jm∂ai = ∂Jm ∂aj and ∂Jm∂wi = ∂Jm ∂wj .
For any m, d ∈ N+, 1 ≤ d ≤ m, suppose that w = (wi) ∈ Rm+ has d different components. Then for any partition P: {1, · · · ,m} = ∪dj=1Ij with Ij1 ∩ Ij2 = ∅ for any j1 6= j2, j1, j2 = 1, · · · , d, define the affine space
MP,(b,v),(m,d) := (a,w) ∈ Rm ⊗ Rm+ : wi = vj for any i ∈ Ij , ∑ i∈Ij ai = bj , j = 1, · · · , d  for some (b, v) ∈ Rd⊗Rd+, where v has exactly d different components. Therefore, for any (a,w) ∈ MP,(b,v),(m,d), we have
Jm(a,w) = ∥∥∥∥∥∥ d∑ j=1 ∑ i∈Ij aie −wit − ρ(t) ∥∥∥∥∥∥ 2
L2[0,∞)
= ∥∥∥∥∥∥ d∑ j=1 bje −vjt − ρ(t) ∥∥∥∥∥∥ 2
L2[0,∞)
= Jd(b, v),
and similarly
∂Jm ∂ak (a,w) = 2 ∫ ∞ 0 e−vst  d∑ j=1 bje −vjt − ρ(t)  dt, k ∈ Is, s = 1, 2, · · · , d, ∂Jm ∂wk (a,w) = 2ak ∫ ∞ 0 (−t)e−vst  d∑ j=1 bje −vjt − ρ(t)
 dt, k ∈ Is, s = 1, 2, · · · , d. Notice that
∂Jd ∂bs (b, v) = 2 ∫ ∞ 0 e−vst  d∑ j=1 bje −vjt − ρ(t)  dt, s = 1, 2, · · · , d, ∂Jd ∂vs (b, v) = 2bs ∫ ∞ 0 (−t)e−vst  d∑ j=1 bje −vjt − ρ(t)
 dt, s = 1, 2, · · · , d, we have
∂Jm ∂ak (a,w) = ∂Jd ∂bs (b, v), bs ∂Jm ∂wk (a,w) = ak ∂Jd ∂vs (b, v), k ∈ Is, s = 1, 2, · · · , d. (126)
11 Since d ≤ min{m,M}, d ∈ N . In fact, for any k ∈ N+, if k /∈ N , there exists i ≤ k such that J(i) has no non-degenerate global minimizers, we have j /∈ N for any j ≥ i, hence M ≤ i − 1 ≤ k − 1. Hence M = ∞ implies N = N+ and M < ∞ implies M ∈ N , and both of them lead to d ∈ N . Therefore, Jd has non-degenerate global minimizers, i.e. there exists (b̂, v̂) ∈ Rd ⊗ Rd+ such that
Jd(b̂, v̂) = inf b∈Rd,v∈Rd+ Jd(b, v), (127)
and (b̂, v̂) takes a non-degenerate form
b̂i 6= 0, v̂i 6= v̂j for any i 6= j, i, j = 1, 2, · · · , d. (128)
By (127), we get ∇Jd(b̂, v̂) = 0. Combining with (126) and (128), we obtain ∇Jm(â, ŵ) = 0 for any (â, ŵ) ∈ MP,(b̂,v̂),(m,d), i.e. (â, ŵ) belongs to a d-coincided critical affine space. Note that the affine space is with the dimension ∑d j=1(|Ij | − 1) = m − d, since there are d linear equality constrains on the m-dimensional vector a. 11By considering the gradient flow dynamic of Jd instead of Jm, a model reduction (from m-dimensional to d-dimensional) is almost completed onMP,(b,v),(m,d), except for some trivial degenerate cases (e.g. ak = 0 or bs = 0).
(ii) Counting. By the structure of affine spaces discussed above, we can identify different affine spaces with respect to the partition P . For counting the number of different partitions P: {1, · · · ,m} = ∪dj=1Ij , it can be decomposed into the following two steps. First, partitioning a set of m labelled objects into d nonempty unlabelled subsets. By definition, the answer is the Stir-
ling number of the second kind { m d } . Second, assign each partition to I1, · · · , Id accordingly.
There are d! ways in total. Therefore, the number of d-coincided critical affine spaces is at least
d! { m d } . The proof is completed.
Combining Lemma D.1, Remark D.1 and Theorem D.1 gives the following theorem, which states that there are much more saddles and degenerate stable points which are not global optimal than global minimizers in the landscape (provided the target being an exponential sum). Theorem D.2. Fix any m ∈ N+ relatively large. Consider the loss Jm associated with the target being a non-degenerate exponential sum, i.e. ρ(t) = ∑m j=1 a ∗ je −w∗j t, where a∗j 6= 0 and w∗i 6= w∗j for any i 6= j, i, j = 1, · · · ,m. Assume that m ∈ N 12 with N defined in (124). Then in the landscape of Jm, the number of coincided critical affine spaces is at least Poly(m) times larger than the number of global minimizers.
Proof. (i) Global minimizers. Since the target is an exponential sum, we have Jm(a,w) ≥ 0 and Jm(ā
∗, w̄∗) = 0, where ā∗ = Pa∗ and w̄∗ = Pw∗ with P ∈ Rm×m to be some permutation matrix. Next we show Jm has no other global minimizers.
Suppose Jm(a,w) = 0, we have m∑ i=1 aie −wit − m∑ j=1 a∗je −w∗j t = 0, t ≥ 0. (129) It is easy to see that for any j = 1, · · · ,m, there exists i(j) such that wi(j) = w∗j . Otherwise, if wi 6= w∗j , i = 1, · · · ,m, by (83) or Lemma C.2, we have a∗j = 0, which is a contradiction. Notice that w∗i 6= w∗j for any i 6= j, different w∗j ’s will correspond to different wi’s, hence the correspondence is one-to-one. Therefore, let wi = w∗j(i), (129) can be rewritten as
0 = m∑ i=1 aie −w∗j(i)t − m∑ i=1 a∗j(i)e −w∗j(i)t = m∑ i=1 (ai − a∗j(i))e −w∗j(i)t, t ≥ 0.
Again by Lemma C.2, we have ai = a∗j(i). That is to say, Jm(a,w) = 0 implies a = Pa ∗ and w = Pw∗ with P ∈ Rm×m to be some permutation matrix. This gives m! global minimizers. (ii) Coincided critical affine spaces. Obviously N 6= ∅, and M = sup N ≥ m. According to
Theorem D.1, for any d ∈ N+, 1 ≤ d ≤ min{m,M} = m, we have at least d! { m d } d-coincided critical affine spaces of Jm. By (i), for any d ≤ m−1, there are no global minimizers in these affine spaces. Counting the total number
m−1∑ d=1 d! { m d } . (130)
(iii) Comparison. To give a bound between (130) and m!, we need an elementary recurrence{ m d } = d { m− 1 d } + { m− 1 d− 1 } .
• For d = m− 1, let pm := {
m m− 1
} , then
pm = (m− 1) { m− 1 m− 1 } + { m− 1 m− 2 } = (m− 1) + pm−1 = · · · = m(m− 1) 2 .
12Although the assumption m ∈ N seems strong, we will provide sufficient conditions to guarantee its validity in Appendix D.1.2. See an complement in Theorem D.7.
• For d = m− 2, let qm := {
m m− 2
} , then
qm = (m− 2) { m− 1 m− 2 } + { m− 1 m− 3 } = (m− 2)pm−1 + qm−1 = · · ·
= 1
24 [2(m− 2)(m− 1)(2m− 3) + 3(m− 2)2(m− 1)2].
Combining above gives
1
m! m−1∑ d=1 d! { m d } > 1 m! [(m− 1)!pm + (m− 2)!qm] = (m+ 1)(3m− 2) 24 ,
which is a quadratic polynominal on m. The proof is completed.
Remark D.2. We only take the last two terms of (130) for a lower bound, which is obviously rather loose. In principle, a Poly(m) bound with higher degrees can be similarly obtained. That is to say, on one hand, there are infinitely many critical points forming affine spaces in the landscape of Jm; on the other hand, we see that even if only counting the number of affine spaces, there are still much less global minimizers (provided the width m relatively large). Remark D.3. When the target ρ is not an exponential sum, it is straightforward to see Theorem D.2 still holds if there is a finite number (with the scale of no more than factorial) of global minimizers.
Now we get down to investigate∇2Jm on the above coincided critical affine spaces. It is shown that ∇2Jm is singular and can have multiple zero eigenvalues. Theorem D.3. Fix any m, d ∈ N+, 1 ≤ d ≤ m. On the d-coincided critical affine spaces (induced by non-degenerate global minimizers of Jd13) of Jm, ∇2Jm is with rank at most m + d, and hence has at least m− d zero eigenvalues.
Proof. A straightforward computation shows that, for k, l = 1, 2, · · · ,m, ∂2Jm ∂ak∂al (a,w) = 2 wk + wl , (131)
∂2Jm ∂ak∂wl (a,w) = −2al (wk + wl)2 , k 6= l, (132)
∂2Jm ∂ak∂wk (a,w) = −ak 2w2k + 2 ∫ ∞ 0 (−t)e−wkt ( m∑ i=1 aie −wit − ρ(t) ) dt. (133)
Let the induced d-coincided critical affine space be MP,(b̂,v̂),(m,d), as is derived in the proof of Theorem D.1. Since (b̂, v̂) is the non-degenerate global minimizer of Jd, we have∫ ∞
0
(−t)e−ŵkt (
m∑ i=1 âie −ŵit − ρ(t)
) dt = ∫ ∞ 0 (−t)e−v̂st  d∑ j=1 b̂je −v̂jt − ρ(t)  dt = 1
2b̂s
∂Jd ∂vs (b̂, v̂) = 0, k ∈ Is, s = 1, 2, · · · , d
for any (â, ŵ) ∈MP,(b̂,v̂),(m,d). This gives
∂2Jm ∂ak∂wk (â, ŵ) = −âk 2ŵ2k . (134)
Now we show that, for any i, j ∈ Is, i 6= j, s = 1, 2, · · · , d, ∇2Jm(â, ŵ)i,: = ∇2Jm(â, ŵ)j,:. In fact, for any k = 1, · · · ,m, let k ∈ Is′ , then by (131),
∂2Jm ∂ai∂ak (â, ŵ) = 2 ŵi + ŵk =
2
v̂s + v̂s′ ,
∂2Jm ∂aj∂ak (â, ŵ) = 2 ŵj + ŵk =
2
v̂s + v̂s′ .
13That is, the affine spaceMP,(b̂,v̂),(m,d). See details in the proof of Theorem D.1
For k 6= i and k 6= j, (132) gives ∂2Jm ∂ai∂wk (â, ŵ) = −2âk (ŵi + ŵk)2 = −2âk (v̂s + v̂s′)2 , ∂2Jm ∂aj∂wk (â, ŵ) = −2âk (ŵj + ŵk)2 = −2âk (v̂s + v̂s′)2 .
By (134), for k = i 6= j, ∂2Jm ∂ai∂wk (â, ŵ) = −âi 2ŵ2i = −âi 2v̂2s , ∂2Jm ∂aj∂wk (â, ŵ) = −2âi (ŵj + ŵi)2 = −âi 2v̂2s ,
and similarly for k = j 6= i, ∂2Jm ∂ai∂wk (â, ŵ) = −2âj (ŵi + ŵj)2 = −âj 2v̂2s , ∂2Jm ∂aj∂wk (â, ŵ) = −âj 2ŵ2j = −âj 2v̂2s .
That is to say, there are at most m + d different rows in the symmetric matrix ∇2Jm(â, ŵ) ∈ R2m×2m, hence rank ( ∇2Jm(â, ŵ) ) ≤ m + d. Therefore, the number of zero eigenvalues of
∇2Jm(â, ŵ) ≥ dim { x ∈ R2m : ∇2Jm(â, ŵ) · x = 0 } = 2m− rank(∇2Jm(â, ŵ)) ≥ m− d. The proof is completed.
Remark D.4. The bound in Theorem D.3 is not sharp, since the estimate on rank ( ∇2Jm ) here is loose as only rows with the same elements are considered. In practice (numerical tests), it is often observed that there are more zero eigenvalues of rank ( ∇2Jm ) on the coincided critical affine space MP,(b̂,v̂),(m,d). Remark D.5. Theorem D.3 shows that, there are local plateaus around the d-coincided critical affine spaces MP,(b̂,v̂),(m,d) for d ≤ m − 1. In addition, the 0-eigenspace of Jm is higherdimensional for smaller d, which may suggest that one can stuck on plateaus more easily.

heading: D.1.2 SUFFICIENT CONDITIONS
text: There is still a gap when connecting Theorem D.1 and Theorem D.2. That is, it is necessary to guarantee sup N relatively large, i.e. J1, J2, · · · , Jd all have non-degenerate global minimizers for d as large as possible. Motivated by Kammler (1979a), we can give some sufficient conditions by limiting the target ρ within a smaller function space, the so-called completely monotonic functions. Definition D.3. F ∈ C[0,∞] ∩ C∞(0,∞) is called completely monotonic, if and only if
(−1)nF (n)(t) ≥ 0, 0 < t <∞, n = 0, 1, · · · , and F (∞) = 0. Remark D.6. Several examples of completely monotonic functions:
• ρ(t) = 1/(1 + t)α for any α > 0;
• The non-degenerate exponential sum with positive coefficients
ρ(t) = m∗∑ j=1 a∗je −w∗j t, 0 ≤ w∗1 < · · · < w∗m∗ , a∗j > 0, j = 1, 2, · · · ,m∗.
Since the space of exponential sums is not close, we turn to consider the problem of finding a best approximation to a given ρ ∈ L2[0,∞) from the set
Vd(R+) := { ρ̂ ∈ Cd[0,∞) : [(D + w1) · · · (D + wd)]ρ̂ = 0 for some w1, · · ·wd ∈ R+ } (135)
with respect to the commonL2-norm, i.e. inf ρ̂∈Vd(R+) ‖ρ̂−ρ‖L2[0,∞),whereD denotes the common differential operator. Obviously Vd(R+) ⊂ L2[0,∞) and Vd(R+) ( Vd+1(R+) for any d ∈ N+. Kammler (1979a) proves the following theorem. Theorem D.4. Assume ρ ∈ L2[0,∞) to be completely monotonic. Then there exists a best approximation ρ̂0 to ρ in Vd(R+), i.e.
‖ρ̂0 − ρ‖L2[0,∞) = inf ρ̂∈Vd(R+) ‖ρ̂− ρ‖L2[0,∞). (136)
When ρ /∈ Vd(R+), any such best approximation admits a non-degenerate form
ρ̂0(t) = d∑ j=1 b̂je −v̂jt, 0 < v̂1 < · · · < v̂d, b̂j > 0, j = 1, 2, · · · , d, (137)
and satisfies the generalized Aigrain-Williams equations
L[ρ̂0](v̂j) = L[ρ](v̂j), j = 1, 2, · · · , d, (138) d
ds L[ρ̂0](s) ∣∣∣ s=v̂j = d ds L[ρ](s) ∣∣∣ s=v̂j , j = 1, 2, · · · , d. (139)
Note that (136) and (137) are pretty similar to Definition D.2, except for a different choice of hypothesis function space. Now we show a connection between these two problems. Theorem D.5. Assume ρ ∈ L2[0,∞) to be completely monotonic, and ρ /∈ Vd(R+) for some d ∈ N+. Then Jd has non-degenerate global minimizers (b̂, v̂) ∈ Rd ⊗ Rd+.
Proof. According to Theorem D.4, there exists a non-degenerate best approximation ρ̂0 to ρ from Vd(R+), i.e.
‖ρ̂0 − ρ‖L2[0,∞) = inf ρ̂∈Vd(R+) ‖ρ̂− ρ‖L2[0,∞), (140)
ρ̂0(t) = d∑ j=1 b̂je −v̂jt, 0 < v̂1 < · · · < v̂d, b̂j > 0, j = 1, 2, · · · , d. (141)
We aim to prove Jd(b̂, v̂) = inf b∈Rd,v∈Rd+ Jd(b, v). Define the following subsets of exponential sums
Vd(R+) := { ρ̂ : ρ̂(t) =
d∑ i=1 aie −wit, ai ∈ R, wi > 0
} ,
Vd,k(R+) := { ρ̂ ∈ Vd(R+) : w = (wi) has k different components } , 1 ≤ k ≤ d,
then we have inf b∈Rd,v∈Rd+ Jd(b, v) = inf ρ̂∈Vd(R+) ‖ρ̂ − ρ‖2L2[0,∞). It is straightforward to verify that
Vd(R+) = ⋃d k=1 Vd,k(R+), and Vd,k(R+) = Vk,k(R+) ( Vk(R+) for k = 1, · · · , d. By (140), we get
‖ρ̂0 − ρ‖2L2[0,∞) = inf ρ̂∈Vd(R+) ‖ρ̂− ρ‖2L2[0,∞) ≤ inf ρ̂∈Vd,d(R+) ‖ρ̂− ρ‖2L2[0,∞).
Since ρ̂0 ∈ Vd,d(R+), we have
Jd(b̂, v̂) = ‖ρ̂0 − ρ‖2L2[0,∞) = inf ρ̂∈Vd,d(R+) ‖ρ̂− ρ‖2L2[0,∞).
The last task is to show inf ρ̂∈Vd(R+) ‖ρ̂ − ρ‖L2[0,∞) = inf ρ̂∈Vd,d(R+) ‖ρ̂ − ρ‖L2[0,∞). In fact, for any ρ̂ ∈ Vk,k, ρ̂(t) = ∑k i=1 aie
−wit, let ã := (a1, · · · , ak, 0), w̃ := (w1, · · · , wk, 1 + max1≤i≤k wi), we get ρ̂(t) := ∑k+1 i=1 ãie −w̃it ∈ Vk+1,k+1, which implies Vk,k ⊂ Vk+1,k+1. Therefore,
inf ρ̂∈Vd(R+) ‖ρ̂− ρ‖L2[0,∞) = inf ρ̂∈ ⋃d k=1 Vd,k(R+) ‖ρ̂− ρ‖L2[0,∞) = min 1≤k≤d
{ inf
ρ̂∈Vd,k(R+) ‖ρ̂− ρ‖L2[0,∞) } = min
1≤k≤d
{ inf
ρ̂∈Vk,k(R+) ‖ρ̂− ρ‖L2[0,∞) } ≥ inf ρ̂∈Vd,d(R+) ‖ρ̂− ρ‖L2[0,∞),
which completes the proof.
Combining Theorem D.1 and Theorem D.5 immediately gives the following result.
Theorem D.6. Assume ρ ∈ L2[0,∞) to be completely monotonic, and ρ /∈ V1(R+). Let D := {d ∈ N+ : ρ /∈ Vd(R+)}, D0 := supD and write m′ := min{m,D0}. Then the total number of
coincided critical affine spaces of Jm is at least ∑m′ d=1 d! { m d } .
Proof. We have 1 ∈ D and thus D 6= ∅, D0 ≥ 1. Since Vd(R+) ( Vd+1(R+) for any d ∈ N+, we have D = {1, 2, · · · , D0} if D0 < ∞, and D = N+ if D0 = ∞.14 Both of them gives {1, 2, · · · ,m′} ⊂ D, i.e. ρ /∈ Vk(R+) for any k ≤ m′. By Theorem D.5, J(k) has non-degenerate global minimizers for any k ≤ m′, i.e. m′ ∈ N . According to Theorem D.1, for any d ∈ N+,
1 ≤ d ≤ m′ = min{m,m′} ≤ min{m,M}, there exists at least d! { m d } d-coincided critical affine
spaces of Jm. Sum over d gives the total number ∑m′ d=1 d! { m d } .
Remark D.7. Examples:
• Suppose the target is ρ(t) = 1/(1+t)α, α > 0, thenD = N+ andD0 =∞. The total number of coincided critical affine spaces of the corresponding Jm is at least ∑m d=1 d! { m d } .
• Suppose the target is an non-degenerate exponential sum with positive coefficients: ρ(t) =∑m j=1 a ∗ je −w∗j t, where a∗j > 0 and w ∗ i 6= w∗j for any i 6= j, i, j = 1, · · · ,m. Then
D = {1, 2, · · · ,m − 1} and D0 = m − 1. The total number of coincided critical affine spaces of the corresponding Jm is at least ∑m−1 d=1 d! { m d } , which is exactly (130).
An complement for Theorem D.2 is as follows.
Theorem D.7. Fix any m ∈ N+ relatively large. Consider the loss Jm associated with the target being a non-degenerate exponential sum with positive coefficients, i.e. ρ(t) = ∑m j=1 a ∗ je −w∗j t, where a∗j > 0 and w ∗ i 6= w∗j for any i 6= j, i, j = 1, · · · ,m. Then in the landscape of Jm, the number of coincided critical affine spaces is at least Poly(m) times larger than the number of global minimizers.
Proof. By Theorem D.2, we only need to show is m ∈ N . Since ρ ∈ L2[0,∞) is completely monotonic, and ρ /∈ Vk(R+) for any k ≤ m − 1, then by Theorem D.5, J(k) has non-degenerate global minimizers for any k ≤ m− 1, i.e. m− 1 ∈ N . The proof is completed by noticing that Jm obviously has non-degenerate global minimizers, e.g. (a∗, w∗).

heading: D.1.3 A LOW-DIMENSIONAL EXAMPLE
text: To further understand the structure of coincided critical affine spaces, we focus on a specific lowdimensional example here. That is
min a∈R2,w∈R2+ J2(a,w) = ∥∥∥∥∥ 2∑ i=1 aie −wit − ρ(t) ∥∥∥∥∥ 2
L2[0,∞)
,
with the target to be a non-degenerate exponential sum ρ(t) = ∑m∗ j=1 a ∗ je −w∗j t, where a∗j 6= 0 and w∗i 6= w∗j for any i 6= j, i, j = 1, · · · ,m∗. As we will show later, the coincided critical affine spaces of J2 contain both saddles and degenerate stable points which are not global optimal.
By Lemma D.1, Remark D.1 and Theorem D.1, we know the 1-coincided critical affine space of J2 exists, and it can be constructed by taking the non-degenerate global minimizer of J1, say (â, ŵ)
14In fact, Vd(R+) ( Vd+1(R+) for any d ∈ N+ implies if ρ /∈ Vd(R+), ρ /∈ Vk(R+) for any k ≤ d, i.e. d ∈ D ⇒ k ∈ D for any k ≤ d; otherwise, if ρ ∈ Vd(R+), ρ ∈ Vl(R+) for any l ≥ d, i.e. d /∈ D ⇒ l /∈ D for any l ≥ d.
with â 6= 0 and ŵ > 0. ThenM(â,ŵ),(2,1) := {(a1, â − a1, ŵ, ŵ) : a1 ∈ R} ∈ R4 is a line 15, and ∇J2(a1, â − a1, ŵ, ŵ) = 0 for any a1 ∈ R. Denote the Hessian of J2 on the lineM(â,ŵ),(2,1) by A(â,ŵ)(a1), i.e. A(â,ŵ)(a1) := ∇2J2(a1, â − a1, ŵ, ŵ). We investigate the landscape of J2 on the lineM(â,ŵ),(2,1) by analyzing the eigenvalue distribution of A(â,ŵ)(a1). Proposition D.1. Suppose m = m∗ = 2, and 0 < w∗1 < w∗2 . Let I1 := [0, â] and I2 := (−∞, 0) ∪ (â,+∞) 16. Then
1. If a∗1a ∗ 2 < 0, the minimal eigenvalue of A(â,ŵ)(a1) is 0 for any a1 ∈ I1, and negative for
any a1 ∈ I2;
2. If a∗1a ∗ 2 > 0 and w ∗ 2/w ∗ 1 < 2 + √ 3, the minimal eigenvalue of A(â,ŵ)(a1) is negative for
any a1 ∈ I1, and 0 for any a1 ∈ I2.
Proof. Write c(w) := ∑m∗ j=1 a ∗ j [ 1
2w(w+w∗j ) 2 − 1(w+w∗j )3
] , and a2 := â − a1. A straightforward
computation shows that
A(â,ŵ)(a1) =  1 ŵ
1 ŵ
−a1 2ŵ2
−a2 2ŵ2
1 ŵ
1 ŵ
−a1 2ŵ2
−a2 2ŵ2
−a1 2ŵ2 −a1 2ŵ2 a21 2ŵ3 + 4c(ŵ)a1
a1a2 2ŵ3
−a2 2ŵ2 −a2 2ŵ2
a1a2 2ŵ3
a22 2ŵ3 + 4c(ŵ)a2  . Considering the congruent transformation of A(â,ŵ)(a1), which does not affect the index of inertia:
A(â,ŵ)(a1) =  1 ŵ
1 ŵ
−a1 2ŵ2
−a2 2ŵ2
1 ŵ
1 ŵ
−a1 2ŵ2
−a2 2ŵ2
−a1 2ŵ2 −a1 2ŵ2 a21 2ŵ3 + 4c(ŵ)a1
a1a2 2ŵ3
−a2 2ŵ2 −a2 2ŵ2
a1a2 2ŵ3
a22 2ŵ3 + 4c(ŵ)a2

→  1 ŵ 0 0 0 0 0 0 0 0 0 a21
4ŵ3 + 4c(ŵ)a1 a1a2 4ŵ3
0 0 a1a24ŵ3 a22 4ŵ3 + 4c(ŵ)a2  , we see thatA(â,ŵ)(a1) has one positive eigenvalue 1/ŵ and one eigenvalue 0. What remains are the
eigenvalues of A′(â,ŵ)(a1) :=
[ a21
4ŵ3 + 4c(ŵ)a1 a1a2 4ŵ3
a1a2 4ŵ3
a22 4ŵ3 + 4c(ŵ)a2
] . To determine their signs, we
compute
det(A′(â,ŵ)(a1)) = a1(â− a1) · 4c(ŵ) ( â
4ŵ3 + 4c(ŵ)
) (142)
= 1
â2ŵ3 a1(â− a1) · âc(ŵ) ·
( â2 + 16ŵ3âc(ŵ) ) . (143)
So we need to analyze the sign of âc(ŵ) and â2 + 16ŵ3âc(ŵ) under different assumptions on (a∗, w∗).
(i) a∗1a ∗ 2 < 0. By the optimality condition of (â, ŵ) for J1, we have
â = 2ŵ m∗∑ j=1 a∗j ŵ + w∗j = 4ŵ2 m∗∑ j=1 a∗j (ŵ + w∗j ) 2 , (144)
and therefore
c(ŵ) = â
8ŵ3 − m∗∑ j=1 a∗j (ŵ + w∗j ) 3 .
15Here we omit the corresponding partition P since it is unique. 16Suppose â > 0 here without loss of generality. If â < 0, we let I1 := [â, 0] and I2 := (−∞, â)∪ (0,+∞)
and the same conclusions hold.
Write vj := w∗j /ŵ, j = 1, 2, we get 0 < v1 < v2, and
â = 2 m∗∑ j=1 a∗j 1 + vj = 4 m∗∑ j=1 a∗j (1 + vj)2 , ŵ3c(ŵ) = â 8 − m∗∑ j=1 a∗j (1 + vj)3 .
Therefore
8ŵ3âc(ŵ) = â2 − 8â m∗∑ j=1 a∗j (1 + vj)3
= 16 m∗∑ j=1 a∗j (1 + vj)2 2 − 16 m∗∑ j=1 a∗j 1 + vj · m∗∑ j=1 a∗j (1 + vj)3
= −16a∗1a∗2(v1 − v2)2
(1 + v1)3(1 + v2)3 (145)
> 0,
which gives âc(ŵ) > 0 and â2 + 16ŵ3âc(ŵ) > 8ŵ3âc(ŵ) > 0.
(ii) a∗1a ∗ 2 > 0, w ∗ 2/w ∗ 1 < 2 +
√ 3. By (145), âc(ŵ) < 0.
â2 + 16ŵ3âc(ŵ) = 3â2 − 16â m∗∑ j=1 a∗j (1 + vj)3
= 16 3 m∗∑ j=1 a∗j (1 + vj)2 2 − 2 m∗∑ j=1 a∗j 1 + vj · m∗∑ j=1 a∗j (1 + vj)3  =
a∗21 u41u 4 2
[ u42 + c 2u41 + 6cu 2 1u 2 2 − 2cu31u2 − 2cu1u32 ] (c = a∗2/a ∗ 1, uj := 1 + vj > 1)
= a∗21 u42
( s4 + c2 + 6cs2 − 2cs− 2cs3 ) (s = u2/u1 > 1)
= a∗21 u42
[ c2 − 2s(s2 − 3s+ 1)c+ s4 ] .
Since 4s2(s2− 3s+ 1)2− 4s4 = 4s2(s− 1)2[(s− 2)2− 3], and 1 < s = u2/u1 = (ŵ+w∗2)/(ŵ+ w∗1) < w ∗ 2/w ∗ 1 < 2 + √ 3, we get ∆c < 0. This implies c2 − 2s(s2 − 3s + 1)c + s4 > 0 and â2 + 16ŵ3âc(ŵ) > 0.
In both (i) and (ii), â2 + 16ŵ3âc(ŵ) > 0, which implies that there is at least one positive diagonal element of A′(â,ŵ)(a1) in a sufficiently small neighborhood of a1 = 0 and a1 = â. By the Rayleigh-Ritz Theorem and Weyl’s Theorem,A′(â,ŵ)(a1) has at least one positive eigenvalue in this neighborhood. However, by (142), det(A′(â,ŵ)(a1)) only changes the sign at a1 = 0 and a1 = â. This implies another eigenvalue of A′(â,ŵ)(a1) changes the sign at a1 = 0 and a1 = â accordingly. By different signs of âc(ŵ) derived in (i) and (ii), and (142), the proof is completed.
Remark D.8. From Proposition D.1, we see that there are both saddles and degenerate stable points of J2 on the critical affine spaces (line)M(â,ŵ),(2,1), and each of them in fact forms affine spaces (lines) respectively, but they are certainly not global minimizers. Therefore, the gradientbased algorithms can get stuck around this affine space, except that it meets saddles with negative eigenvalues of large magnitude.

heading: E MOMENTUM HELPS TRAINING: QUADRATIC EXAMPLES
text: In practice, it is often the case that training is trapped in some very flat regions (plateaus), where the loss function has rather small gradients and negative eigenvalues of Hessian. Now we illustrate the escape dynamics (from plateaus) via a simple quadratic example.
Consider the loss function f(x) = (x21 − x22)/2 with 0 < 1. We check the escaping performance for continuous-in-time analogs of two optimization algorithms: gradient decent (GD) and momentum (heavy ball) method.
(1) Gradient decent
Consider the gradient flow of f(x) with an initial value x0 = (δ, 1)>, where 0 < δ 1 and δ = O( ). Thus ‖∇f(x0)‖2 = O( ), and{
x′1(τ) = −x1(τ), x1(0) = δ x′2(τ) = x2(τ), x2(0) = 1
⇒ { x1(τ) = δe −τ
x2(τ) = e τ
⇒ f(x(τ)) = (δ2e−2τ − e2 τ )/2 =: `1(τ). It is easy to show that there are different timescales of `1(τ). In fact, when τ = O(1/ ), `1(τ) = O( 2)e−|O(1/ )| − e|O(1)| = O( ). However, when τ continuous to increase, say
τ ≥ 1 2 ln δ0 =: τ 1 , (146)
where δ0 > 0 denotes the gap satisfying = o(δ0), we get `1(τ) ≤ `1(τ 1) = O( 2) − e2 · 1 2 ln
δ0 /2 = O( 2)− δ0/2 < −δ0/4 for any τ ≥ τ 1 .
(2) Momentum
The momentum algorithm has the update rule
xk+1 = xk − η∇f(xk) + ρ(xk − xk−1), (147) where ρ ∈ R, η > 0 is the learning rate. The continuous-in-time analog can be derived as (similar to the arguments in Su et al. (2014))
0 = ρ xk+1 − 2xk + xk−1
η + (1− ρ) √ η xk+1 − xk√ η +∇f(xk)
≈ ρx′′(t) + (1− ρ)√ η x′(t) +∇f(x(t)),
with xk := x(k √ η) and the step size √ η of the simple finite differences 17. Let x1 = x0−η∇f(x0), we also get x′(0) = −√η∇f(x(0)).
To facilitate a comparison to GD, we take η = 1 18 and ρ = 1 19. Plugging the expression of f , we can solve the ODE
x′′(t) +∇f(x(t)) = 0⇔ { x′′1(τ) + x1(τ) = 0, x1(0) = δ, x ′ 1(0) = −δ
x′′2(τ)− x2(τ) = 0, x2(0) = 1, x′2(0) = ⇒ { x1(τ) = δ(cos τ − sin τ) x2(τ) = 1+ √ 2 e √ τ + 1− √ 2 e − √ τ
⇒ f(x(τ)) = 1 2
[ δ2(cos τ − sin τ)2 − ( 1 + √
2 e √ τ +
1− √
2 e− √ τ
)2] .
Write `2(τ) := f(x(τ)). It is not hard to show that there are still different timescales of `2(τ). In fact, when τ = O(1/ √ ), `2(τ) = O( 2)|O(1)| − |O(1)|(e|O(1)| + e−|O(1)|)2 = O( ). However, when τ continuous to increase, say
τ ≥ 1 2 √ ln 4δ0 =: τ 2 , (148)
we get `2(τ 2) = O( 2) − ( 1+
√
2 ) 2e2 √ τ/2 + O( ) = O( ) − δ0 < −δ0/2, hence `2(τ) ≤
O( )− δ0 < −δ0/2 for any τ ≥ τ 2 . Combining (1) and (2), we have the following conclusions.
17It is easy to check that the error of discretization is of order O( √ η). 18In the continuous-in-time analog of GD, i.e. the gradient flow, the step size is taken as 1. 19As is seen later, ρ = 1 not only simplifies the analysis, but also helps to obtain the best acceleration.
• For both training dynamics, there are different timescales in the loss function. That is to say, relatively long time is needed to escape from the plateaus; • Comparing (146) and (148), we get different timescales of escaping: O (1/ · ln(1/ )) for
GD and O (1/ √ · ln(1/ )) for momentum. Just like the convex case, where momentum improves the convergence rate by weakening the dependence on condition number, we see momentum can also help to escape rather flat saddles.

Review summaries:
This paper studies approximation and optimization of linear RNNs for learning linear functions, from the perspective of the memory-properties of the temporal sequence. It shows that linear functionals can be approximated by a linear RNN, with the rate of approximation depending on the long-term memory of the process.  It also shows that the training dynamics slow down for certain linear functionals with long-term memory. 

The paper provides a theoretical examination of the challenge of fitting recurrent neural networks (RNNs) to fit processes with long memory (or long-range dependence). Dubbed the “curse of memory”, the author(s) restrict to the case of linear activation functions, and show that for processes with increased spatial dependence: (a) the width of the layers must increase exponentially to _guarantee_ accurate approximations under a provided bound, and (b) a gradient-based optimization algorithm will take exponentially more time to converge. Sufficient details for reproducing the experiments are provided. 

This paper reports a mathematical study of approximation properties of linear RNNs. The first part  reports a universal approximation theorem, and presents an analysis of how efficient the approximation is. In particular, it is shown that approximating a slowly decaying, power-law temporal filter requires a large number of neurons, a property the authors refer to as the "curse of memory".