Year: 2018

Abstract:
We study discrete time dynamical systems governed by the state equation ht+1 = φ(Aht +But). HereA,B are weight matrices, φ is an activation function, and ut is the input data. This relation is the backbone of recurrent neural networks (e.g. LSTMs) which have broad applications in sequential learning tasks. We utilize stochastic gradient descent to learn the weight matrices from a finite input/state trajectory {ut,ht}t=0. We prove that SGD estimate linearly converges to the ground truth weights while using near-optimal sample size. Our results apply to increasing activations whose derivatives are bounded away from zero. The analysis is based on i) a novel SGD convergence result with nonlinear activations and ii) careful statistical characterization of the state vector. Numerical experiments verify the fast convergence of SGD on ReLU and leaky ReLU in consistence with our theory.

heading: 1 INTRODUCTION
text: A wide range of problems involve sequential data with a natural temporal ordering. Examples include natural language processing, time series prediction, system identification, and control design, among others. State-of-the-art algorithms for sequential problems often stem from dynamical systems theory and are tailored to learn from temporally dependent data. Linear models and algorithms; such as Kalman filter, PID controller, and linear dynamical systems, have a long history and are utilized in control theory since 1960’s with great success (Brown et al. (1992); Ho & Kalman (1966); Åström & Hägglund (1995)). More recently, nonlinear models such as recurrent neural networks (RNN) found applications in complex tasks such as machine translation and speech recognition (Bahdanau et al. (2014); Graves et al. (2013); Hochreiter & Schmidhuber (1997)). Unlike feedforward neural networks, RNNs are dynamical systems that use their internal state to process inputs. The goal of this work is to shed light on the inner workings of RNNs from a theoretical point of view. In particular, we focus on the RNN state equation which is characterized by a nonlinear activation function φ, state weight matrixA, and input weight matrixB as follows
ht+1 = φ(Aht +But), (1.1)
Here ht is the state vector and ut is the input data at timestamp t. This equation is the source of dynamic behavior of RNNs and distinguishes RNN from feedforward networks. The weight matrices A and B govern the dynamics of the state equation and are inferred from data. We will explore the statistical and computational efficiency of stochastic gradient descent (SGD) for learning these weight matrices.
Contributions: Suppose we are given a finite trajectory of input/state pairs (ut,ht)Nt=0 generated from the state equation (1.1). We consider a least-squares regression obtained from N equations; with inputs (ut,ht)Nt=1 and outputs (ht+1) N t=1. For a class of activation functions including leaky ReLU and for stable systems1, we show that SGD linearly converges to the ground truth weight matrices while requiring near-optimal trajectory length N . In particular, the required sample size is O(n+ p) where n and p are the dimensions of the state and input vectors respectively. The results are extended to unstable systems when the samples are collected from multiple independent RNN trajectories rather than a single trajectory. Our theory applies to increasing activation functions
1Throughout this work, a system is called stable if the spectral norm of the state matrixA is less than 1.
whose derivatives are bounded away from zero, which includes leaky ReLU, and Gaussian input data. Numerical experiments on ReLU and leaky ReLU corroborate our theory and demonstrate that SGD converges faster as the activation slope increases. To obtain our results, we i) characterize the statistical properties of the state vector (e.g. well-conditioned covariance) and ii) derive a novel SGD convergence result with nonlinear activations; which may be of independent interest. As a whole, this paper provides a step towards foundational understanding of RNN training via SGD.

heading: 1.1 RELATED WORK
text: Our work is related to the recent optimization and statistics literature on linear dynamical systems (LDS) and neural networks.
Linear dynamical systems: The state-equation (1.1) reduces to a LDS when φ is the linear activation (φ(x) = x). Identifying the weight matrices is a core problem in linear system identification and is related to the optimal control problem (e.g. linear quadratic regulator) with unknown system dynamics. While these problems are studied since 1950’s (Ljung (1998; 1987); Åström & Eykhoff (1971)), our work is closer to the recent literature that provides data dependent bounds and characterize the non-asymptotic learning performance. Recht and coauthors have a series of papers exploring optimal control problem (Simchowitz et al. (2018); Tu et al. (2018; 2017)). In particular, Hardt et al. (2016) shows gradient descent learns single-input-single-output (SISO) LDS with polynomial guarantees. Oymak & Ozay (2018) and Faradonbeh et al. (2018) provide sample complexity bounds for learning LDS. Sanandaji et al. (2011b;a); Pereira et al. (2010) study the identification of sparse systems.
Neural networks: There is a growing literature on the theoretical aspects of deep learning and provable algorithms for training neural networks. Most of the existing results are concerned with feedforward networks. Ge et al. (2017); Li & Yuan (2017); Mei et al. (2018b); Soltanolkotabi (2017); Janzamin et al. (2015); Zhong et al. (2017b) consider learning fully-connected shallow networks with gradient descent. Mei et al. (2018a); Soltanolkotabi et al. (2017); Foster et al. (2018) analyze empirical landscape of related nonlinear learning problems. Brutzkus & Globerson (2017); Zhong et al. (2017a); Du et al. (2017); Goel et al. (2018) address convolutional neural networks; which is an efficient weight-sharing architecture. Brutzkus et al. (2017); Wang et al. (2018) studies over-parameterized networks when data is linearly separable. Janzamin et al. (2015); Oymak & Soltanolkotabi (2018) utilize tensor decomposition techniques for learning feedforward neural nets. For recurrent networks, Sedghi & Anandkumar (2016) proposed tensor algorithms with polynomial guarantees and Khrulkov et al. (2017) studied their expressive power. More recently, Miller & Hardt (2018) showed that stable RNNs can be approximated by feed-forward networks.

