Gaussian Approximation and Multiplier Bootstrap for Polyak-Ruppert Averaged Linear Stochastic Approximation with Applications to TD Learning

Authors: sergey  samsonov
eric  moulines
ecole  polytechnique
qi-man  shao
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 10:35:15

Abstract: In this paper, we obtain the Berry–Esseen bound for multivariate normal approximation for the Polyak-Ruppert averaged iterates of the linear stochastic approximation (LSA) algorithm with decreasing step size. Moreover, we prove the non-asymptotic validity of the confidence intervals for parameter estimation with LSA based on multiplier bootstrap. This procedure updates the LSA estimate together with a set of randomly perturbed LSA estimates upon the arrival of subsequent observations. We illustrate our findings in the setting of temporal difference learning with linear function approximation.

Main paper: 
 introduction: 
 Stochastic approximation (SA) methods are a central component for solving various optimization problems that arise in machine learning [32, 26] , empirical risk minimization [72] and reinforcement learning [42, 67] . There is a vast number of contributions in the literature, which cover both asymptotic [48, 53] and non-asymptotic [45, 15, 36] properties of the SA estimates. The primarily important property among the asymptotic ones of the SA estimates is their asymptotic normality [53] , which is important due to its role in constructing (asymptotic) confidence intervals and hypothesis testing [71] . However, a natural question of the rate of convergence in the appropriate central limit theorems (CLT) is not well addressed in literature even in the relatively simple setting of the linear stochastic approximation (LSA) [21] , [34] , [9] .Alternatively, confidence sets for SA algorithms can be constructed in a non-asymptotic manner based on concentration inequalities [4] . These bounds are often regarded as loose [60] , yielding suboptimal performance of the statistical procedures based on the latter estimates [28] . In contrast, for statistical inference procedures based on independent and identically distributed (i.i.d.) observations, such as M -estimators [71] , there is a machinery of non-parametric methods for constructing confidence sets with the bootstrap [19, 58] . This approach is accompanied with theoretical guarantees, showing the non-asymptotic validity of the bootstrap-based confidence intervals for parameters in linear regression [64] and statistical tests [12] . Extending theoretical guarantees to a non-classical situation with online learning algorithms encounters serious problems, essentially related to the problem of obtaining rate of convergence in the corresponding CLTs. At the same time, many phenomena arising in the analysis of nonlinear SA algorithms already appear in the analysis of LSA problems.The LSA procedure aims to find an approximate solution for the linear system Āθ = b with a unique solution θ based on a sequence of observations {(A(Z k ), b(Z k ))} k∈N . Here A : Z → R d×d and b : Z → R d are measurable functions and (Z k ) k∈N is a sequence of noise variables taking values in some measurable space (Z, Z) with a distribution π satisfying E[A(Z k )] = Ā and E[b(Z k )] = b. We focus on the setting of independent and identically distributed (i.i.d.) observations {Z k } k∈N . With a sequence of decreasing step sizes (α k ) k∈N and the starting point θ 0 ∈ R d , we consider the estimates { θn } n∈N given byEQUATIONHere, we have fixed the size of the burn-in period (see, e.g., [16, 44] ) to n 0 = n. Provided that n is large enough, the burn-in size affects only a constant factor in the subsequent bounds. The sequence {θ k } k∈N corresponds to the standard LSA iterates, while { θn } n∈N corresponds to the Polyak-Ruppert (PR) averaged iterates [59, 53] . It is known that θn is asymptotically normal with a minimax-optimal covariance matrix (see [53] and [23] for discussion). Specifically, under appropriate technical conditions on the step sizes {α k } and noisy observations {A(Z k )}, it holds that √ n( θn -θ )d → N (0, Σ ∞ ) ,where Σ ∞ is the asymptotic covariance matrix defined later in Section 3.1. There is a long list of contributions to the non-asymptotic analysis of θn , particularly [43] and [16] , which study moment and Bernstein-type concentration bounds for √ n( θn -θ ). Unfortunately, such bounds do not imply Berry-Esseen type inequalities for √ n( θn -θ ), that is, they do not allow us to control the quantityEQUATIONwhere Conv(R d ) refers to the set of convex sets in R d . While the Berry-Esseen bounds are a popular subject of study in probability theory, starting from the classical work [20] , most results are obtained for sums of random variables or martingale difference sequences [52, 8] . We can only mention a few results for SA algorithms, see Section 2 for more details. This paper aims to provide the latter bounds for the specific setting of the LSA procedure. Our primary contribution is twofold:• We establish a BerryEsseen bound for accuracy of normal approximation of the distribution of Polyak-Ruppert averaged LSA iterates with a polynomially decreasing step size. Our results suggest that the best rate of normal approximation, in the sense of (2) , is of order n -1/4 up to logarithmic factors in n, where n denotes the number of samples. Interestingly, this rate is achieved with an aggressive step size, α k = c 0 / √ k. Our proof technique follows the Berry-Esseen bounds for nonlinear statistics provided in [63] .• We provide non-asymptotic confidence bounds for the distribution of the PR-averaged statistic √ n( θn -θ ) using the multiplier bootstrap procedure. In particular, our bounds imply that the quantiles of the exact distribution of √ n( θn -θ ) can be approximated at a rate of n -1/4 , where n is the number of samples used in the procedure, provided that n is sufficiently large (see A4 for exact conditions). To the best of our knowledge, this is the first non-asymptotic bound on the accuracy of bootstrap approximation in SA algorithms. We apply the proposed methodology to the temporal difference learning (TD) algorithm for policy evaluation in reinforcement learning.The rest of the paper is organized as follows. In Section 2, we provide a literature review on the non-asymptotic analysis of the LSA algorithm and bootstrap methods. Next, in Section 3, we analyze the convergence rate of Polyak-Ruppert averaged LSA iterates to the normal distribution. In Section 4, we discuss the multiplier bootstrap approach for LSA and establish bounds on the accuracy of approximating the quantiles of the true distribution. Finally, we apply our findings to TD learning and present numerical illustrations in Section 5.Notations. For matrix A ∈ R d×d we denote by A its operator norm. For symmetric matrix Q = Q 0 , Q ∈ R d×d and x ∈ R d we define the corresponding norm x Q = x Qx, and define the respective matrix Q-norm of the matrix B ∈ R d×d by B Q = sup x =0 Bx Q / x Q . For sequences a n and b n , we write a n b n if there exist a constant c > 0 such that a n ≤ cb n for c > 0. For simplicity we state the main results of the paper up to constant factors. 
related works: 
 Among contributions to the analysis of the LSA algorithm, we should mention the papers [53, 34, 9, 6] . These works investigate the asymptotic properties of the LSA estimates (such as asymptotic normality and almost sure convergence) under i.i.d. and Markov noise. Non-asymptotic results for the LSA and PR-averaged LSA estimates were obtained in [55, 47, 7, 35, 44] , where MSE bounds were established, and in [43, 17, 16] , which provided high-probability error bounds. The latter results enable the construction of Bernstein-type confidence intervals for the error θn -θ . Unfortunately, the corresponding bounds typically depend on unknown problem properties of (1), related to the design matrix Ā and the noise variables A(Z k ), b(Z k ). For this reason, applying these error bounds in practice is complicated. Furthermore, concentration bounds for the LSA error [43, 17, 16] do not imply convergence rates of the rescaled error √ n( θn -θ ) to the normal distribution in Wasserstein or Kolmogorov distance. Non-asymptotic convergence rates were previously studied in [2] using the Stein method, but the resulting rate corresponds to a smoothed Wasserstein distance. Recent work [65] investigates convergence rates to the normal distribution in Wasserstein distance for LSA with Markovian observations. Both papers yield bounds that are less tight with respect to their dependence on trajectory length n than those presented in the present work, see a detailed comparison after Theorem 2.A popular method for constructing confidence intervals in the context of parametric estimation is based on the bootstrap approach ( [19] ). Its analysis has attracted many contributions, in particular a series of papers [12] and [13] that validate a bootstrap procedure for a test based on the maximum of a large number of statistics. Their study shows a close relationship between bootstrap validity results, Gaussian comparison and anticoncentration bounds for rectangular sets. The papers [64] and [27] investigate the applicability of likelihood-based statistics for finite samples and large parameter dimensions under possible model misspecification. The important step in proving bootstrap validity is again based on Gaussian comparison and anticoncentration bounds, but now for spherical sets. The bootstrap procedure for spectral projectors of covariance matrices is discussed in [46] and [31] . The authors follow the same steps to prove the validity of the bootstrap.Extending the classical bootstrap approach to online learning algorithms is a challenge. For example, the iterates {θ k } k∈N determined by (1) are not necessarily stored in memory, which makes the classical bootstrap inapplicable. This problem can be solved by performing randomly perturbed updates of the online procedure, as proposed in [22] for the iterates of the Stochastic Gradient Descent (SGD) algorithm. The authors in [56] used the same procedure for the case of Markov noise and policy evaluation algorithms in reinforcement learning, but in both papers the authors only consider the asymptotic validity. In our paper we use the same multiplier bootstrap approach (see Section 4), but we provide an explicit error bound for the bootstrap approximation of the distribution of the statistics √ n( θn -θ ).In addition to the bootstrap approach, one can also use the pivotal statistics [37, 40, 41] or various estimates of the asymptotic covariance matrix [73] to construct the confidence intervals for θ . The latter approach can be based on the plug-in estimators [39] , batch mean estimators [11] or in combination with the multiplier bootstrap approach [74] . However, the theoretical guarantees for mentioned methods remain purely asymptotic. 
accuracy of normal approximation for lsa: 
 We first study the rate of normal approximation for the tail-averaged LSA procedure. When there is no risk of ambiguity, we shall use simply the notationsA n = A(Z n ) and b n = b(Z n ).Starting from the definition (1), we get with elementary transformations thatEQUATIONwhere we have setε n = ε(Z n ) with ε(z) = Ã(z)θ -b(z) , Ã(z) = A(z) -Ā , b(z) = b(z) -b .Here the random variable ε(Z n ) can be viewed as a noise, measured at the optimal point θ . We now assume the following technical conditions: A 1. Sequence {Z n } n∈N is a sequence of i.i.d. random variables defined on a probability space (Ω, F, P) with distribution π.A 2. Z A(z)dπ(z) = Ā and Z b(z)dπ(z) = b, with the matrix -Ā being Hurwitz. Moreover, ε ∞ = sup z∈Z ε(z) < +∞, and the mapping z → A(z) is bounded, that is,EQUATIONMoreover, for the noise covariance matrixEQUATIONit holds that its smallest eigenvalue is bounded away from 0, that is,EQUATIONIt is possible to change (4) to the moment-type bound as it was previously considered in [43] and [16] , see the detailed discussion after Theorem 2. The fact that the matrix -Ā is Hurwitz implies that the linear system Āθ = b has a unique solution θ . Moreover, this fact is sufficient to show that the matrix I -α Ā is a contraction in an appropriate matrix Q-norm for small enough α > 0. Precisely, the following result holds: Proposition 1. Let -Ā be a Hurwitz matrix. Then for any P = P I, there exists a unique matrixQ = Q I, satisfying the Lyapunov equation Ā Q + Q Ā = P . Moreover, setting a = λmin(P ) 2 Q, andEQUATIONwhereEQUATIONEQUATIONThe proof of Proposition 1 is provided in Appendix D.1. Note that it is possible to set P = I as in [18] , yet we will observe that other choices of P could be more beneficial. Now consider an assumption on the step sizes α k and number of observations n: A3. The step sizes {α k } k∈N has a form α k = c 0 /k γ , where γ ∈ [1/2; 1) and c 0 ∈ (0; α ∞ ]. Moreover, we assume that n ≥ d, andEQUATIONThe main aim of lower bounding n is to ensure that the number of observations is large enough in order that the LSA error related to the choice of initial condition θ 0 -θ becomes small.3.1 Central limit theorem for Polyak-Ruppert averaged LSA iterates.It is known that the assumptions A1-A3 guarantee that the CLT applies to the iterates of θn , namely,EQUATIONwhere the asymptotic covariance matrix Σ ∞ has a formEQUATIONand Σ ε is defined in (5) . This result can be found for example in [53] and [23] . We are interested in the Berry-Esseen type bound for the rate of convergence in (10) , that is, we aim to bound ρ Conv n defined in (2) w.r.t. the available sample size n. We control ρ Conv n using a method from [63] based on randomized multivariate concentration inequality. Below we briefly state its setting and required definitions. Let X 1 , . . . , X n be independent random variables taking values on X and T = T (X 1 , . . . , X n ) be a general d-dimensional statistics such that T = W + D, whereEQUATIONξ = h (X ) and h : X → R d is a Borel measurable function. Here the statistics D can be non-linear and is treated as an error term, which is "small" compared to W in an appropriate sense. Assume thatEQUATIONwhere D ( ) = D(X 1 , . . . , X -1 , X , X +1 , . . . , X n ) and X is an independent copy of X . This result is due to [63, Theorem 2.1] . One can modify the bound (13) for the setting whenn =1 E[ξ ξ ] = Σ 0.This result due to [63, Corollary 2.3] . Following the construction (12), we set T = √ n Ā( θn -θ ) and consider it as a nonlinear statistic of i.i.d. random variables Z 1 , . . . , Z 2n , which drive the LSA dynamics (1) . We can exactly represent T as a sum of linear (W ) and non-linear parts (D), whereW = - 1 √ n 2n-1 k=n ε k+1 , D = 1 √ n θ n -θ α n - 1 √ n θ 2n -θ α 2n - 1 √ n 2n k=n+1 (A k -Ā)(θ k-1 -θ ) + 1 √ n 2n k=n+1 θ k-1 -θ 1 α k - 1 α k-1 .The proof of this result can be bound in Proposition 3. To obtain a bound for the approximation accuracy in (2) using the bound (13), we need to upper boundE 1/2 [ D(Z 1 , . . . , Z 2n ) 2 ] and E[ D - D (i) ].The first result below provides a second moment bound on D: Theorem 1. Assume A1, A2, and A3. Then we obtain the following error bound:E 1/2 D(Z 1 , . . . , Z n ) 2 √ κ Q ε ∞ √ a 1 n (1-γ)/2 + C A n γ/2 + √ κ Q C A n (1-γ)/2 √ a + √ κ Q n 2γ-1 a exp - c 0 an 1-γ 2(1 -γ) θ 0 -θ ,where stands for inequality up to an absolute constant.The proof of Theorem 1 is provided in Appendix A.3. Now it remains to upper bound the term E[ D -D (i) ], which is done in Appendix B.1 using the synchronous coupling methods [10] .Combining these bounds, we obtain the following theorem: Theorem 2. Assume A1, A2, and A3. Then the following bound holds:EQUATIONwhere∆ 1 = ∆ 1 (n, a, C A , Tr Σ ε ) is a polynomial function defined in (31), and constants C 1 , C 2 , depending upon a, C A , κ Q , Tr Σ ε , are defined in (32).The proof of Theorem 2 is provided in Appendix B. Note that the assumption A 2 requires that ε(Z 1 ) is almost sure bounded. It is a strong assumption, but can be partially relaxed. Following the stability of matrix products technique, used in [17, Proposition 3] , it is possible to consider the setting when the random variable Ã(Z 1 ) has only finite number of moments. In particular, a finite third moment of Ã(Z 1 ) implies that ε(Z 1 ) also has a finite third moment, which is sufficient to obtain a counterpart to Theorem 1. However, this generalization requires non-trivial technical work on generalizing the stability of matrix products result (see Corollary 4 in Appendix D ).Note that the bound of Theorem 2 predicts the optimal error of normal approximation for Polyak-Ruppert averaged estimates of order n -1/4 , which is achieved with the aggressive step size α k = c 0 / √ k, that is, when setting γ = 1/2 in (14) . In this case we obtain the optimized boundρ Conv n C 3 λ min n 1/4 + d 1/2 ε 3 ∞ λ 3/2 min √ n + ∆ 1 exp -c 0 a √ n λ min θ 0 -θ ,where (32) .C 3 = C 3 (a, C A , κ Q , Tr Σ ε , ε ∞ ) is provided inDiscussion. Our proof technique of Theorem 2 reveals an interesting feature: fastest rate of convergence in the convex distance ρ Conv n corresponds to the learning rate schedule that admits the fastest decay of the second-order term in the MSE bound for remainder statistics D (see Theorem 1). Results similar to the one of Theorem 2 have been recently obtained in the literature in [65] and [2] . The author in [65] considers the LSA problem specified to the temporal-difference learning (see Section 5) with Markov noise and obtains convergence rate in Wasserstein distance of order n -1/4 , which corresponds to the "optimal" step size schedule α k = c 0 /k 3/4 . Using the bound of [49, eq. ( 3)] (see also section 2 in [57] ), this result yield a suboptimal bound of order n -1/8 for the convex distance ρ Conv n . Such an upper bound may be loose for some classes of distributions, but it is not clear if in particular setting of LSA the bound of [65] could imply scaling of order n -1/4 for ρ Conv n . At the same time, in case of X 1 , . . . , X n forming a Markov chain in (12) there is no available counterpart of the bound (13) . Generalizing ( 13) is an interesting research direction that would allow to obtain a counterpart of Theorem 2 in case of Markovian dynamics. Similarly, the result of [2] holds for much stronger metrics, which controls the convergence of moments of twice differentiable functions. We provide additional details about connections between this metric and ρ Conv n in Appendix B.2. At the same time, the authors in [2] cover the non-linear setting of PR-averaged iterates of stochastic gradient descent algorithm under strong convexity. Remark 1. The leading (with respect to n) terms of the bound from Theorem 1 have an implicit dependence on the problem dimension d due to the presence of λ min . Yet the result of Theorem 1 can be improved in a sense of dependence in dimension if one is interested not in the rates of convergence for √ n( θn -θ ), but in the projected iterated √ nΠ ( θn -θ ) for some Π ∈ R d×m , m ≤ d. If this is the case, one may apply (13) for the class Conv m = Conv(R m ) of convex sets in R m and obtain, setting step sizeEQUATIONΣ (Π) ε = ΠΣ ε Π , that ρ Conv n C 4 λ min n 1/4 + m 1/2 ε 3 ∞ λ 3/2 min √ n + ∆ 1 e -c0a √ n λ min θ 0 -θ ,and the constantC 4 = C 4 (a, C A , κ Q , Tr Σ (Π) ε , ε ∞ ) is provided in (32). Remark 2.Results similar to Theorem 1 can be obtained not only for the Polyak-Ruppert averaged estimator θn , but also for the last iterate θ n . In particular, it is known (see e.g. [23] ), that the last iterate error θ n -θ is also asymptotically normal:θn-θ * √ αn → N (0, Σ last ) ,where the covariance matrix Σ last is different from Σ ∞ . In such a case Σ last can be found as a solution to appropriate Lyapunov equation, see [23] . Then, we can use the perturbation-expansion technique from [1] , and, applying similar technique of randomized concentration inequalities [63] (see (13) ), we can obtain the Berry-Esseen boundsup B∈Conv(R d ) P θn-θ * √ αn ∈ B -P(Σ 1/2 last η ∈ B) √ α n . 
