Streaming Linear System Identification with Reverse Experience Replay
Year: 2021

Abstract:
We consider the problem of estimating a linear time-invariant (LTI) dynamical system from a single trajectory via streaming algorithms, which is encountered in several applications including reinforcement learning (RL) and time-series analysis. While the LTI system estimation problem is well-studied in the offline setting, the practically important streaming/online setting has received little attention. Standard streaming methods like stochastic gradient descent (SGD) are unlikely to work since streaming points can be highly correlated. In this work, we propose a novel streaming algorithm, SGD with Reverse Experience Replay (SGD − RER), that is inspired by the experience replay (ER) technique popular in the RL literature. SGD− RER divides data into small buffers and runs SGD backwards on the data stored in the individual buffers. We show that this algorithm exactly deconstructs the dependency structure and obtains information theoretically optimal guarantees for both parameter error and prediction error. Thus, we provide the first – to the best of our knowledge – optimal SGD-style algorithm for the classical problem of linear system identification with a first order oracle. Furthermore, SGD− RER can be applied to more general settings like sparse LTI identification with known sparsity pattern, and non-linear dynamical systems. Our work demonstrates that the knowledge of data dependency structure can aid us in designing statistically and computationally efficient algorithms which can “decorrelate” streaming samples.

heading: 1 Introduction
text: In this paper, we study the problem of learning linear-time invariant (LTI) systems, where the goal is to estimate the matrix A∗ ∈ Rd×d from the given samples (X0, . . . , XT ) that obey:
Xτ+1 = A ∗Xτ + ητ , Xτ ∈ Rd, ητ i.i.d.∼ µ, (1) where µ is an unbiased noise distribution. The problem is central in control theory and reinforcement learning (RL) literature [1, 2]. It is also equivalent to estimating Vector Autoregressive (VAR) model popular in the time-series analysis literature [3], where it has been used in several applications like finding gene regulatory information network [4].
Despite a long line of classical literature for the problem, most of the existing results focus on the offline setting, where all the samples (X0, . . . , XT ) are available apriori. In this setting, ordinary 35th Conference on Neural Information Processing Systems (NeurIPS 2021).
least squares (OLS) method that estimates A as, Â = argminA ∑T−1
τ=0 ∥Xτ+1 −AXτ∥2 is known to be nearly optimal [5, 6]. However, such offline solutions do not apply to the streaming setting – where A∗ needs to be estimated online – that has applications in several domains like RL, large-scale forecasting systems, recommendation systems [7, 8].
In this paper, we study the above mentioned problem of learning LTI systems via first order gradient oracle with streaming data. The goal is to design an estimator that provides accurate estimation while ensuring nearly optimal time complexity and space complexity that is nearly independent of T . Note that due to specific form arising in linear regression, the optimal solution to OLS can be estimated in online fashion using Sherman-Morrison-Woodbury formula. But such a solution is limited and does not apply to practically important settings like generalized non-linear dynamical system or when A∗ is high-dimensional and has special structure like low-rank or sparsity [9, 10].
So, in this work, we focus on designing Stochastic Gradient Descent (SGD) style methods that can work directly with first order gradient oracle, and hence is more widely applicable to the settings mentioned above. In fact, after the first appearance of this manuscript, the algorithm (SGD− RER) and the techniques introduced in this paper were used to obtain near-optimal guarantees for learning certain classes of non-linear dynamical systems [11] as well as in Q-learning tabular MDPs in RL [12]. We note that prior to [11], even optimal offline algorithms were unknown for such non-linear systems.
SGD is a popular method for general streaming settings, and has been shown to be optimal for problems like streaming linear regression [13]. However, when the data has temporal dependencies, as in the estimation of linear dynamical systems, such a naive implementation of SGD may not perform well as observed in [14, 15]. In fact, for linear system identification, our experiments suggest that SGD suffers from a non-zero bias (Section 6). In order to address temporal dependencies in data, practitioners use a heuristic called experience replay, which maintains a buffer of points, and samples points randomly from the buffer. However, for linear system identification, experience replay does not seem to provide an accurate unbiased estimator for reasonable buffer sizes (see Section 6).
In this work, we propose reverse experience replay for linear system identification. Our method maintains a small buffer of points, but instead of random ordering, we replay the points in a reverse order. We show that this algorithm exactly unravels the temporal correlations to obtain a consistent estimator for A∗. Similar to the standard linear regression problem with i.i.d. samples, we can break the error in two parts: a) bias: that depends on the initial error ∥A0 −A∗∥, b) variance: the steady state error due to noise η. We show that our proposed method, under fairly standard assumptions and with a small buffer size, is able to decrease the bias at fast rate, while the variance error is nearly optimal (see Theorem 1), matching the information theoretic lower bounds [5, Theorem 2.3]. To the best of our knowledge, we provide first non-trivial analysis for a purely streaming SGD-style algorithm with optimal computation complexity and nearly bounded space complexity that is dependent logarithmically on T . We note here that the idea of reverse experience replay was independently discovered in experimental reinforcement learning by [16] based on reverse replay observed in Hippocampal place cells [17] in Neurobiology. We also refer to [18] for more on this connection.
In addition to the transition matrix estimation error ∥A−A∗∥, we also provide analysis of prediction error, i.e., E[∥AX −A∗X∥2] (see Theorem 2). Here again, we bound the bias and the variance part of the error separately. We further derive new lower bounds for prediction error (see Theorem 4) and show that our algorithm is minimax optimal, under standard assumptions on the model. As mentioned earlier, our method work with general first order oracles, hence applies to more general problems like sparse LTI estimation with known sparsity structure and unlike online OLS methods, SGD− RER has nearly optimal time complexity. Finally, we also provide empirical validation of our method on simulated data, and demonstrate that the proposed method is indeed able to provide error rate similar to the OLS method while methods like SGD and standard experience replay, lead to biased estimates.
Related Work. Due to applications in RL, recently LTI system identification has been widely studied. In particular, [19] studied the problem in offline setting under the “stability" condition, i.e., the spectral radius (ρ(A∗)) of A∗ is a constant bounded away from 1. The sequence of papers [5, 6, 20, 21] provide optimal analyses of the offline OLS estimator beyond assumptions of stability. That is, they show that OLS recovers A∗ near optimally even the process defined by (1) is stable but does not mix within time T (when ρ(A∗) is 1− O(1/T )) or is unstable (when ρ(A∗) is larger
than 1). Further [5, 22] provide information theoretic lower bounds for the LTI system identification problem. [11, 23, 24] consider the problem of identifying non-linear dynamical systems of the form Xt+1 = ϕ(A
∗Xt) + ηt where ϕ is a one dimensional link function which acts co-ordinate wise. In this setting, however, there is no closed for expressions for the estimator of A∗. [23, 24] give offline algorithms whose error guarantees are worse off by factors of mixing time whereas [11] obtains near optimal offline and streaming algorithms for this setting. In fact, [11] uses SGD− RER which was first introduced in this work in order to obtain the streaming algorithm.
LTI identification problem has been studied in time series forecasting literature as well. For example, [25] obtains asymptotic consistency results for system identification problem and [26, 27] consider the problem of finite time recovery. Both consider a certain parameterized predictor for a linear system with empirical risk minimization for the parameter and analyzes the deviation from population risk. Similarly, [28] also studies generalization error guarantees. In contrast, our work is able to provide precise bias and variance (similar to generalization error) of the estimator in the streaming setting, and show that the asymptotic error is minimax optimal.
[29] studied SISO systems with observations (xτ , yτ ) ∈ R2 and a hidden state hτ which is high dimensional, thus their model and applications are significantly different than the LTI system we study. For the SISO system, [29] analyzes SGD to provide error bounds contain (a large) polynomial in the hidden state dimension. Here, the hidden state has an evolution similar to Equation 1 whereas x1, . . . , xT are drawn i.i.d from some distribution.
System identification has been studied in the context of partially observed LTI systems as well. Recent works [19, 30–34] focus on identifying a certain Hankel-like matrix of the system. These are not directly comparable to the fully observed setting in this work since the model parameters are identifiable only upto a similarity transformation in the partially observed setting.
Recently, there has been an exciting line of work in the related domain of online control (see [35–38] and references therein). The state equation studied in these papers also contain an additive term of Buτ for some unknown matrix B and a control signal uτ and the noise ητ is either stochastic (as in [35]) or adversarial (as in [36–38]). The goal is to output control signals uτ after observing X1, . . . , Xτ , such that the cost ∑ τ cτ (Xτ , uτ ) is minimized for some sequence of convex costs cτ . We focus on the LTI system identification(or estimation) problem while the goal of the above mentioned line of work is to design an online controller.
We also note here another line of works [32, 39–44] focused on online prediction of both fully observed and partially observed LTI systems, and the similar problem of time series forecasting by regret minimization [28, 45]. In particular, the main goal there is to design online prediction algorithms minimizing regret against a certain class (for instance, against a Kalman filter with knowledge of the system parameters in the case of partially observed LTI systems). The situation considered in our work is different in atleast two aspects: 1) we focus significantly on parameter recovery or system identification and 2) our notion of prediction is prediction at stationarity which can be thought of as one-step regret (compared to T–step regret for instance in [39, 40]).
Finally, [9] considers offline sparse linear regression with ℓ1 penalty where the feature vector is derived from an auto regressive model. Similarly, [14] considers the problem of linear regression where the feature vectors come from a Markov chain. This line of work is different from ours in that we try to estimate the parameters of the Markov process itself.
Paper Organization. We provide the problem definition and introduce the notations in the next section. We then present our algorithm and the key intuition behind it in Section 3. We then present our main result in Section 4 and provide a proof sketch in Section 5. Finally, we present simulation results in Section 6.

