Near-optimal Offline and Streaming Algorithms for Learning Non-Linear Dynamical Systems
Year: 2021

Abstract:
We consider the setting of vector valued non-linear dynamical systems Xt+1 = φ(AXt) + ηt, where ηt is unbiased noise and φ : R→ R is a known link function that satisfies certain expansivity property. The goal is to learn A∗ from a single trajectory X1, · · · , XT of dependent or correlated samples. While the problem is well-studied in the linear case, where φ is identity, with optimal error rates even for non-mixing systems, existing results in the non-linear case hold only for mixing systems. In this work, we improve existing results for learning nonlinear systems in a number of ways: a) we provide the first offline algorithm that can learn non-linear dynamical systems without the mixing assumption, b) we significantly improve upon the sample complexity of existing results for mixing systems, c) in the much harder one-pass, streaming setting we study a SGD with Reverse Experience Replay (SGD− RER) method, and demonstrate that for mixing systems, it achieves the same sample complexity as our offline algorithm, d) we justify the expansivity assumption by showing that for the popular ReLU link function — a non-expansive but easy to learn link function with i.i.d. samples — any method would require exponentially many samples (with respect to dimension of Xt) from the dynamical system. We validate our results via simulations and demonstrate that a naive application of SGD can be highly sub-optimal. Indeed, our work demonstrates that for correlated data, specialized methods designed for the dependency structure in data can significantly outperform standard SGD based methods.