multiplier bootstrap for lsa: 
 In order to perform statistical inference with the Polyak-Ruppert estimator θn , we propose an online bootstrap resampling procedure, which recursively updates the LSA estimate as well as a large number of randomly perturbed LSA estimates, upon the arrival of each data point. The suggested procedure follows the one outlined in [22] . It has the following advantages: it does not rely on the asymptotic distribution of the error √ n( θn -θ ), does not require to know the moments of √ n( θn -θ ) or its asymptotic covariance matrix Σ ∞ , and does not involve any data splitting.We state the suggested procedure as follows. Let W 2n = {W } 1≤ ≤2n be a set of i.i.d. random variables, independent of Z 2n = {Z } 1≤ ≤2n , with E[W 1 ] = 1 and Var[W 1 ] = 1. We write, respectively, P b = P(•|Z 2n ) and E b = E(•|Z 2n ) for the corresponding conditional probability and expectation. In parallel with procedure (1) that generates {θ k } 1≤k≤2n and θn , we generate M independent samples (w n , . . . , w 2n ), 1 ≤ ≤ M distributed as W 2n , and recursively update M randomly perturbed LSA estimates, that is,EQUATIONWe use a short notation θb n for θb,1 n . The key idea of the procedure (15) is that the "Bootstrap-world" distribution (that is, the one conditional on Z 2n ) of the perturbed samples √ n( θb n -θn ) is close to the distribution of the quantity of interest, that is, √ n( θn -θ ). Precisely, the main result of this section will show that the quantity supEQUATIONis small. Although an analytic expression for P b ( √ n( θb n -θn ) ∈ B) is not available, one can approximate it from numerical simulations according to (15) by generating sufficiently large number M of perturbed trajectories. Standard arguments, see e.g. [62, Section 5.1] suggest that the accuracy of Monte-Carlo approximation is of order M -1/2 . To analyze the suggested procedure, we shall impose an additional assumption on the trajectory length n: A4. Assumption A3 holds with γ = 1/2, and c 0 ≤ 1/(C 2 A κ Q e). Moreover, settingEQUATION)it holds thatEQUATIONMoreover, we assume that for λ min defined in (6) it holds thatEQUATIONNote that the new bound (18) simply states that √ n/ log 2 (n) is sufficiently large, since h(n) scales as log 2 n. We discuss the assumption A4 in more details in the proof scheme. Now we formulate the main result of this section. We analyze only the setting of polynomially decaying step size with γ = 1/2, since decay rate of (16) essentially depends on the approximation rate of Theorem 2, with the fastest rate achieved when γ = 1/2. For other learning rates the decay rate of right-hand side in Theorem 3 will be slower. Theorem 3. Assume A1, A2, A3 with γ = 1/2, and A4. Then with P -probability at least 1 -6/n it holds that supB∈Conv(R d ) |P b ( √ n( θb n -θn ) ∈ B) -P( √ n( θn -θ ) ∈ B)| κ 2 Q (C 4 A ∨1)(1 + ε 2 ∞ ) log n a 5/2 λ min n 1/4 + √ d √ n ε 3 ∞ λ 3/2 min + κ Q ε ∞ √ log n √ λ min + κ Q (1 + ε 2 ∞ /λ min ) log n √ n + ∆ 2 e -(c0/2)a √ n λ min θ 0 -θ ,where (42) .∆ 2 = ∆ 2 (n, a, C A , ε ∞ ) is a polynomial function defined inThe proof of Theorem 3 is based on the Gaussian approximation performed both in the "real" world and bootstrap world together with an appropriate Gaussian comparison inequality. The main steps of the proof are illustrated by the following scheme:Real world:√ n Ā( θn -θ ) ξ ∼ N (0, Σ ε ) Bootstrap world: √ n Ā( θb n -θn ) ξ b ∼ N (0, Σ b ε )Gaussian approximation, Th. 2Gaussian comparison, Theorem 5Gaussian approx. in Bootstrap world, Th. 4In the above scheme we have denoted by Σ b ε = n -1 2n-1 =n ε ε the sample covariance matrix approximating Σ ε . Gaussian approximation for the true distribution of √ n Ā( θn -θ ) follows from Theorem 2. Proof of Gaussian approximation in the Bootstrap world Theorem 4 is also based on the inequality (13) , but is more complicated and involves the expansion analysis of the LSA error from [1] . This technique allows to separate the LSA error into different scales with respect to the step sizes {α k }, see Appendix C.4 for details. However, this technique requires to impose additional assumption A4 -eq. (18) . Proof of the Gaussian comparison part of Theorem 5 is based on Pinsker's inequality and matrix Bernstein inequality. The latter result requires that n is large enough to ensure that minimal eigenvalue of Σ b ε is close to λ min , justifying the assumption A4 -eq. (19) . Detailed proof if provided in Appendix C.Discussion. We emphasize that the Gaussian approximation result of Theorem 2 (with Bootstrap world generalization in Theorem 4) is a key result to prove the above bootstrap validity. This argument was missing in the earlier works studying confidence intervals for stochastic optimization algorithms [11, 73, 74] , where the authors considered procedures to estimate Σ ∞ in (11) . They combine nonasymptotic bounds on the accuracy of recovering Σ ∞ with only asymptotic validity of the resulting confidence intervals. We expect that our proof technique for Theorem 2 can be used to provide similar non-asymptotic validity results for outlined approaches for constructing confidence intervals based on the estimation of the asymptotic covariance matrix. Corollary 1. (Set of Euclidean balls or ellipsoids) Suppose that we are interested in estimating quantile of a given order α ∈ (0, 1) and some matrix B ∈ R d×d , that is, the quantityt α = inf{t > 0 : P( √ n B( θn -θ ) ≥ t) ≤ α}.We define its counterpart in the Bootstrap world, t b α , ast b α = inf{t > 0 : P b ( √ n B( θb n -θn ) ≥ t) ≤ α}.Note that t b α is defines with respect to the bootstrap measure, therefore, it depends on the data Z 2n . This bootstrap critical value t b α is applied in the Bootstrap world to build the confidence setE(α) = {θ ∈ R d : √ n B(θ -θn ) ≤ t b α }. Theorem 3 justifies this construction and evaluate the coverage probability of the true value θ by this set. It states thatP(θ / ∈ E(α)) = P( √ n B( θn -θ ) > t b α ) ≈ α ,with the error of order n -1/4 in the right-hand side. Although an analytic expression for t b α is not available, one can approximate it by generating a large number M of independent samples of W n and computing from them the empirical distribution function of √ n B( θb n -θn ) , following (15). Remark 3. A natural question that arises after Theorem 3 is whether it is possible to prove similar bounds for the iterates of first-order stochastic optimization algorithms. There are several MSE bounds for corresponding algorithms with explicit dependence on the step size α k ; see, for example, [45, 5] . Therefore, we expect that it is possible to obtain a counterpart to Theorem 2. At the same time, for general first-order stochastic optimization algorithms, unlike LSA, there are no counterparts to the precise error expansions of [1] . Thus, proving the counterpart of Theorem 3 in this setting is more challenging. Similarly, we emphasize that generalizations of the procedure (15) to cases where {Z k } k∈N are dependent, for example, form a Markov chain, are complicated. The approach of [63] is not directly applicable in this setting, and appropriate generalization of (13) is a separate and challenging research direction.Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)? Answer: [Yes] Justification: Numerical results are stated with a complete description of the environments that are used, as well as the precise sets of hyperparameters that we used. The code (in Python) is provided as supplementary with the paper, making it easy for one to reproduce our numerical experiments. At the same time, tracing the second-order terms in the normal approximation is computationally involved and can take sufficiently large amount of time.Guidelines:• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways.For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