heading: 2 Problem Setting and Notation
text: In this section, we first introduce the data generation model, the required assumptions and then provide the precision problem definition. Throughout the paper, we use ∥A∥ to denote the operator norm of A unless otherwise specified. ∥A∥F denotes the Frobenius norm of A. σi(A) denotes the i-th largest singular value of A, i.e., σmax(A) = σ1(A). κ(A) := σmax(A)/σmin(A) denotes the condition number of A. ρ(A) denotes the spectral radius of A. For two symmetric matrices
A,B ∈ Rd×d we say A ⪯ B if B −A is positive semidefinite (psd). For notational simplicity, we use C to denote a constant, and it’s value can be different in different equations.
Linear Dynamical System/VAR(1) model. Given an initial (possibly random) data point X0 which is independent of the noise sequence, we generate the (X0, . . . , XT ) from the VAR model as:
Xτ+1 = A ∗Xτ + ητ , 0 ≤ τ ≤ T − 1, (2)
where A∗ ∈ Rd×d be the transition matrix. Let η1, . . . , ηT ∈ Rd be an i.i.d noise sequence with 0 mean and finite second moment with probability measure µ. We will denote this model by VAR(A∗, µ). We also make the following assumptions about A∗, µ, and X0: Assumption 1. External Stability. ∥A∗∥ < 1 Assumption 2. Sub-Gaussian Noise. µ has co-variance Σ and for all x ∈ Rd, ⟨x, ητ ⟩ is Cµ⟨x,Σ ·x⟩ sub-Gaussian. Further, Σ is full rank. Also, let µ4 := E [ ∥ητ∥4 ] be the fourth moment of the noise.
Assumption 3. Stationarity. X0 ∼ π, the stationary distribution corresponding to (A∗, µ). Let M4 := E [ ∥X0∥4 ] .
Due to Assumption 1, we can show that the law of the iterate XT from the VAR model defined above converges to a stationary distribution π as T → ∞ for arbitrary choice of X0 and has a mixing time of the order τmix = O ( 1
1−∥A∗∥
) . For simplicity, we will absorb Cµ into other
constants. Finally, we will use (Z0, . . . , ZT ) ∼ VAR(A∗, µ) to mean that Z0, . . . , ZT is a stationary sequence corresponding to VAR(A∗, µ). We also note that the covariance matrix under stationarity, G := EX∼πXX⊤ = ∑∞ s=0 A
∗sΣ(A∗⊤)s ⪰ Σ. Remark. It is indeed possible to replace Assumption 1 with the weaker condition on the spectral radius of A∗: ρ(A∗) < 1. While our results still hold in this case, the bound might have additional condition number factors. See Section A.1 for more details.
Remark. The full rank assumption on Σ is needed for polynomial sample complexity [46].
Problem Statement. Let (X0, X1, · · · , XT ) be sampled from VAR(A∗, µ) model for a fixed horizon T . Then, the goal is to design and analyze an online algorithm that uses only first order gradient oracle to estimate the system matrix A∗. That is, at each time-step τ , we obtain gradient for the transition (Xτ , Xτ+1) and output estimate Aτ . The goal is to ensure that each Aτ has small estimation error wrt A∗; naturally, we would expect better estimation error with increasing τ . We quantify estimation error using the following two loss functions:

heading: 1. Parameter error: Lop(A;A∗, µ) = ∥A−A∗∥
text: 2. Prediction error at stationarity: Lpred(A;A∗, µ) := EXτ∼π∥Xτ+1 −AXτ∥2
Note that the problem is equivalent to d linear regression problems, but with dependent samples, making it significantly more challenging. Whenever Assumption 1 holds, stationary distribution π exists, so the prediction error Lpred is meaningful. Furthermore: Lpred(A) − Lpred(A∗) = Tr [ (A−A∗)⊤(A−A∗)G ] where G := EX∼πXX⊤.

heading: 3 Algorithm
text: As mentioned in related works, the standard OLS estimator that minimizes the empirical loss is known to be nearly optimal in the offline setting [5]:
ÂOLS = argmin A T−1∑ τ=0 ∥AXτ −Xτ+1∥2 . (3)
Note that for least squares loss, one can indeed maintain covariance matrix and residual vector to compute the OLS solution online. But such a solution does not work if we have access to only gradients and breaks down even for generalized linear models, whereas as the techniques introduced in this work has been extended to non-linear systems [11].
On the other hand, using standard SGD we can obtain update to A efficiently by using gradient at the current point. That is, assuming A0 = 0, we get the following SGD update (for all τ ≥ 0):
Aτ+1 = Aτ − 2γ(AτXτ −Xτ+1)X⊤τ , (4) where γ is the stepsize. While SGD is known to be an optimal estimator in certain streaming problems with i.i.d. data, for the VAR(A∗, µ) problem the standard SGD does not apply, as samples (Xτ , Xτ+1) and (Xτ+1, Xτ+2) are highly correlated. To see why this is the case, let us unroll the recursion for two steps and using Equation (2):
A2 −A∗ = (A0 −A∗)(I − 2γX0X⊤0 )(I − 2γX1X⊤1 ) + 2γη1X⊤1 + 2γη0X⊤0 (I − 2γX1X⊤1 ).
Note that the last term does not have 0 mean because X1 depends on η0 by Equation (2). Even in the case when A0 = A∗, this means that EA2 ̸= A∗ in general. In fact, in Section 6, we show empirically that SGD with constant step-size converges to a significantly larger error than OLS, even when T is very large. This shows that we cannot naively treat this problem as a collection of d linear regressions. This is consistent with the results in [14, 15] which show a similar behavior for constant step-size SGD with dependent data. Now, one can use techniques like data drop that drops a large fraction of points (either explicitly or during the mathematical analysis) from the stream to obtain nearly independent samples [14, 47], but such methods waste a lot of samples and have significantly suboptimal error rate than OLS.
So, the goal is to design a streaming method for the problem of learning dynamical systems that at each time-step t provides an accurate estimate of A∗, while also ensuring small space+time complexity.We now present a novel algorithm that addresses the above mentioned problem.