heading: 1 Introduction
text: Non-linear dynamical systems (NLDS) are commonly used to model the data in a variety of domains like control theory, time-series analysis, and reinforement learning (RL) [1–4]. Standard NLDS models the data points (X0, X1, . . . , XT ) as:
Xt+1 = φ(A ∗Xt) + ηt, (1)
where Xt ∈ Rd are the states, ηt ∈ Rd are i.i.d. noise vectors, A∗ ∈ Rd×d and φ : R → R is an increasing function called the ‘link function’. Here, φ is supposed to act component wise over Rd. System identification problem is a foundational problem for NLDS, i.e., given (X0, X1, . . . , XT ) generated from (1), the goal is to estimate A∗ accurately from a single trajectory (X0, X1, . . . , XT ). 35th Conference on Neural Information Processing Systems (NeurIPS 2021).
The system identification problem is heavily studied in control theory [5–8] as well as time-series analysis [9]. For instance, the non-linear dynamical system considered here has an application in modeling non-linear distortions in power amplifiers [10]. The problem is challenging as data points X0, X1, . . . , XT are not i.i.d. as usually encountered in machine learning, but form a Markov process. If the mixing time τmix of the process is finite (τmix <∞), then we can make the data approximately i.i.d. by considering only the points separated by Õ(τmix) time. While this allows using standard techniques for i.i.d. data, it reduces the effective number of samples to O( Tτmix ), which typically gives an error of the order O( τmixT ). In fact, even the state-of-the-art results have error bounds which are sub-optimal by a factor of τmix.
Interestingly, for the special case of linear systems, i.e., when φ(x) = x, the results are significantly stronger. For example, [11, 12] showed that the matrix A∗ can be estimated with an error O(1/T ) even when the mixing time τmix > T . But these results rely on the fact that for linear systems, the estimation problem reduces to an ordinary least squares (OLS) problem for which a closed form expression is available and can be analyzed effectively.
On the other hand, NLDS do not admit such closed form expressions. In fact the existing techniques mostly rely on mixing time arguments to induce i.i.d. like behavior in a subset of the points which leads to sub-optimal rates by τmix factor. Similarly, a direct application of uniform convergence results [13] to show that the minimizer of the empirical risk is close to the population minimizer still gives sub-optimal rates as off-the-shelf concentration inequalities (cf. [14]) incur an additional factor of mixing time. Finally, existing results are mostly focused on offline setting, and do not apply to the case where the data points are streaming which is critical in several practical problems like reinforcement learning (RL) and control theory.
In this work, we provide algorithms and their corresponding error rates for the NLDS system identification problem in both offline and online setting, assuming the link function to be expansive (Assumption 1). The main highlight of our results is that the error rates are independent of the mixing time τmix, which to the best of our knowledge is first such result for any non-linear system identification in any setting. In fact, for offline setting, our analysis holds even for systems which do not mix within time T and even for marginally stable systems which do not mix at all. Furthermore, we analyze SGD-Reverse Experience Replay (SGD-RER) method, we provide the first streaming method for NLDS identification with error rate that is independent of τmix (in the leading order term) while still ensuring small space and time complexity. This algorithm was first discovered in the experimental RL setting in [15] based on Hippocampal reverse replay observed in biological networks [16–18]. It was introduced independently in [19] for the case of linear systems and efficiently unravels the complex dependency structure present in the problem. Finally, through a lower bound for ReLU— a non-expansive function—we provide strong justification for why expansivity might be necessary for a non-trivial result.
Instead of mixing time arguments, our proofs for learning NLDS without mixing use a natural exponential martingale of the kind considered in the analysis of self normalized process ([20, 21]). For streaming setting, while we do use mixing time arguments (proof of Theorems 2 and 3), we combine them with a delicate stability analysis of the specific algorithm and the machinery developed in [19] to obtain strong error bounds. See Section 6 for a description of these techniques.
Our Contributions. Key contributions of the paper are summarized below:
1. Assuming expansive and monotonic link function φ and sub-Gaussian noise, we show that the offline Quasi Newton Method (Algorithm 1) estimates the parameter A∗ with near optimal errors of the order O(1/T ), even when the dynamics does not mix within time T .
2. Assuming mixing NLDS, finite fourth moment on the noise, and expansive monotonic link function, we show that offline Quasi Newton Method again estimates the parameter A∗ with near-optimal error of O(1/T ), independent of mixing time τmix.
3. We give a one-pass, streaming algorithm inspired by SGD− RER method by [19], and show that it achieves near-optimal error rates under the assumption of sub-Gaussian noise, NLDS stability (see section 2.1 for the definition), uniform expansivity and second differentiability of the link function.
4. We then show that learning with ReLU link function, which is non-expansive but is known to be easy to learn with if data points are all i.i.d. [22], requires exponential (in d) many samples.
We believe that the techniques developed in this work can be extended to provide efficient algorithms for learning with dependent data in more general settings.
Related Works. NLDS has been studied in a variety of domains like time-series and recurrent neural networks (RNN). [9] studies specific NLDS models from a time series perspective and establishes non-asymptotic convergence bounds for natural estimators; their error rates suffer from mixing time factor τmix. [23] considers asymptotic learning of NLDS via neural networks trained using SGD, whereas [24] shows that overparametrized LSTMs trained with SGD learn to memorize the given data. [25–27] consider learning dynamical systems of the form ht+1 = φ(A∗ht +B∗ut) for states ht and inputs ut; this setting is different from standard NLDS model we study. [28] considers the non linear dynamical systems of the form xt+1 = Aφ(xt, ut) + ηt which φ is a known non-linearity and matrix A is to be estimated. [29, 30] consider essentially linear dynamics but allow for certain non-linearities that can be modeled as process noise. All these again differ from the model we consider.
Standard NLDS identification (1) has received a lot of attention recently, with results by [31, 32] being the most relevant. [31] uses uniform convergence results via. mixing time arguments to obtain parameter estimation error for offline SGD. [32] obtains similar bounds via. the analysis of the GLMtron algorithm [33]. However, both these works suffer from sub-optimal dependence on the mixing time. We refer to Table 1 for a comparison of the results. [34], which appeared after the initial manuscript of this work, considers the question from a perspective of time series forecasting. This work considers sparsity in A∗ and an unknown link function which is estimated with isotonic regression. Their recovery guarantees eschew the mixing time dependence. However, the setting, assumptions and the error rates are incomparable to our setting.
[32] also obtains within sample prediction error in the case when φ is not uniformly expansive along with parameter recovery bounds when φ is the ReLU function and the driving noise is Gaussian. However, the parameter estimation bounds for ReLU suffer from an exponential dependence on the dimension d and mixing time τmix. In Theorem 4 we establish that indeed we cannot improve the exponential dependence in the dimension d for the case of parameter estimation. We note that the exponential dependence arises due to the dynamics present in the system since ReLU regression with isotropic i.i.d. data in well specified case has only a polynomial dependence in d [22].
Linear system identification (LSI) literature has been well studied with strong minimax optimal bounds [11, 35, 36]. These results primarily consider the (convex) empirical square loss which has a closed form solution. However, the square loss in the non-linear case is non-convex. Under the assumption that the link function is increasing, we consider a convex proxy loss which is widely used in generalized linear regression literature [22, 33, 37]. Similarly, GLMtron algorithm for learning NLDS in[32] (see Equation (3)) also considers a similar proxy loss. In [38], the authors consider a family of GLMtron-like algorithms call Reflectron under the i.i.d. data setting. But they compare the performance of these algorithms experimentally on an NLDS similar to one considered in this work under low rank assumption on the system matrix.
Finally, streaming setting for LSI has been recently studied in different model settings [19, 39]. These methods observe that by exploiting techniques like experience replay ([40]) along with squared loss error, one can obtain strong error rates. SGD− RER method studied in this work is inspired by a similar method by [19] which was primarily studied for the linear case.