applications to the td learning and numerical results: 
 We illustrate our findings for the setting of temporal difference (TD) learning algorithm [66, 67] for policy evaluation in RL. Non-asymptotic error bounds for this algorithm attracted lot of contributions [43, 16, 30, 51, 38] . At the same time, confidence intervals for TD were studied in [22, 56] only in terms of their asymptotic validity. In the TD algorithm we consider a discounted MDP (Markov Decision Process) given by a tuple (S, A, P, r, γ). Here S and A stand for state and action spaces, and γ ∈ (0, 1) is a discount factor. Assume that S is a complete metric space with metric d S and Borel σ-algebra B(S). P stands for the transition kernel P(B|s, a), which determines the probability of moving from state s to a set B ∈ B(S) when action a is performed. Reward function r : S × A → [0, 1] is assumed to be deterministic. Policy π(•|s) is the distribution over action space A corresponding to agent's action preferences in state s ∈ S. We aim to estimate value functionV π (s) = E ∞ k=0 γ k r(s k , a k )|s 0 = s ,where a k ∼ π(•|s k ), and s k+1 ∼ P(•|s k , a k ) for any k ∈ N. Define the transition kernel under π, P π (B|s) = A P(B|s, a)π(da|s) , (20) which corresponds to the 1-step transition probability from state s to a set B ∈ B(S). The state space S here can be arbitrary. It is a common option to consider the linear function approximation for V π (s), defined for s ∈ S, θ ∈ R d , and a feature mapping ϕ :S → R d as V π θ (s) = ϕ (s)θ.Here d is the dimension of feature space. Our goal is to find a parameter θ which is defined as a unique solution to the projected Bellman equation, see [70] , which defines the best linear approximation of V π . We denote by µ the invariant distribution over the state space S induced by P π (•|s) in (20) . We define the design matrix Σ ϕ asEQUATION)Consider the following assumptions on the generative mechanism and on the feature mapping ϕ(•):TD 1. Tuples (s, a, s ) are generated i.i.d.with s ∼ µ, a ∼ π(•|s), s ∼ P(•|s, a) . TD 2. Matrix Σ ϕ is non-degenerate with the minimal eigenvalue λ min (Σ ϕ ) > 0. Moreover, the feature mapping ϕ(•) satisfies sup s∈S ϕ(s) ≤ 1.In the setting of linear function approximation the estimation of V π (s) reduces to estimating θ ∈ R d , which can be done via the LSA procedure. Here, the k-th step randomness is given by the tuple Z k = (s k , a k , s k ). Then, the corresponding LSA update can be written asEQUATIONwhere A k and b k are given, respectively, byA k = ϕ(s k ){ϕ(s k ) -γϕ(s k )} , b k = ϕ(s k )r(s k , a k ) .We provide the expressions for the corresponding system matrix Ā = E[A k ] and the right-hand side b in Appendix E. We verify that assumption A2 holds and, furthermore, we provide a tighter counterpart to the result of Proposition 1. This result closely follows [51] and [61] . Proposition 2. Let {θ} k∈N be a sequence of TD updates generated by (22) under TD 1 and TD 2.Then this update scheme satisfies assumption A2 withC A = 2(1 + γ) , ε ∞ = 2(1 + γ)( θ + 1) , moreover, one can check that I -α Ā 2 ≤ 1 -αa with a = (1 -γ)λ min (Σ ϕ ) , α ∞ = (1 -γ)/(1 + γ) 2 , that is, Proposition 1 holds with Q = I.Proof of Proposition 2 is provided in Appendix E. Since all the assumptions in A2 are fulfilled, we can verify tightness of the bound Theorem 2 for different learning rate schedules α k in (22) .Numerical results. Efficiency of the multiplier bootstrap approach (15) to the problems of constructing confidence sets in online algorithms has been demonstrated in the works [22] and [56] . We aim to illustrate the tightness of our bounds for normal approximation outlined in Theorem 2 in the setting of TD learning with linear function approximation. To this end, we consider the classical Garnet problem [3] , in the simplified version proposed by [25] . This problem is characterized by the number of states N s , number of actions a, and branching factor b (i.e. the number of neighbors of each state in the MDP). We set these values to N s = 10, a = 2 and b = 3, and aim to evaluate the value function of the randomly generated policy π(•|s). Details on the way the policy π is set can be found in Appendix F. We consider the problem of policy evaluation in this MDP using the TD learning algorithm with identity feature mapping, that is, φ(s) = e s (that is, s-th coordinate vector) for s ∈ {1, . . . , N s }. We run the procedure ( 22) with the learning rates α k = c 0 /k γ and different powers γ ∈ {0.5, 0.7, 0.75}. For each of the experiments we aim to estimate the supremumEQUATIONη ∼ N (0, I Ns ), and show that this supremum scales as n -1/4 when γ = 1/2 and admits slower decay for other powers of γ. We approximate true probability P( Σ1/2∞ η ≤ x) by the corresponding empirical probabilities based on sample of size M n. Second, for n ∈ {1600, . . . , 1638400}, where next sample size is twice larger than the previous one, we generate N = 1638400 trajectories of TD algorithm and approximate the distribution of √ n θn -θ based on the corresponding empirical distribution. We report our results in Figure 1 , showing that the smallest values of ∆ n correspond to the step size schedule γ = 1/2, moreover, the decay rate n -1/4 seems to be tight, otherwise one should expect further decay of ∆ n n 1/4 . Additional simulations are provided in Appendix F.  
conclusion: 
 In this paper, we have established, to the best of our knowledge, the first fully non-asymptotic confidence bounds for parameter estimation in the LSA algorithm using the multiplier bootstrap. This result is based on a novel Berry-Esseen bound for the Polyak-Ruppert averaged LSA iterates, which is of independent interest. Our paper suggests several interesting directions for further research. First, our Berry-Esseen bounds are obtained using the randomized concentration inequality [63] , and it would be valuable to generalize this approach to the setting of Markov chains. Second, it is natural to extend our results to the first-order gradient methods, both for stochastic optimization and variational inequalities. Third, it becomes possible to prove the fully non-asymptotic validity of confidence intervals obtained with plug-in techniques or other estimators of the asymptotic covariance matrix of θn . These could then be compared with the multiplier bootstrap confidence intervals in terms of their dependence on problem dimension d and other instance-dependent quantities. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes]Justification: All necessary information to reproduce experiments is provided in Appendix F.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?Answer: [NA]Justification: Not applicable: paper does not release new assets.Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?Answer: [NA] Justification: Not applicable: paper does not involve crowdsourcing nor research on human subjects.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