heading: 3.1 SGD with Reverse Experience Replay
text: We now discuss a novel algorithm called SGD with Reverse Experience Replay (SGD− RER) that addresses the problem of learning stationary auto-regressive models (or linear dynamical systems) in the streaming setting. Our method is inspired by the experience replay technique [48], used extensively in RL to break temporal correlations between dependent data. We make the following crucial observation. Suppose in Equation (4), instead of processing the samples in the order (X1, X2) → (X2, X3) → · · · → (XT−1, XT ), we process it in the reverse order. That is: (XT−1, XT ) → (XT−2, XT−1)→ · · · → (X1, X2). Then,
A2 −A∗ = (A0 −A∗)(I − 2γXT−1X⊤T−1)(I − 2γXT−2X⊤T−2) + 2γηT−2X⊤T−2 +2γηT−1X ⊤ T−1(I − 2γXT−2X⊤T−2) (5)
Now, observe that (XT−2, XT−1) are independent of ηT−1. Therefore the problematic last term, 2γηT−1X ⊤ T−1(I − 2γXT−2X⊤T−2), now has expectation 0. So the updates for reverse order SGD would be unbiased. This, however, requires us to know all the data points beforehand which is infeasible in the streaming setting. We alleviate this issue by designing SGD− RER, which is the online variant of the above algorithm. SGD− RER uses a buffer of large enough size to store values of consecutive data points and then performs reverse SGD in each of these buffers and then discards this buffer. Experience replay methods also use such (small) buffers of data, but typically samples point randomly from the buffer instead of the reverse order that we propose. We refer to Figure 1 for an illustration of the proposed data processing order.
We present a pseudocode of SGD − RER in Algorithm 1. Note that the algorithm forms nonoverlapping buffers of size S = B + u. Here B is the actual size of the buffer while u samples are used to interleave between two buffers so that the buffers are almost independent of each other. Now
Algorithm 1: SGD− RER Input :Streaming data {Xτ}, horizon T , buffer size B, buffer gap u, bound R, tail average start: a Output :Estimate Âa,t, for all a < t ≤ N − 1; N = T/(B + u)
1 begin 2 Step-size: γ ← 18RB , 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̂a,t = 0 7 for i← 0 to B − 1 do 9 At−1i+1 ← A t−1 i − 2γ(A t−1 i X t−1 S−1−i −X t−1 S−i) ( Xt−1S−1−i
)⊤ 10 end 11 At0 = A t−1 B
12 If t > a, then Âa,t ← 1t−a ∑t τ=a+1 A τ−1 B 13 end 14 end
within a buffer, we perform the usual SGD but with samples read in reverse order. Formally, suppose we index our buffers by t = 0, 1, 2, · · · and let S = B + u be the total samples (including those that were dropped) in the buffers. Let N denote the total number of buffers in horizon T . Within each buffer t, we index the samples as Xti where i = 0, 1, 2, · · · , S − 1. That is Xti ≡ XtS+i is the i-th sample in buffer t. Similarly ηti ≡ ηtS+i. Further let Xt−i ≡ Xt(S−1)−i. Similarly we set ηt−i ≡ ηt(S−1)−i Then, the algorithm performs the recursion stated in Line 1 of Algorithm 1. Note that the recursion can also be written as,
At−1i+1 −A ∗ = ( At−1i −A ∗) (I − 2γXt−1−i Xt−1−i ⊤)+ 2γηt−1−i Xt−1−i . (6) for 1 ≤ t ≤ N and 0 ≤ i ≤ B − 1 with At0 = At−1B and A00 = A0. We then ignore the iterates corresponding to first a buffers as part of the burn-in period, and output average of the remaining iterates (t > a) at each step as that step’s estimator (see Line 2 of Algorithm 1). That is, we have the tail-averaged iterate:
Âa,t = 1
t− a t∑ τ=a+1 Aτ−1B . (7)
We output the new iterate Âa,t only at the end of each buffer t. At intermediate steps, (t− 1)B+1 ≤ τ ≤ tB, we output Âa,t−1. Also, note that the tail average can be computed in small space and time complexity, by using a running sum of the tail iterates. The update for each point is rank-one, so can be computed in time linear in number of parameters (O(d2)). In the next section, we show that despite using small buffer size S = B + u (that depends logarithmically on T ), and by throwing away a small constant–independent of any problem parameter–fraction of points u in each buffer, we are still able to provide error bound similar to that of OLS.