heading: 2 Problem Statement
text: Let φ : R → R be an increasing, 1-Lipschitz function such that φ(0) = 0. Suppose X0 ∈ Rd is a random variable and A∗ ∈ Rd×d. We consider the following non-linear dynamical system (NLDS):
Xt+1 = φ(A ∗Xt) + ηt, (2)
where the noise sequence η0, . . . , ηT is i.i.d random vectors independent of X0. The noise ηt is such that Eηt = 0, Eηtη>t = σ2I for some σ > 0. We will also assume that M4 := E‖ηt‖4 < ∞. Let µ be the law of noise η. We denote the model above as NLDS(A∗, µ, φ). Whenever a stationary distribution exists for the process, we will denote it by π(A∗, µ, φ) or just π when the process is clear from context. We will call the trajectory X0, X1, . . . , XT ‘stationary’ if X0 is distributed according to the measure π(A∗, µ, φ). Unless specified otherwise, we take X0 = 0 almost surely.
t=0 ‖φ(AXt) −Xt+1‖2. However,
when the link function φ is not linear, then this would be non-convex and hard to optimize. Instead, we use a convex proxy loss given by:
Lprox(A;X) = 1
T T−1∑ t=0 d∑ i=1 φ̄(〈ai, Xt〉)− 〈ei, Xt+1〉〈ai, Xt〉 , (3)
where φ̄ is the indefinite integral of the link function φ and ai is the i-th row of A. Note that the gradient of Lprox(A;X) with respect to A is given by:
∇Lprox(A;X) = 1
T T−1∑ t=0 (φ(AXt)−Xt+1)X>t . (4)
When the model is clear from context and the stationary distribution exists, we will denote the second moment matrix under the stationary distribution by G := E[XtX>t ]. Note that G E[ηtη>t ] = σ2I . Also, the empirical second moment matrix is denoted by Ĝ := 1T ∑T−1 t=0 XtX > t .

heading: 2.1 Assumptions
text: We now state the assumptions below and use only a subset of the assumptions for each result.
Assumption 1 (Lipschitzness and Uniform Expansivity). φ is 1-Lipschitz and |φ(x) − φ(y)| ≥ ζ|x− y|, for some ζ > 0.
Note that when φ is only weakly differentiable but satisfy Assumption 1, with a slight abuse of notation, we will write down φ(x)− φ(y) = φ′(β)(x− y) for some φ′(β) ∈ [ζ, 1]. Assumption 2 (Bounded 2nd Derivative). φ is twice continuously differentiable and |φ′′| is bounded. Assumption 3 (Noise Sub-Gaussianity). For any unit norm vector x ∈ Rd, we have 〈ηt, x〉 to be sub-Gaussian with variance proxy Cησ2.
Next, we extend the definition of exponential stability in [31] to ‘exponential regularity’ to allow unstable systems.
Assumption 4 (Exponential Regularity). Let XT = hT−1(X0, η0, . . . , ηT ) be the function representation of XT . We say that NLDS(A∗, µ, φ) is (Cρ, ρ) exponentially regular if for any choice of T ∈ N and X0, X ′0, η0, . . . , ηT ∈ Rd:
‖hT (X0, η0, . . . , ηT )− hT (X ′0, η0, . . . , ηT )‖2 ≤ CρρT−1‖X0 −X ′0‖2 .
When ρ < 1, we will call the system stable. When ρ = 1 we will call it ‘possibly marginally stable’ and when ρ > 1, we will call it ‘possibly unstable’.
Algorithm 1: Quasi Newton Method Input :Offline data {X0, . . . , XT }, horizon T , no. of iterations m, link function φ, step size γ Output :Estimate Am
1 begin 2 A00 = 0 /*Initialization*/ 3 Ĝ← 1T ∑T−1 t=0 XtX > t ; If Ĝ is not invertible, then return Am = 0 4 for i← 0 to m− 1 do 6 Ai+1 ← Ai − 2γ (∇Lprox(Ai;X)) Ĝ−1
Note that when Assumption 4 holds with ρ < 1, the system necessarily mixes and converges to a stationary distribution as T →∞. Such systems forget their initial conditions in time scales of the order τmix = O( 1+logCρ
log 1 ρ
) = O (
1+logCρ 1−ρ
) , and hence we use this as a proxy for the mixing time. In
what follows, when we say ‘the system does not mix’ we either mean that it does not mix within time T or it does not converge to a stationary distribution (ex: ρ ≥ 1). Assumption 5 (Norm Boundedness). ‖A∗‖op = ρ < 1
That is, if A∗ satisfies Assumption 5, we have for arbitrary X,X ′ ∈ Rd: ‖φ(A∗X)− φ(A∗X ′)‖ ≤ ρ ‖X −X ′‖ and ∥∥(φ ◦A∗)k(X)∥∥ ≤ ρk ‖X‖. Hence, for such A∗, NLDS is necessarily stable.