institutional review board (irb) approvals or equivalent for research with human subjects: 
 Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?Answer: [NA]  
 
 Appendix: A.1 Expansion of the LSA-PR error Proposition 3. The following expansion holds:√ n Ā( θn -θ ) = - 1 √ n 2n k=n+1 ε k W + 1 √ n θ n -θ α n D1 - 1 √ n θ 2n -θ α 2n D2 - 1 √ n 2n k=n+1 (A k -Ā)(θ k-1 -θ ) D3 + 1 √ n 2n k=n+1 θ k-1 -θ 1 α k - 1 α k-1 D4 (24)Proof. We use expansion (3) and rewrite it asθ k -θ = (I -α k Ā)(θ k-1 -θ ) -α k (A k -Ā)(θ k-1 -θ ) -α k ε k .The previous equation implies, after algebraic manipulation and division by α k , thatĀ(θ k-1 -θ ) = θ k-1 -θ α k - θ k -θ α k -(A k -Ā)(θ k-1 -θ ) -ε k .Taking average for k from n + 1 to 2n and multiplying by √ n we obtain the statement.We begin with of technical lemma on the behaviour of the last iterate of the LSA procedure θ n given in (1). This result is classical and can be traced back to [53] , with the modern exposition re-appearing in [7, 14, 17] . We provide the proof here for completeness. In the proof we use conventions for the product of random matrices Γ 1:k outlined in (67) . Proposition 4. Assume A1, A2, and A3. Then for any k ≥ n, where n satisfies (9), it holds for 2 ≤ p ≤ log n 2 , thatE 1/p [ θ k -θ p ] ≤ √ κ Q e exp -(a/2) k =1 α θ 0 -θ + 4e √ κ Q ε ∞ p √ a √ α k .Proof. Expanding the decomposition (66), we obtain thatEQUATIONand we bound the both term separately. Since the sample size n satisfies (9), we get applying Corollary 4 (see (69) ), that for 2 ≤ p ≤ log n it holdsE 1/p [ Γ 1:k {θ 0 -θ } p ] ≤ √ κ Q e exp -(a/2) k =1 α θ 0 -θ .Now we proceed with the second term in (25) . Applying Burholder's inequality [50] , we obtain thatE 1/p [ k j=1 α j Γ j+1:k ε j p ] ≤ p E 2/p k j=1 α 2 j Γ j+1:k ε j 2 p/2 1/2 ≤ p k j=1 α 2 j E 2/p Γ j+1:k ε j p 1/2 ≤ p √ κ Q e ε ∞ k j=1 α 2 j k =j+1 1 - aα 4 1/2 ≤ 4e √ κ Q ε ∞ p √ a √ α k .Corollary 2. Under assumptions of Proposition 4, it holds thatP (∀k ∈ [n, 2n -1] : θ k -θ ≥ g(k, θ 0 -θ , n)) ≤ 1 n ,where we have definedg(k, θ 0 -θ , n) = √ κ Q e 2 exp -(a/2) k =1 α θ 0 -θ + 8e 2 √ κ Q ε ∞ log n √ a √ α k .Proof. We first note that Lemma 1 implies, setting δ = 1/n 2 , that for every fixed k ∈ [n; 2n -1],P θ k -θ ≥ √ κ Q e 2 exp -(a/2) k =1 α θ 0 -θ + 8e 2 √ κ Q ε ∞ log n √ a √ α k ≤ 1 n 2 .Application of the union bound concludes the proof.We conclude this part with a simple consequence of Markov's inequality. Lemma 1. Fix δ ∈ (0, 1) and let Y be a positive random variable, such thatE 1/p [Y p ] ≤ C 1 + C 2 p for any 2 ≤ p ≤ log (1/δ). Then it holds with probability at least 1 -δ, that Y ≤ eC 1 + eC 2 log (1/δ) .Proof. Applying Markov's inequality, for any t ≥ 0 we get thatP(Y ≥ t) ≤ E[Y p ] t p ≤ (C 1 + C 2 p) p t p .Now we set p = log (1/δ), t = eC 1 + eC 2 log (1/δ), and aim to check that(C 1 + C 2 log (1/δ)) log (1/δ) (eC 1 + eC 2 log (1/δ)) log (1/δ) ≤ δ .Taking logarithms from both sides, the latter inequality is equivalent tolog (1/δ) log C 1 + C 2 log (1/δ) e(C 1 + C 2 log (1/δ))≤ log δ , which turns into exact equality.Proof. Since both terms in the right-hand side of the error bound of Proposition 4 scales linearly with √ κ Q , for simplicity we do not trace it in the subsequent bounds (i.e. assume κ Q = 1), and then keep the required scaling with κ Q only in the final bounds. The decomposition ( 24) is a key element of our proof and allows to treat different error sources D 1 -D 4 separately. For the last iterate we have, using Proposition 4, thatE 1/2 θ n -θ 2 ε ∞ √ a √ α n + exp -(a/2) n =1 α θ 0 -θ E 1/2 θ 2n -θ 2 ε ∞ √ a √ α 2n + exp -(a/2) 2n =1 α θ 0 -θ .Thus, using thatn k=1 α k ≥ c0n 1-γ 1-γ , we obtain that E 1/2 D 1 2 ε ∞ √ an (1-γ)/2 + exp - c 0 an 1-γ 2(1 -γ) θ 0 -θ , E 1/2 D 2 2 ε ∞ √ an (1-γ)/2 + exp - c 0 a(2n) 1-γ 1 -γ θ 0 -θ .Now we proceed with D 3 . Since it is a sum of a martingale-difference sequence w.r.t. F k = σ(Z , ≤ k), we get using Proposition 4, thatE D 3 2 C 2 A n 2n-1 k=n E[ θ k-1 -θ 2 ] C 2 A n 2n k=n+1 ε 2 ∞ α k a + C 2 A n 2n k=n+1 exp -a k =1 α θ 0 -θ 2 C 2 A n 2n k=n+1 ε 2 ∞ α k a + C 2 A nα 2n exp -a n =1 α 2n k=n+1 α k exp -a k =n+1 α S1 θ 0 -θ 2 C 2 A ε 2 ∞ an γ + C 2 A n 1-γ a exp - c 0 an 1-γ 1 -γ θ 0 -θ 2 ,where we additionally used that S 1 1/a due to Lemma 3. Now it remains to bound the term D 4 from the representation (24) . Using Minkowski's inequality and Proposition 4, we get thatE 1/2 D 4 2 1 √ n 2n k=n+1 E 1/2 θ k -θ 2 1 α k - 1 α k-1 1 √ n 2n k=n+1 ε ∞ (k γ -(k -1) γ ) √ a √ α k + 1 √ n 2n k=n+1 (k γ -(k -1) γ ) exp -(a/2) k =1 α θ 0 -θ ε ∞ √ a √ n 2n k=n+1 1 k 1-γ/2 + n 2γ-1 a exp - c 0 an 1-γ 2(1 -γ) θ 0 -θ ε ∞ √ an (1-γ)/2 + n 2γ-1 a exp - c 0 an 1-γ 2(1 -γ) θ 0 -θ .Combining the estimates above yields the result of Theorem 1.We conclude this section with some technical lemmas. Lemma 2 (Lemma 24 in [18] ). Let b > 0 and (α k ) k≥0 be a non-increasing sequence such thatα 0 ≤ 1/b. Then n+1 j=1 α j n+1 l=j+1 (1 -α l b) = 1 b 1 - n+1 l=1 (1 -α l b)Proof. The proof of this statement is given in [18] , we provide it here for completeness. Let us denoteu j:n+1 = n+1 l=j (1 -α l b). Then, for j ∈ {1, . . . , n + 1}, u j+1:n+1 -u j:n+1 = bα j u j+1:n+1 . Hence, n+1 j=1 α j n+1 l=j+1 (1 -α l b) = 1 b n+1 j=1 (u j+1:n+1 -u j:n+1 ) = b -1 (1 -u 1:n+1 ) ,and the statement follows.Lemma 3 (Modified Lemma 25 in [18] ). Let b > 0 and let α = c 0 / γ , γ ∈ [1/2; 1), such that c 0 ≤ 1/b. Then for any n satisfyingn ≥ 2 + 2 2γ c 0 b 1/(1-γ), andEQUATION)any k ≥ n, it holds that k+1 j=1 α 2 j k+1 =j+1 (1 -α b) ≤ (4/b)α k+1 .Proof. From elementary algebra, we obtain thatEQUATIONwhere we used the fact that(1 + x) γ ≤ 1 + γx for γ ∈ [1/2; 1) and x ∈ [0, 1]. Hence, α α +1 ≤ 1 + γ .Thus we obtain that, since k ≥ nk+1 j=1 α 2 j k+1 =j+1 (1 -α b) = α k+1 k+1 j=1 α j k+1 =j+1 α -1 α (1 -α b) ≤ α k+1 k+1 j=1 α j k+1 =j+1 1 + γ -1 (1 -α b) ≤ α k+1 k+1 j=1 α j exp    n =j+1 γ -1    exp    - n =j+1 α b    exp k+1 =n+1 γ -1 exp - k+1 =n+1 α b ≤ α k+1 k+1 j=1 α j exp    n =j+1 γ -1    exp    - n =j+1 α b    exp - b 2 k+1 =n+1 α b ,where we have used in the last identity that, since n satisfies (26) , that forEQUATION)We will now prove that for j ≤ n -1, it holdsEQUATIONFor j ≥ n/2, the bound (29) directly follows from (28) . We now turn to the proof of (29) for j ≤ n/2. Note first thatn =j+1 γ -1 ≤ n k=2 γ -1 ≤ γ(log(n) + 1) .On the other hand, we getn =j+1 1 γ ≥ n+1 j+1 dx x γ = (n + 1) 1-γ -(j + 1) 1-γ 1 -γ .Comparing the above bounds, to ensure that (29) holds, it is enough to check thatEQUATIONNote that (30) is guaranteed by (26) . Thus we obtain thatk+1 j=1 α 2 j k+1 =j+1 (1 -α b) ≤ α k+1 k+1 j=1 α j exp -(b/2) k+1 =j+1 α ≤ α k+1 k+1 j=1 α j k+1 =j+1 1 -(b/4)α ≤ (4/b) α k+1 ,where the last inequality follows from Lemma 2.We first define explicitly the remainder terms outlined in the statement of Theorem 2:EQUATIONand constants C 1 , C 2 are given, respectively, byEQUATIONTo complete the proof we only need to combine (13) with the bounds of Theorem 1. Applying the Cauchy-Schwartz inequality, we getE[ D W ] ≤ E 1/2 [ D 2 ]E 1/2 [ W 2 ] √ κ Q ε ∞ √ Tr Σ ε √ a 1 n (1-γ)/2 + C A n γ/2 + √ κ Q C A √ Tr Σ ε n (1-γ)/2 √ a + √ Tr Σ ε n 2γ-1 a exp - c 0 an 1-γ 2(1 -γ) θ 0 -θ .Now it remains to bound the last term in (13) . Using the Cauchy-Schwartz inequality and Lemma 4, we obtain thatE[ 2n-1 i=n ε i D -D (i) ] ≤ E 1/2 [ ε 1 2 ] 2n-1 i=n E 1/2 [ D -D (i) 2 ] κ Q (Tr Σ ε + C A √ Tr Σ ε ) a n γ-1/2 + κ Q C A √ Tr Σ ε √ a ε ∞ + Tr Σ ε + C A n 1/2-γ/2 + κ Q √ Tr Σ ε a + C A √ Tr Σ ε n γ-1/2 a 2 exp - c 0 an 1-γ 2(1 -γ) θ 0 -θ ,and the statement follows.B.1 Proof of auxiliary lemmas for Theorem 2.Our proof of Theorem 2 is based on the key lemma below, which allows us to boundE 1/2 [ D - D (i) 2 ] for i ∈ {n + 1, . . . , 2n}.Lemma 4. Assume A1, A2, and A3. Then2n i=n+1 E 1/2 [ D -D (i) 2 ] κ Q ( √ Tr Σ ε + C A ) a n γ-1/2 + κ Q C A √ a ε ∞ + Tr Σ ε + C A n 1/2-γ/2 +κ Q 1 a + C A n γ-1/2 a 2 exp - c 0 an 1-γ 2(1 -γ) θ 0 -θ .Proof. Since both terms in the right-hand side of the error bound of Proposition 4 scales linearly with √ κ Q , for simplicity we do not trace it in the subsequent bounds (i.e. assume κ Q = 1), and then keep the required scaling with κ Q only in the final bounds. Consider the sequences of noise variables (Z 1 , . . . , Z i-1 , Z i , Z i+1 , . . . , Z 2n ) and (Z 1 , . . . , Z i-1 , Z i , Z i+1 , . . . , Z 2n ) , which differs only in position i, n + 1 ≤ i ≤ 2n, with Z i being an independent copy of Z i . Consider the associated SA processesEQUATIONwhereY k = Z k for k = i and Y i = Z i .From the above representations we easily observe thatEQUATIONwhereε i = ε(Z i ) and ε i = ε(Z i ).Representation (34) implies, together with Proposition 4, thatEQUATIONMoreover, for any j > i one observes, expanding (33) , thatEQUATIONNow we note that, using Minkowski's inequality,EQUATIONand bound the respective differences separately. Recall that here D 1 -D 4 are defined in (24) , andD (i) 1 -D (i)4 are their respective counterparts with Z i substituted with Z i . First we note that the termD 1 = D (i)1 for any n + 1 ≤ i ≤ 2n. Next, using (36) and (35) , we getE 1/2 [ D 2 -D (i) 2 2 ] = 1 √ nα 2n E 1/2 [ θ 2n -θ (i) 2n 2 ] ≤ 1 √ nα 2n E 1/2 2n k=i+1 (I -α k A k ) 2 E 1/2 [ θ i -θ (i) i 2 ] (a) α i ( √ Tr Σ ε + C A ε ∞ ) √ nα 2n exp - a 2 2n k=i+1 α k + α i √ nα 2n exp - a 2 2n k=1 α k θ 0 -θ .In the inequality (a) above we additionally used the stability of matrix product introduced from Corollary 4. Summing the above inequality for i = n + 1 to 2n and applying Lemma 2, we getEQUATION+ 1 a √ nα 2n exp - a 2 n k=1 α k θ 0 -θ √ Tr Σ ε + C A ε ∞ a n γ-1/2 + n γ-1/2 a exp - c 0 an 1-γ 2(1 -γ) θ 0 -θ .Now we proceed with the difference D 3 -D (i)3 . Using (24) , we getD 3 -D (i) 3 = 1 √ n (A i -A i )(θ i-1 -θ ) + 1 √ n 2n k=i+1 (A k -Ā)(θ k-1 -θ (i) k-1 ) .The expression above is a sum of martingale-difference terms w.r.t. filtrationF k = σ(Z i , Z , ≤ k).Hence, we get, using (36) , thatEQUATION)C 2 A ε 2 ∞ α i na + C 2 A θ 0 -θ 2 n exp -a i-1 j=1 α j + C 2 A n E[ θ i -θ (i) i 2 ] 2n k=i+1 exp -a k-1 j=i+1 α j .Using now the bound (35) , we obtain thatE[ θ i -θ (i) i 2 ] 2n k=i+1 exp -a k-1 j=i+1 α j α 2 i Tr Σ ε + C 2 A α 3 i a 2n k=i+1 exp -a k-1 j=i+1 α j + α 2 i θ 0 -θ 2 2n k=i+1 exp -a k-1 j=1 α i α 2 i Tr Σ ε α 2n + C 2 A α 3 i aα 2n 2n k=i+1 α k exp -a k-1 j=i+1 α j + α 2 i α 2n θ 0 -θ 2 2n k=i+1 α k exp -a k-1 j=1 α i (a) α 2 i Tr Σ ε aα 2n + C 2 A α 3 i a 2 α 2n + α 2 i aα 2n θ 0 -θ 2 exp -a i-1 j=1 α j (b) α 2 i (Tr Σ ε + C 2 A ) aα 2n + α 2 i aα 2n θ 0 -θ 2 exp -a i-1 j=1 α j .In the above formula in (a) we additionally used thatEQUATIONand in (b) we used that α i ≤ a. Hence, combining everything in (39) , and using additionally thatα i ≤ a, κ Q ≥ 1, we get E 1/2 [ D 3 -D (i) 3 2 ] C A √ na ε ∞ √ α i + α i ( √ Tr Σ ε + C A ) √ α 2n + C A √ n 1 + α i √ aα 2n exp - a 2 i-1 j=1 α i θ 0 -θ .Summing the above inequality for i = n + 1 to 2n, and using that α k = c 0 /k γ , we get2n i=n+1 √ α i n 1-γ/2 , 2n i=n+1 α i √ α 2n n 1-γ/2 ,and, hence,2n i=n+1 E 1/2 [ D 3 -D (i) 3 2 ] C A √ a ε ∞ + Tr Σ ε + C A n 1/2-γ/2 + C A √ nα 2n θ 0 -θ 2n i=n+1 α i exp - a 2 i-1 j=1 α j (b) C A √ a ε ∞ + Tr Σ ε + C A n 1/2-γ/2 + C A n γ-1/2 a θ 0 -θ exp - a 2 n j=1 α j C A √ a ε ∞ + Tr Σ ε + C A n 1/2-γ/2 + C A n γ-1/2 a 2 θ 0 -θ exp - c 0 an 1-γ 2(1 -γ) ,where for the last identity we used the fact that α k = c 0 /k γ , and for the bound (b) we used (40) . It remains to upper bound the difference D 4 -D4 . Note first that, proceeding as in (27) , we getα k-1 -α k ≤ γ k -1 1 α k-1 1 (k -1) 1-γ .Using now the definition of D 4 in ( 24), we have thatE 1/2 [ D 4 -D (i) 4 2 ] = 1 √ n E 1/2 [ 2n k=i+1 θ k-1 -θ (i) k-1 1 α k - 1 α k-1 2 ] ≤ 1 √ n E 1/2 [ θ i -θ (i) i 2 ] 2n k=i+1 1 α k - 1 α k-1 exp - a 2 k-1 j=i+1 α j .Hence, using the bound (35) and taking sum for i = n + 1 to 2n,2n i=n+1 E 1/2 [ D 4 -D (i) 4 2 ] √ Tr Σ ε + C A √ n 2n i=n+1 α i 2n k=i+1 1 α k - 1 α k-1 exp -a k-1 j=i+1 α j + θ 0 -θ √ n 2n i=n+1 α i 2n k=i+1 1 α k - 1 α k-1 exp - a 2 k-1 j=1 α j .Changing now the summation order, we obtain that2n i=n+1 α i 2n k=i+1 1 α k - 1 α k-1 exp -a k-1 j=i+1 α j 1 a 2n k=n+2 1 α k - 1 α k-1 = 1 a 1 α n+1 - 1 α 2n 1 aα 2n .Hence, combining the above bounds, we getEQUATION√ Tr Σ ε + C A a n γ-1/2 + θ 0 -θ √ na exp - c 0 an 1-γ 2(1 -γ) .It remains now to combine (38) , (39) , and (41) in (37) .In this section we closely follow the exposition outlined in [24] . Consider two R d -valued random variables X and Y . Then the integral probability metric [75] , associated with the class of functionsH = {h : R d → R, E[|h(X)|] < ∞, E[|h(Y )|] < ∞}, is defined as d H (X, Y ) = sup h∈H E[h(X)] -E[h(Y )] .Different choices of H induce different metrics, in particular, we mention the following:H K = {1 x≤u , u = (u 1 , . . . , u d ) ∈ R d } H Conv = {1 x∈B , B ∈ Conv(R d )} H W = {h : R d → R, h Lip ≤ 1} H [m] = {h : R d → R, h m-1 is Lipschitz with |h| j ≤ 1 , 1 ≤ j ≤ m} ,where Conv(R d ) refers to the set of convex sets in R d , h Lip = sup x =y h(x)-h(y)x-y, and the quantity |h| j is defined as|h| j = max i1,...,ij ∈{1,...,d} ∂ j h(u) ∂u i1 . . . ∂u ij ∞ .In other words, for m ∈ N, the class H d K (X, Y ) ≤ ρ Conv n (X, Y )for any random vectors X and Y . Other relations are more involved. When Y is a multivariate normal vector, it is known (see e.g. [49] ) thatρ Conv n (X, Y ) ≤ C d W (X, Y ) ,where the constant C in the above inequality depends on the covariance matrix of vector Y . This inequality justifies comparison of our bounds of Theorem 2 with the result of [65] . The authors in [2] considered integral probability metric d [2] and obtained rate of convergenced [2] ( √ n( θn -θ ), Y ) ≤ C 1 √ n ,where Y ∼ N (0, Σ ∞ ), and C 1 in the above inequality stands for a constant depending upon problem dimension d and other instance-dependent parameters from A2. Applying the result of [24, Proposition 2.6] yieldsd K ( √ n( θn -θ ), Y ) d [2] ( √ n( θn -θ ), Y ) 1/3 1 n 1/6 .Thus, the result of [2] implies rate of convergence of √ n( θn -θ ) to normal law N (0, Σ ∞ ) of order n -1/6 in a sense of Kolmogorov distance d K . Our result of Theorem 2 implies the respective rate of order n -1/4 . At the same time, it is not clear if ρ Conv n can be directly related to d [2] .We first define explicitly the remainder term ∆ 2 outlined in the statement of Theorem 3, that is,EQUATION)We now define the following sets, with the convention α = c 0 / √ :EQUATION+ 8e 2 √ κ Q ε ∞ log n √ a √ α k , Ω 2 =    n + 1 ≤ m ≤ k ≤ 2n : Γ m:k ≤ √ κ Q e 2 k j=m 1 - aα j 4    , Ω 3 = Σ -1/2 ε Σ b ε Σ -1/2 ε -I ≤ 4 ε ∞ log(n) σn + 4(1 + ε 2 ∞ /σ 2 ) log(n) n , Ω 4 = ∀ ∈ [n, 2n -1] : 2n k= +1 (A k -Ā)Γ +1:k-1 ≤ 8 C A √ κ Q e 2 √ log n √ aα + 6 C A √ κ Q e log n ,Ω 5 =    ∀h ∈ [1; n] , ∀m ∈ [n, 2n -h] : m+h =m+1 α (A -Ā) Q ≤ 2 C A √ κ Q m+h =m+1 α 2 log(2n 4 )    , Ω 6 = Σ b ε -Σ ε ≤ 4 ε ∞ Σ ε log(n) n + 4( Σ ε + ε 2 ∞ ) log(n) n ,Then, due to Corollary 2, we have that P(Ω 1 ) ≥ 1 -1 n . Similarly, due to Corollary 5, P(Ω 2 ) ≥ 1 -1 n . The bounds on P(Ω 3 ) and P(Ω 4 ) follows from Lemma 5 and Lemma 6, respectively. Similarly, Proposition 7 implies that P(Ω 5 ) ≥ 1 -1 n . Hence, based on the sets above, we can constructΩ 0 = Ω 1 ∩ Ω 2 ∩ Ω 3 ∩ Ω 4 ∩ Ω 5 ∩ Ω 6 , such that P(Ω 0 ) ≥ 1 -6n . All further on, we restrict ourselves to the event Ω 0 . Restricting to this event, we obtain that, with Minkowski's inequality,sup B∈Conv(R d ) |P b ( √ n( θb n -θn ) ∈ B) -P( √ n( θn -θ ) ∈ B)| ≤ sup B∈Conv(R d ) P b ( √ n( θb n -θn ) ∈ B) -P b (ξ b ∈ B) + sup B∈Conv(R d ) P(ξ ∈ B) -P b (ξ b ∈ B) + sup B∈Conv(R d ) P √ n( θn -θ ) ∈ B -P(ξ ∈ B) ,where we setξ b ∼ N (0, Ā-1 Σ b ε Ā-), Σ b ε = n -1 n =1 ε ε , and ξ ∼ N (0, Σ ∞ ), where Σ ∞ = Ā-1 Σ εĀ-. Now we control the first supremum using Theorem 4, second one using Theorem 5, and third with Theorem 2. Lemma 5. Assume A1, A2, A3 with γ = 1/2, and A4. ThenP(Ω 3 ) ≥ 1 -1/n.Proof. The proof follows directly from the matrix Bernstein inequality, e.g. [69] . We note thatΣ -1/2 ε ε ε Σ -1/2 ε -I ≤ 1 + Σ -1/2 ε ε 2 ∞ . and 2n k=n+1 E[(Σ -1/2 ε ε ε Σ -1/2 ε -I) 2 ] ≤ nE Σ -1/2 ε ε 2 ∞ .Lemma 6. Assume A1, A2, A3 with γ = 1/2, and A4. ThenP (Ω 4 ∩ Ω 2 ) ≥ 1 - 1 n . Proof. Denote X k = (A k -Ā)Γ +1:k-1 .and letF k,l+1 = σ{Z j , + 1 ≤ j ≤ k}, + 1 ≤ k ≤ 2n. Then E[X k |F k-1, +1 ] = 0. Let S = n k= +1 X k .Note that on Ω 2 , quadratic variation of S can be controlled asVar 2 := max( 2n k= +1 E[X k X k |F k-1,l+1 ] , 2n k= +1 E[X k X k |F k-1,l+1 ] ) ≤ κ Q e 4 C 2 A 2n k= +1 k-1 j= +1 (1 -aα j /4) 2 ≤ 4κ Q e 4 C 2 A aα Furthermore, on Ω 2 X k ≤ √ κ Q e 2 C A k-1 j= +1 (1 -aα j /4) ≤ √ κ Q e 2 C A .It remains to apply the Freedman inequality for matrix-values martingales [68] and use the union bound over ∈ [n, 2n -1].Lemma 7. Assume A1, A2, A3 with γ = 1/2, and A4. ThenP(Ω 5 ) ≥ 1 -1/n. Proof. We first fix h ∈ [1; n], m ∈ [n, 2n -h],and consider the random variableT n = m+h =m+1 α (A -Ā) .Then we control its variance as max(m+h =m+1 α 2 E[(A -Ā)(A -Ā) ] , m+h =m+1 α 2 E[(A -Ā) (A -Ā)] ) ≤ C A m+h =m+1 α 2 , moreover, (A -Ā)(A -Ā) ≤ C 2 A .Applying now the matrix Bernstein inequality [69] , we obtain that with probability at least 1 -1/n 3 , we haveT n ≤ C A 2 m+h =m+1 α 2 log (2n 3 d) + α m+1 C A 3 log (2n 3 d) ≤ 2 C A m+h =m+1 α 2 log (2n 4 ) .In the last line here we used that d ≤ n. Rest of the proof follows by taking union bound over h and m together withB 2 Q ≤ κ Q B 2 valid for any matrix B ∈ R d×d .Lemma 8. Assume A1, A2, A3 with γ = 1/2, and A4. ThenP(Ω 6 ) ≥ 1 -1/n. Proof. It is easy to check that ε ε -Σ ε ≤ Σ ε + ε 2 ∞ . Moreover, Var 2 := n =1 (ε ε -Σ ε ) 2 ≤ n ε 2 ∞ Σ ε .It remains to apply the matrix Bernstein inequality together with the bound n ≥ d from A3.The main result of this section is the following theorem. Theorem 4. Assume A1, A2, A3 with γ = 1/2, and A4. Then, conditionally on the event Ω 0 , the following error bound holds:sup B∈Conv(R d ) P b ( √ n( θb n -θn ) B) -P b (ξ b ∈ B) d 1/2 ε 3 ∞ λ 3/2 min √ n + κ 3/2 Q (C 3 A ∨1) ε 2 ∞ log n a 3/2 λ min n 1/4 + κ 3/2 Q (C 3 A ∨1) ε ∞ n 1/4 √ log n a 3/2 λ min exp - a 4 n j=1 α j θ 0 -θ ,whereξ b ∼ N (0, Ā-1 Σ b ε Ā-) and Σ b ε = n -1 n =1 ε ε .Proof. Since both terms in the right-hand side of the error bound of Proposition 4 scales linearly with √ κ Q , for simplicity we do not trace it in the subsequent bounds (i.e. assume κ Q = 1), and then keep the required scaling with κ Q only in the final bounds. Recall first that the quantities θ b k and θ k are defined in (15) . We start from the following decomposition:θ b k -θ k = (I-α k Ā)(θ b k-1 -θ k-1 )-α k (w k -1)ε k -α k (A k -Ā)(θ b k-1 -θ k-1 )-α k (w k -1)A k (θ b k-1 -θ ) .Taking average for k from n + 1 to 2n, we get after multiplying byEQUATIONThe formula (44) resembles the key representationEQUATIONand D b 1 -D b 5 are defined in (44) . Now we aim to apply the result of [63] :EQUATIONwhere ξ = 1 √ n (w -1)ε . We finish the proof by the application of the formula (46) . In order to bound the quantities E b D b 2 and E b D (b,i) 2 , we apply the respective results of Proposition 5 and Proposition 6, respectively. Namely, applying the Cauchy-Schwartz inequality together with Proposition 5, we get that on the event Ω 0 it holdsE b [ D b W b ] ≤ E b [ D b 2 ] 1/2 {E b [ W b 2 ]} 1/2 κ 2 Q (C 4 A ∨1) ε ∞ Tr Σ b ε log n n 1/4 a 5/2 + κ 3/2 Q ε ∞ (C 3 A ∨1)n 1/4 √ a + (C 5 A ∨1) √ log n √ na exp - c 0 a √ n 2 θ 0 -θ .Similarly, applying Minkowski's inequality and Proposition 6, we obtain thatE b [ 2n-1 i=n ξ i D b -D (b,i) ] ≤ E b [ ξ 1 2 ] 1/2 2n-1 i=n E b [ D b -D (b,i) 2 ] 1/2 κ 3/2 Q (C 3 A ∨1) ε 2 ∞ log n a 3/2 n + κ Q ε 2 ∞ C 2 A a 2 √ n + κ 3/2 Q (C 3 A ∨1) ε ∞ n 1/4 √ log n a 3/2 exp - a 4 n j=1 α j θ 0 -θ .Now it remains to combine the bounds above in (46) .Proposition 5. Assume A1, A2, A3 with γ = 1/2, and A4. Then, conditionally on the event Ω 0 , the following error bound holds:E b D b (w 1 , . . . , w 2n , Z 1 , . . . , Z 2n ) 2 1/2 κ 2 Q (C 4 A ∨1) ε ∞ log n n 1/4 a 5/2 + κ 3/2 Q (C 3 A ∨1)n 1/4 √ a + (C 5 A ∨1) √ log n √ na exp - c 0 a √ n 2 θ 0 -θ ,where stands for inequality up to an absolute constant.Proof of Proposition 5 is provided below in Appendix C.5. The lemma below is a direct counterpart of Lemma 4. Proposition 6. Assume A1, A2, A3 with γ = 1/2, and A4. Then, conditionally on the event Ω 0 , the following error bound holds:2n i=n+1 E b [ D b -D (b,i) 2 ] 1/2 (C 3 A ∨1) ε ∞ a 3/2 1/4 log n + ε ∞ C 2 A a 2 + (C 3 A ∨1)n 3/4 √ log n a 3/2 exp - a 4 n j=1 α j θ 0 -θ ,Proof of Proposition 6 is provided below in Appendix C.6. Lemma 9. For any k ≥ n on the set Ω 0 the following inequality holds:E b [ θ b k -θ k 2 ] α k ε 2 ∞ C 2 A a 3 + C 2 A k k j=1 (1 -aα j /4) 2 θ 0 -θ 2 .Proof. A direct application of Lemma 11 with L = 0 yields thatE b [ θ b k -θ k 2 ] α k ε 2 ∞ a 1 + C 2 A a 2 + C 2 A k k j=1 (1 -aα j /4) 2 θ 0 -θ 2 .Now to complete the proof it remains to notice that C A ≥ a.Lemma 10. For any matrix-valued sequences (U n ) n∈N , (V n ) n∈N and for any M ∈ N, it holds that:M k=1 U k - M k=1 V k = M k=1 { M j=k+1 V j }(U k -V k ){ k-1 j=1 U j } .Theorem 5. Assume A1 and A2. Then on the setΩ 3 sup B∈Conv(R d ) |P(ξ ∈ B) -P b (ξ b ∈ B)| ≤ 4 Σ -1/2 ε ε ∞ d log n n + 4 √ d(1 + Σ -1/2 ε ε 2 ∞ ) log n nProof. We will use the following inequalityEQUATIONApplying ( 47) we obtain supB∈Conv(R d ) |P(ξ ∈ B) -P b (ξ b ∈ B) ≤ √ d 2 Σ -1/2 ε Σ b ε Σ -1/2 ε -I .It remains apply definition of Ω 3 .For the analysis of the difference term θ b k -θ k we use the perturbation expansion technique introduced in [1] , see also [16] . Within this approach, we represent the fluctuation component of the error θ(fl) n defined in (66) as θ(fl) n = J (0) n + H (0) n , where the latter terms are defined by the following pair of recursionsEQUATIONH (0) n = (I -α n A(Z n )) H (0) n-1 -α n Ã(Z n )J (0) n-1 , H (0) 0 = 0 . Moreover, it is known that for L ≥ 1 the term H (0)n can be further decomposed as follows:H (0) n = L =1 J ( ) n + H (L) n .Here the terms J ( )n and H ( )n are given by the following recurrences:EQUATIONThe expansion depth L here controls the desired approximation accuracy. Informally, one can show that E . Using the outlined expansion, we prove the following lemma: Lemma 11. Assume A1, A2, A3 with γ = 1/2, and A4. Then for any k ≥ n and L ∈ N the following decomposition holds:EQUATIONwhereEQUATIONand the quantities ε are defined asEQUATION)Proof. We start from the decompositionEQUATION)Expanding the recurrence above till k = n, and using the fact that θ b n = θ n , we get running the recurrence (54) ,θ b k -θ k = - k =n+1 α (w -1)Γ b +1:k ε .Hence, proceeding as in (48) , we obtain the representationJ (b,0) k = (I -α k A k ) J (b,0) k-1 -α k (w k -1) εk , J (b,0) 0 = 0 , H (b,0) k = (I -α k w k A k ) H (b,0) k-1 -α k (w k -1)A k J (b,0) k-1 , H (b,0) 0 = 0 .It is easy to check that J(b,0) k + H (b,0) k = θ b k -θ k .Similarly, with further expansion of H (b,0) k along the lines of ( 49), we arrive at the decomposition (50) . Since w k for k = n + 1, . . . , 2n are i.i.d., we get using the definition of the events Ω 1 and Ω 2 , that on the event Ω 0 :EQUATIONwhere for the last bound we have additionally used that α log 2 n/a ≤ 1 for ≥ n. The latter bound is guaranteed by A4. Hence, using the bound (55) together with the definition of J b,0 k , we obtain thatE b [ J b,0 k 2 ] = k =n+1 α 2 Γ +1:k ε 2 = k =n+1 α 2 Γ +1:k A (θ -1 -θ ) + ε 2 ε 2 ∞ k =n+1 α 2 k j= +1 (1 -aα j /4) 2 + C 2 A k =n+1 α 2 k j=1 (1 -aα j /4) 2 θ 0 -θ 2 + ε 2 ∞ C 2 A log 2 n a k =n+1 α 3 k j= +1 (1 -aα j /4) 2 ε 2 ∞ α k a + C 2 A log k n k j=1 (1 -aα j /4) 2 θ 0 -θ 2 + α 2 k ε 2 ∞ C 2 A log 2 n a 2 ε 2 ∞ α k a + C 2 A k j=1 (1 -aα j /4) 2 θ 0 -θ 2 ,where we additionally used the fact that k ∈ [n; 2n] and n satisfies A4. Assume now that the bound on J b,j-1k has a form E b [ J b,j-1 k 2 ] ε 2 ∞ α j k a j + C 2j A k =1 (1 -aα /4) 2 θ 0 -θ 2 .Then, using the martingale property of J b,j k , we write thatE b [ J b,j k 2 ] = k =n+1 α 2 E b [ Γ +1:k A J b,j-1 -1 2 ] k =n+1 α j+2 ε 2 ∞ C 2j A a j k j= +1 (1 -aα j /4) 2 + C 2j+2 A k =n+1 α 2 k j=1 (1 aα j /4) 2 θ 0 -θ 2 α j+1 k ε 2 ∞ C 2j A a j+1 + C 2j+2 A k j=1 (1 -aα j /4) 2 θ 0 -θ 2 ,and thus the bound ( 52) is proved. Moreover, using ( 51) and Minkowski's inequality, we obtain that(E b [ H b,L k 2 ]) 1/2 ≤ C A k =n+1 α (E b [ Γ b +1:k 2 ]) 1/2 (E b [ J b,L -1 2 ]) 1/2 C A k =n+1 α (L+3)/2 ε ∞ C L A a (L+1)/2 k j= +1(1 -aα j /4)EQUATIONRecall that the quantity D b is defined in (45) . Since θ b n = θ n , we conclude that D b 1 = 0. To estimate other terms we will use the main error decomposition outlined in Lemma 11, that is, the expansionθ b k -θ k = L =0 J b, k + H b,L k ,applied with different L ≥ 0. To bound D b 2 we take L = 0 and obtainE b [ D b 2 2 ] E b [ J b,j 2n 2 ] + E b [ H b,L 2n 2 ] ε 2 ∞ nα 2n a 1 + C 2 A a 2 + C 2 A α 2 2n 2n j=1 (1 -aα j /4) 2 θ 0 -θ 2 ε 2 ∞ a √ n 1 + C 2 A a 2 + n C 2 A 2n j=1 (1 -aα j /4) 2 θ 0 -θ 2 .To estimate D b 3 we note thatD b 3 = 1 √ n 2n k=n+1 (w k -1)A k (θ b k-1 -θ ) = D b 3,1 + D b 3,2 ,where we have set, respectively,D b 3,1 = 1 √ n 2n k=n+1 (w k -1)A k (θ b k-1 -θ k-1 ), D b 3,2 = 1 √ n 2n k=n+1 (w k -1)A k (θ k-1 -θ ) .It follows from Lemma 9 that on the event Ω 0 it holdsE b [ D b 3,1 2 ] ≤ C 2 A n 2n k=n+1 E b [ θ b k-1 -θ k-1 2 ] C 2 A n 2n k=n+1 α k ε 2 ∞ a 1 + C 2 A a 2 + C 4 A 2n k=n+1 k j=1 (1 -aα j /4) 2 θ 0 -θ 2 C A ε 2 ∞ a √ n 1 + C 2 A a 2 + C 4 A √ n a exp -c 0 a √ n θ 0 -θ 2Moreover, on the set Ω 0 it holds (sinceΩ 0 ⊆ Ω 1 ), that E b [ D b 3,2 2 ] = 1 n 2n k=n+1 A k (θ k-1 -θ ) 2 C 2 A n 2n k=n+1 exp -a k =1 α θ 0 -θ 2 + α k ε 2 ∞ log 2 n a C 2 A naα 2n exp -a n =1 α θ 0 -θ 2 + C 2 A ε 2 ∞ log 2 n na 2n k=n+1 α k C 2 A a √ n exp -c 0 a √ n θ 0 -θ 2 + C 2 A ε 2 ∞ log 2 n a √ n .Combining the above bounds, we getE b [ D b 3 2 ] E b [ D b 3,1 2 ] + E b [ D b 3,2 2 ] C 2 A ε 2 ∞ log 2 n a √ n 1 + C 2 A a 2 + C 4 A √ n a exp -c 0 a √ n θ 0 -θ 2 .Now we proceed with the term D b 4 . Applying Minkowski's inequality, we get{E b [ D b 4 2 ]} 1/2 ≤ 1 √ n 2n k=n+1 1 α k - 1 α k-1 {E b [ θ b k-1 -θ k-1 2 ]} 1/2 1 √ n 2n k=n+1 1 √ k   √ α k ε ∞ √ a 1 + C A a + C A √ k k j=1 (1 -aα j /4) θ 0 -θ   ε ∞ n 1/4 √ a 1 + C A a + C A √ n 2n k=n+1 k j=1 (1 -aα j /4) θ 0 -θ ε ∞ n 1/4 √ a 1 + C A a + C A exp - c 0 a √ n 2 θ 0 -θ .It remains to upper bound the term D b 5 . Using the decomposition, suggested by Lemma 11 with L = 2, we get thatD b 5 = 1 √ n 2n k=n+1 (A k -Ā)(θ b k-1 -θ k-1 ) = 1 √ n 2n k=n+1 (A k -Ā)J b,0 k-1 D b 5,1 + 1 √ n 2n k=n+1 (A k -Ā)J b,1 k-1 D b 5,2 + 1 √ n 2n k=n+1 (A k -Ā)J b,2 k-1 D b 5,3 + 1 √ n 2n k=n+1 (A k -Ā)H b,2 k-1 D b 5,4.Here we have to consider expansion until H b,2 , since dealing with the latter term (outlined as D b 5,4 in the above expansion) is possible only with Minkowski's inequality. Now we consider the summands D b 5,1 -D b 5,4 separately. Consider first the term D b 5,1 . Changing the summation order, we obtainD b 5,1 = - 1 √ n 2n k=n+1 (A k -Ā) k-1 =n+1 α (w -1)Γ +1:k-1 ε = - 1 √ n 2n-1 =n+1 α (w -1) 2n k= +1 (A k -Ā)Γ +1:k-1 ε .Then on the event Ω 0 we get, since Ω 0 ⊆ Ω 4 , and using that n satisfies A4,EQUATION)Combining the above bound together with the one provided by ( 55), we obtain thatE b [ D b 5,1 2 ] ε 2 ∞ n 2n-1 =n+1 α log n a + C 2 A n 2n-1 =n+1 α log n a -1 j=1 1 -aα j /2 2 θ 0 -θ 2 ε 2 ∞ log n √ na + C 2 A log n a 2 n exp -a n j=1 α j θ 0 -θ 2 .Similarly, for the term D b 5,2 we get, changing the order of summation, thatD b 5,2 = 1 √ n 2n-1 =n+1 α (w -1) 2n k= +1 (A k -Ā)Γ +1:k-1 A J b,0 -1 .Hence, using the bound (56) together with (52), we getE b [ D b 5,2 2 ] 1 n 2n-1 =n+1 α log n a C 2 A α ε 2 ∞ a + C 2 A -1 j=1 (1 -aα j /4) 2 θ 0 -θ 2 C 2 A ε 2 ∞ log n na 2 2n-1 =n+1 α 2 + C 4 A log n na 2n-1 =n+1 α -1 j=1 (1 -aα j /4) 2 θ 0 -θ 2 C 2 A ε 2 ∞ log n na 2 + C 4 A log n na exp -(a/2) n j=1 α j θ 0 -θ 2 .We proceed with D b 5,3 . We change the summation order and proceed exactly as with D b 5,2 . Indeed,D b 5,3 = 1 √ n 2n-1 =n+1 α (w -1) 2n k= +1 (A k -Ā)Γ +1:k-1 A J b,1 -1 ,andE b [ D b 5,3 2 ] 1 n 2n-1 =n+1 α log n a C 2 A α 2 ε 2 ∞ C 2 A a 2 + C 4 A -1 j=1 (1 -aα j /4) 2 θ 0 -θ 2 C 4 A ε 2 ∞ log n n 3/2 a 3 + C 6 A log n na exp -(a/2) n j=1 α j θ 0 -θ 2 .It remains to upper bound D b 5,4 . Proceeding as above, we change the summation order, and obtainD b 5,4 = 1 √ n 2n-1 =n+1 α (w -1) 2n k= +1 (A k -Ā)Γ b +1:k-1 A J b,2 -1 .Applying Minkowski's inequality, we get(E b [ D b 5,4 2 ]) 1/2 C 2 A √ n 2n-1 =n+1 α 2n k= +1 (E b [ Γ b +1:k-1 2 ]) 1/2 (E b [ J b,2 -1 2 ]) 1/2 C 2 A √ n 2n-1 =n+1 α 5/2 2n k= +1 exp - a 4 k-1 j= +1 α j ε ∞ C 2 A a 3/2 + C 5 A √ n 2n-1 =n+1 α 2n-1 k= +1 exp - a 4 k-1 j=1 α j θ 0 -θ C 4 A ε ∞ √ na 5/2 2n-1 =n+1 α 3/2 + C 5 A √ na exp - a 4 n j=1 α j θ 0 -θ C 4 A ε ∞ n 1/4 a 5/2 + C 5 A √ na exp - a 4 n j=1 α j θ 0 -θ .Now the result follows from the representation (45) and combinations of the above bounds forD b 1 - D b 5 .Consider the sequences of weights (w 1 , . . . , w i-1 , w i , w i+1 , . . . , w 2n ) and (w 1 , . . . , w i-1 , w i , w i+1 , . . . , w 2n ) ,which differs only in position i, n + 1 ≤ i ≤ 2n, with w i being an independent copy of w i . Consider the associated SA processesEQUATIONwhere w(i) k = w k for k = i and w (i) i = w i .Respective random variables D b and D (b,i) are based on the first and second sequences from (57) , respectively, and are constructed according to the equation (44) . From the above representations we easily observe thatθ b k = θ (b,i) k for k < i, moreover, θ b i -θ (b,i) i = -α i (w i -w i ) A(Z i ))θ b i-1 -b(Z i ) = -α i (w i -w i ) A(Z i ))(θ b i-1 -θ i-1 ) -b(Z i ) = -α i (w i -w i ) A(Z i ))(θ b i-1 -θ i-1 ) + εi .whereε i = ε(Z i ) and ε i = ε(Z i ).From the above representation we get, applying Lemma 9 and (55), thatEQUATIONα i C A ε ∞ + α 3/2 i ε ∞ C 2 A a 3/2 + (C 2 A ∨1) √ i i j=1 (1 -aα j /4) θ 0 -θ α i ε ∞ C 2 A a + (C 2 A ∨1) √ i i j=1 (1 -aα j /4) θ 0 -θ ,where for the last line we have additionally assumed that α i a for i ≥ n. Moreover, for any j > i one observes, expanding (58) , thatθ b j -θ (b,i) j = j k=i+1 (I -α k w k A(Z k )) (θ b i -θ (b,i) i ) = Γ b i+1:j (θ b i -θ (b,i) i) .Thus, similarly to (37) , we obtain that{E b [ D b -D (b,i) 2 ]} 1/2 ≤ 5 j=1 {E b [ D b j -D (b,i) j 2 ]} 1/2 ,and bound the respective differences separately. By the construction of the process above, we note thatD b 1 = D (b,i) 1. Proceeding further, and using the equation ( 44), we obtain that{E b [ D b 2 -D (b,i) 2 2 ]} 1/2 = 1 √ nα 2n E b [ θ b 2n -θ (b,i) 2n 2 ] 1/2 ≤ 1 √ nα 2n E b [ Γ b i+1:2n 2 ] 1/2 E b [ θ b i -θ (b,i) i 2 ] 1/2 α i ε ∞ C 2 A √ nα 2n a exp - a 4 2n j=i+1 α j + (C 2 A ∨1) √ i √ nα 2n exp - a 4 2n j=1 α j θ 0 -θ .Thus, taking sum for i from n + 1 to 2n, and applying Lemma 2, we get thatEQUATION)Here in the last line above we used a particular form α k = c 0 / √ k, and relied on the boundn 3/4 exp - a 4 2n j=n+1 α j ≤ 1 ,which is guaranteed by the lower bound on the trajectory length n of the form√ n log n ≥ 3 2( √ 2 -1)ac 0 .The latter condition is guaranteed by A4. Now we proceed withD b 3 -D (b,i) 3. Using its definition in (44), we getD b 3 -D (b,i) 3 = 1 √ n (w i -w i )A i (θ b i-1 -θ ) + 1 √ n 2n k=i+1 (w k -1)A k (θ b k-1 -θ (b,i) k-1 ) .Since the latter term is a martingale-difference, we obtain thatE b [ D b 3 -D (b,i) 3 2 ] C 2 A n E b [ θ b i-1 -θ 2 ] + C 2 A n 2n k=i+1 E b [ θ b k-1 -θ (b,i) k-1 2 ] C 2 A n E b [ θ b i-1 -θ i-1 2 ] + C 2 A n θ i-1 -θ 2 + C 2 A n 2n k=i+1 E b [ θ b k-1 -θ (b,i) k-1 2 ] C 2 A n E b [ θ b i-1 -θ i-1 2 ] + C 2 A n θ i-1 -θ 2 + C 2 A n 2n k=i+1 E b [ Γ b i+1:k-1 2 θ b i -θ (b,i) i 2 ] .Hence we obtain, using (59) together with the definition of Ω 1 , thatEQUATIONConsidering the latter term in the sum, we obtainT 1 α 2 i ε 2 ∞ C 6 A nα 2n a 2 2n k=i+1 α k exp - a 2 k-1 j=i+1 α j + (C 6 A ∨1)i n 2n k=i+1 exp - a 2 k-1 j=1 α j θ 0 -θ 2 α i ε 2 ∞ C 6 A na 3 + (C 6 A ∨1)i nα 2n exp - a 2 i j=1 α j θ 0 -θ 2 .Thus, summing the equations ( 61) for i from n + 1 to 2n, we obtain thatEQUATIONUsing now the definition of D b 4 in (44) and Minkowski's inequality, we write{E b [ D b 4 -D (b,i) 4 2 ]} 1/2 ≤ 1 √ n {E b [ θ b i -θ (b,i) i 2 ]} 1/2 2n k=i+1 1 α k - 1 α k-1 exp - a 4 k-1 j=i+1 α j α i ε ∞ C 2 A √ na 2n k=i+1 α k exp - a 4 k-1 j=i+1 α j + (C 2 A ∨1) √ i √ n 2n k=i+1 α k exp - a 4 k-1 j=1 α j θ 0 -θ α i ε ∞ C 2 A √ na 2 + (C 2 A ∨1) √ i √ na exp - a 4 i j=1 α j θ 0 -θ .Thus, taking sum for i from n + 1 to 2n, we getEQUATIONSimilarly, with the definition of D b 5 in (44), we writeD b 5 -D (b,i) 5 = 1 √ n 2n k=i+1 (A k -Ā)(θ b k-1 -θ (b,i) k-1 ) = 1 √ n 2n k=i+1 (A k -Ā)Γ b i+1:k-1 (θ b i -θ (b,i) i ) = 1 √ n 2n k=i+1 (A k -Ā)Γ i+1:k-1 (θ b i -θ (b,i) i ) T2 + 1 √ n 2n k=i+1 (A k -Ā)(Γ b i+1:k-1 -Γ i+1:k-1 ) (θ b i -θ (b,i) i ) T3.Now we bound the terms T 2 and T 3 separately. Indeed, for the term T 2 we get, applying the definition of the set Ω 4 , thatE b [ T 2 2 ] 1 n C 2 A log n aα i + C 2 A log 2 n E b [ θ b i -θ (b,i) i 2 ] C 2 A log n naα i E b [ θ b i -θ (b,i) i 2 ] .In the above bounds we have used that α ≤ 1 a log n . For the term T 3 we get, applying Lemma 10, that for any vectorv ∈ R d , 2n k=i+1 (A k -Ā)(Γ b i+1:k-1 -Γ i+1:k-1 )v = 2n k=i+1 k-1 =i+1 (A k -Ā)Γ +1:k-1 α (w -1)A Γ b i+1: -1 v = 2n-1 =i+1 α (w -1) 2n k= +1 (A k -Ā)Γ +1:k-1 A Γ b i+1: -1 v .From the above representation we obtain, using the definition of the set Ω 4 , thatE b [ T 3 2 ] C 2 A n 2n-1 =i+1 α 2 C 2 A log n aα + C 2 A log 2 n exp - a 4 -1 j=i+1 α j E b [ θ b i -θ (b,i) i 2 ] C 4 A n 2n-1 =i+1 α log n a exp - a 4 -1 j=i+1 α j E b [ θ b i -θ (b,i) i 2 ] C 4 A log n na 2 E b [ θ b i -θ (b,i) i 2 ] .Combining the above bounds, we obtain thatE b [ D b 5 -D (b,i) 5 2 ] C 2 A log n na 1 α i + C 2 A a E b [ θ b i -θ (b,i) i 2 ] C 2 A log n naα i E b [ θ b i -θ (b,i) i 2 ] ,where we have additionally used that α i ≤ a/ C 2 A . Thus, using the upper bound (59), we obtain thatEQUATIONNow it remains to combine the bounds outlined above in ( 60), ( 62), (63) , and ( 64), and the statement follows.The fact that there exists a unique matrix Q, such that the following Lyapunov equation holds:EQUATIONfollows directly from [54, Lemma 9.1, p. 140]. In order to show the second part of the statement, we note that for any non-zero vector x ∈ R d , we havex (I -α Ā) Q(I -α Ā)x x Qx = 1 -α x ( Ā Q + Q Ā)x x Qx + α 2 x Ā Q Āx x Qx = 1 -α x P x x Qx + α 2 x Ā Q Āx x Qx ≤ 1 -α λ min (P ) Q + α 2 Ā 2 Q λ min (Q) ≤ 1 -αa ,where we set a = 1 2λ min (P ) λ max (Q) ,and used the fact that α ≤ α ∞ , where α ∞ is defined in (7) .Our analysis of the bootstrap procedure and the last iterate error of LSA procedure is based on the error expansion technique from [1] , see also [16] . Namely, to perform the expansion, we decompose the LSA iterates θ n defined in (1) into a transient and fluctuation terms:θ k -θ = θ(tr) n + θ(fl)n , where we have defined the quantitiesEQUATIONsettingEQUATION)The dependence of Γ m:k upon the stepsizes (α j ) is implicit in (67) . Here the quantity θ(tr) k is the transient component of the error, which determines the rate at which the initial error θ 0 -θ is forgotten. The term θ(fl) k corresponds to the fluctuation component of the error and is determined by the oscillations of the last iterate θ k around θ .In order to bound the moment E[ θ k -θ p ], we first prove a stability results on the products of random matrices Γ m:k arising in the LSA recursion. Towards this aim we first introduce some notations and definitions. For a matrix B ∈ R d×d we denote by (σ (B)) d =1 its singular values. For q ≥ 1, the Shatten q-norm of B is denoted by B q = { d =1 σ q (B)} 1/q . For q, p ≥ 1 and a random matrix X we write X q,p = {E[ X p q ]} 1/p . Our proof technique is based on the stability results arising in [29] , see also [16] . Lemma 12 (Proposition 15 in [16] ). Let {Y } ∈N be an independent sequence and P be a positive definite matrix. Assume that for each ∈ N there exist m ∈ (0, 1) and σ > 0such that E[Y ] 2 P ≤ 1 -m and Y -E[Y ] P ≤ σ almost surely. Define Z k = k =0 Y = Y k Z k-1, for k ≥ 1 and starting from Z 0 . Then, for any 2 ≤ q ≤ p and k ≥ 1,EQUATIONwhere we recall that κ P = λ -1 min (P )λ max (P ).Now we aim to bound Γ m:k defined in (67) using Lemma 12. We identify the latter with k =m Y , where Y = I -α A , ≥ 1, and Y 0 = I. Applying the bound (8), we getE[Y ] 2 Q = I -α Ā 2 Q ≤ 1 -aα . Further, assumption A2 implies that almost surely, Y -E[Y ] Q = α A -Ā Q ≤ α √ κ Q C A = b Q α .Therefore, (68) holds with m = aα and σ = b Q α . As I p = d 1/p , we obtain the following corollary. Corollary 3. Assume A1 and A2. Then, for any α ∈ [0, α ∞ ], 2 ≤ q ≤ p, and 1 ≤ m ≤ k, it holdsE 1/q [ Γ m:k q ] ≤ Γ m:k p,q ≤ √ κ Q d 1/p k =m (1 -aα + (p -1)b 2 Q α 2 ) ,where α ∞ was defined in (7) ,and b Q = √ κ Q C A .Corollary 4. Assume A1, A2, and A3. Then for any 2 ≤ q ≤ log n, and any k ≥ n, 1 ≤ m ≤ k, it holds thatEQUATIONwhere α ∞ is defined in (7) . Moreover,EQUATIONProof. We first apply the result of Corollary 3. Indeed, for k ≥ n, and any 2 ≤ q ≤ p, it holds,setting b Q = √ κ Q C A , that E 1/q [ Γ m:k q ] ≤ √ κ Q d 1/p k =m (1 -aα + (p -1)b 2 Q α 2 ) ≤ √ κ Q d 1/p exp -a k =m α + (p -1)b 2 Q k =m α 2 .Note that, setting p = log n, and provided that n satisfies (9), we easily obtain that, for ≥ n/2,EQUATIONHence, for m ≥ n/2, we haveE 1/q [ Γ m:k q ] ≤ √ κ Q e exp - a 2 k =m α ,and the statement follows. Suppose now that m < n/2. In such a case we have, applying (71) , thatEQUATIONand we need to bound the first term in the product. We first consider α = c 0 -1/2 , and use the inequalitiesEQUATIONandEQUATIONThus, it is enough to satisfy the constraint(log n)b 2 Q c 2 0 (1 + log n -log m) ≤ ac 0 ( √ n - √ m) .Since m < n/2, it is enough to ensure that(1 + log n)(log n)b 2 Q c 2 0 ≤ ac 0 ( √ n -n/2) , or, equivalently, √ n (1 + log n) log n ≥ c 0 b 2 Q a(1 -1/ √ 2) ,which is granted by A3. Combining the above bounds in (72) , we obtain that the lemma's statement (69) holds for the step size α = c 0 / 1/2 . Similarly, for α = c 0 / γ with γ ∈ (1/2; 1), we get for m ≥ n/2 thatE 1/q [ Γ m:k q ] ≤ √ κ Q e exp - a 2 k =m α ,since the relation (71) holds. Similarly, for m < n/2, the desired upper bound would follow from the inequalityn =m 1 2γ ≤ n m-1 dx x 2γ = (m -1) 1-2γ -n 1-2γ 2γ -1 ≤ 1 2γ -1 ,together with an inequality(log n)b 2 Q c 2 0 2γ -1 ≤ (a/2)c 0 (n 1-γ -(n/2) 1-γ ) .The latter inequality can be re-written asn 1-γ log n ≥ 2c 0 b 2 Q a(2γ -1)(1 -(1/2) 1-γ ,which is also granted by A3. Combining the above inequalities implies that (69) holds for α = c 0 / γ . The bound (70) can be immediately obtained from (69) using the fact that e -x ≤ 1 -x/2 for x ∈ [0; 1].Corollary 5. Under conditions of Corollary 4 it holds with P -probability at least 1 -1/n 2 thatΓ m:k ≤ √ κ Q e 2 exp -(a/2) k =m α ,andEQUATIONProof. It is sufficient to choose q = 2 log n and use Markov's inequality together with the union bound.Proposition 7. Assume A1, A2, A3 with γ = 1/2, and A4. Then on the set Ω 5 defined in (43) , it holds for any n ≤ m ≤ k ≤ 2n, thatE b [ Γ b m+1:k 2 ] 1/2 ≤ κ 3/2 Q e 9/8 exp - a 4 k =m+1 α .Proof. Our proof relies on the auxiliary result of Lemma 13 below together with the blocking technique. Indeed, let us represent k -m = N h + r , where r < h and h = h(n) is a block size defined in (17) . Then we obtain, using the independence of bootstrap weights w m+1 , . . . , w k , that{E b [ Γ b m+1:k 2 ]} 1/2 ≤ √ κ Q {E b [ Γ b m+1:k 2 Q ]} 1/2 = √ κ Q N j=1 E b [ Γ b m+1+(j-1)h:m+jh 2 Q ] 1/2 E b [ Γ b m+1+N h:k 2 Q ] 1/2 ≤ √ κ Q exp - a 4 k =m+1 α E b [ Γ b m+1+N h:k 2 Q ] 1/2 exp a 4 k =m+1+N h:k α .In the last inequality we applied Lemma 13 to each of the blocks of length h in the first bound. It remains to upper bound the residual terms. Since the remainder block has length less then h, we have due to (79) (which holds according to A4), thatexp a 4 k =m+1+N h:k α ≤ exp α ∞ a 4 ≤ e 1/8 ,where the last inequality is due to Proposition 1. Next,E b [ Γ b m+1+N h:k 2 Q ] 1/2 ≤ κ Q k =m+1+N h:k {E b [ (I -α w A ) 2 ]} 1/2 ≤ κ Q k =m+1+N h:k {E b [(1 + α |w | C A ) 2 ]} 1/2 ≤ κ Q k =m+1+N h:k {E b [1 + 2α |w | C A +α 2 w 2 C 2 A ]} 1/2 . Since E[|w |] ≤ E[w 2 ] ≤ (E[w ]) 2 + Var w = √ 2 , we get from previous bound E b [ Γ b m+1+N h:k 2 Q ] 1/2 ≤ κ Q k =m+1+N h:k (1 + 2 √ 2α C A +2α 2 C 2 A ) 1/2 ≤ κ Q exp √ 2 C A k =m+1+N h:k α ≤ κ Q e √ 2 C A c0h/ √ n ≤ κ Q e ,where in the last line we additionally used (18) .Lemma 13. Assume A1, A2, A3 with γ = 1/2, and A4. On the set Ω 5 defined in (43) , it holds for h = h(n) defined in (17) and any m ∈ [n; 2n -h], thatE b [ Γ b m+1:m+h 2 Q ] 1/2 ≤ exp - a 4 m+h =m+1 α .Proof. Recall that we use the notationE b [•] = E[•|Z 2n ],where Z 2n = (Z 1 , . . . , Z 2n ) are the random variables used in the construction of the iterates {θ k } 1≤k≤n in (1) .Let h ∈ N be a block length, which value will be determined later, and consider a productEQUATIONExpanding the product of matrices (75) , we obtainEQUATIONwhere S = m+h =m+1 α (w -1)A is a linear statistics in {w } m+h =m+1 , and the remainder R collects the higher-order terms in the productsR = h r=2 (-1) r (i1,...,ir)∈I r r u=1 α iu w iu A iu . with I r = {(i 1 , . . . , i r ) ∈ {m + 1, . . . , m + h} r : i 1 < • • • < i r }.We first consider the contracting part in matrix Q-norm. Indeed, applying (8) , we obtain thatI - m+h =m+1 α Ā 2 Q ≤ 1 -a m+h =m+1 α ,provided that h is set in such a manner that m+h =m+1 α ≤ α ∞ , where α ∞ is defined in (7) . Hence, we get from the above inequality that for any u ∈ R d , it holds thatI - m+h =m+1 α Ā Q ≤ 1 -(a/2) m+h =m+1 α .Now we need to estimate the remainders in the representation (76). On the set Ω 5 , it holds thatm+h =m+1 α (A -Ā) Q ≤ 2 C A √ κ Q m+h =m+1 α 2 log(2n 4 ) .Moreover, it is straightforward to check thatE b [ S 2 Q ] ≤ C 2 A κ Q m+h =m+1 α 2 .In order to bound the remainder term R, we note thatE b [ R Q ] ≤ h r=2 h r α r m+1 (2 C A ) r κ r/2 Q ≤ α 2 m+1 (2 C A ) 2 κ Q h-2 r=0 h r + 2 α r m+1 (2 C A ) r κ r/2 Q ≤ α 2 m+1 h 2 (2 C A ) 2 κ Q 2 exp 2α m+1 C A κ 1/2 Q ≤ α 2 m+1 h 2 (2 C A ) 2 κ Q e 2 .To complete the proof it remains to set the parameter h in such a way that we can guaranteeEQUATIONkeeping at the same time the constraintEQUATIONRecall that α = c 0 / √ . Thus, using the bounds ( 73) and ( 74), we obtain thatEQUATIONandEQUATIONHence, taking into account (79) and (80), and 1 m+1 ≤ 1 m , the inequality (77) would follow from the boundEQUATIONSince log (1 + x) ≤ x for x ≥ 0 and c 0 C 2 A κ Q e ≤ 1, the latter inequality is satisfied ifC A √ κ Q √ h √ m 1 + 2 log(2n 4 ) + h 2 2m ≤ a 2 ( √ m + h - √ m) .√ m + h - √ m = √ m( 1 + h/m -1) ≥ √ m( √ 2 -1)h m = ( √ 2 -1)h √ m ,which follows from an elementary inequality√ 1 + x ≥ 1 + ( √ 2 -1)x, valid for 0 ≤ x ≤ 1. Hence, (81) would from the inequalityEQUATION) Setting h is such a manner that h √ m ≤ a( √ 2 -1) 2 ,inequality (82) would follow fromC A √ κ Q √ h √ m 1 + 2 log(2n 4 ) ≤ a( √ 2 -1)h 4 √ m .The latter inequality is satisfied, if the block size h satisfiesh ≥ 4 C A κ 1/2 Q ( √ 2 -1)a 2 (1 + 2 log (2n 4 )) 2 .Thus, setting h(n) as in (17) , all previous inequalities will be fulfilled, provided thath(n) √ n ≤ a( √ 2-1) 2 c0h(n) √ n ≤ α ∞ .Here last inequality follows from (78) and the following simple bounds, where we use that m ≥ n and √1 + x ≤ 1 + x/2: m+h =m+1 α ≤ n+h =n+1 α = c 0 n+h =n+1 1 √ ≤ c 0 n+h n dx √ x = 2c 0 ( √ n + h - √ n) ≤ c 0 h √ n . Now (76) implies that E b [ Γ b m+1:m+h 2 Q ] 1/2 ≤ 1 -(a/4) m+h =m+1 α ,and the statement follows from an elementary inequality 1 + x ≤ e x .Recall that the temporal difference learning algorithm in the LSA's setting can be written asEQUATIONwhere A k and b k are given byEQUATIONRecall that our aim is to estimate the agent's value function V π (s) = E[ ∞ k=0 k r(s k , a k )|s 0 = s] , where a k ∼ π(•|s k ),We first need to check that the matrix Ā + Ā , where Ā is defined in (86), is positive-definite. In order to show this fact we closely follow the exposition of [61, Lemma 18] and [51, Lemma 5] . Define a random matrix A as an independent copy of A k from (84), that is,A = ϕ(s){ϕ(s) -γϕ(s )} ,where s ∼ µ, and s ∼ P π (•|s). With the definition of A, we get thatA + A = ϕ(s){ϕ(s) -γϕ(s )} + {ϕ(s) -γϕ(s )}ϕ(s) = 2ϕ(s)ϕ(s) -γ{ϕ(s)ϕ(s ) + ϕ(s )ϕ(s) } (2 -γ)ϕ(s)ϕ(s) -γϕ(s )ϕ(s ) ,where we used an elementary inequality uv + vu (uu + vv ) valid for any u, v ∈ R d . Hence, with the definition of Σ ϕ in (21), we getEQUATIONHence, Ā + Ā is positive-definite, and we can set P = Ā + Ā in the right-hand side of the Lyapunov equation (65) . Obviously, Q = I is a solution to the corresponding Lyapunov equationĀ Q + Q Ā = Ā + Ā .Moreover, applying [61, Lemma 18] , we obtainEQUATIONHence, we get for α ≤ (1 -γ)/(1 + γ) 2 , and applying (87) and (88), that(I -α Ā) (I -α Ā) = I -α( Ā + Ā) + α 2 Ā Ā I -2α(1 -γ)Σ ϕ + α 2 (1 + γ) 2 Σ ϕ I -α(1 -γ)Σ ϕ (1 -α(1 -γ)λ min (Σ ϕ ))I .Hence, the bound (8) holds with a = (1 -γ)λ min (Σ ϕ ) and α ∞ = (1 -γ)/(1 + γ) 2 .Here we provide some details on numerical experiments. Code to run experiments is provided in https://github.com/svsamsonov/BootstrapLSA. For the considered Garnet problem we choose the policy π in the following way. For any a ∈ A, we setπ(a|s) = U (s) a |A| i=1 U (s) i ,where the U (s) i are independent random variables following uniform distribution U[0, 1]. Here we assume that each action a ∈ A can be selected at any state s ∈ {1, . . . , N s }. We generate an instance of Garnet problem with mentioned parameters, and find analytically the true parameter θ . In order to estimate the supremum∆ n := sup x∈R P( √ n θn -θ ≤ x) -P( Σ 1/2 ∞ η ≤ x) ,η ∼ N (0, I Ns ), and show that this supremum scales as n -1/4 when γ = 1/2 and admits slower decay for other powers of γ. We first approximate true probability P( Σ 1/2 ∞ η ≤ x) by the corresponding empirical probabilities based on sample of size M n. We fix M = 2 • 10 7 . We choose trajectory lengths n ∈ {1600, 3200, 6400, 12800, 25600, 51200, 102400, 204800, 409600, 819200, 1638400} , fix the length of burn-in period n 0 = 102400, and generate N = 1638400 independent trajectories starting in the fixed point θ 0 ∈ R Ns . We set the learning rate schedule as α k = c 0 /k γ and try different values γ ∈ {0.5, 0.65, 0.7, 0.75}. In Figure 2 we provide respective plots for the rescaled RMSE √ n θn -θ and its non-rescaled version θn -θ , averaged over N trajectories. Then in Figure 3 we provide the expressions for ∆ n and ∆ n n 1/4 , both in log-scale for y axis and without log-scaling. Unfortunately, even the chosen order of trajectory length n seems to be insufficient in order to significantly distinguish, for example, between γ = 0.5 and γ = 0.65. However, learning rate schedule with faster decay performs consistently worse in terms of ∆ n . The order of lines in Figure 3-(b ) is remarkable. However, note that the current experiment is already rather computationally intense for artificial problem and takes about 12 hours of compute on a Core i9 -10920x processor with 12 cores with 3.7 GHz (for a bunch of N = 1638400 trajectories of length n = 1638400). Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Answer: [Yes] Justification: All our theoretical results are provided with references to assumptions, that are stated in Section 3 and Section 4. All results are given with proofs, that are correctly referenced for each theorem and corollary. Guidelines:• The answer NA means that the paper does not include theoretical results.• All the theorems, formulas, and proofs in the paper should be numbered and crossreferenced. • All assumptions should be clearly stated or referenced in the statement of any theorems.• The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition. • Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material. • Theorems and Lemmas that the proof relies upon should be properly referenced.Justification: This paper is of purely theoretical nature, and the proposed methods do not deal with sensitive attributes that could induce unfairness or privacy issues.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 