heading: 4 Main Results
text: We now state our main results with leading order terms. For simplicity, we only state the results for the tail average ÂN
2 ,N but a similar result holds for any Âa,t when a = Ω(dBκ(G) log2 T ). We
refer to Section A for complete statements. Recall the problem setting, and the covariance matrix G := EX∼πXX⊺. Before stating the results, we choose the parameters B,R, α and u as follows, which can be estimated using upper bounds on ∥A∗∥: 1. d ≤ Poly(T ). We use this to bound the norm of covariates in the next item.

heading: 2. α ≥ 22 ; R ≥ C(α)Tr(Σ) log T
text: 1−∥A∗∥2 = O(dτmix log T ) s.t. P [ ∥Xτ∥2 ≤ R, τ ≤ T ] ≥ 1− 1Tα . See
lemma 9 in appendix.
3. u ≥ α log T log
( 1 ∥A∗∥ ) = O(τmix log T ); B = 10u
For all the results below, we suppose that Assumptions 1, 2 and 3 hold, the stream of samples Xτ is sampled from VAR(A∗, µ) model described in Section 2 and that R,B, α and u are chosen as above. Further we hide some mild conditions on N and T .
Theorem 1 (Informal version of Theorem 5). Let the step size γ < min (
C Bσmin(G) , 18BR
) for some
constant C depending only on Cµ. Then, with probability at least 1− 1T 100 , we have:
Lop(ÂN 2 ,N , A∗, µ) ≤ C
√ (d+ log T )σmax(Σ)
Tσmin(G) + Lower Order Terms .
Theorem 2 (Informal version of Theorem 6). Consider the setting of Theorem 1 but where the step size γ = min ( 1 2R , c BR ) for some constant 0 < c < 1. Then, the following holds:
E [ Lpred(ÂN
2 ,N ;A
∗, µ) ] − Tr(Σ) ≤ C dTr(Σ)
T + Lower Order Terms
where “lower order” is with respect to dT .
See Section F.1, Section F.3 for a detailed proof of the parameter error bound and see Section G.1, Section G.2 for a detailed proof of the prediction error bound.
We now make the following observations:
(1) The dominant term in our bound on Lop (Theorem 1) matches the information theoretically optimal bound (up to logarithmic factors) for the VAR(A∗, µ) estimation problem [5] as long as ∥A∗∥ ≤ 1− 1
T ξ for ξ ∈ (0, 1/2). Note that despite working with dependent data, leading term
in our error bound is nearly independent of mixing time τmix. In contrast, most of the existing streaming/SGD style methods for dependent data have strong dependence on τmix [14]. (2) SGD for linear regression with independent data [13, 49], but with similar problem setting incurs error O(dTr(Σ)T ) for Lpred. So our bound for SGD− RER matches the independent data setting bound in the minimax sense. (3) The space complexity of our method is O(Bd+ d2) where B = O(τmix log T ) is independent of d and only logarithmically dependent on T . (4) Sparse matrices with known support: Suppose A∗ is known to be sparse and we know the support (say by running L1 regularized OLS on a small set of samples). Let sj denote the sparsity of row j of A∗. Then the SGD− RER algorithm can be modified to run row by row such that it operates only on the support of row j. That is the covariates can be projected onto the support of each row. Then it can be shown that the prediction error is bounded as O (∑d j=1 σ 2 j sj/T
) where σ2j is the j-th diagonal entry of Σ. Note that SGD − RER requires only O(|supp(A∗)|) operations per iteration while applying online version of standard OLS would require O(d2) operations. In the simple case of Σ = σ2I , we note that G ⪰ σ2I and hence the bound for Lpred becomes O ( |supp(A∗)|
T
) . We refer to Section O for a sketch of this extension.
Next, we show that our error bounds are nearly information theoretically optimal. For the lower bound on Lop we directly use [5, Theorem 2.3]. Theorem 3. Let ρ < 1 and δ ∈ (0, 1/4). Let µ be the distribution N (0, σ2I). For any estimator Â ∈ F , there exists an matrix A∗ ∈ Rd×d where A∗ = ρO for some orthogonal matrix O such that |σmax(A∗)| = ρ and we have that with probability at least δ:
∥Â−A∗∥ = Ω √
(d+ log(1/δ))(1− ρ) T . (8)
Notice that in the setting of Theorem 3, we have G = ∑∞
i=0 σ 2(A∗)i(A∗)i,⊤ = σ
2
1−ρ2 I . Therefore, σmin(G) = 1 1−ρ2 ∼ 1 1−ρ . The bound in Theorem 1 matches the above minimax bound up to logarithmic factors.
Next we consider the prediction loss. We fix dimension d and horizon T and consider the class of VAR modelsM such that Assumptions 1, 2, and 3 hold such that Tr(Σ(µ)) = β ∈ R+ be fixed. Let
F be the class of all estimators for parameter A∗ given data (Z0, . . . , ZT ). We want to lower bound the minimax error:
Lminmax(M) := inf f∈F sup (A∗,µ)∈M E(Zt)∼VAR(A∗,µ)Lpred(f(Z0, . . . , ZT );A ∗, µ)− Lpred(A∗;A∗, µ).
Theorem 4. For some universal constant c, we have: Lminmax(M) ≥ cβ(d− 1)min ( 1
T , 1 d2
) , where β = Tr(Σ(µ)).
Note that the theorem shows that our algorithm is minimax optimal with respect to the prediction loss at stationarity, Lpred. See Section M for a detailed proof of the above lower bound.