heading: 3 Offline Learning with Quasi Newton Method
text: In this section we consider estimating A∗ using a single trajectory (X1, . . . , XT ) from NLDS(A∗, µ, φ). To this end, we study an offline Quasi Newton Method (Algorithm 1) where the iterates descend in the directions of the gradient of Lprox normalized by the inverse of the empirical second moment matrix Ĝ := 1T ∑T−1 t=0 XtX > t . That is, the iterates follow an approximation of the standard Newton update.
We now present analysis of Algorithm 1 in two settings: a) Theorem 1 provides estimation error for possibly unstable systems with sub-Gaussian noise that is close to the minimax optimal error incurred in the linear system identification case, b) Theorem 2 provides similarly tight estimation error for mixing systems but with heavy-tailed noise. Theorem 1 (Learning Without Mixing). Suppose Assumptions 1, 3 and 4 hold with expansivity factor ζ and regularity parameters (Cρ, ρ). Let C̄, C̄3 be constants depending only on Cη, and let
δ ∈ (0, 12 ). Let R ∗ := C2ρCηdσ
2 (∑T−1 t=1 ρ t )2 log( 4Tdδ ), and assume
1. The number samples T ≥ C̄3 ( d log ( R∗
σ2 ) + log 1δ ) 2. Step size γ = 14
3. m ≥ 10ζ · log ( ‖A0−A∗‖2F ·TR ∗ σ2d2 ) Then, the output Am of Algorithm 1 after m iterations and λmin ( Ĝ )
satisfy with probability at-least ≥ 1− δ:
‖Am −A∗‖2F ≤ C̄σ 2
Tζ2λmin(Ĝ)
[ d2 log ( 1 + R ∗
σ2
) + d log ( 2d δ )] ,
λmin
( Ĝ ) ≥ σ 2
2 .
Note that as λmin(Ĝ) & σ2, the error rate scales as ≈ d2/T , independent of τmix ≈ 1/(1− ρ). The theorem also holds for non-mixing or possibly unstable systems as long as ρ < 1 + CT . Furthermore, the error bound above is similar to the minimax optimal bound by [11] for the linear setting, i.e., when
φ(x) = x. As the link function φ tends to decrease the information in x, intuitively lower bound for linear setting should apply for NLDS as well, which would imply our error rate to be optimal; we leave further investigation into lower bound of NLDS identification for future work. Interestingly, in the linear case whenever the smallest singular value σmin(A∗) > 1 + , it can be show than λmin(Ĝ) grows exponentially with T , leading to an exponentially small error. It is not clear how to arrive at such a growth lower bound in the non-linear case.
The computational complexity of the algorithm scales as m · T which depends only logarithmically on τmix. Interestingly, the algorithm is almost hyperparameter free, and does not require knoweldge of parameters σ, τmix, ζ.
Also note that the stationary points of Algorithm 1 and GLMtron ( [32]) are the same. So, the stronger error rate in the result above compared to the result by [32] is due to a sharper analysis. However, in dynamical systems of the form 1, the squared norm of the iterates grow as d1−ρ even in the stable case. Hence, the GLMtron algorithm requires step sizes to be ≈ 1−ρd which implies significantly slower convergence rate for large τmix = 1/(1− ρ). In contrast, convergence rate for Algorithm 1 depends at most logarithmically on τmix.
Theorem 2 (Learning with Heavy Tail Noise). Suppose Assumptions 1 and 4 hold. In Assumption 4, let ρ < 1. Let X0, . . . , XT be a stationary trajectory drawn from NLDS(A∗, µ, φ) and Am be the m-th iterate of Algorithm 1. For some universal constants C,C1, C0 > 0, whenever δ ∈ (0, 12 ), R∗ := 4TdC2ρσ 2
(1−ρ)2δ and
1. T ≥ Cd log( 1δ ) log( R∗ σ2 ) max  4C6ρM4 (1−ρ)4(1−ρ2)σ4 , log ( R∗C1Cρ σ2 ) log ( 1 ρ ) 
2. Step size γ = 14 3. m ≥ 10ζ · log ( ‖A0−A∗‖2F ·TR ∗ σ2d2 ) There exists an eventW ∈ σ(X0, η0, . . . , ηT−1) with P(W) ≥ 1− δ and :
E [ ‖Am −A∗‖2F1(W) ] ≤ C0d 2σ2
ζ2Tλmin (G) .
Obtaining High Probability Bounds: The bound above shows that the expectation of the error restricted to a high probability set is small. This, along with Markov inequality, shows that we can have an error of at-most Cd 2σ2
ζ2Tλmin(G) with probability at-least 23 . This can be boosted to a high
probability bound by splitting the horizon T into K contiguous segments with a gap of O(τmix log T ) to maintain approximate independence (see Section C.2). We then run the Quasi Newton method on each of these ‘split’ data sets to obtain nearly independent estimates Â1, . . . , ÂK , which each have error at-most Cd
2σ2K ζ2Tλmin(G) with probability at-least 23 . Using a standard high-dimensional median
of means estimator (see [41, Algorithm 3]) for Â1, . . . , ÂK , we obtain error bounds of the order Cd2σ2K
ζ2Tλmin(G) with probability at least 1− e−Ω(K).
We refer to Section 6 for a high-level exposition of the key ideas in the analysis and Section B for the full proof of Theorems 1 and 2.