Decision: The present paper proves normal approximation (CLT) and validity of a bootstrap procedure for the law of the Polyak-Ruppert averaged iterates of linear stochastic approximation (LSA). Their main results have n^{-1/4} convergence rates (up to polylogs and dimension-dependent terms) under reasonable assumptions. Intriguingly, the authors suggest that the magnitude of the error in their CLT might be optimal. 

The paper was generally well-received by the reviewers, with scores 5, 6 and 7. The issues raised in the revision period do not seem significant, and can be resolved in the final version. In that connection, I would encourage the authors to add their remarks on weaker assumptions to the text (cf. the response to reviewer uNZU). 

Overall, I believe the paper deserves to be accepted: the quantitative CLT is quite interesting, and the bootstrap results seem to be new.

Reviews: review: paper_summary: paper_summary: ## Overview

Let $Z, Z_1, \dots, Z_n$ be i.i.d. random elements with a common distribution $\pi$ over $\mathbf{Z}$. Given $A : \mathbf{Z} \to \mathbb{R}^{d\times d}$ and $b : \mathbf{Z} \to \mathbb{R}^d$, the goal of the LSA procedure is to find the unique solution $\theta^\star$ of
$$\mathbb{E}\left(A(Z)\theta^\star - b(Z)\right) = 0. $$