heading: 5 Idea Behind Proofs
text: In this section, we provide an overview of the key techniques to prove our results. As observed in the discussion following Equation (5), when the data is processed in the reverse order within a buffer, it behaves similar to SGD for linear regression with i.i.d. data. Due to the gaps of size u, we can take the buffers to be approximately independent. Therefore, we analyze the algorithm as follows:
1. Analyze reverse order within a buffer using the property noted in Equation (5). 2. Treat different buffers to be i.i.d. due to gap and present an i.i.d data type analysis.
To execute the proposed proof strategy, we introduce the following technical notions:
Coupled Process. For the real data points (Xτ ), the points in different buffers are weakly dependent. In order to make the analysis straight forward, we introduce the fictitious coupled process X̃τ such that ∥∥∥X̃τ −Xτ∥∥∥ ≲ 1Tα for large enough α, for every data point Xτ used by SGD− RER. We have the additional property that the successive buffers are actually independent for this coupled process. We refer to Definition 1 in the appendix for the construction of the coupled process X̃τ .
Suppose we run SGD−RER with the coupled process X̃τ instead of Xτ to obtain the coupled iterates Ãti. We can then show that Ã t i ≈ Ati. Thus it suffices analyze the coupled iterates Ãti. We refer to Sections B and C for the details.
Bias Variance Decomposition. We consider the standard bias variance decomposition with individual buffers as the basic unit as opposed to individual data points. We refer to Section D for the details. We decompose the error in the iterates into the bias part ( Ãt−1,bB −A∗ ) = (A0 − A∗) ∏t−1 s=0 H̃
s 0,B−1 and the variance part(
Ãt−1,vB
) = 2γ ∑t r=1 ∑B−1 j=0 η t−r −j X̃ t−r,⊤ −j H̃ t−r j+1,B−1 ∏1 s=r−1 H̃ t−s 0,B−1 where the matrices
H̃s0,B−1 = ∏B−1 i=0 ( I − 2γX̃s−iX̃ s,⊤ −i ) are the independent ’contraction’ matrices associated with each buffer s. This result in the geometric decay of the initial distance between (A0 − A∗). The variance part is due to the inherent noise present in the data. In Section F.1 we first establish the exponential decay of the ‘bias’. We then consider the second moment of the variance term. Observe that the distinct terms in the expression for ( Ãt−1,vB ) are uncorrelated either due to reverse order within a buffer as noted in Equation (5) or due to independence between the data in distinct buffers (due to coupling). This allows us to split the second moment into diagonal terms with non-zero mean and cross terms with zero mean. Diagonal terms are analyzed via a recursive argument in Claim 1 and the following discussion in order to remove dependence on mixing time factors. The analysis for parameter recovery (the result of Theorem 2) is similar but we bound the relevant exponential moments using sub-Gaussianity of the noise sequence ηt to obtain high-probability bounds which when combined with standard ϵ-net arguments give us guarantees for the operator norm error Lop.
Averaged Iterates. We then combine the bias and variance bounds obtained for individual iterates in Section F.1 to analyze the tail averaged output. Using techniques standard in the analysis of SGD for linear regression, we finally show that this averaging leads error rates of the order d 2
T . We refer to Sections E (for parameter recover) and G (for prediction error) for the detailed results.
Picking the Step Sizes and Conditioning. Due to the auto-regressive nature of the data generation, the iterates can grow to be of the size O( d1−ρ ). The step sizes need to be set small enough so that the γ∥XτX⊤τ ∥ ≤ 1 in order for the SGD− RER iterations to not diverge to infinity. In the statement of Theorem 2, we condition on the event where ∥Xτ∥2 are all bounded by a sufficiently large number R for every τ in order to ensure this property. The relevant events where the norm is bounded are defined in Section B. Conditioning on these events results in previously zero mean terms to be not zero mean. Routine calculations using triangle inequality and Cauchy-Schwarz inequality ensure that the means are still of the order 1Tα for any fixed constant α > 0. Furthermore, we actually require step sizes such that γ ∥∥∑ τ∈Buffer XτX ⊤ τ
∥∥ ≤ 1 to show exponential contraction of H̃s0,B−1 matrices due to the Grammian G as described next.
Probabilistic Results. We establish some properties of H̃s0,B−1, which are products of dependent random matrices in Section L. Specifically we refer to Lemmas 28, 29, 30, and 31 which establish that ∥∥∥∏t−1s=0 H̃s0,B−1∥∥∥ ≲ (1− γBσmin(G))t with high probability.