heading: 4 Streaming Learning with SGD-RER
text: In this section, we consider the one-pass, streaming setting, where the data points are presented in a streaming fashion. The goal is to continuously produce better estimates ofA∗ while also ensuring that the space and the time complexity of the algorithm is small. This disallows approaches that would just store all the observed points and then apply offline Algorithm 1 to produce strong estimation error. Such one-pass streaming algorithms are critical in a variety of settings like large-scale and online time-series analysis [42, 43], TD learning in RL [44], econometrics.
Algorithm 2: SGD− RER Input :Streaming data {Xτ}, horizon T , buffer size B, buffer gap u, bound R, tail start: t0 ≤ N/2, link function φ, step size γ Output :Estimate Ât0,t, for all t0 ≤ t ≤ N − 1; N = T/(B + u)
1 begin 2 Total buffer size: S ← B + u, Number of buffers: N ← T/S 3 A00 = 0 /*Initialization*/ 4 for t← 1 to N do 5 Form buffer Buft−1 = {Xt−10 , . . . , X t−1 S−1}, where, X t−1 i ← X(t−1)·S+i
6 If ∃i, s.t., ∥∥Xt−1i ∥∥2 > R, then return Ât0,t = 0 7 for i← 0 to B − 1 do 9 At−1i+1 ← A t−1 i − 2γ [ φ(At−1i X t−1 S−i−1)−X t−1 S−i ] Xt−1,>S−i−1
10 At0 = A t−1 B 11 If t ≥ t0 + 1, then Ât0,t ← 1t−t0 ∑t τ=t0+1 Aτ−1B
To address this problem, we consider SGD− RER (Algorithm 2) which was introduced in [19] in the context of linear system identification (LSI). We apply the method for NLDS identification as well. SGD− RER uses SGD like updates, but the data is processed in a different order than it is received from the dynamical system. This algorithm is based on the observation made in [19] that for LSI, when SGD is run on the least squares loss in the forward order, there are spurious correlations which prevent the algorithm’s convergence to the optimum parameter A∗. Surprisingly, considering the data in the reverse order exactly unravels these correlations to resolve the problem. Reverse order traversal of data, even though one pass, does not give a streaming algorithm. Hence, we divide the data into multiple buffers of size B and leave of size u between the buffers (See Figure 1). The data within each buffer is processed in the reverse order whereas the buffers themselves are processed in the order received. See Figure 1 for an illustration of the processing order. The gaps u are set large enough so that the buffers behave approximately independently. Setting B ≥ 10u we note that this simple strategy improves the sample efficiency compared to naive data dropping since we use most of the samples for estimating A∗. We now present the main result for streaming setting.
Theorem 3 (Streaming Algorithm). Suppose Assumptions 1, 2, 3 and 5 hold and that the data points are stationary. Set α = 100, R = 16(α+2)dCησ
2 log T 1−ρ , u ≥ 2α log T
log( 1 ρ ) , B ≥( C̄1 d (1−ρ)(1−ρ2) log ( d 1−ρ ) , 10u ) for a global constant C̄1 dependent only Cη and α. Let N = T/(B + u) be the number of buffers. Finally, set step-size γ = CT ν where ν = 6.5/7 and let T be
large enough such that γ ≤ min (
ζ 4BR(1+ζ) , 1 2R ) . If N/2 > t0 > c1 log TζγBλmin(G) = Θ(T ν log T ) for
some large enough constant c1 > 0, then output Ât0,N of Algorithm 2 satisfies:
E [ ‖Ât0,N −A∗‖2F ] ≤ C d 2σ2 log T
Tλmin (G) ζ2 + Lower Order Terms (5)
where C is a constant dependent on Cη , α.
Remark 1. The lower order terms are of the order Poly(R,B, β, 1/ζ, 1/λmin, ‖φ′′‖)γ7/2T 2 + γ2Rσ2d 1
Tα/2−2 + ‖A0 −A∗‖2F
[ e−c2ζγBλmint0
Tζγλmin
] . We refer to the proof in Section C.13 for details.
Remark 2. Although the bound in Theorem 3 is given for the algorithmic iterate at the end of the horizon, the proof shows that in fact we can bound the error of the iterates at the end of each buffer after (1 + c)t0 i.e. if t ≥ (1 + c)t0 for some c > 0 then we obtain
E [ ‖Ât0,t −A∗‖2F ] ≤ C d 2σ2 log T
(tB)λmin (G) ζ2 + Lower Order Terms
Note that the estimation error above matches the error by offline method up to log factors (see Theorem 2). Furthermore, while the method requires NLDS to be mixing, i.e., ρ < 1, but the leading term in error rate does not have an explicit dependence on it. Moreover, the space complexity of the method is only B ·d which scales as d2/((1−ρ)(1−ρ2)), i.e., it is 1/((1−ρ)(1−ρ2)) ∼ τ2mix factor worse than the obvious lower bound of O(d2) to store A. We leave further investigation into space complexity optimization or tightening the lower bound for future work. Also, note that u ≤ B/10, so SGD− RER wastes only about 10% of the samples. Finally, the algorithm requires a reasonable upper bound on ρ to set up various hyperparameters like R, u,B. However, it is not clear how to estimate such an upper bound only using the data, and seems like an interesting open question.
See Section 6 for an explanation of the elements involved in the analysis of the algorithm and to Section C.1 for a detailed overview of the proof.