heading: 2 PROBLEM SETUP
text: We first introduce the notation. ‖ · ‖ returns the spectral norm of a matrix and smin(·) returns the minimum singular value. The activation φ : R → R applies entry-wise if its input is a vector. Throughout, φ is assumed to be a 1-Lipschitz function. With proper scaling of its parameters, the system (1.1) with a Lipschitz activation can be transformed into a system with 1-Lipschitz activation. The functions Σ[·] and var[·] return the covariance of a random vector and variance of a random variable respectively. In is the identity matrix of size n × n. Normal distribution with mean µ and covariance Σ is denoted by N (µ,Σ). Throughout, c, C, c0, c1, . . . denote positive absolute constants.
Setup: We consider the dynamical system parametrized by an activation function φ(·) and weight matricesA ∈ Rn×n,B ∈ Rn×p as described in (1.1). Here, ht is the n dimensional state-vector and ut is the p dimensional input to the system at time t. As mentioned previously, (1.1) corresponds to the state equation of a recurrent neural network. For most RNNs of interest, the state ht is hidden and we only get to interact with ht via an additional output equation. For Elman networks Elman (1990), this equation is characterized by some output activation φy and output weights C,D as follows
yt = φy(Cht +Dut). (2.1)
In this work, our attention is restricted to the state equation (1.1); which corresponds to setting yt = ht+1 in the output equation. To analyze (1.1) in a non-asymptotic data-dependent setup, we assume a finite input/state trajectory of {ut,ht}Nt=0 generated by some ground truth weight matrices
Algorithm 1 Learning state equations with nonlinear activations
1: Inputs: (yt,ht,ut)Nt=1 sampled from a trajectory. Scaling µ, learning rate η. Initialization A0,B0. 2: Outputs: Estimates Â, B̂ of the weight matricesA,B. 3: xt ← [µhTt uTt ]T for 1 ≤ t ≤ N . 4: Θ0 ← [µ−1A0 B0] 5: for τ from 1 to END do 6: Pick γτ from {1, 2, . . . , N} uniformly at random. 7: Θτ ← Θτ−1 − η∇Lγτ (Θτ−1) 8: end for 9: return [Â B̂]← ΘEND [ µIn 0
0 Ip
] .
(A,B). Our goal is learning the unknown weightsA andB in a data and computationally efficient way. In essence, we will show that, if the trajectory length satisfies N & n+ p, SGD can quickly and provably accomplish this goal using a constant step size.
Appoach: Our approach is described in Algorithm 1. It takes two hyperparameters; the scaling factor µ and learning rate η. Using the RNN trajectory, we constructN triples of the form {ut,ht,ht+1}Nt=1. We formulate a regression problem by defining the output vector yt, input vector xt, and the target parameter C as follows
yt = ht+1 , xt = [ µht ut ] ∈ Rn+p , C = [µ−1AB] ∈ Rn×(n+p). (2.2)
With this reparameterization, we find the input/output identity yt = φ(Cxt). We will consider the least-squares regression given by
L(Θ) = 1 N N∑ t=1 Lt(Θ) where Lt(Θ) = 1 2 ‖yt − φ(Θxt)‖2`2 . (2.3)
For learning the ground truth parameter C, we utilize SGD on the loss function (2.3) with a constant learning rate η. Starting from an initial point Θ0, after END SGD iterations, Algrorithm 1 returns an estimate Ĉ = ΘEND. Estimates of A and B are decoded from the left and right submatrices of Ĉ respectively.

heading: 3 MAIN RESULTS

heading: 3.1 PRELIMINARIES
text: The analysis of the state equation naturally depends on the choice of the activation function; which is the source of nonlinearity. We first define a class of Lipschitz and increasing activation functions.
Definition 3.1 (β-increasing activation). Given 1 ≥ β ≥ 0, the activation function φ satisfies φ(0) = 0 and 1 ≥ φ′(x) ≥ β for all x ∈ R.
Our results will apply to strictly increasing activations where φ is β-increasing for some β > 0. Observe that, this excludes ReLU activation which has zero derivative for negative values. However, it includes Leaky ReLU which is a generalization of ReLU. Parameterized by 1 ≥ β ≥ 0, Leaky ReLU is a β-increasing function given by
LReLU(x) = max(βx, x). (3.1)
In general, given an increasing and 1-Lipschitz activation φ, a β-increasing function φβ can be obtained by blending φ with the linear activation, i.e. φβ(x) = (1− β)φ(x) + βx. A critical property that enables SGD is that the state-vector covariance Σ[ht] is well-conditioned under proper assumptions. The lemma below provides upper and lower bounds on this covariance matrix in terms of problem variables.
Lemma 3.2 (State vector covariance). Consider the state equation (1.1) where h0 = 0 and ut i.i.d.∼ N (0, Ip). Define the upper bound term Bt as
Bt = ‖B‖ √ 1− ‖A‖2t 1− ‖A‖2 . (3.2)
• Suppose φ is 1-Lipschitz and φ(0) = 0. Then, for all t ≥ 0, Σ[ht] B2t In.
• Suppose φ is a β-increasing function and p ≥ n. Then, Σ[ht] β2smin(B)2In.
As a natural extension from linear dynamical systems, we will say the system is stable if ‖A‖ < 1 and unstable otherwise. For activations we consider, stability implies that if the input is set to 0, state vector ht will exponentially converge to 0 i.e. the system forgets the past states quickly. This is also the reason (Bt)t≥0 sequence converges for stable systems and diverges otherwise. The condition number of the covariance will play a critical role in our analysis. Using Lemma 3.2, this number can be upper bounded by ρ defined as
ρ =
( B∞
βsmin(B)
)2 = ( ‖B‖
smin(B)
)2 1
β2(1− ‖A‖2) . (3.3)
Observe that, the condition number ofB appears inside the ρ term.

heading: 3.2 LEARNING FROM SINGLE TRAJECTORY
text: Our main result applies to stable systems (‖A‖ < 1) and provides a non-asymptotic convergence guarantee for SGD in terms of the upper bound on the state vector covariance. This result characterizes the sample complexity and the rate of convergence of SGD; and also provides insights into the role of activation function and the spectral norm ofA. Theorem 3.3 (Main result). Let {ut,ht+1}Nt=1 be a finite trajectory generated from the state equation (1.1). Suppose ‖A‖ < 1, φ is β-increasing, h0 = 0, p ≥ n, and ut
i.i.d.∼ N (0, Ip). Let ρ be same as (3.3) and c, C, c0 be properly chosen absolute constants. Pick the trajectory length N to satisfy N ≥ CLρ2(n+ p), where L = 1− log(cnρ)log ‖A‖ . Pick scaling µ = 1/B∞, learning rate η = c0 β2
ρn(n+p) , and consider the loss function (2.3). With probability 1 − 4N exp(−100n) − 8L exp(−O( NLρ2 )), starting from an initial point Θ0, for all τ ≥ 0, the SGD iterations described in Algorithm 1 satisfies
E[‖Θτ −C‖2F ] ≤ (1− c0 β4
2ρ2n(n+ p) )τ‖Θ0 −C‖2F . (3.4)
Here the expectation is over the randomness of the SGD updates.
Sample complexity: Theorem 3.3 essentially requires N & (n+ p)/β4 samples for learning. This can be seen by unpacking (3.3) and ignoring the logarithmic L term and the condition number of B. Observe that O(n + p) growth achieves near-optimal sample size for our problem. Each state equation (1.1) consists of n sub-equations (one for each entry of ht+1). We collect N state equations to obtain a system of Nn equations. On the other hand, the total number of unknown parameters inA andB are n(n+ p). This implies Theorem 3.3 is applicable as soon as the problem is mildly overdetermined i.e. Nn & n(n+ p).
Computational complexity: Theorem 3.3 requires O(n(n + p) log 1ε ) iterations to reach εneighborhood of the ground truth. Our analysis reveals that, this rate can be accelerated if the state vector is zero-mean. This happens for odd activation functions satisfying φ(−x) = −φ(x) (e.g. linear activation). The result below is a corollary and requires ×n less iterations. Theorem 3.4 (Faster learning for odd activations). Consider the same setup provided in Theorem 3.3. Additionally, assume that φ is an odd function. Pick scaling µ = 1/B∞, learning rate η = c0 β 2
ρ(n+p) , and consider the loss function (2.3). With probability 1 − 4N exp(−100n) − 8L exp(−O( NLρ2 )), starting from an initial point Θ0, for all τ ≥ 0, the SGD iterations described in Algorithm 1 satisfies
E[‖Θτ −C‖2F ] ≤ (1− c0 β4
2ρ2(n+ p) )τ‖Θ0 −C‖2F , (3.5)
where the expectation is over the randomness of the SGD updates.
Another aspect of the convergence rate is the dependence on β. In terms of β, the SGD error (3.4) decays as (1−O(β8))τ . While it is not clear how optimal is the exponent 8, numerical experiments in Section 6 demonstrate that larger β indeed results in drastically faster convergence.

heading: 4 MAIN IDEAS AND PROOF STRATEGY
text: We first outline our high-level proof strategy for Theorem 3.3; which brings together ideas from statistics and optimization.
1. We first show that input data is well-behaved by proving that state-vector ht has a wellconditioned covariance as discussed in Lemma 3.2 and shown in Appendix B. The key idea is if φ is β-increasing, then the random input data ut provides sufficient excitation for the output state ht+1.
2. Even if individual samples are well-behaved, analyzing (2.3) is still challenging due to temporal dependencies between the samples. These dependencies prevent us from directly using statistical learning results that typically assume i.i.d. samples. We show that the dependency between samples at time t and t+ T decay exponentially fast in separation T (for stable systems). This is outlined in Appendix C.
3. This observation allows us to obtain nearly independent data by subsampling the original trajectory to get (hiT ,uiT )i≥0. Thanks to exponential decay, a logarithmically small T can be chosen to generate large subtrajectories of size N/T . Appendix D uses additional perturbation arguments to establish the well-behavedness of the overall data matrix.
4. To conclude, we obtain a deterministic result which establishes fast convergence result for β-increasing activations and well-behaved dataset. This is provided in Theorem 4.1 and proved in Appendix A.
The first three steps are related to the statistical nature of the problem which can be decoupled from the last step. Specifically, the last step derives a deterministic result that establishes the linear convergence of SGD for β-increasing functions. For linear convergence proofs, a typical strategy is showing the strong convexity of the loss function i.e. showing that, for some α > 0 and all points v,u, the gradient satisfies
〈∇L(v)−∇L(u),v − u〉 ≥ α‖v − u‖2`2 .
The core idea of our convergence result is that the strong convexity parameter of the loss function with β-increasing activations can be connected to the loss function with linear activations. In particular, recalling (2.3), set ylint = Cxt and define the linear loss to be
Llin(Θ) = 1 2N N∑ i=1 ‖ylint −Θxt‖2`2 .
Denoting the strong convexity parameter of the original loss by αφ and that of linear loss by αlin, we argue that αφ ≥ β2αlin; which allows us to establish a convergence result as soon as αlin is strictly positive. Next result is our SGD convergence theorem which follows from this discussion.
Theorem 4.1 (Deterministic convergence). Suppose a data set {xi,yi}Ni=1 is given; where output yi is related to input xi via yi = φ(〈xi,θ〉) for some θ ∈ Rn. Suppose β > 0 and φ is a β-increasing. Let γ+ ≥ γ− > 0 be scalars. Assume that input samples satisfy the bounds
γ+In 1
N N∑ i=1 xix T i γ−In , ‖xi‖2`2 ≤ B for all i.
Let {rτ}∞τ=0 be a sequence of i.i.d. integers uniformly distributed between 1 to N . Then, starting from an arbitrary point θ0, setting learning rate η = β2γ− γ+B
, for all τ ≥ 0, the SGD iterations for quadratic loss
θτ+1 = θτ − η(φ(xTrτθτ )− yrτ )φ ′(xTrτθτ )xrτ , (4.1)
satisfies the error bound
E[‖θτ − θ‖2`2 ] ≤ ‖θ0 − θ‖ 2 `2(1− β4γ2− γ+B )τ , (4.2)
where the expectation is over the random selection of the SGD iterations {rτ}∞τ=0.
This theorem provides a clean convergence rate for SGD for β-increasing activations and naturally generalizes standard results on linear regression which corresponds to β = 1. We remark that related results appear in the literature on generalized linear models. Kakade et al. (2011); Foster et al. (2018); Mei et al. (2018a) provide learning theoretic loss/gradient/hessian convergence results for isotonic regression, robust regression, and β-increasing activations. Goel et al. (2018) establishes a similar result for leaky ReLU activations under the assumption of symmetric input distribution and infinitely many samples (i.e. in population limit). Compared to these, we establish a deterministic linear convergence guarantee for SGD that works whenever the data matrix is full rank. We believe extensions to proximal gradient methods might be beneficial for high-dimensional nonlinear problems (e.g. sparse/low-rank approximation, manifold constraints Cai et al. (2010); Beck & Teboulle (2009); Oymak et al. (2018); Agarwal et al. (2010); Pereira et al. (2010)) and is left as a future work.
To derive our main results in Section 3, we need to address the first three steps outlined earlier and determine the conditions under which Theorem 4.1 is applicable to the data obtained from RNN state equation with high probability. Below we provide desirable characteristics of the state vector; which enables our statistical results. Assumption 1 (Well-behaved state vector). Let L > 1 be an integer. There exists positive scalars γ+, γ−, θ and an absolute constant C > 0 such that θ ≤ 3 √ n and the following holds
• Lower bound: Σ[hL−1] γ−In,
• Upper bound: for all t, the state vector satisfies
Σ[ht] γ+In , ‖ht − E[ht]‖ψ2 ≤ C √ γ+ and ‖E[ht]‖`2 ≤ θ √ γ+. (4.3)
Here ‖·‖ψ2 returns the subgaussian norm of a vector (see Def. 5.22 of Vershynin (2010)).
Assumption 1 ensures that covariance is well-conditioned, state vector is well-concentrated, and it has a reasonably small expectation. Our next theorem establishes statistical guarantees for learning the RNN state equation based on this assumption. Theorem 4.2 (General result). Let {ut,ht+1}Nt=1 be a lengthN trajectory of the state equation (1.1). Suppose ‖A‖ < 1, φ is β-increasing, h0 = 0, and ut
i.i.d.∼ N (0, Ip). Given scalars γ+ ≥ γ− > 0, set the condition number as ρ = γ+/γ−. For absolute constants C, c, c0 > 0, choose trajectory length N to satisfy
N ≥ CLρ2(n+ p) where L = d1− log (cnρ) log ‖A‖ e.
Suppose Assumption 1 holds with L, γ+, γ−, θ. Pick scaling to be µ = 1/ √ γ+ and learning rate to be η = c0 β 2
ρ(θ+ √ 2)2(n+p) . With probability 1− 4N exp(−100n)− 8L exp(−O( NLρ2 )), starting from Θ0, for all τ ≥ 0, the SGD iterations on loss (2.3) as described in Algorithm 1 satisfies
E[‖Θτ −C‖2F ] ≤ (1− c0 β4
2ρ2(θ + √ 2)2(n+ p) )τ‖Θ0 −C‖2F , (4.4)
where the expectation is over the randomness of SGD updates.
The advantage of this theorem is that, it isolates the optimization problem from the statistical properties of state vector. If one can prove tighter bounds on achievable (γ+, γ−, θ), it will immediately imply improved performance for SGD. In particular, Theorems 3.3 and 3.4 are simple corollaries of Theorem 4.2 with proper choices.
• Theorem 3.3 follows by setting γ+ = B2∞, γ− = β2smin(B)2, and θ = √ n. • Theorem 3.4 follows by setting γ+ = B2∞, γ− = β2smin(B)2, and θ = 0.

heading: 5 LEARNING UNSTABLE SYSTEMS
text: So far, we considered learning from a single RNN trajectory for stable systems (‖A‖ < 1). For such systems, as the time goes on, the impact of the earlier states disappear. In our analysis, this allows us to split a single trajectory into multiple nearly-independent trajectories. This approach will not work for unstable systems (A is arbitrary) where the impact of older states may be amplified over time. To address this, we consider a model where the data is sampled from multiple independent trajectories.
Suppose N independent trajectories of the state-equation (1.1) are available. Pick some integer T0 ≥ 1. Denoting the ith trajectory by the triple (h(i)t+1,h (i) t ,u (i) t )t≥0, we collect a single sample from each trajectory at time T0 to obtain the triple (h (i) T0+1 ,h (i) T0 ,u (i) T0
). To utilize the existing optimization framework (2.3); for 1 ≤ i ≤ N , we set,
(yi,hi,ui) = (h (i) T0+1 ,h (i) T0 ,u (i) T0 ). (5.1)
With this setup, we can again use the SGD Algorithm 1 to learn the weightsA andB. The crucial difference compared to Section 3 is that, the samples (yi,hi,ui)Ni=1 are now independent of each other; hence, the analysis is simplified. As previously, having an upper bound on the condition number of the state-vector covariance is critical. This upper bound can be shown to be ρ defined as
ρ = { ρ̄ if n > 1 ρ̄ 1−β
2|A|2 1−(β|A|)2T0 if n = 1
where ρ̄ = B2T0
β2smin(B)2 . (5.2)
The ρ̄ term is similar to the earlier definition (3.3); however it involves BT0 rather than B∞. This modification is indeed necessary since B∞ =∞ when ‖A‖ > 1. On the other hand, note that, B2T0 grows proportional to ‖A‖2T0 ; which results in exponentially bad condition number in T0. Our ρ definition remedies this issue for single-output systems; where n = 1 andA is a scalar. In particular, when β = 1 (e.g. φ is linear) ρ becomes equal to the correct value 12. The next theorem provides our result on unstable systems in terms of this condition number and other model parameters.
Theorem 5.1 (Unstable systems). Suppose we are given N independent trajectories (h(i)t ,u (i) t )t≥0 for 1 ≤ i ≤ N . Each trajectory is sampled at time T0 to obtain N samples (yi,hi,ui)Ni=1 where the ith sample is given by (5.1). Suppose the sample size satisfies
N ≥ Cρ2(n+ p)
where ρ is given by (5.2). Assume the initial states are 0, φ is β-increasing, p ≥ n, and ut i.i.d.∼ N (0, Ip). Set scaling µ = 1/ √ BT0 , learning rate η = c0 β2
ρn(n+p) , and run SGD over the equations described in (2.2) and (2.3). Starting from Θ0, with probability 1 − 2N exp(−100(n + p)) − 4 exp(−O(Nρ2 )), all SGD iterations satisfy
E[‖Θτ −C‖2F ] ≤ (1− c0 β4
2ρ2n(n+ p) )τ‖Θ0 −C‖2F ,
where the expectation is over the randomness of the SGD updates.

heading: 6 NUMERICAL EXPERIMENTS
text: We conducted experiments on ReLU and Leaky ReLU activations. Let us first describe the experimental setup. We pick the state dimension n = 50 and the input dimension p = 100. We choose the ground truth matrixA to be a scaled random unitary matrix; which ensures that all singular values of A are equal. B is generated with i.i.d.N (0, 1) entries. Instead of using the theoretical scaling choice, we determine the scaling µ from empirical covariance matrices outlined in Algorithm 2. Similar to our proof strategy, this algorithm equalizes the spectral norms of the input and state covariances to speed up convergence. We also empirically determined the learning rate and used η = 1/100 in all experiments.
2Clearly, any nonzero 1× 1 covariance matrix has condition number 1. However, due to subtleties in the proof strategy, we don’t use ρ = 1 for β < 1. Obtaining tighter bounds on the subgaussian norm of the state-vector would help resolve this issue.
Algorithm 2 Empirical hyperparameter selection.
1: Inputs: (ht,ut)Nt=1 sampled from a trajectory. 2: Outputs: Scaling µ. 3: Form the empirical covariance matrix Σh from {ht}Nt=1. 4: Form the empirical covariance matrix Σu from {ut}Nt=1. 5: return √ ‖Σu‖/‖Σh‖.
Evaluation: We consider two performance measures in the experiments. Let Ĉ be an estimate of the ground truth parameter C = [µ−1A B]. The first measure is the normalized error defined as ‖Ĉ −C‖2F /‖C‖2F . The second measure is the normalized loss defined as∑N
i=1 ‖yt − φ(Ĉxt)‖2`2∑N i=1 ‖yt‖2`2 .
In all experiments, we run Algorithm 1 for 50000 SGD iterations and plot these measures as a function of τ ; by using the estimate available at the end of the τ th SGD iteration for 0 ≤ τ ≤ 50000. Each curve is obtained by averaging the outcomes of 20 independent realizations.Our first experiments use N = 500; which is mildly larger than the total dimension n + p = 150. In Figure 1, we plot the Leaky ReLU errors with varying slopes as described in (3.1). Here β = 0 corresponds to ReLU and β = 1 is the linear model. In consistence with our theory, SGD achieves linear convergence and as β increases, the rate of convergence drastically improves3. The improvement is more visible for less stable systems driven byA with a larger spectral norm. In particular, while ReLU converges for small ‖A‖, SGD gets stuck before reaching the ground truth when ‖A‖ = 0.8. To understand, how well SGD fits the training data, in Figure 2a, we plotted the normalized loss for ReLU activation. For more unstable system (‖A‖ = 0.9), training loss stagnates in a similar fashion to the parameter error. We also verified that the norm of the overall gradient ‖∇L(Θτ )‖F continues to decay (where Θτ is the τ th SGD iterate); which implies that SGD converges before reaching a global minima. AsA becomes more stable, rate of convergence improves and linear rate is visible. Finally, to better understand the population landscape of the quadratic loss with ReLU activations, Figure 2b repeats the same ReLU experiments while increasing the sample size five times to N = 2500. For this more overdetermined problem, SGD converges even for ‖A‖ = 0.9; indicating that
• population landscape of loss with ReLU activation is well-behaved, • however ReLU problem requires more data compared to the Leaky ReLU for finding global
minima.
Overall, as predicted by our theory, experiments verify that SGD indeed quickly finds the optimal weight matrices of the state equation (1.1) and as the activation slope β increases, the convergence rate improves. 3Note that convergence becomes faster for larger β under the realizable model i.e. there exists a ground truth state equation with activation slope β that can fit the observed trajectory. This is consistent with the technical setup our results are proven. Also note that the data distribution in the experiments changes with the activation slope β. If the dataset is fixed and not realizable, the results may be different as we vary the slope β.

heading: 7 CONCLUSIONS
text: This work showed that SGD can learn the nonlinear dynamical system (1.1); which is characterized by weight matrices and an activation function. This problem is of interest for recurrent neural networks as well as nonlinear system identification. We showed that efficient learning is possible with optimal sample complexity and good computational performance. Our results apply to strictly increasing activations such as Leaky ReLU. We empirically showed that Leaky ReLU converges faster than ReLU and requires less samples; in consistence with our theory. We list a few unanswered problems that would provide further insights into recurrent neural networks. • Covariance of the state-vector: Our results depend on the covariance of the state-vector and requires it to be positive definite. One might be able to improve the current bounds on the condition number and relax the assumptions on the activation function. Deriving similar performance bounds for ReLU is particularly interesting. • Hidden state: For RNNs, the state vector is hidden and is observed through an additional equation (2.1); which further complicates the optimization landscape. Even for linear dynamical systems, learning the (A,B,C,D) system ((1.1), (2.1)) is a non-trivial task Ho & Kalman (1966); Hardt et al. (2016). What can be said when we add the nonlinear activations? • Classification task: In this work, we used normally distributed input and least-squares regression for our theoretical guarantees. More realistic input distributions might provide better insight into contemporary problems, such as natural language processing; where the goal is closer to classification (e.g. finding the best translation from another language).

heading: A DETERMINISTIC CONVERGENCE RESULT FOR SGD
text: Proof of Theorem 4.1. Given two distinct scalars a, b; define φ′(a, b) = φ(a)−φ(b) a−b . φ ′(a, b) ≥ β since φ is β-increasing. Definewτ to be the residualwτ = θτ − θ. Observing
φ(xTrτ θτ )− yrτ = φ ′(xTrτ θτ ,x T rτ θ)x T rτwτ ,
the SGD recursion obeys
‖wτ+1‖2`2 = ‖wτ − η(φ(x T rτ θτ )− yrτ )φ ′(xTrτ θτ )xrτ ‖ 2 `2 .
= ‖wτ − ηxrτφ ′(xTrτ θτ )φ ′(xTrτ θτ ,x T rτ θ)x T rτwτ‖ 2 `2 = ‖(I − ηGrτ )wτ‖ 2 `2
whereGrτ = xrτφ ′(xTrτ θτ )φ ′(xTrτ θτ ,x T rτ θ)x T rτ . Since φ is 1-Lipschitz and β-increasing,Grτ is a positivesemidefinite matrix satisfying
xrτx T rτ Grτ β 2xrτx T rτ ,
GTrτGrτ xrτx T rτxrτx T rτ Bxrτx T rτ .
Consequently, we find the following bounds in expectation
γ+In E[Grτ ] β 2γ−In, (A.1)
E[GTrτGrτ ] Bγ+In.
Observe that (A.1) essentially lower bounds the strong convexity parameter of the problem with β2γ−; which is the strong convexity of the identical problem with the linear activation (i.e. β = 1). However, we only consider strong convexity around the ground truth parameter θ i.e. we restricted our attention to (θ,θτ ) pairs. With this, wτ+1 can be controlled as,
E[‖wτ+1‖2`2 ] = E[‖(I − ηGrτ )wτ‖ 2 `2 ]
= ‖wτ‖2`2 − 2η E[w T τ Grτwτ ] + η 2 E[wTτ GTrτGrτwτ ] ≤ ‖wτ‖2`2(1− 2ηβ 2γ− + η 2Bγ+).
Setting η = β 2γ− γ+B , we find the advertised bound
E[‖wτ+1‖2`2 ] ≤ E[‖wτ‖ 2 `2 ](1− β4γ2− γ+B ).
Applying induction over the iterations τ , we find the advertised bound (4.2)
E[‖wτ‖2`2 ] ≤ ‖w0‖ 2 `2(1− β4γ2− γ+B )τ .
Lemma A.1 (Merging L splits). Assume matrices X(i) ∈ RNi×q are given for 1 ≤ i ≤ L. Suppose for all 1 ≤ i ≤ L, rows ofX(i) has `2 norm at most √ B and eachX(i) satisfies
γ+In X(i)
T X(i)
Ni γ−In.
Set N = ∑L i=1 Ni and form the concatenated matrix X =  X(1) X(2)
... X(L) . Denote ith row of X by xi. Then, for each i, ‖xi‖2`2 ≤ B and
γ+In XTX
N =
1
N N∑ i=1 xix T i γ−In.
Proof. The bound on the rows ‖xi‖`2 directly follows by assumption. For the remaining result, first observe thatXTX = ∑L i=1X (i)TX(i). Next, we have
Nγ+In = L∑ i=1 Niγ+In L∑ i=1 X(i) T X(i) L∑ i=1 Niγ−In = Nγ−In.
Combining these two yields the desired upper/lower bounds onXTX/N .