Given a decreasing sequence of step sizes $\alpha_k$ and a starting point $\theta_0$, the standard LSA is given by
$$\theta_k = \theta_{k-1} - \alpha_k(A(Z_k)\theta_{k-1} - b(Z_k)) $$
and the Polyak-Ruppert averaged LSA is given by
$$\overline{\theta_n} = \frac{1}{n} \sum_{k=n}^{2n-1} \theta_k. $$

The authors provide Berry-Essen-type bounds for the Gaussian approximation of $\sqrt{n} \left( \overline{\theta_n} - \theta^\star \right)$ and for the corresponding multiplicative bootstrap process. Namely, for the Gaussian approximation result they upper bound the quantity
$$\rho_n = \sup_\text{B convex} \left| \mathbb{P}\left( \sqrt{n} \left( \overline{\theta_n} - \theta^\star \right) \in B \right) - \mathbb{P}\left( \Sigma_\infty^\frac{1}{2} \eta \in B \right)\right| $$
where $\eta \sim \mathcal{N}(0,I_d)$. And for the bootstrap approximation result they upper bound
$$\rho_n^b = \sup_\text{B convex} \left| \mathbb{P}\left( \left.\sqrt{n} \left( \overline{\theta_n^b} - \overline{\theta_n} \right) \in B \right| Z_1, \dots, Z_{2n} \right) - \mathbb{P}\left( \sqrt{n} \left( \overline{\theta_n} - \theta^\star \right) \in B \right)\right| $$
where $\overline{\theta_n^b}$ are obtained from a multiplier bootstrap process. It is interesting to notice that this multiplier process can be evaluated online without keeping a history in memory.

