Title: STOCHASTIC GRADIENT DESCENT LEARNS STATE EQUATIONS WITH NONLINEAR ACTIVATIONS

Abstract: We study discrete time dynamical systems governed by the state equation h t+1 = φ(Ah t + Bu t ). Here A, B are weight matrices, φ is an activation function, and u t 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 {u t , h t } N 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.
(2.1)
In this work, our attention is restricted to the state equation (1.1); which corresponds to setting y t = h t+1 in the output equation. To analyze (1.1) in a non-asymptotic data-dependent setup, we assume a finite input/state trajectory of {u t , h t } N t=0 generated by some ground truth weight matrices

1. INTRODUCTION:

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 matrix A, and input weight matrix B as follows

h t+1 = φ(Ah t + Bu t ), (1.1)

Here h t is the state vector and u t 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 (u t , h t ) N t=0 generated from the state equation (1.1). We consider a least-squares regression obtained from N equations; with inputs (u t , h t ) N t=1 and outputs (h t+1 ) N t=1 . For a class of activation functions including leaky ReLU and for stable systems 1 , 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 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.

1.1. RELATED WORK:

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. 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.

2. PROBLEM SETUP:

We first introduce the notation. • returns the spectral norm of a matrix and s min (•) 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. I n is the identity matrix of size n × n. Normal distribution with mean µ and covariance Σ is denoted by N (µ, Σ). Throughout, c, C, c 0 , c 1 , . . . denote positive absolute constants.

Setup: We consider the dynamical system parametrized by an activation function φ(•) and weight matrices A ∈ R n×n , B ∈ R n×p as described in (1.1). Here, h t is the n dimensional state-vector and u t 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 h t is hidden and we only get to interact with h t 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 Algorithm 1 Learning state equations with nonlinear activations 1: Inputs: (y t , h t , u t ) N t=1 sampled from a trajectory. Scaling µ, learning rate η. Initialization A 0 , B 0 . 2: Outputs: Estimates Â, B of the weight matrices A, B.

3: x t ← [µh T t u T t ] T for 1 ≤ t ≤ N . 4: Θ 0 ← [µ -1 A 0 B 0 ] 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 [ Â B] ← Θ END µI n 0 0 I p .

(A, B). Our goal is learning the unknown weights A and B 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 construct N triples of the form {u t , h t , h t+1 } N t=1 . We formulate a regression problem by defining the output vector y t , input vector x t , and the target parameter C as follows

y t = h t+1 , x t = µh t u t ∈ R n+p , C = [µ -1 A B] ∈ R n×(n+p) . (2.2)

With this reparameterization, we find the input/output identity y t = φ(Cx t ). We will consider the least-squares regression given by

L(Θ) = 1 N N t=1 L t (Θ) where L t (Θ) = 1 2 y t -φ(Θx t ) 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 Ĉ = Θ END . Estimates of A and B are decoded from the left and right submatrices of Ĉ respectively.

3. MAIN RESULTS:



3.1. PRELIMINARIES:

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 Σ[h t ] 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 h 0 = 0 and u t i.i.d.

∼ N (0, I p ). Define the upper bound term B t as

B t = B 1 -A 2t 1 -A 2 . (3.2)

• Suppose φ is 1-Lipschitz and φ(0) = 0. Then, for all t ≥ 0,

Σ[h t ] B 2 t I n . • Suppose φ is a β-increasing function and p ≥ n. Then, Σ[h t ] β 2 s min (B) 2 I n .

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 h t will exponentially converge to 0 i.e. the system forgets the past states quickly. This is also the reason (B t ) 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 ∞ βs min (B) 2 = B s min (B) 2 1 β 2 (1 -A 2 ) . (3.3)

Observe that, the condition number of B appears inside the ρ term.

3.2. LEARNING FROM SINGLE TRAJECTORY:

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 of A. Theorem 3.3 (Main result). Let {u t , h t+1 } N t=1 be a finite trajectory generated from the state equation (1.1). Suppose A < 1, φ is β-increasing, h 0 = 0, p ≥ n, and u t i.i.d.