heading: 6 Experiments
text: In this section, we compare performance of our SGD− RER method on synthetic data against the performance of standard baselines OLS and SGD, along with SGD−ER method that applies standard experience replay technique, but where points from a buffer are sampled randomly. Synthetic data: We sample data from VAR(A∗, µ) with X0 = 0, µ ∼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 ⌈d/2⌉ entries on diagonal being ρ and the remaining diagonal entries are set to ρ/3. We set d = 5, ρ = 0.9 and σ2 = 1. We fix a horizon T = 107 and set the buffer size as B = 100 and u = 10. To estimate R from the data, we use the first ⌊2 log T ⌋ = 32 samples and set R as the sum of the norms of these samples. We let the stepsize to be γ = 12R which is aggressive compared to our theorems. We start the SGD−RER and other SGD-like algorithms from the second buffer onward.
For tail averaging, as described in algorithm 1, we ignore the first ⌊log T ⌋ = 16 buffers, and maintain a running tail average at the end of each of the subsequent buffers. In figure 2, we plot the parameter errors ∥∥∥Âlog T,t −A∗∥∥∥ and ∥∥∥Â0,t −A∗∥∥∥ versus the buffer index t as the algorithm runs for horizon T . For OLS, we include samples in the first buffer as well (which were used for estimating R). Clearly, SGD−RER has very similar performance as that of OLS whereas SGD−ER and SGD seem to display residual bias for the chosen step-size (which is logarithmic in the horizon T ) and buffer lengths. We also observe a similar behavior when we choose A∗ = ρI .