heading: 5 Exponential Lower Bounds for Non-Expansive Link Functions
text: The previous results showed that we can efficiently recover the matrix A∗ given that the link function is uniformly expansive. We now consider non-expansive functions and show that parameter recovery is hard in this case. In particular, we show that even for the case of φ = ReLU, the noise being N (0, I), and ‖A∗‖ ≤ 12 , the error has an information theoretic lower bound which is exponential in the dimension. We note that this is consistent with Theorem 3 in [32] which too has an exponential dependence on the dimension (since the matrix K I). Before stating the results, we introduce some notation. Consider any algorithm A, with accepts input (X0, . . . , XT ) and outputs an estimate Â ∈ Rd×d. For simplicity of calculation, we will assume that X0 = 0 and Xt+1 = ReLU(A∗Xt) + ηt. Since the mixing time is O(1), similar results should hold for stationary sequences. We define the loss L(A, T, A∗) = E‖Â − A∗‖2F , where the expectation is over the randomness in the data and the algorithm. By Θ( 12 ), we denote all the the elements of B ∈ Rd×d such that ‖B‖ ≤ 12 . The minimax loss is defined as:
L(Θ( 12 ), T ) := infA sup A∗∈Θ( 12 ) L(A, T, A∗) .
Theorem 4 (ReLU Lower Bound). For universal constants c0, c1 > 0, we have: L(Θ( 12 ), T ) ≥ c0 min ( 1, exp(c1d)
T
) .
We prove the theorem above using the two point method. We find a family of A∗ in Θ( 12 , T ) such that 〈Xt, ed〉 = 〈ηt−1, ed〉 with probability at-least 1− exp(−Ω(d)). Therefore, with a large probability, we only observe noise in the last co-ordinate and hence do not obtain any information regarding the last row of A∗ (i.e, a∗d). We refer to Section G for a full proof.