∼ N (0, I p ). Let ρ be same as (3.3) and c, C, c 0 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 η = c 0 β 2 ρn(n+p) , and consider the loss function (2.3). With probability 1 -4N exp(-100n) -8L exp(-O( N Lρ 2 )), starting from an initial point Θ 0 , for all τ ≥ 0, the SGD iterations described in Algorithm 1 satisfies

E[ Θ τ -C 2 F ] ≤ (1 -c 0 β 4 2ρ 2 n(n + p) ) τ Θ 0 -C 2 F . (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 h t+1 ). We collect N state equations to obtain a system of N n equations. On the other hand, the total number of unknown parameters in A and B are n(n + p). This implies Theorem 3.3 is applicable as soon as the problem is mildly overdetermined i.e. N n 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 η = c 0 β 2 ρ(n+p) , and consider the loss function (2.3). With probability 1 -4N exp(-100n) -8L exp(-O( N Lρ 2 )), starting from an initial point Θ 0 , for all τ ≥ 0, the SGD iterations described in Algorithm 1 satisfies

β 4 2ρ 2 (n + p) ) τ Θ 0 -C 2 F , (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.

4. MAIN IDEAS AND PROOF STRATEGY:

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 h t 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 u t provides sufficient excitation for the output state h t+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 (h iT , u iT ) 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 y lin t = Cx t and define the linear loss to be

L lin (Θ) = 1 2N N i=1 y lin t -Θx t 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 {x i , y i } N i=1 is given; where output y i is related to input x i via y i = φ( x i , θ ) for some θ ∈ R n . Suppose β > 0 and φ is a β-increasing. Let γ + ≥ γ -> 0 be scalars. Assume that input samples satisfy the bounds

γ + I n 1 N N i=1 x i x T i γ -I n , x i 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 = θ τ -η(φ(x T rτ θ τ ) -y rτ )φ (x T rτ θ τ )x rτ , (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. 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: Σ[h L-1 ] γ -I n ,

• Upper bound: for all t, the state vector satisfies

Σ[h t ] γ + I n , h t -E[h t ] ψ2 ≤ C √ γ + and E[h t ] 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 {u t , h t+1 } N t=1 be a length N trajectory of the state equation (1.1). Suppose A < 1, φ is β-increasing, h 0 = 0, and u t i.i.d.

∼ N (0, I p ). Given scalars γ + ≥ γ -> 0, set the condition number as ρ = γ + /γ -. For absolute constants C, c, c 0 > 0, choose trajectory length N to satisfy

N ≥ CLρ 2 (n + p) where L = 1 - log (cnρ) log A .

Suppose Assumption 1 holds with L, γ + , γ -, θ. Pick scaling to be µ = 1/ √ γ + and learning rate to be η = c 0

β 2 ρ(θ+ √ 2) 2 (n+p) . With probability 1 -4N exp(-100n) -8L exp(-O( N Lρ 2 ))

, starting from Θ 0 , for all τ ≥ 0, the SGD iterations on loss (2.3) as described in Algorithm 1 satisfies

E[ Θ τ -C 2 F ] ≤ (1 -c 0 β 4 2ρ 2 (θ + √ 2) 2 (n + p) ) τ Θ 0 -C 2 F , (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

γ + = B 2 ∞ , γ -= β 2 s min (B) 2 , and θ = √ n. • Theorem 3.4 follows by setting γ + = B 2 ∞ , γ -= β 2 s min (B) 2

, and θ = 0.

5. LEARNING UNSTABLE SYSTEMS:

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 T 0 ≥ 1. Denoting the ith trajectory by the triple (h

EQUATION

t ) t≥0 , we collect a single sample from each trajectory at time T 0 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,

(y i , h i , u i ) = (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 weights A and B. The crucial difference compared to Section 3 is that, the samples (y i , h i , u i ) N i=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|) 2T 0 if n = 1 where ρ = B 2 T0 β 2 s min (B) 2 .

(5.2)

The ρ term is similar to the earlier definition (3.3); however it involves B T0 rather than B ∞ . This modification is indeed necessary since B ∞ = ∞ when A > 1. On the other hand, note that, B2 T0 grows proportional to A 2T0 ; which results in exponentially bad condition number in T 0 . Our ρ definition remedies this issue for single-output systems; where n = 1 and A is a scalar. In particular, when β = 1 (e.g. φ is linear) ρ becomes equal to the correct value 1 2 . 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 T 0 to obtain N samples (y i , h i , u i ) N i=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 u t i.i.d. ∼ N (0, I p ). Set scaling µ = 1/ B T0 , learning rate η = c 0 β 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 2 F ] ≤ (1 -c 0 β 4 2ρ 2 n(n + p) ) τ Θ 0 -C 2 F ,

where the expectation is over the randomness of the SGD updates.

6. NUMERICAL EXPERIMENTS:

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 matrix A 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.

Under review as a conference paper at ICLR 2019 Algorithm 2 Empirical hyperparameter selection.

1: Inputs: (h t , u t ) N t=1 sampled from a trajectory. 2: Outputs: Scaling µ. Evaluation: We consider two performance measures in the experiments. Let Ĉ be an estimate of the ground truth parameter

C = [µ -1 A B].

The first measure is the normalized error defined as

Ĉ -C 2 F / C 2 F .

The second measure is the normalized loss defined as

N i=1 y t -φ( Ĉx t ) 2 2 N i=1 y t 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 by A 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. As A 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. 

7. CONCLUSIONS:

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).

A . φ (a, b) ≥ β since φ is β-increasing. Define wτ to be the residual wτ = θτθ. Observing φ(x T rτ θτ ) -yr τ = φ (x T rτ θτ , x T rτ θ)x T rτ wτ , the SGD recursion obeys

wτ+1 2 2 = wτ -η(φ(x T rτ θτ ) -yr τ )φ (x T rτ θτ )xr τ 2 2 . = wτ -ηxr τ φ (x T rτ θτ )φ (x T rτ θτ , x T rτ θ)x T rτ wτ 2 2 = (I -ηGr τ )wτ 2 2

where Gr τ = xr τ φ (x T rτ θτ )φ (x T rτ θτ , x T rτ θ)x T rτ . Since φ is 1-Lipschitz and β-increasing, Gr τ is a positivesemidefinite matrix satisfying

xr τ x T rτ Gr τ β 2 xr τ x T rτ , G T rτ 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[G T rτ 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[w T τ G T rτ Gr τ wτ ] ≤ wτ 2 2 (1 -2ηβ 2 γ-+ η 2 Bγ+).

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) ∈ R N i ×q are given for 1 ≤ i ≤ L. Suppose for all 1 ≤ i ≤ L, rows of X (i) has 2 norm at most √ B and each

X (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 X T X 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 that X T X = L i=1 X (i) T X (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 on X T X/N .

B PROPERTIES OF THE NONLINEAR STATE EQUATIONS:

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 ∈ R p is defined as x ψa = sup v∈B p v T x ψa where B p 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} i =τ (i.e. all except uτ ), ht+1 is A t-τ B Lipschitz function of uτ for 0 ≤ τ ≤ t.

Proof. Fixing {ui} i =τ , 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τ . ∼ 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 = [u T 0 . . . u T t-1 ] T ∈ R tp .

• There exists an absolute constant c > 0 such that ht -E[ht] ψ 2 ≤ cBt and Σ[ht] B 2 t 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 qt, we find a scalar L f satisfying,

f (qt) -f ( qt) 2 ≤ L f qt -qt 2 . (B.1)

Define the vectors, {ai} t i=0 , as follows

ai = [ ûT 0 . . . ûT i-1 u T i . . . u T t-1 ] T .

Observing that a0 = qt, at = qt, we write the telescopic sum,

f (qt) -f ( qt) 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, û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 -ûi 2 .

To bound the sum, we apply the Cauchy-Schwarz inequality; which yields

|f (qt) -f ( qt)| ≤ t-1 i=0 A t-1-i B ui -ûi 2 ≤ ( t-1 i=0 A 2(t-1-i) B 2 ) 1/2 ( t-1 i=0 ui -ûi 2 2 ) 1/2 q t -qt 2 ≤ B 2 (1 -A 2t ) 1 -A 2 qt -qt 2 = Bt qt -qt 2 . (B.

2)

The final line achieves the inequality (B.1) with L f = 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 as N (0, Itp). Since ht a Bt Lipschitz function of qt, for any fixed unit length vector v, αv := v T ht = v T 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(- t 2 2B 2 t ).

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] ≤ B 2 t (page 49 of Ledoux ( 2001)), which implies the covariance bound

Σ[ht] B 2 t 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 2i tr(BB T ) ≤ tr(BB T ) 1 -A 2(t+1) 1 -A 2 .

Now using the fact that tr(BB T ) ≤ 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(- t 2 2B 2 t

).

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 ⊂ R n 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 ∼ N (0, Ip). Suppose φ is a β-increasing function for some constant β > 0. For any t ≥ 1, the state vector obeys

P(ht+1 ⊂ S) = P(Aht + But ⊂ S ) = P(A(-ht) + B(-ut) ⊂ S ) (B.

Σ[ht] β 2 smin(BB T )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(Σ)I d ) and g2 ∼ N (0, Σ -smin(Σ)I d ) are independent multivariate Gaussians.

Since But ∼ N (0, BB T ), setting Σ = BB T 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 ≥ 0, via

Σ[ht+1] = Σ[φ(g1 + g2 + Aht)] β 2 sminIn.

The next theorem applies to multiple-input-single-output (MISO) systems where A is a scalar and B 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)] ≥ β 2 var[X]. Recursively, this yields