## Overall proof arguments

To provide the Gaussian approximation result the authors write
$$\theta_n - \theta^\star = (I - \alpha_n A(Z_n))(\theta_{n-1} - \theta^\star) - \alpha_n \varepsilon(Z_n) $$
where
$$\varepsilon(z) = A(z)\theta^\star - b(z). $$

The authors observe that $\varepsilon(Z_n)$ can be viewed as a noise, which is assumed to be bounded. Meanwhile, the operator $ I - \alpha_n A(Z_n) $ is a random perturbation around $ I - \alpha_n \mathbb{E}A(Z) $, which is shown to act as a contraction in an appropriate norm, provided that $-\mathbb{E}A(Z)$ is Hurwitz. Thus, when one take the PR average the noise is expected to behave as the sum of i.i.d. noises and the contraction term must shrink. This is formally done writing the PR average in the form of Theorem 2.1 [reference 60 of the paper] and bounding the terms given by the latter theorem.

The proof the bootstrap approximation result follows the standard practice in this literature:
- First, conditionally on the sample, a Gaussian approximation is obtained relating the bootstrap to a Gaussian with its covariance.
- Second, a Gaussian comparison theorem relates the latter Gaussian with the desired Gaussian, with some concentration results being used to bound their difference in high probability.

## Main claims and observations