heading: 7 Conclusion
text: In this paper, we studied the problem of linear system identification in streaming setting and provided an efficient algorithm (SGD− RER). We proved that SGD− RER achieves nearly minimax optimal error rate, both in terms of parameter error as well as prediction error. Furthermore, using experiments, we validated that standard SGD as well as SGD with experience replay can have large bias error. Our algorithm and analysis demonstrates that the knowledge of dependency structure can aid us in designing accurate algorithms for dependent data.
This work opens up a myriad of open questions about learning from dependent data in general and Markov processes in particular. Our work currently assumes a specific Markovian dependency structure – extending the intuition and techniques to handle more general data dependencies is an interesting open question. Further, our work does not address the question of recovering a sparse system matrix with unknown sparsity pattern. So online learning of such linear dynamical systems with (unknown) sparsity pattern or low-rank structure is an exciting question with applications to domains like bioinformatics. Moreover, even in our linear setting, extending SGD − RER to the situation of partially observed states with or without control inputs would be another direction to pursue. Finally, it would be interesting to understand how the techniques introduced in this work perform in practical RL settings where learning with data from Markov processes is essential.
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. Part of this work was done when S.S.K was visiting Microsoft Research Lab India Pvt Ltd during summer 2020.

Review summaries:
This paper considers systems identification in vectorized auto-regressive model with stochastic noise. The paper focuses on the practicality of system identification with nearly optimal computation and sample complexity instead of using offline-style methods such as OLS. The authors propose a new SI algorithm based on SGD that uses reverse experience replay to address the challenges due to the inherent correlation between samples. Performance bounds on the accuracy of estimating the transition matrix as well as prediction error are given, both of which nearly matching the lower bounds. Numerical simulations show that SGD-RER is competitive to OLS, which enjoying a smaller time complexity. 

The paper proposes an SGD-style algorithm for linear dynamical system identification. The authors claim that the proposed algorithm is the first of its kind to have near-optimal error rates. The key idea behind the algorithm is to run SGD in reverse on streaming segments of the data. 

The authors consider a rather special case of a linear dynamical system, where there is no hidden state and the transition matrix is "externally stable" (less than 1 in the operator norm), process noise is sub-gaussian, and the stationary distribution has bounded fourth moment. In this case, they suggest running the SGD "backwards" in small buffers, which they call "reverse experience replay". This leads to the independence of the observations, and the unbiased nature of the SGD updates.