var[ht+1] = var[φ(Aht + But)] ≥ β 2 var[Aht + But] = β 2 (|A| 2 var[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)] β 2 var[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 

= H(τ) ∈ R N ×n , H = H(τ) ∈ R N ×n , Ũ = Ũ (τ ) ∈ R N ×p , and X = X(τ) ∈ R N ×(n+p) as H = [ h(1) . . . h( N ) ] T , H = [h (1) . . . h ( N ) ] T , Ũ = [u (1) . . . u ( N ) ] T , X = [µ H Ũ ].

Q = [γ -1/2 + H Ũ ] ∈ R N ×(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 perturbed Q matrices given by,

Q = [γ -1/2 + M Ũ ], (D.2) satisfy (Θ + √ 2) 2 QT 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[ HT Ũ ] = 0 due to independence. Also, for i > 1, the truncated state vector h(i) is statistically identical to hL-

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 . . . q N ] T and note that Q = [q1 QT ] T . Applying Theorem F.1 on Q and Corollary F.2 on Q, we find that, with the desired probability,

EQUATION

.

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 = 1 - log (cnρ) log A , N0 = N L ≥ 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 of X has 2 norm at most c0 √ p + n where c0 is an absolute constant.

• X T X obeys the bound

(Θ + √ 2) 2 In+p X T X N ρ -1 In+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(τ) , Ũ (τ ) , X(τ) . Observe that X 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 matrix X.

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(τ) , Ũ (τ ) , set Q = [µ H(τ) Ũ (τ ) ]. 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 matrices M satisfying M -H(τ) ≤ √ γ-N0/10 and Q as in (D.2) obeys

(Θ + √ 2) 2 QT 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 with M = H(τ) choice; which sets Q = X(τ) in (D.5). Applying Lemma C.6, we find that, with probability 1 -2 Nτ 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{200c 2 0 , 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 = 1 - log (cnρ) log A . (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 )-2 Nτ 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 of X 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 matrix X. Proof. The tth row of X is equal to xt = [

h T t √ γ + u T t ] 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} N t=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.

E PROOFS OF MAIN RESULTS:

E.1 PROOF OF LEMMA 3.2

The statement follows from upper bound Lemma B.3 and lower bound Lemma B.5.

E.2 PROOF OF THEOREM 4.2:

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) ∈ R n+p denote the ith rows of Θ, C respectively. Observe that the square-loss is separable along the rows of C via Θ -C 2 F = 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