The authors show that taking a step size of order $\frac{1}{\sqrt{k}}$ yields the best possible convergence rate in their bounds. They provide empirical evidence that this rate is optimal.
They also claim to be the first ones to fully provide a non-asymptotic bootstrap approximation result.
strength_weakness: strength_weakness: The main claim of the paper is that their bounds suggest an optimal convergence rate of $n^{-\frac{1}{4}}$ when taking $\alpha_k = \frac{c}{\sqrt{k}}$. The bottleneck of this convergence rate comes from an application of Cauchy-Schwarz inequality, the boundness of $\varepsilon$ and the MSE bound on $D$ given by Theorem 1. Meanwhile, the authors also provide experimental evidence that this convergence rate is indeed optimal (although not considering gamma values below 0.5 in Figures 1 or 3). A counter-example or a deeper discussion on this convergence rate would be beneficial to the paper: is it an artifact of the proof or is there reason to believe it cannot be improved?
questions: questions: Some observations:
- Line 54 has a typo in "Berry-Essee".
- In Equation (15) the $\ell$ can be removed from $\theta_{k}^{b,\ell}$ to enhance clarity. It is only used to explain that to evaluate the probability in practice one must run several samples of the bootstrap process.
- Line 263 has a typo in "date".
- The Equation after line 303 is using $\phi$ instead of $\varphi$.
- Line 307 asks for a sequence of "TD(0) updates", what does the "(0)" stands for?
- Figure 1 lacks x and y labels. The legend can also be improved for clarity.
limitations: limitations: The authors adequately addressed the limitations and the potential impact of their work.
paper_summary: paper_summary: The present paper studies linear stochastic approximation with martingale difference noise and diminishing step-sizes. The authors obtain Berry-Esseen bounds for the parameter sequence with Polyak-Ruppert averaging as well as a generalization of finite-time bounds for estimation confidence intervals for parameters in LSA. The obtained Berry-Esseen bounds are illustrated by a TD learning numerical example.
strength_weakness: strength_weakness: Although the authors did a good job in outlining the scope and contributions of the paper, the analysis and main text are hard to follow given the number of symbols and equations. It is easy for a reader to get lost/distracted midway and it is very hard to keep track of the definitions of each of the terms 

I understand that this is an issue with theory papers like this, but would encourage the authors to move unnecessary terms or inequalities to the appendix (e.g. (9) in A3. The exact lower bound adds very little in the main text in my opinion. Its definition could have been postponed to the Appendix.)

The numerical experiments are also weak since many of the plots are not related to the contributions of the paper itself. Plot c seems to be the only one directly related to the theorems of the paper, but it still does not illustrate the theory that well. Maybe including a plot of C k^{-1/4}  to plot (b) where C is a constant for comparison could help in inferring convergence rates.

Also, there is a mistake in the label of Figure 1. Subfigure (b) is mentioned twice.
questions: questions: - Could the authors clarify if using Polyak Ruppert averaging is necessary to obtain such a Berry-Esseen bound? It would be exciting to see bounds for unaveraged estimates as well.

-Could the authors  run  the experiment  supporting figure (c) for longer to see if the curve with \gamma = 1/2 willnot continue to go upwards eventually?
limitations: limitations: The authors clearly identified the limitations of their results through a clear list of assumptions
paper_summary: paper_summary: The paper presents advancements in the theoretical understanding of the linear stochastic approximation (LSA) algorithm. 

It establishes the Berry–Esseen bound for the normal approximation of Polyak-Ruppert averaged iterates, achieving an optimal rate with an aggressive step size of $\alpha_k \approx k^{-1/2}$. Additionally, it demonstrates the non-asymptotic validity of confidence intervals using a novel multiplier bootstrap procedure, marking a first in this domain. 

The practical utility of these theoretical results is showcased through applications in temporal difference (TD) learning for reinforcement learning.
strength_weakness: strength_weakness: 1. **Strong Assumptions**: The requirement that $\epsilon(z)$ is uniformly bounded is a strong assumption that might limit the applicability of the results. Relaxing this condition to weaker ones, such as finite moments, could enhance the paper's generality.

2. **Missing References**: There are some missing references, see the limitations.

3. **Discussion on Lower Bounds**: The paper focuses on deriving upper bounds but lacks a discussion on the tightness and potential lower bounds. Addressing this aspect could provide a more comprehensive understanding of the bounds' efficacy and limitations.
questions: questions: See the Limitations.
limitations: limitations: 1. **Missing References on Statistical Inference**: Statistical inference for nonlinear stochastic approximation has been considered by [1*] and [2*], which is highly related to this manuscript but hasn’t been cited. Note that [1*] provides a Berry-Esseen-like bound for the whole trajectory rather than the averaged iterates. These references could be cited after [36] in line 116.
     - [1*] Li, Xiang, Jiadong Liang, and Zhihua Zhang. "Online statistical inference for nonlinear stochastic approximation with Markovian data." arXiv preprint arXiv:2302.07690 (2023).
     - [2*] Li, Xiang, et al. "A statistical analysis of Polyak-Ruppert averaged Q-learning." International Conference on Artificial Intelligence and Statistics. PMLR, 2023.

2. **Assumption on Bounded $\epsilon(z)$**: Theorem 1 requires that $\epsilon(z)$ is uniformly bounded, which is a strong condition. Is it possible to relax this condition to a weaker one, such as $\epsilon(z)$ only having a finite order of moments (such as the fourth order moment or smaller)?

3. **Tightness of Derived Upper Bounds**: The paper provides upper bounds, but it would be insightful to discuss the tightness of these bounds. Are there any thoughts or conjectures regarding potential lower bounds?

#### Minor Corrections:

1. **Figure 1**: The last subfigure should be labeled as (c).

Comments: comment: Dear referee,

Please kindly let us know if you have any follow-up questions that need further clarification. Your insights are valuable to us, and we stand ready to provide any additional information that might be helpful.
We thank the referee for their comments. We will include a discussion on the Berry-Esseen result for last iterate as well as a longer run for simulations.
Yes, we will include the corresponding discussion on the optimal rates in martingale CLT in the revised version of the text, as well as a comment on TD(0). Fetching the example provided in [Bolthausen, 1982] into the LSA paradigm is not immediate, but we work in this direction, and, of course, if we succeed to construct such a lower bound, we will include it in the final text.
I appreciate the author's responses. 

I believe that it would be beneficial to include some discussion/remarks on the final version about the extension of their results for estimates without PR as in the response provided by the authors. 

I apologize for not taking enough time to fully grasp the experiments in the paper, but I understand them now. Thank you for providing a longer run.
I thank the authors for their rebuttal.
The comment made in the rebuttal on the optimal $n^{1/4}$-rate obtained in [Bolthausen, 1982] for the martingale Berry-Essen should be on the main text. I wonder if the authors can find an example matching their upper bound (thus showing it is optimal), maybe drawing inspiration from the example in Section 6 of [Bolthausen, 1982].
I thank the authors for the explanation of the meaning of TD(0), it should also be included on the revised text.