heading: 6 Proof Sketch
text: Quasi Newton Method. Let a∗i be the i-th row of A∗ and ai(l) be the i-th row of Al both in column vector form. The proofs of Theorems 1 and 2 follow once we consider the lyapunov function ∆l,i = ‖Ĝ1/2(ai(l)− a∗i )‖ and show that
∆l+1,i ≤ (1− 2γζ)∆l + γ‖Ĝ−1/2N̂i‖ (6)
Where N̂i := 1T ∑T−1 t=0 〈ei, ηt〉Xt. In the case of Theorem 1, we use the sub-Gaussianity of the noise
sequence and a martingale argument to obtain a high probability upper bound on ∑d i=1 ‖Ĝ−1/2N̂i‖2
(see Lemma 1). In the heavy tailed case considered in Theorem 2, we use mixing time arguments along with Payley-Zygmund inequality to show the high probability lower isometry Ĝ c0G for some universal constant c0 (see Lemma 2). Using this lower isometry, we can replace Ĝ in Equation (6) with G. The upper bounds follow once we note that EN̂>i G−1N̂i = σ 2d T .
SGD-RER. Due to the observations made in Section 4, we can split the analysis into the following parts, which are explained in detail below.
1. Analyze the reverse order SGD within the buffers. 2. Treat successive buffer as independent samples. 3. Give a bias-variance decomposition similar to the case of linear regression. 4. Use algorithmic stability to control ’spurious’ coupling introduced by non-linearity in the bias-
variance decomposition.
Coupled Process. We deal with the dependence between buffers using a fictitious coupled process, constructed just for the sake of analysis (see Definition 2). Leveraging the gap u, this process (X̃τ ) is constructed such that X̃τ ≈ Xτ with high probability and the ‘coupled buffers’ containing data X̃ instead of X are exactly independent. Since X̃τ ≈ Xτ , the output of SGD− RER run with the fictitious coupled process should be close output of SGD− RER run with the actual data points. We then use the strategy outlined above to analyze SGD− RER with the coupled process. In the analysis given for SGD− RER, all the quantities with ·̃ involve the coupled process X̃ instead of the real process X .
Non-Linear Bias Variance Decomposition. We use the mean value theorem to linearize the nonlinear problem. This works effectively when the step size γ is a vanishing function of the horizon T . Observe that the a single SGD/ SGD− RER step for a single row can be written as:
a′i − a∗i = ai − a∗i − 2γ(φ(〈ai, Xτ 〉)− φ(〈a∗i , Xτ 〉))Xτ + 2γ〈ητ , ei〉Xτ = ( I − 2γφ′(βτ )XτX>τ ) (ai − a∗i ) + 2γ〈ητ , ei〉Xτ (7)
In the second step, we have used the mean value theorem. Equation (7) can be interpreted as follows: the matrix ( I − 2γφ′(βτ )XτX>τ ) ’contracts’ the distance between ai and a∗i whereas the noise 2γ〈ητ , ei〉Xτ is due to the inherent uncertainty. This gives us a bias-variance decomposition similar to the case of SGD with linear regression. We refer to Section C.5 for details on unrolling the recursion in Equation (7) to obtain the exact bias-variance decomposition.
Algorithmic Stability: Unfortunately, non-linearities result in a ‘coupling’ between the contraction matrices through the iterates via the first derivative φ′(βτ ) due to reverse order traversal. This is an important issue since unrolling the recursion in (7), we encounter terms such as 〈ητ , ei〉(I − 2γφ′(βτ−1)Xτ−1X > τ−1)Xτ , which have zero mean in the linear case. However, in the non-linear case, βτ−1 depends on ητ due to reverse order traversal. We show that such dependencies are ‘weak’ using the idea of algorithmic stability ([45, 46]). In particular, we establish that the output of the algorithm is not affected too much if we re-sample the entire data trajectory by independently re-sampling a single noise co-ordinate (ητ becomes η′τ and βτ−1 becomes β ′ τ−1) when the step size γ is small enough (in other words, the output is stable under small perturbations). Via second derivative arguments, we show that βτ−1 ≈ β′τ−1. Now observe that resampling noise ητ does not affect the past value of data i.e, Xτ , Xτ−1 and is independent of β′τ−1 by construction. Therefore
0 = E〈ητ , ei〉 ( I − 2γφ′(β′τ−1)Xτ−1X>τ−1 ) Xτ ≈ E〈ητ , ei〉 ( I − 2γφ′(βτ−1)Xτ−1X>τ−1 ) Xτ
Such a resampling procedure is also explored in [47] for the analysis of SGD with random reshuffling.
We put together all the ingredients above in order to prove the error bounds given in Theorem 3.

heading: 7 Experiments
text: In this section, we compare performance of our methods SGD− RER and Quasi Newton method on synthetic data against the performance of standard baselines SGD (called ‘Forward SGD’ here),
GLMtron, along with the SGD− ER method that applies standard experience replay technique i.e, the points from a buffer are sampled randomly instead of the reverse order. Since GLMtron and Quasi Newton Method are offline and SGD− RER, SGD and SGD− ER are streaming, we compare the algorithms by plotting parameter error measured by the Frobenius norm with respect to the compute time. We also compare error vs. number of iterations for the streaming algorithms. We show the results of additional experiments by considering various buffer sizes and heavy tailed noise in Section A.
Synthetic data: We sample data from NLDS(A∗, µ, φ) where µ ∼ N (0, σ2I) and A∗ ∈ Rd×d is generated from the "RandBiMod" distribution. That is, A∗ = UΛU> with random orthogonal U , and Λ is diagonal with dd/2e entries on diagonal being ρ and the remaining diagonal entries are set to ρ/3. φ is the leaky ReLU function given by φ(x) = 0.5x1(x < 0) + x1(x ≥ 0). We set d = 5, ρ = 0.98 and σ2 = 1. We set a horizon of T = 105.
Algorithm Parameters We set B = 240 and u = 10 for the buffer size and gap size respectively for both SGD− RER and SGD− ER and use full averaging (i.e, θ = 0 in Algorithm 2 ). We set the step size γ = 5 log TT for SGD, SGD− RER, and SGD− ER and γnewton = 0.2 and γGLMtron = 0.017. From Figure 2 observe that SGD− ER and SGD obtain sub-optimal results compared SGD− RER, Quasi Newton Method and GLMtron. After a single pass, the performance of SGD− RER almost matches that of the offline algorithms. The step sizes for GLMtron have to be chosen to be small in-order to ensure that the algorithm does not diverge as noted in Section 3, which slows down its convergence time compared to the Quasi Newton method. We set the step size to be as large as possible without obtaining divergence to infinity.