EQUATION

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) τ , . . . , θ 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.

E.4 LEARNING UNSTABLE SYSTEMS:

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 Σ[hT 0 ] γ-In ,

hT 0 -E[hT 0 ] ψ 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

EQUATION

t ) t≥0 for 1 ≤ i ≤ N . Sample each trajectory at time T0 to obtain N samples (yi, hi, ui) N i=1 where ith sample is (yi, hi, ui) = (h

(i) T 0 +1 , h (i) T 0 , u (i) T 0 ).

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

EQUATION

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 2 F ] ≤ (1 -c0 β 4 2ρ 2 (θ + √ 2) 2 (n + p) ) τ Θ0 -C 2 F ,

where the expectation is over the randomness of the SGD updates.

:

Throughout this work, a system is called stable if the spectral norm of the state matrix A is less than 1.

Clearly, 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.

Note 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 β.

Figures and Tables:

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.

Figure 1: SGD convergence behavior for Leaky ReLUs with varying minimum slope β. Figures a) and b) repeat the same experiments. The difference is the spectral norm of the ground truth state matrix A.

Figure 2: SGD convergence behavior for ReLU with varying spectral norm of the state matrix A. Figures a) and b) repeats the same experiments. The difference is that a) uses N = 500 trajectory length whereas b) uses N = 2500 (i.e. ×5 more data). Shaded regions highlight the one standard deviation around the mean.