heading: B PROPERTIES OF THE NONLINEAR STATE EQUATIONS
text: This section characterizes the properties of the state vector ht when input sequence is normally distributed. These bounds will be crucial for obtaining upper and lower bounds for the singular values of the data matrix X = [x1 . . . xN ]
T described in (2.2). For probabilistic arguments, we will use the properties of subgaussian random variables. Orlicz norm provides a general framework that subsumes subgaussianity.
Definition B.1 (Orlicz norms). For a scalar random variable Orlicz-a norm is defined as
‖X‖ψa = sup k≥1
k−1/a(E[|X|k])1/k
Orlicz-a norm of a vector x ∈ Rp is defined as ‖x‖ψa = supv∈Bp ‖vTx‖ψa where Bp is the unit `2 ball. The subexponential norm is the Orlicz-1 norm ‖·‖ψ1 and the subgaussian norm is the Orlicz-2 norm ‖·‖ψ2 . Lemma B.2 (Lipschitz properties of the state vector). Consider the state equation (1.1). Suppose activation φ is 1-Lipschitz. Observe that ht+1 is a deterministic function of the input sequence {uτ}tτ=0. Fixing all vectors {ui}i6=τ (i.e. all except uτ ), ht+1 is ‖A‖t−τ‖B‖ Lipschitz function of uτ for 0 ≤ τ ≤ t.
Proof. Fixing {ui}i 6=τ , denote ht+1 as a function of uτ by ht+1(uτ ). Given a pair of vectors uτ ,u′τ using 1-Lipschitzness of φ, for any t > τ , we have
‖ht+1(uτ )− ht+1(u′τ )‖`2 ≤ ‖φ(Aht(uτ ) +But)− φ(Aht(u ′ τ ) +But)‖`2
≤ ‖A(ht(uτ )− ht(u′τ ))‖`2 ≤ ‖A‖‖ht(uτ )− ht(u′τ )‖`2 .
Proceeding with this recursion until t = τ , we find
‖ht+1(uτ )− ht+1(u′τ )‖`2 ≤ ‖A‖ t−τ‖hτ+1(uτ )− hτ+1(u′τ )‖`2
≤ ‖A‖t−τ‖φ(Ahτ +Buτ )− φ(Ahτ +Bu′τ )‖`2 ≤ ‖A‖t−τ‖B‖‖uτ − u′τ‖`2 .
This bound implies ht+1(uτ ) is ‖A‖t−τ‖B‖ Lipschitz function of uτ .
Lemma B.3 (Upper bound). Consider the state equation governed by equation (1.1). Suppose ut i.i.d.∼ N (0, Ip), φ is 1-Lipschitz, φ(0) = 0 and h0 = 0. Recall the definition (3.2) of Bt. We have the following properties
• ht is a Bt-Lipschitz function of the vector qt = [uT0 . . . uTt−1]T ∈ Rtp.
• There exists an absolute constant c > 0 such that ‖ht − E[ht]‖ψ2 ≤ cBt and Σ[ht] B2t In.
• ht satisfies
E[‖ht‖2`2 ] ≤ tr(BB T )
1− ‖A‖2t 1− ‖A‖2 ≤ min{n, p}B 2 t .
Also, there exists an absolute constant c > 0 such that for any m ≥ n, with probability 1 − 2 exp(−100m), ‖ht‖`2 ≤ c √ mBt.
Proof. i) Bounding Lipschitz constant: Observe that ht is a deterministic function of qt i.e. ht = f(qt) for some function f . To bound Lipschitz constant of f , for all (deterministic) vector pairs qt and q̂t, we find a scalar Lf satisfying,
‖f(qt)− f(q̂t)‖`2 ≤ Lf‖qt − q̂t‖`2 . (B.1)
Define the vectors, {ai}ti=0, as follows
ai = [û T 0 . . . û T i−1 u T i . . . u T t−1] T .
Observing that a0 = qt, at = q̂t, we write the telescopic sum,
‖f(qt)− f(q̂t)‖`2 ≤ t−1∑ i=0 ‖f(ai+1)− f(ai)‖`2 .
Focusing on the individual terms f(ai+1)− f(ai), observe that the only difference is the ui, ûi terms. Viewing ht as a function of ui and applying Lemma B.2,
‖f(ai+1)− f(ai)‖`2 ≤ ‖A‖ t−1−i‖B‖‖ui − ûi‖`2 .
To bound the sum, we apply the Cauchy-Schwarz inequality; which yields
|f(qt)− f(q̂t)| ≤ t−1∑ i=0 ‖A‖t−1−i‖B‖‖ui − ûi‖`2
≤ ( t−1∑ i=0 ‖A‖2(t−1−i)‖B‖2)1/2 ( t−1∑ i=0 ‖ui − ûi‖2`2) 1/2
︸ ︷︷ ︸ ‖qt−q̂t‖`2
≤
√ ‖B‖2(1− ‖A‖2t)
1− ‖A‖2 ‖qt − q̂t‖`2
= Bt‖qt − q̂t‖`2 . (B.2)
The final line achieves the inequality (B.1) with Lf = Bt hence ht is Bt Lipschitz function of qt.
ii) Bounding subgaussian norm: When ut i.i.d.∼ N (0, Ip), the vector qt is distributed asN (0, Itp). Since ht a Bt Lipschitz function of qt, for any fixed unit length vector v, αv := vTht = vT f(qt) is still Bt-Lipschitz function of qt. Hence, using Gaussian concentration of Lipschitz functions, αv satisfies
P(|αv − E[αv]| ≥ t) ≤ 2 exp(− t2
2B2t ).
This implies that for any v, αv − E[αv] is O(Bt) subgaussian Vershynin (2010). This is true for all unit v, hence using Definition B.1, the vector ht satisfies ‖ht − E[ht]‖ψ2 ≤ O(Bt) as well. Secondly, Bt-Lipschitz function of a Gaussian vector obeys the variance inequality var[αv] ≤ B2t (page 49 of Ledoux (2001)), which implies the covariance bound Σ[ht] B2t In. iii) Bounding `2-norm: To obtain this result, we first bound E[‖ht‖2`2 ]. Since φ is 1-Lipschitz and φ(0) = 0, we have the deterministic relation ‖ht+1‖`2 ≤ ‖Aht +But‖`2 . Taking squares of both sides, expanding the right hand side, and using the independence of ht,ut and the covariance information of ut, we obtain
E[‖ht+1‖2`2 ] ≤ E[‖Aht +But‖ 2 `2 ] = E[‖Aht‖ 2 `2 ] + E[‖But‖ 2 `2 ] (B.3)
≤ ‖A‖2 E[‖ht‖2`2 ] + tr(BB T ). (B.4)
Now that the recursion is established, expanding ht on the right hand side until h0 = 0, we obtain
E[‖ht+1‖2`2 ] ≤ t∑ i=0 ‖A‖2itr(BBT ) ≤ tr(BBT )1− ‖A‖ 2(t+1) 1− ‖A‖2 .
Now using the fact that tr(BBT ) ≤ rank(B)‖B‖2 ≤ min{n, p}‖B‖2, we find
E[‖ht+1‖`2 ] 2 ≤ E[‖ht+1‖2`2 ] ≤ min{n, p}B 2 t+1.
Finally, using the fact that ht is Bt-Lipschitz function and utilizing Gaussian concentration of qt ∼ N (0, Itp), we find
P(‖ht+1‖`2 − E[‖ht+1‖`2 ] ≥ t) ≤ exp(− t2
2B2t ).
Setting t = (c − 1) √ mBt for sufficiently large c > 0, we find P(‖ht‖`2 ≥ √ nBt + (c − 1) √ mBt) ≤ exp(−100m).
Lemma B.4 (Odd activations). Suppose φ is strictly increasing and obeys φ(x) = −φ(−x) for all x and h0 = 0. Consider the state equation (1.1) driven ut i.i.d.∼ N (0, Ip). We have that E[ht] = 0.
Proof. We will inductively show that {ht}t≥0 has a symmetric distribution around 0. Suppose the vector ht satisfies this assumption. Let S ⊂ Rn be a set. We will argue that P(ht+1 ⊂ S) = P(ht+1 ⊂ −S). Since φ is strictly increasing, it is bijective on vectors, and we can define the unique inverse set S′ = φ−1(S). Also since φ is odd, φ(−S′) = −S. Since ht,ut are independent and symmetric, we reach the desired conclusion as follows
P(ht+1 ⊂ S) = P(Aht +But ⊂ S′) = P(A(−ht) +B(−ut) ⊂ S′) (B.5) = P(Aht +But ⊂ −S′) = P(φ(Aht +But) ⊂ φ(−S′)) = P(ht+1 ⊂ −S). (B.6)
Theorem B.5 (State-vector lower bound). Consider the nonlinear state equation (1.1) with {ut}t≥0 i.i.d.∼ N (0, Ip). Suppose φ is a β-increasing function for some constant β > 0. For any t ≥ 1, the state vector obeys
Σ[ht] β2smin(BBT )In.
Proof. The proof is an application of Lemma B.7. The main idea is to write ht as sum of two independent vectors, one of which has independent entries. Consider a multivariate Gaussian vector g ∼ N (0,Σ). g is statistically identical to g1 + g2 where g1 ∼ N (0, smin(Σ)Id) and g2 ∼ N (0,Σ − smin(Σ)Id) are independent multivariate Gaussians.
Since But ∼ N (0,BBT ), setting Σ = BBT and smin = smin(Σ), we have that But ∼ g1 + g2 where g1, g2 are independent and g1 ∼ N (0, sminIn) and g2 ∼ N (0,Σ− sminIn). Consequently, we may write
But +Aht ∼ g1 + g2 +Aht.
For lower bound, the crucial component will be the g1 term; which has i.i.d. entries. Applying Lemma B.7 by setting x = g1 and y = g2 +Aht, and using the fact that ht, g1, g2 are all independent of each other, we find the advertised bound, for all t ≥ 0, via
Σ[ht+1] = Σ[φ(g1 + g2 +Aht)] β2sminIn.
The next theorem applies to multiple-input-single-output (MISO) systems whereA is a scalar andB is a row vector. The goal is refining the lower bound of Theorem B.5.
Theorem B.6 (MISO lower bound). Consider the setup of Theorem B.5 with single output i.e. n = 1. For any t ≥ 1, the state vector obeys
var[ht] ≥ β2‖B‖2`2 1− (β|A|)2t
1− β2|A|2 .
Proof. For any random variable X , applying Lemma B.7, we have var[φ(X)] ≥ β2var[X]. Recursively, this yields
var[ht+1] = var[φ(Aht +But)] ≥ β2var[Aht +But] = β2(|A|2var[ht] + ‖B‖2`2).
Expanding these inequalities till h0, we obtain the desired bound
var[ht] ≥ t∑ i=1 (βi|A|i−1‖B‖`2) 2.
Lemma B.7 (Vector lower bound). Suppose φ is a β-increasing function. Let x = [x1 . . . xn]T be a vector with i.i.d. entries distributed as xi ∼ X . Let y be a random vector independent of x. Then,
Σ[φ(x+ y)] β2var[X]In.
Proof. We first apply law of total covariance (e.g. Lemma B.8) to simplify the problem using the following lower bound based on the independence of x and y,
Σ[φ(x+ y)] Ey[Σ[φ(x+ y) ∣∣ y]] (B.7)
= Ey[Σx[φ(x+ y)]]. (B.8)
Now, focusing on the covariance Σx[φ(x + y)], fixing a realization of y, and using the fact that x has i.i.d. entries; φ(x+ y) has independent entries as φ applies entry-wise. This implies that Σx[φ(x+ y)] is a diagonal matrix. Consequently, its lowest eigenvalue is the minimum variance over all entries,
Σx[φ(x+ y)] min 1≤i≤n var[φ(xi + yi)]In.
Fortunately, Lemma B.9 provides the lower bound var[φ(xi + yi)] ≥ β2var[X]. Since this lower bound holds for any fixed realization of y, it still holds after taking expectation over y; which concludes the proof.
The next two lemmas are helper results for Lemma B.7 and are provided for the sake of completeness.
Lemma B.8 (Law of total covariance). Let x,y be two random vectors and assume y has finite covariance. Then Σ[y] = E[Σ[y ∣∣ x]] + Σ[E[y ∣∣ x]].
Proof. First, write Σ[y] = E[yyT ]− E[y]E[yT ]. Then, applying the law of total expectation to each term, Σ[y] = E[E[yyT ∣∣ x]]− E[E[y ∣∣ x]]E[E[yT ∣∣ x]]. Next, we can write the conditional expectation as E[E[yyT
∣∣ x]] = E[Σ[y ∣∣ x]] + E[E[y ∣∣ x]E[y ∣∣ x]]T . To conclude, we obtain the covariance of E[y
∣∣ x] via the difference, E[E[y
∣∣ x]E[y ∣∣ x]]T − E[E[y ∣∣ x]]E[E[yT ∣∣ x]] = Σ[E[y ∣∣ x]], which yields the desired bound.
Lemma B.9 (Scalar lower bound). Suppose φ is a β-increasing function with β > 0 as defined in Definition 3.1. Given a random variable X and a scalar y, we have
var[φ(X + y)] ≥ β2var[X].
Proof. Since φ is β-increasing, it is invertible and φ−1 is strictly increasing. Additionally, φ−1 is 1/β Lipschitz since, |φ(a)− φ(b)| ≥ β|a− b| =⇒ |a− b| ≥ β|φ−1(a)− φ−1(b)|. Using this observation and the fact that E[X] minimizes E(X − α)2 over α, var[φ(X + y)] can be lower bounded as follows
var[φ(X + y)] = E(φ(X + y)− E[φ(X + y)])2
≥ β2 E((X + y)− φ−1(E[φ(X + y)]))2
≥ β2 E(X + y − E[X + y])2
= β2 E(X − EX)2 = β2var[X]. Note that, the final line is the desired conclusion.

heading: C TRUNCATING STABLE SYSTEMS
text: One of the challenges in analyzing dynamical systems is the fact that samples from the same trajectory have temporal dependence. This section shows that, for stable systems, the impact of the past states decay exponentially fast and the system can be approximated by using the recent inputs only. We first define the truncation of the state vector. Definition C.1 (Truncated state vector). Suppose φ(0) = 0, initial condition h0 = 0, and consider the state equation (1.1). Given a timestamp t, L-truncation of the state vector ht is denoted by h̄t,L and is equal to qt where
qτ+1 = φ(Aqτ +Bu ′ τ ) , q0 = 0 (C.1)
is the state vector generated by the inputs u′τ satisfying
u′τ = { 0 if τ < t− L uτ else .
In words, L truncated state vector h̄t,L is obtained by unrolling ht until time t− L and setting the contribution of the state vector ht−L to 0. This way, h̄t,L depends only on the variables {uτ}t−1τ=t−L.
The following lemma states that impact of truncation can be made fairly small for stable systems (‖A‖ < 1). Lemma C.2 (Truncation impact – deterministic). Consider the state vector ht and its L-truncation h̄t,L from Definition C.1. Suppose φ is 1-Lipschitz. We have that
‖ht − h̄t,L‖`2 ≤ { 0 if t ≤ L ‖A‖L‖ht−L‖`2 else .
Proof. When t ≤ L, Definition C.1 implies u′τ = uτ hence ht = qt = h̄t,L. When t > L, we again use Definition C.1 and recall that u′τ = 0 until time τ = t− L− 1. For all t− L < τ ≤ t, using 1-Lipschitzness of φ, we have that
‖hτ − qτ‖`2 = ‖φ(Ahτ−1 +Buτ−1)− φ(Aqτ−1 +Buτ−1)‖`2 ≤ ‖(Ahτ−1 +Buτ−1)− (Aqτ−1 +Buτ−1)‖`2 ≤ ‖A(hτ−1 − qτ−1)‖`2 ≤ ‖A‖‖hτ−1 − qτ−1‖`2 .
Applying this recursion between t− L < τ ≤ t and using the fact that qt−L = 0 implies the advertised result ‖ht − qt‖`2 ≤ ‖A‖
L‖ht−L − qt−L‖`2 ≤ ‖A‖L‖ht−L‖`2 .

heading: C.1 NEAR INDEPENDENCE OF SUB-TRAJECTORIES
text: We will now argue that, for stable systems, a single trajectory can be split into multiple nearly independent trajectories. First, we describe how the sub-trajectories are constructed.
Definition C.3 (Sub-trajectory). Let sampling rate L ≥ 1 and offset 1 ≤ τ̄ ≤ L be two integers. Let N̄ = N̄τ̄ be the largest integer obeying (N̄ − 1)L + τ̄ ≤ N . We sample the trajectory {ht,ut}Nt=0 at the points τ̄ , τ̄ + L, . . . , τ̄ + (N̄ − 1)L+ τ̄ and define the τ̄ th sub-trajectory as
(h(i),u(i)) := (h(i,τ̄),u(i,τ̄)) = (h(i−1)L+τ̄ ,u(i−1)L+τ̄ ).
Definition C.4 (Truncated sub-trajectory). Consider the state equation (1.1) and recall Definition C.1. Given offset τ̄ and sampling rate L, for 1 ≤ i ≤ N̄ , the ith truncated sub-trajectory states are {h̄(i)}N̄i=1 where the ith state is defined as
h̄(i) = h̄L(i−1)+τ̄ ,L−1.
The truncated samples are independent of each other as shown in the next lemma.
Lemma C.5. Consider the truncated states of Definition C.4. If (1.1) is generated by independent vectors {ut}t≥0, for any offset τ̄ and sampling rate L, the vectors {h̄(i)}N̄i=1, {u(i)}N̄i=1 are all independent of each other.
Proof. By construction h̄(i) only depends on the vectors {uτ}L(i−1)+τ̄−1τ=L(i−2)+τ̄+1. Note that the dependence ranges [L(i− 2) + τ̄ + 1, L(i− 1) + τ̄ − 1] are disjoint intervals for different i’s; hence (h̄(i))N̄i=1 are independent of each other. To show the independence of u(i) and h̄(i); observe that inputs u(i) = uL(i−1)+τ̄ have timestamp τ̄ modulo L; which is not covered by the dependence range of (h̄(i))N̄i=1.
If the input is randomly generated, Lemma C.2 can be combined with a probabilistic bound on ht, to show that truncated states h̄(i) are fairly close to the actual states h(i).
Lemma C.6 (Truncation impact – random). Given offset τ̄ and sampling rate L, consider the state vectors of the sub-trajectory {h(i)}N̄i=1 and L− 1-truncations (h̄(i))N̄i=1. Suppose {ut}t≥0
i.i.d.∼ N (0, Ip), ‖A‖ < 1, h0 = 0, φ is 1-Lipschitz, and φ(0) = 0. Also suppose upper bound (4.3) of Assumption 1 holds for some θ ≤√ n, γ+ > 0. There exists an absolute constant c > 0 such that with probability at least 1− 2N̄ exp(−100n), for all 1 ≤ i ≤ N̄ , the following bound holds
‖h(i) − h̄(i)‖`2 ≤ c √ n‖A‖L−1√γ+.
In particular, we can always pick γ+ = B2∞ (via Lemma B.3).
Proof. Using Assumption 1, we can apply Lemma F.3 on vectors {h(i−2)L+τ̄+1}N̄i=1. Using a union bound, with desired probability, all vectors obey
‖h(i−2)L+τ̄+1 − E[h(i−2)L+τ̄+1]‖`2 ≤ (c− 1) √ nγ+,
for sufficiently large c. Since θ ≤ √ n, triangle inequality implies ‖h(i−2)L+τ̄+1‖`2 ≤ c √ nγ+. Now, applying Lemma C.2, for all 1 ≤ i ≤ N̄ , we find
‖h(i) − h̄(i)‖`2 = ‖h(i−1)L+τ̄ − h̄(i−1)L+τ̄ ,L−1‖`2 ≤ ‖A‖L−1‖h(i−2)L+τ̄+1‖`2 ≤ c‖A‖L−1√nγ+.

heading: D PROPERTIES OF THE DATA MATRIX
text: This section utilizes the probabilistic estimates from Section B to provide bounds on the condition number of data matrices obtained from the RNN trajectory (1.1). Following (2.2), these matricesH,U andX are defined as
H = [h1 . . . hN ] T , U = H = [u1 . . . uN ] T , X = [x1 . . . xN ] T . (D.1)
The challenge is that, the state matrix H has dependent rows; which will be addressed by carefully splitting the trajectory {ut,ht}Nt=0 into multiple sub-trajectories which are internally weakly dependent as discussed in Section C. We first define the matrices obtained from these sub-trajectories.
Definition D.1. Given sampling rate L and offset τ̄ , consider the L-subsampled trajectory {h(i),u(i)}N̄i=1 as described in Definitions C.3 and C.4. Define the matrices H̄ = H̄(τ̄) ∈ RN̄×n, H̃ = H̃(τ̄) ∈ RN̄×n, Ũ = Ũ (τ̄) ∈ RN̄×p, and X̃ = X̃(τ̄) ∈ RN̄×(n+p) as
H̄ = [h̄(1) . . . h̄(N̄)]T , H̃ = [h(1) . . . h(N̄)]T , Ũ = [u(1) . . . u(N̄)]T , X̃ = [µH̃ Ũ ].
Lemma D.2 (Handling perturbation). Consider the nonlinear state equation (1.1). Given sampling rate L > 0 and offset τ̄ , consider the matrices H̄, H̃, X̃ of Definition D.1 and let Q = [γ−1/2+ H̄ Ũ ] ∈ RN̄×(n+p). Suppose Assumption 1 holds, φ is β-increasing, and ut
i.i.d.∼ N (0, Ip). There exists an absolute constant C > 0 such that if N̄ ≥ C γ 2 +
γ2− (n + p), with probability 1 − 8 exp(−c γ 2 − γ2+ N̄), for all matrices M obeying
‖M − H̄‖ ≤ √ γ−N̄
10 , the perturbedQ matrices given by,
Q̃ = [γ −1/2 + M Ũ ], (D.2)
satisfy
(Θ + √ 2)2 Q̃ T Q̃
N̄ γ− 2γ+ . (D.3)
Proof. This result is a direct application of Theorem F.1 after determining minimum/maximum eigenvalues of population covariance. The cross covariance obeys E[H̄T Ũ ] = 0 due to independence. Also, for i > 1, the truncated state vector h̄(i) is statistically identical tohL−1 hence Σ[h̄(i)] γ−In. Consequently, Σ[u(i)] = Ip, 1 γ+ Σ[h̄(i)] In for all i and γ−γ+ In 1 γ+ Σ[h̄(i)] for all i > 1. Hence, setting qi = [ 1√ γ+ h̄(i)
u(i)
] , for all i > 1
γ− γ+ In Σ[qi] In.
Set the matrix Q̄ = [q2 . . . qN̄ ] T and note thatQ = [q1 Q̄T ]T . Applying Theorem F.1 on Q̄ and Corollary F.2 onQ, we find that, with the desired probability,
θ + √
3/2 ≥ 1√ N̄ ‖Q‖ ≥ 1√ N̄ smin(Q) ≥ 1√ N̄ smin(Q̄) ≥ √ N − 1 N √ 2γ− 3γ+ ≥ 0.99× √ 2γ− 3γ+ .
Setting E = M − H̄ and observing Q̃ = Q+ [γ−1/2+ E 0], the impact of the perturbation E can be bounded naively via smin(Q) − γ−1/2+ ‖E‖ ≤ smin(Q̃) ≤ ‖Q̃‖ ≤ ‖Q‖ + γ −1/2 + ‖E‖. Using the assumed bound on ‖E‖, this yields
θ + √
2 ≥ 1√ N̄ ‖Q̃‖ ≥ 1√ N̄ smin(Q̃) ≥ √ γ− 2γ+ .
This final inequality is identical to the desired bound (D.3).
Theorem D.3 (Data matrix condition). Consider the nonlinear state-equation (1.1). Given γ+ ≥ γ− > 0, define the condition number ρ = γ+
γ− . For some absolute constants c, C > 0, pick a trajectory length N where
L = d1− log (cnρ) log ‖A‖ e , N0 = b N L c ≥ Cρ2(n+ p),
and pick scaling µ = 1√ γ+
. Suppose ‖A‖ < 1, φ is β-increasing, ut i.i.d.∼ N (0, Ip), and Assumption
1 holds with γ+, γ−, θ, L. Matrix X = [x1 . . . xN ]T of (D.1) satisfies the following with probability 1− 4N exp(−100n)− 8L exp(−O(N0/ρ2)).
• Each row ofX has `2 norm at most c0 √ p+ n where c0 is an absolute constant.
• XTX obeys the bound
(Θ + √ 2)2In+p XTX
N ρ−1In+p/2. (D.4)
Proof. The first statement on `2-norm bound can be concluded from Lemma D.4 and holds with probability 1−2N exp(−100(n+p)). To show the second statement, for a fixed offset 1 ≤ τ̄ ≤ L, consider Definition D.1 and the matrices H̃(τ̄), Ũ (τ̄), X̃(τ̄). Observe thatX is obtained by merging multiple sub-trajectory matrices
{X̃(τ̄)}Lτ̄=1. We will first show the advertised bound for an individual X̃(τ̄) by applying Lemma D.2 and then apply Lemma A.1 to obtain the bound on the combined matrixX .
Recall that N̄τ̄ is the length of the τ̄ th sub-trajectory i.e. number of rows of X̃(τ̄). By construction 2N0 ≥ N̄τ̄ ≥ N0 for all 1 ≤ τ̄ ≤ L. Given 1 ≤ τ̄ ≤ L and triple H̄(τ̄), H̃(τ̄), Ũ (τ̄), setQ = [µH̄(τ̄) Ũ (τ̄)]. Since N0 is chosen to be large enough, applying Theorem D.2 with µ = 1/ √ γ+ choice, and noting ρ = γ+/γ−, we find that, with probability 1− 4 exp(−c1N0/ρ2), all matricesM satisfying ‖M − H̄(τ̄)‖ ≤ √ γ−N0/10 and Q̃ as in (D.2) obeys
(Θ + √ 2)2 Q̃ T Q̃
N ρ−1/2. (D.5)
Let us call this Event 1. To proceed, we will argue that with high probability ‖H̃(τ̄) − H̄(τ̄)‖ is small so that the bound above is applicable withM = H̃(τ̄) choice; which sets Q̃ = X̃(τ̄) in (D.5). Applying Lemma C.6, we find that, with probability 1− 2N̄τ̄ exp(−100n),
‖H̄(τ̄) − H̃(τ̄)‖ ≤ √ 2N0 max{‖h(i) − h̄(i)‖`2} ≤ c0 √ 2N0 √ nγ+‖A‖L−1.
Let us call this Event 2. We will show that our choice of L ensures right hand side is small enough and guarantees ‖H̄(τ̄) − H̃(τ̄)‖ ≤ √ γ−N0/10. Set c = max{200c20, 1}. Desired claim follows by taking logarithms of upper/lower bounds and cancelling out √ N0 terms as follows
c0 √ n‖A‖L−1√γ+ ≤ √ γ−/10 √ 2 ⇐⇒ (L− 1) log ‖A‖+ log√cnρ ≤ 0 (D.6)
⇐⇒ − log cnρ 2 log ‖A‖ ≤ L− 1 (D.7)
⇐= L = d1− log (cnρ) log ‖A‖ e. (D.8)
Here we use the fact that log ‖A‖ < 0 since ‖A‖ < 1 and cnρ ≥ 0. Consequently, both Event 1 and Event 2 hold with probability 1−4 exp(−c1N0/ρ2)−2N̄τ̄ exp(−100n), implying (D.5) holds with Q̃ = X̃(τ̄). Union bounding this over 1 ≤ τ̄ ≤ L, (D.5) uniformly holds with Q̃ = X̃(τ̄) and all rows ofX are `2-bounded with probability 1− 4N exp(−100n)− 8L exp(−c1N0/ρ2). Applying Lemma A.1 on (X̃(τ̄))Lτ̄=1, we conclude with the bound (D.4) on the merged matrixX .
Lemma D.4 (`2-bound on rows). Consider the setup of Theorem D.3. With probability 1− 2N exp(−100(n+ p)), each row ofX has `2-norm at most c √ p+ n for some constant c > 0.
Proof. The tth row of X is equal to xt = [ hTt√ γ+ uTt ] T . Since ‖ht − E[ht]‖ψ2 ≤ O( √ γ+) and ‖ut‖ψ2 ≤ O(1), we have that ‖xt − E[xt]‖ψ2 ≤ O(1). Now, applying Lemma F.3 on all rows {xt}Nt=1, and using a union bound, with probability at least 1− 2N exp(−100(n+ p)), we have that ‖xt − E[xt]‖`2 ≤ c √ n+ p for all t. To conclude, note that ‖E[xt]‖`2 = ‖E[ht]‖`2/ √ γ+ ≤ θ ≤ 3 √ n via Assumption 1.

heading: E PROOFS OF MAIN RESULTS

heading: E.1 PROOF OF LEMMA 3.2
text: The statement follows from upper bound Lemma B.3 and lower bound Lemma B.5.

heading: E.2 PROOF OF THEOREM 4.2
text: Proof. To prove this theorem, we combine Theorem D.3 with deterministic SGD convergence result of Theorem 4.1. Applying Theorem D.3, with the desired probability, inequality (D.4) holds and for all i, input data satisfies the bound ‖xi‖`2 ≤ √ (n+ p)/(2c0) for a sufficiently small constant c0 > 0. As the next step, we will argue that these two events imply the convergence of SGD.
Let θ(i), c(i) ∈ Rn+p denote the ith rows of Θ,C respectively. Observe that the square-loss is separable along the rows of C via ‖Θ−C‖2F = ∑n i=1 ‖θ (i) − c(i)‖2`2 . Hence, SGD updates each row c (i) via its own state equation yt,i = φ( 〈 c(i),xt 〉 ),
where yt,i is the ith entry of yt. Consequently, we can establish the convergence result for an individual row of C. Convergence of all individual rows will imply the convergence of the overall matrix Θτ to the ground
truth C. Pick a row index i (1 ≤ i ≤ n), set c = c(i) and denote ith row of Θτ by θτ . Also denote the label corresponding to ith row by yt = yt,i. With this notation, SGD over (2.3) runs SGD over the ith row with equations yt = φ(〈c,xt〉) and with loss functions
L(θ) = N−1 N∑ t=1 Lt(θ), Lt(θ) = 1 2 (yt − φ(〈θ,xt〉))2.
Substituting our high-probability bounds on xt (e.g. (D.4)) into Theorem 4.1, we can set B = (n+ p)/(2c0), γ+ = (θ + √ 2)2, and γ− = ρ−1/2. Consequently, using the learning rate η = c0 β 2ρ−1
(θ+ √ 2)2(n+p) , for all τ ≥ 0,
the τ th SGD iteration θτ obeys
E[‖θτ − c‖2`2 ] ≤ ‖θ0 − c‖ 2 `2(1− c0
β4ρ−2
2(θ + √ 2)2(n+ p) )τ , (E.1)
where the expectation is over the random selection of SGD updates. This establishes the convergence for a particular row of C. Summing up these inequalities (E.1) over all rows θ(1)τ , . . . ,θ (n) τ (which converge to c(1), . . . , c(n) respectively) yields the targeted bound (4.4).

heading: E.3 PROOFS OF MAIN RESULTS ON STABLE SYSTEMS

heading: E.3.1 PROOF OF THEOREM 3.3
text: Proof. Applying Lemmas B.3 and 3.2, independent of L, Assumption 1 holds with parameters
γ+ = B 2 ∞ , γ− = β 2smin(B) 2 , θ = √ 6n− √ 2 ≥ √ n.
This yields (θ + √ 2)2 = 6n. Hence, we can apply Theorem 4.2 with the learning rate η = c0 β 2
6ρn(n+p) where
ρ = B2∞ β2smin(B)2 = γ+ γ− , (E.2)
and convergence rate 1− β 2η 2ρ . To conclude with the stated result, we use the change of variable c0/6→ c0.

heading: E.3.2 PROOF OF THEOREM 3.4
text: Proof. The proof is similar to that of Theorem 3.3. Applying Lemmas B.3, B.4, and 3.2, independent of L, Assumption 1 holds with parameters
γ+ = B 2 ∞ , γ− = smin(B) 2 , θ = 0.
Hence, we again apply Theorem 4.2 with the learning rate η = c0 β 2
2ρ(n+p) where ρ is given by (E.2). Use the
change of variable c0/2→ c0 to conclude with the stated result.

heading: E.4 LEARNING UNSTABLE SYSTEMS
text: In a similar fashion to Section 4, we provide a more general result on unstable systems that makes a parametric assumption on the statistical properties of the state vector. Assumption 2 (Well-behaved state vector – single timestamp). Given timestamp T0 > 0, there exists positive scalars γ+, γ−, θ and an absolute constant C > 0 such that θ ≤ 3 √ n and the following holds
γ+In Σ[hT0 ] γ−In , ‖hT0 − E[hT0 ]‖ψ2 ≤ C √ γ+ and ‖E[ht]‖`2 ≤ θ √ γ+. (E.3)
The next theorem provides the parametrized result on unstable systems based on this assumption.
Theorem E.1 (Unstable system - general). Suppose we are given N independent trajectories (h(i)t ,u (i) t )t≥0 for 1 ≤ i ≤ N . Sample each trajectory at time T0 to obtain N samples (yi,hi,ui)Ni=1 where ith sample is
(yi,hi,ui) = (h (i) T0+1 ,h (i) T0 ,u (i) T0 ).
Let C, c0 > 0 be absolute constants. Suppose Assumption 1 holds with T0 and sample size satisfies N ≥ Cρ2(n+ p) where ρ = γ+/γ−. Assume φ is β-increasing, zero initial state conditions, and ut
i.i.d.∼ N (0, Ip). Set scaling to be µ = 1/ √ γ+ and learning rate to be η = c0 β 2
ρ(θ+ √ 2)2(n+p) . Starting from Θ0, we run SGD
over the equations described in (2.2) and (2.3). With probability 1−2N exp(−100(n+p))−4 exp(−O( N ρ2 )), all iterates satisfy
E[‖Θi −C‖2F ] ≤ (1− c0 β4
2ρ2(θ + √ 2)2(n+ p) )τ‖Θ0 −C‖2F ,
where the expectation is over the randomness of the SGD updates.
Proof. Set xi = [γ −1/2 + h T i u T i ] T and X = [x1 . . . xN ]T . Since X has i.i.d. rows, we can apply Theorem F.1 and Lemma F.3 to find with the desired probability that
• Rows of xi satisfy ‖xi − E[xi]‖ψ2 ≤ O(1) and E[‖xi‖`2 ] ≤ 3 √ n, hence all rows of X obeys ‖xi‖`2 ≤ √ (n+ p)/(2c0),
• X satisfies
(θ + √ 2)2 X TX
N ρ−1/2.
To proceed, using γ− = ρ−1/2, B = (n+ p)/(2c0), and γ+ = (θ + √
2)2, we apply Theorem 4.1 on the loss function (2.3); which yields the desired result.

heading: E.5 PROOF OF THEOREM 5.1
text: Proof. The proof is a corollary of Theorem E.1. We need to substitute the proper values in Assumption 2. Applying Lemma B.3, we can substitute γ+ = B2T0 and θ = √ 6n− √ 2 ≥ √ n. Next, we need to find a lower bound. Applying Lemma 3.2 for n > 1 and Lemma B.6 for n = 1, we can substitute γ− = γ+/ρ with the ρ definition of (5.2). With these, the result follows as an immediate corollary of Theorem E.1.

heading: F SUPPLEMENTARY STATISTICAL RESULTS
text: The following theorem bounds the empirical covariance of matrices with independent subgaussian rows. Given a random vector x, define the de-biasing operation as zm(x) = x− E[x].
Theorem F.1. Let A ∈ Rn×d be a matrix with independent subgaussian rows {ai}ni=1 satisfying ‖zm(ai)‖ψ2 ≤ O(K) and Σ[ai] K2Id for some K > 0 and ‖E[ai]‖`2 ≤ θ. Suppose Σ[ai] λId. Suppose n ≥ O(K4d/λ2). Then, each of the following happens with probability at least 1− 2 exp(−cK−4λ2n),
• θ + √
3/2K ≥ 1√ n ‖A‖.
• Suppose all rows ofA have equal expectations. Then 1√ n smin(A) ≥
√ 2λ/3.
Proof. Let E = E[A], Ā = A− E[A], āi = zm(ai). We will decomposeA = Ā+E hence we will first focus on bounding the upper and lower singular values of Ā by studying the random processes Xv = ‖Āv‖2`2 and Yv = Xv − E[Xv] over the unit sphere Sd−1. First, we provide a deviation bound for the quantity supv∈Sd−1 |Yv|. To achieve this, we will utilize Talagrand’s mixed tail bound and show that increments of Yv are subexpoential. Pick two unit vectors v,u ∈ Rd. Write x = u+ v,y = u− v. We have that
Xu−Xv = ‖Āu‖2`2−‖Āv‖ 2 `2 = ‖Ā(x+ y)/2‖ 2 `2−‖Ā(x− y)/2‖ 2 `2 = x T ĀT Āy = n∑ i=1 (āTi x)(ā T i y).
Letting x̂ = x/‖x‖`2 , ŷ = y/‖y‖`2 , observe that, multiplication of subgaussians xT āi,yT āi obey
‖(xT āi)(yT āi)‖ψ1 ≤ O(‖x‖`2‖y‖`2K 2) ≤ O(K2‖y‖`2).
Centering this subexponential variable around zero introduces a factor of 2 when bounding subexponential norm and yields ‖(xT āi)(yT āi)− E[(xT āi)(yT āi)]‖ψ1 ≤ O(K2‖y‖`2). Now, using the fact that Yu − Yv is sum of n independent zero-mean subexponential random variables, we have the tail bound
P(n−1|Yu − Yv| ≥ t) ≤ 2 exp(−c′nmin{ t2
K4‖y‖2`2 ,
t
K2‖y‖`2 }).
Applying Talagrand’s chaining bound for mixed tail processes with distance metrics ρ2 = K2‖·‖`2√
n , ρ1 =
K2‖·‖`2 n
, (Theorem 3.5 of Dirksen (2013) or Theorem 2.2.23 of Talagrand (2014)) and using the fact that for unit sphere Sd−1, Talagrand’s γ functionals (see Talagrand (2014)) obey γ1(Sd−1), γ22(Sd−1) ≤ O(d),
n−1 sup v∈Sd−1
|Yv| ≤ cK2( √ d/n+ d/n+ t/ √ n), (F.1)
with probability 1 − 2 exp(−min{t2, √ nt}). Since n ≥ Cλ−2K4d for sufficiently large C > 0, picking t = 1 16c K−2λ √ n, with probability 1− 2 exp(−O(K−4λ2n)), we ensure that right hand side of (F.1) is less than λ/8. This leads to the following inequalities
1 n ‖ĀT Ā− E[ĀT Ā]‖ ≤ λ 8 =⇒ 9K
2
8 Id
1 n ĀT Ā 7λ 8 Id. (F.2)
=⇒ 9 8 K ≥ 1√ n ‖Ā‖ ≥ smin(Ā) ≥
√ 7
8 λ.
Upper bound on spectral norm: For spectral norm ofA, we use the triangle inequality
1√ n ‖A‖ ≤ 1√ n (‖E‖+ ‖Ā‖) ≤ max 1≤i≤n ‖E[ai]‖`2 + 9K/8 ≤ θ +
√ 3/2K.
Lower bound on minimum singular value: This part assumes that all row expectations are same. Denote the size n all ones vector by 1n and define the process Zv = 1√n1 T n Āv. Observe that ĀT1n = ∑n i=1 āi ∈ R d is a vector satisfying ‖ĀT1n/ √ n‖ψ2 ≤ O(K). Hence, again using n ≥ CK4λ−2d for sufficiently large C > 0, applying Lemma F.3 with m = c0K−4λ2n > d by picking a sufficiently small constant c0 > 1/C, with probability at least 1− 2 exp(−100c0K−4λ2n)
1√ n
sup ‖v‖`2=1
|Zv| = 1
n ‖ĀT1n‖`2 ≤
1
12 KK−2λ ≤
√ λ 12 .
Let P = In − 1n1n1 T n be the projection onto the orthogonal complement of the all ones vector. Note that PEv = 0 as the rows of E are equal. With this observation, with desired probability, for any unit length v,
‖Av‖`2 ≥ ‖PAv‖`2 = ‖PĀv‖`2 ≥ ‖Āv‖`2 − |Zv| (F.3)
≥ smin(Ā)− sup v∈Sd−1
|Zv| ≥ ( √ 7/8− 1/12) √ λn, (F.4)
which implies smin(A)/ √ n ≥ √ 2λ/3.
The corollary below is obtained by slightly modifying the proof above by using 1 n ‖ĀT Ā− E[ĀT Ā]‖ ≤ K
2
8 in line (F.2) and only focusing on the spectral norm bound.
Corollary F.2. Let A ∈ Rn×d be a matrix with independent {ai}ni=1 subgaussian rows satisfying ‖zm(ai)‖ψ2 ≤ O(K) and Σ[ai] K2Id for some K > 0 and ‖E[ai]‖`2 ≤ θ. Suppose Σ[ai] λId. Suppose n ≥ O(K2d). Then, with probability at least 1− 4 exp(−cK−2n),
θ + √
3/2K ≥ 1√ n ‖A‖.
The following lemma is fairly standard and is proved for the sake of completeness.
Lemma F.3 (Subgaussian vector length). Let a ∈ Rn be a zero-mean subgaussian vector with ‖a‖ψ2 ≤ L. Then, for any m ≥ n, there exists C > 0 such that
P(‖a‖`2 ≤ CL √ m) ≥ 1− 2 exp(−100m).
Proof. We can pick a 1/2 cover C of the unit `2-sphere with size log |C| ≤ 2n. For any v ∈ C, subgaussianity implies, P(|vTa| ≥ t) ≤ 2 exp(− ct 2
2L2 ). Setting t = CL
√ m for sufficiently large constant C > 0, and union
bounding over all v ∈ C, we find
P( ⋂ v∈C ‖v‖`2 ≤ CL √ m) ≥ 1− 2 exp(2n− cC 2L2m 2L2 ) ≤ 1− 2 exp(−100m).
To conclude, let v(a) ∈ C be a’s neighbor satisfying ‖v − a‖a‖`2 ‖`2 ≤ 1/2. Hence, we have
‖a‖`2 ≤ ‖(a− v(a)) Ta‖`2 + ‖v Ta‖`2 ≤ ‖a‖`2/2 + CL √ m =⇒ ‖a‖`2 ≤ 2CL √ m.
To conclude, use the change of variable C → C/2.

Review summaries:
The paper studies discrete time dynamical systems with a non-linear state equation.  They assume the non-linear function is assumed to be \beta-increasing like leaky ReLU. Under this setting, the authors prove that for the given state equation for stable systems with random gaussian input at each time step, running SGD on a fixed length trajectory gives logarithmic convergence.

This work considers the problem of learning a non-linear dynamical system in which the output equals the state.  Under several assumptions (input is Gaussian, non-linear activation is strictly increasing, stable system) it is shown that SGD converges linearly to the ground truth system with near-optimal sample complexity. The proof idea is to reduce this problem to the problem of learning a single non-linear neuron in the case that the covariance matrix of the data is well-conditioned. The main challenge is to show the covariance is well-conditioned under the reduction. In a nutshell, this is done by splitting the trajectory to sub-trajectories with independent states and using results from random matrix theory on matrices with independent rows.