heading: 8 Conclusion
text: In this work, we studied the problem of learning non-linear dynamical systems of the form (1) from a single trajectory and analyzed offline and online algorithms to obtain near-optimal error guarantees. In particular we showed that mixing time based arguments are not necessary for learning certain classes of non-linear dynamical systems. Even though we show that one cannot hope for efficient parameter recovery with non-expansive link functions like ReLU, we do not deal with the problem of minimizing the ‘prediction’ error - where we output a good predictor based on samples, without parameter recovery. We believe that this problem would require significantly different set of techniques than the ones established in this work and hope to investigate this in future work. Presently our work only deals with specific kinds of Markovian time evolution. It would also be interesting to understand in general, the kind of structures which allow for learning without mixing based arguments. Another related direction is to design simple and efficient algorithms like SGD− RER for learning with various models of dependent data by unraveling the dependency structure present in the data.
Acknowledgments and Disclosure of Funding
D.N. was supported in part by NSF grant DMS-2022448. S.S.K was supported in part by Teaching Assistantship (TA) from EECS, MIT.

Review summaries:
The authors consider non-linear dynamical systems of the form $X_{t+1} = \phi(A^* X_t) + \eta_t$, where $\eta_t$ is a noise term, $A$ is a $d \times d$ matrix and $\phi$ acts component-wise, and provide some results relating to the estimation of $A^* $ under certain assumptions. A key assumption is that $\phi$ is "expansive", which I understand to be a sort of converse Lipschitz condition. (Since 1-Lipschitzness is also assumed, this means that $\zeta |x - y| \le |\phi(x) - \phi(y)| \le |x - y|$ for some $\zeta \in (0, 1)$. I wasn't sure why this was called "expansive", because it seems to me more to be a limit on contraction than expansion per se, but this seemed understandable.) Under these conditions, they provide:  1. An algorithm that, given sub-Gaussian noise, can recover $A^*$ with squared-norm error scaling in $O(d^2/T)$ (though I have questions about this, see main review). 2. An algorithm that, given system stability and finite-fourth-moment noise, can recover $A^*$ with squared-norm error scaling in $O(d^2/T)$ (though I have questions about this too). 3. A _streaming_ algorithm that, given sub-Gaussian noise, system stability and twice-differentiable $\phi$, can recover $A^*$ with squared-norm error scaling in $O(d^2/T)$. 4. A negative result with the non-expansive ReLU function, showing a lower bound on squared-norm error in this case of $\Theta(e^{cd}/T)$, implying a need for exponentially many samples to recover $A^*$ satisfactorily. 

This work has four major contributions: a) authors provide the first offline algorithm that can learn non-linear dynamical systems without the mixing assumption, b) authors significantly improve upon the sample complexity of existing results for mixing systems, c) in the much harder one-pass, streaming setting authors study a SGD with Reverse Experience Replay (SGD − RER) method, and demonstrate that for mixing systems, it achieves the same sample complexity as our offline algorithm, d) authors justify the expansivity assumption by showing that for the popular ReLU link function — a non-expansive but easy to learn link function with i.i.d. samples — any method would require exponentially many samples (with respect to dimension of Xt) from the dynamical system.   This paper is comprehensive with both solid theoretical analysis and simulations.  

The authors study a very special case of "non-linear dynamical systems" where there is no hidden state and there is some expansive, component-wise "link function" \phi in the evolution X(t+1) = \phi(AX(t)) + noise. For this very special case, they show both off-line and on-line algorithms. They also show lower bounds on the sample complexity of learning systems with with non-expansive link functions (such as ReLU).  

This work considers learning NLDS. Previous works can only learn mixing systems with sub-gaussian noise in the offline setting, with sample complexity depending on the mixing time. This work improves existing results in several ways, assuming the link function is expansive: 1, it's shown that the offline Quasi Newton Method achieves near-optimal convergence rate for non-mixing systems. 2, the offline Quasi Newton Method still achieves near-optimal convergence rate with heavy-tail noise, assuming mixing systems. 3, a new algorithm SGD-RER is introduced, which achieves near-optimal convergence rate in the online setting. 4, an exponential lower bound is proven for the ReLU link function, showing the expansivity assumption is necessary.