DETERMINISTIC CONVERGENCE RESULT FOR SGD Proof of Theorem 4.1. Given two distinct scalars a, b; define φ (a, b) = φ(a)-φ(b) a-b

Lemma B.3 (Upper bound). Consider the state equation governed by equation (1.1). Suppose ut i.i.d.

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.

1≤i≤n var[φ(xi + yi)]In.Fortunately, Lemma B.9 provides the lower bound var[φ(xi + yi)] ≥ β 2 var[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]].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

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

Lemma D.4 ( 2-bound on rows). Consider the setup of Theorem D.3. With probability 1 -2N exp(-100(n + p)), each row of X has 2-norm at most c √ p + n for some constant c > 0.

to c (1) , . . . , c (n) respectively) yields the targeted bound (4.4).E.3 PROOFS OF MAIN RESULTSON STABLE SYSTEMS E.3.1 PROOF OF THEOREM 3.3 Proof. Applying Lemmas B.3 and 3.2, independent of L, Assumption 1 holds with parameters γ+ = B 2 ∞ , γ-= β 2 smin(B) 2 , θ = = 6n. Hence, we can apply Theorem 4.2 with the learning rate η = c0 β 2 6ρn(n+p) 1 -β 2 η 2ρ . To conclude with the stated result, we use the change of variable c0/6 → c0. E.3.2 PROOF OF THEOREM 3.4 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.

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 constraintsCai et al. (2010);Beck & Teboulle (2009);

Table content: <table/>

Peer review (iclr2019):
 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. This work tackles a very challenging problem and the results are interesting. The guarantees are strong – linear convergence to the ground truth parameters and near-optimal sample size. Given that not much is known on deep non-linear networks, I think that the result is significant. The main weakness of the paper is the assumption that the state equals the output. Another minor weakness is the clarity and presentation of results: 1. The proof outline of the main result is hard to follow. There is no proof outline of Theorem 4.2 in the main text. The proof is highly technical and there are many technical ideas that were moved to the appendix. For instance, the proofs in sections C and D are not mentioned in the main text. I suggest to write a summary of the steps required to prove the main result and how all of the technical ideas are combined together. 2. There is no reference and comparison to the paper of Mei et al. [1] that study single neuron models. 3. It is claimed that by increasing beta the convergence is faster. However, I am not sure why this is meaningful. By changing beta the ground truth changes as well. For beta = 0 the ground truth dynamical system is linear and for beta = 1 the ground truth is a non-linear dynamical system with ReLU. Since a ReLU network is more expressive, generally in the case of beta = 1 the ground truth is more difficult to learn than beta = 0. Therefore, we should expect convergence to be slower than beta=0 or not occur at all. Am I missing something? 4. The Gaussian assumption is not stated clearly in the text. It can be deduced only from the statements of the theorems and the conclusion section. 5. In Theorem F.1, it is claimed that all rows of E are equal. However, in the statement of the theorem it is not mentioned that the rows of A are identically distributed. Should this assumption be included in the statement? [1] Mei, Song, Yu Bai, and Andrea Montanari. "The landscape of empirical risk for non-convex losses." arXiv preprint arXiv:1607.06534 (2016).‏ -----------Revision------------------------ I am not changing the score. I disagree with AnonReviewer2 regarding the significance of the results. The assumption that the states are observed is indeed a weakness of the paper. However, understanding non-linear dynamical systems is extremely challenging and this paper provides strong convergence guarantees. Furthermore, there are several insights in the analysis that may be useful in future work.