The Collusion of Memory and Nonlinearity in Stochastic Approximation With Constant Stepsize

Authors: dongyan lucy huo
yixuan  zhang
yudong  chen
qiaomin  xie
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 00:08:56

Abstract: In this work, we investigate stochastic approximation (SA) with Markovian data and nonlinear updates under constant stepsize $\alpha>0$. Existing work has primarily focused on either i.i.d. data or linear update rules. We take a new perspective and carefully examine the simultaneous presence of Markovian dependency of data and nonlinear update rules, delineating how the interplay between these two structures leads to complications that are not captured by prior techniques. By leveraging the smoothness and recurrence properties of the SA updates, we develop a fine-grained analysis of the correlation between the SA iterates $\theta_k$ and Markovian data $x_k$. This enables us to overcome the obstacles in existing analysis and establish for the first time the weak convergence of the joint process $(x_k, \theta_k)$. Furthermore, we present a precise characterization of the asymptotic bias of the SA iterates, given by $\mathbb{E}[\theta_\infty]-\theta^\ast=\alpha(b_\textup{m}+b_\textup{n}+b_\textup{c})+\mathcal{O}(\alpha^{3/2})$. Here, $b_\textup{m}$ is associated with the Markovian noise, $b_\textup{n}$ is tied to the nonlinearity of the SA operator, and notably, $b_\textup{c}$ represents a multiplicative interaction between the Markovian noise and the nonlinearity of the operator, which is absent in previous works. As a by-product of our analysis, we derive finite-time bounds on higher moment $\mathbb{E}[||\theta_k-\theta^\ast||^{2p}]$ and present non-asymptotic geometric convergence rates for the iterates, along with a Central Limit Theorem.

Main paper: 
 introduction: 
 Stochastic Approximation (SA) is an iterative scheme for solving fixed-point equations using noisy observations. Its application spans various domains including stochastic control [9, 40] , reinforcement learning (RL) [4, 62] and stochastic optimization [43] . A typical SA algorithm takes the form θ k+1 = θ k +αg(θ k , x k ), where (x k ) k≥0 represents the underlying noisy data sequence and α > 0 is a constant stepsize. The goal of SA is to approximate the target solution θ * that solves E x∼π [g(θ * , x)] = 0, with π being the stationary distribution of the stochastic process (x k ) k≥0 .SA subsumes many important algorithms. A prime example is stochastic gradient descent (SGD) for minimizing a function J(θ) given noisy estimates g(θ, x) of its gradient. Linear SA schemes include SGD for quadratic objective functions, as well as various RL algorithms such as linear TD-Learning (in which g is not the gradient of any function and standard SGD results do not apply).Of particular interest to us are SA updates given by a nonlinear function g(θ, x) of θ. One motivating example is learning a Generalized Linear Model (GLM) y ≈ σ(z ⊤ θ) with a nonlinear mean function σ : R → R. A power approach, developed in [19, 36, 37, 64] , uses a surrogate loss function, where the corresponding SGD update takes the form θ k+1 = θ k + α(σ(w ⊤ k θ k ) -y k )w k , where x k = (w k , y k ) is the observed covariate-response pair. Common choices of σ include the identity map for linear regression, the sigmoid function for logistic regression, as well as Rectified Linear Unit (ReLU) and its various smoothed versions (e.g., ELU and SoftPlus) for ReLU regression [7, 18, 19, 30, 36] . Furthermore, we are interested in the setting where the data sequence (x k ) k≥0 forms a Markov chain, going beyond the common i.i.d. data setting. The Markovian model captures a wide range of SA problems in machine learning where stochastic data exhibit serial dependence [5, 33, 38, 53] .Classical work on SA focuses on diminishing stepsizes [10, 58] . Constant stepsize schemes have recently gained popularity due to easy parameter tuning, fast initial convergence, and robust empirical performance. Non-asymptotic error bounds have been obtained for constant stepsize SA [17, 60] .Recent work further provides fine-grained characterization of the distributional and steady-state behaviors of the iterates [20, 33, 45, 67, 69] . Two recurring themes in these results are weak convergence of the distribution of θ t and the presence of an asymptotic bias E[θ ∞ ] -θ * ∝ α, both having important implications for iterate averaging, bias reduction and statistical inference [34] .Note that most previous work studied the nonlinear update setting and Markovian data setting separately-e.g., in [20, 67] for nonlinear SGD with i.i.d. data, and in [33, 34] for Markovian linear SA. The linearity or i.i.d. assumptions imposed in these prior works are restrictive, especially in the face of modern machine/reinforcement learning paradigms where nonlinear models are the norm and dependent data is common. Moreover, the absence of prior work dealing with Markovian nonlinear SA is not merely an overlook-as argued below, this setting is significantly more challenging.Our Contributions In this work, we study constant-stepsize SA with both Markovian data and nonlinear update. In Section 3, we elucidate the new challenges that arise from the simultaneous presence of these two structures, which break key steps in previous analyses of the i.i.d. or linear setting. Due to the interaction between these two structures, establishing weak convergence is far from obvious, and the asymptotic bias exhibits new behaviors. Consequently, analyzing the nonlinear Markovian setting requires more than simply combining previous techniques.To address the above confounding complication, we exploit the smoothness and recurrence structures of the SA update, thereby developing a fine-grained analysis of the correlation of the parameter θ k and data x k . This allows us to establish for the first time the weak convergence of the joint process (x k , θ k ) k≥0 to a unique invariant distribution, represented by the limiting random variable (x ∞ , θ ∞ ). As a by-product of our analysis, we derive finite-time bounds on E[∥θ k -θ * ∥ 2p ], the 2p-th moments of the errors, generalizing the results in [17, 20, 60] to higher moments and to the nonlinear Markovian setting. In addition, we prove a Central Limit Theorem (CLT) for averaged iterates.Moreover, we show that nonlinearity and Markovian structure contribute in a multiplicative way to the asymptotic bias of the SA iterates. We obtain the following bias characterization: E[θ Importantly, we identify the compound term b c , which is absent in both nonlinear SA with i.i.d. data and linear SA with Markovian data. We explore the algorithmic implications of the above results on Polyak-Ruppert (PR) averaging [35, 56, 59 ] and Richardson-Romberg (RR) extrapolation [31] . We show that PR averaging reduces the variance but not the bias, whereas RR extrapolation eliminates the leading bias term α(b m + b n + b c ), reducing the asymptotic bias to a higher order of α.Postponing a detailed literature review to Section 6, here we remark on the recent works most relevant to ours [1, 44, 45, 46] , all studying Markovian nonlinear SA. The authors of [1] present an upper bound for the PR-averaged iterates and, similar to our results, demonstrate the effectiveness of RR-extrapolation in reducing bias, but lack weak convergence result for last iterates. In [44] , the authors suggest adopting the ordinary differential equation framework to prove weak convergence of iterates and derive an upper bound for the asymptotic bias, which contrasts with our equality characterization with a closed-form solution for the leading-order bias. In [45] , the authors prove weak convergence of (x t , θ t ) using coupling, but only in the linear setting, not for nonlinear SA. In the latter setting, their weak convergence analysis is thwarted by challenges similar to what we elucidate in Section 3, the interplay between nonlinearity and Markovian data leading to "double recursions." The coupling technique differs as well: we couple two processes by sharing data x t = x ′ t , while in [45] , they initialize two processes with different x 0 and x ′ 0 , and analyze the stopping time τ when θ τ = θ ′ τ . Moreover, they only present an upper bound for asymptotic bias, while ours presents a fine-grained characterization in Theorem 4.6 necessary for justifying RR-extrapolation. Lastly, the paper [46] discusses stepsize selection and its impact on the asymptotic statistics of PR-averaged SA with both constant and diminishing stepsizes.Notations The Euclidean norm is denoted by ∥ • ∥. The notation "u ⊗ v" represents the tensor product between vectors u and v, and "u ⊗k " denotes the k-th tensor power of vector u. The ball with radius β is B(β) := {θ ∈ R d : ∥θ∥ ≤ β}. L(z) denotes the distribution of a random vector z and Var(z) its covariance matrix. Let P 2 (R d ) be the space of square-integrable distributions on R d and P 2 (X × R d ) be the space of distributions ν on X × R d with square-integrable second marginal on R d . The Wasserstein-2 between two probability measures µ and ν in P 2 (R d ) is defined as W 2 (µ, ν) = inf ψ∈Π(µ,ν) E[∥θ -θ ′ ∥ 2 ] 1 2 : L(θ) = µ, L(θ ′ ) = ν , where Π(µ, ν) is the set of all couplings between µ and ν. Extending to X × R d , we define the metric d (x, θ), (x ′ , θ ′ ) := 1{x ̸ = x ′ } + ∥θ -θ ′ ∥ 2 , and denote by W2 the extended Wasserstein-2 distance w.r.t. d.The lowercase letter c and its derivatives c ′ , c 0 , etc. denote universal numerical constants, whose value may change from line to line. We use s ≡ s(θ 0 , θ * , µ, L, R) and its derivatives to denote quantities (scalars, vectors, or matrices) that are independent of the stepsize α and the iteration index k, but may depend on the initialization θ 0 , SA primitives θ * , µ and L, and the coefficient R for the geometric mixing rate of (x k ) in Assumption 1. As we are primarily interested in dependence on α and k, we adopt the following big-O notation: ∥f ∥ = O(h(α, k)) if it holds that ∥f ∥ ≤ s • ∥h(α, k)∥. 
problem setup and preliminaries: 
 Let (x k ) k≥0 be a Markov chain on a general state space X . Consider the following projected stochastic approximation (SA) iteration:EQUATIONwhere g : R d × X → R d is a deterministic function, {ξ k } k≥1 are i.i.d. zero-mean random fields, α > 0 is a constant stepsize, and Π B(β) (θ) := arg min z:∥z∥≤β ∥z -θ∥ is the projection operator. We shall omit the superscript (α) in θ k when the dependence on α is clear from the context. In this work, we also consider the projection-free variant of the iteration (2.1) with β = ∞.We denote by π the stationary distribution of the Markov chain (x k ) k≥1 and define the shorthand ḡ(θ) := E π [g(θ, x)], where E π [•] denotes the expectation with respect to x ∼ π. The algorithm (2.1) computes an estimation of the target vector θ * that solves the steady-state equation E π [g(θ, x)] = 0. Our general goal is to characterize the relationship between the iterate θ k and the target solution θ * .In the following, we state the assumptions needed for our main results. For a more detailed discussion of the assumptions, we refer readers to Appendix B. Assumption 1 (Uniform Ergodicity). (x k ) k≥0 is a uniformly ergodic Markov chain on a countable state space (X , B(X )) with transition kernel P and a unique stationary distribution π. That is, there exist constants r ∈ [0, 1) and R > 0 such that ∥P k (x, •) -π∥ TV ≤ Rr k , ∀x ∈ X .The countable state space ensures separability under the 1{x ̸ = x ′ } metric, necessary for constructing a valid coupling in the invariance proof and establishing a well-defined P * for bias characterization.We keep the notation general to allow future extensions to general state space and broader applicability of our results. All irreducible, aperiodic, and finite state space Markov chains are uniformly ergodic.The uniform ergodicity assumption is common in prior work on SA with Markovian noise [6, 21, 24, 33, 48] . Relaxing this uniform ergodicity assumption, in the style of [45, 52, 60] is possible but orthogonal to our focus, and thus we do not pursue this direction in this work.Assumption 2 (Differentiability and Linear Growth). For each x ∈ X , the function g(θ, x) is three times continuously differentiable in θ with uniformly bounded first to third derivatives, i.e., sup θ∈R d ∥g (i) (θ, x)∥ < +∞ for i = 1, 2, 3, x ∈ X . Moreover, there exists a constant L 1 > 0 such that (1)∥g (i) (θ, x) -g (i) (θ ′ , x)∥ ≤ L 1 , for all θ, θ ′ ∈ R d , i = 0, 1, 2 and x ∈ X , and (2) ∥g(0, x)∥ ≤ L 1 for all x ∈ X .The linear growth condition in Assumption 2 implies that g(θ, x) is L 1 -Lipschitz w.r.t. θ uniformly in x. When g is a linear function, i.e., g(θ, x) = A(x)θ + b(x), this assumption is satisfied with sup x∈X ∥A(x)∥ < ∞ and sup x∈X ∥b(x)∥ < ∞, which are commonly assumed for linear SA. The above assumption immediately implies that the growth rate of ∥g∥ and ∥ḡ∥ will be at most linear in θ, i.e., ∥g(θ, x)∥ ≤ L 1 (∥θ -θ * ∥ + 1) and ∥ḡ(θ)∥ ≤ L 1 (∥θ -θ * ∥ + 1).). There exists µ > 0 such that ⟨θ -θ ′ , ḡ(θ) -ḡ(θ ′ )⟩ ≤ -µ∥θ -θ ′ ∥ 2 , ∀θ, θ ′ ∈ R d . Consequently, the target equation ḡ(θ) = 0 has a unique solution θ * .When g is a gradient field, Assumption 3 is equivalent to strong convexity. For notational simplicity, we assume the strong monotonicity parameter satisfies µ ≤ 1 -r, where r is the convergence factor in Assumption 1. For general µ, our results remain valid with µ replaced by min{µ, 1 -r}.We next consider the noise. Denote by F k the filtration generated by {x t , θ t , ξ t+1 } k-1 t=0 ∪ {x k , θ k }. Assumption 4 (Noise Sequence). Let p ∈ Z + be given. The noise sequence (ξ k ) k≥1 is a collection of i.i.d. random fields satisfying the following conditions with L 2,p > 0:E[ξ k+1 (θ)|F k ] = 0 and E 1/(2p) [∥ξ 1 (θ)∥ 2p ] ≤ L 2,p (∥θ -θ * ∥ + 1), ∀θ ∈ R d .(2.2)Define C(θ) = E[ξ 1 (θ) ⊗2] and assume that C(θ) is at least twice differentiable. There also existM ϵ , k ϵ ≥ 0 such that for θ ∈ R d , we have max i=1,2 C (i) (θ) ≤ M ϵ 1 + ∥θ -θ * ∥ kϵ .In the sequel, we set L := L 1 + L 2 , and without loss of generality, we assume L ≥ 1.When p = 1, the second inequality in (2.2) only requires linear growth in expectation, which relaxes the almost sure linear growth condition in [17] . The constraint on the covariance matrix C(θ) is lenient and satisfied in most regular enough settings, as shown in [20] . 
analytical challenges and techniques: 
 In this section, we elaborate on the challenges and techniques in proving the above results.Previous work has established weak convergence of (x k , θ k ) separately for nonlinear SA with i.i.d. data, and for Markovian linear SA. The high-level approaches used in two representative prior works can be summarized as follows. The work [20] on nonlinear SGD leverages local linearization of g through Taylor expansion. The work [33] on Markovian linear SA exploits the mixing property of the Markovian noise to regain approximate independence, particularly between x k and θ k-τ for sufficiently large τ. It is tempting to expect that nonlinear SA can be analyzed by combining these two approaches. Perhaps surprisingly, such a simple combination would not work due to the interplay between nonlinearity and Markovian structures.To demonstrate this challenge, let us seek to establish weak convergence in the Wasserstein distance W 2 via forward coupling [29] , an approach employed by both [20, 33] as well as others [25] . Specifically, we consider two SA iterate sequences (θ[1] k ) k≥0 and (θ [2]k ) k≥0 from different initializations θ [1] 0 and θ [2] 0 coupled by sharing the data sequence (x k ) k≥0 : θ[1] k+1 = θ [1] k + αg(θ [1] k , x k ) and θ [2] k+1 = θ [2] k + αg(θ [2]k , x k ). To establish convergence in W 2 , we consider the difference sequenceEQUATIONand it suffices to prove w k converges to 0 in mean square:E[∥w k+1 ∥ 2 ] ≲ ρ k E[∥w 0 ∥ 2 ] for ρ < 1.With this goal in mind and following the idea from [20] , one may first linearize the right-hand side of the difference dynamic (3.1) and obtain the approximationw k+1 ≈ w k + αg ′ (θ [2] k , x k )w k . (3.2) +2α E g ′ (θ [2] k , x k ) -g ′ (θ [2] k-τ , x k ) | F k-τ ≲ ρE[∥w k-τ ∥ 2 ] + αE E ⟨w k-τ , g(θ [1] k , x k ) -g(θ [2] k , x k ) -g(θ [1] k-τ , x k ) + g(θ [2] k-τ , x k )⟩ ♠ | F k-τ ,where we approximatew k ≈ w k-τ , w t g ′ (θ [2] t , x k ) ≈ (g(θ [1] t , x k ) -g(θ [2]k , x k )) for t = k, k -τ and obtain the second term in the last inequality. Next, we propose employing two different Taylor expansions to prove that ♠ is of higher orders of α. We first apply the Taylor expansion to g(θ[1] k , x k )- g(θ [2] k , x k ) and g(θ [1] k-τ , x k ) -g(θ [2] k-τ , x k ). However, this only achieves ♠ ≲ ∥w k ∥ 2 ∥w k ∥ + ατ T 1 , where T 1 = min(∥θ [1] k ∥, ∥θ [2] k ∥, ∥θ [1] k-τ ∥, ∥θ [2] k-τ ∥) + 1. When θ [1]k and θ [2] k are not close to each order, i.e., when ∥w k ∥ is large, ♠ is not necessarily of higher order. Therefore, we consider a second type of Taylor expansion on g(θ[1] k , x k ) -g(θ [1] k-τ , x k ) and g(θ [2] k , x k ) -g(θ [2]k-τ , x k ). The intuition for the second type of Taylor expansion is to analyze and bound ♠ by the small distance between θ[j] k and θ [j] k-τ for j ∈ {1, 2}, even when ∥w k ∥ is large. This achieves ♠ ≲ ∥w k ∥ατ T 1 ∥w k ∥ + ατ T 1 .Simultaneously applying the two Taylor expansions will yield ♠ ≲ ατ ∥w k ∥ 2 T 1 . Finally, we overcome this challenge by carefully analyzing the boundness of T 1 ; see Theorem 4.1 and its proof.In parallel to the above coupling approach, we also explore an alternative approach by verifying the joint Markov chain (x k , θ k ) satisfies certain irreducibility and Lyapunov drift conditions, which in turn imply the chain is ergodic. To apply this approach, we exploit additional properties of the SA noise, namely miniorization, which is satisfied in many applications where additional randomness is injected to the SA update. While the high level strategy of this approach is well developed [23, 50] , carrying out the analysis of each step is technically involved. In particular, we need to translate the minorization property of the noise to the irreducibility of the joint chain (x k , θ k ), which is nontrivial in the presence of Markovian noise and nonlinearity. 
weak convergence of projected sa: 
 Our first main result proves the ergodicity of the joint process (x k , θ k ) k≥0 of the projected SA (2.1). Theorem 4.1 (Ergodicity of Projected SA). Suppose that Assumption 1-4 (p = 1) hold. The projected SA (2.1) is applied with radius parameter 2∥θ * ∥ ≤ β < ∞. For stepsize α > 0 that satisfies the constraintατ α ≤ µ (940+96β)L 2 , the Markov chain (x k , θ k ) k≥0 converges to a unique stationary distribution να ∈ P 2 (X × R d ). Let ν α := L(θ ∞ ) be the second marginal of να . For k ≥ 2τ α , it holds that W 2 (L(θ k ), ν α ) ≤ W2 (L(x k , θ k ), να ) ≤ (1 -αµ) k/2 • s(θ 0 , θ * , µ, L, R).Theorem 4.1 generalizes prior weak convergence results for constant stepsize SA/SGD either under i.i.d. noise [20, 67] or linear update [33, 45] . Our stepsize condition ατ α ≲ µ/L 2 coincides with [33, 60] on linear SA, a special case of our setting.The proof of Theorem 4.1 highlights the stabilizing effect of the projection operation in (2.1) . This effect, together with the smoothness of update function g, controls how the Markovian correlation propagates through the nonlinear update, allowing us to overcome the challenges discussed in Section 3. It is unclear whether our proof, which is based on Markov chain coupling, can be fully generalized to SA without projection. Nevertheless, we show that such a generalization is possible for a sub-family of nonlinear SA where g possesses the additional structure termed "asymptotic linearity", which is satisfied by, e.g., SGD applied to certain settings of logistic regression. For a formal statement of this result and proof, we refer the readers to Appendix E.As a by-product of our analysis, we establish the following non-asymptotic 2p-th moment bound on the error θ k -θ * . Let θ t+1/2 := θ t + α(g(θ t , x t ) + ξ t+1 (θ t )) denote the pre-projection iterate. EQUATIONατ ) p for sufficiently large k, generalizing the results of [17, 20, 60] to higher moments and the nonlinear Markovian setting. Notably, this result holds even without the projection operation in the SA update (2.1), i.e., β = ∞. Furthermore, Proposition 4.2 can be used to derive high-probability tail bounds using the Markov inequality.Parallel to the coupling approach, we consider an alternative approach for establishing weak convergence via verifying irreducibility, positive Harris recurrence, and V -uniform ergodicity [50] of the Markov chain (x k , θ k ). This approach applies to nonlinear SA even without projection. To verify irreducibility, we exploit the following additional noise structure.Assumption 5 (Noise Minorization). For each θ ∈ R d , the distribution of the random variable ξ 1 (θ), denoted by ζ θ , can be decomposed asζ θ = ζ 1,θ + ζ 2,θ, where the measure ζ 1,θ has a density, denoted by p θ , which satisfies inf θ∈C p θ (t) > 0 for any bounded set C and any t ∈ R d .A similar assumption is considered in [5, 67] . This assumption is mild and satisfied by any continuous random field supported on R d . Introducing such (small) continuous noise is often part of the algorithm design for inducing privacy [2, 22] or exploration [28, 55] . Without Assumption 5, the chain may fail to be irreducible even when the other assumptions are satisfied; see [33] for a counterexample.Under Assumption 5, we obtain the following ergodicity result paralleling Theorem 4.1.Theorem 4.3 (Ergodicity of SA -Minorization). Suppose that Assumption 1-3, Assumption 4(p = 1), and Assumption 5 hold. For stepsize α > 0 that satisfies the constraint ατ α L 2 < c 2 µ, the Markov chain (x k , θ k ) k≥0 of (2.1) with β = ∞ is V -uniformly ergodic with Lyapunov function V (x, θ) = ∥θ -θ * ∥ 2 + 1 and a unique stationary distribution να ∈ P 2 (X × R d ). Moreover, defining the V -norm ∥ν∥ V := |ν(dx)|V (x), we haveEQUATIONwhere the constants ρ ∈ (0, 1) and κ ∈ (0, ∞) may depend on α.In the sequel, letEQUATION∞ ) denote the random vector whose law is the stationary distribution να given in Theorem 4.1. As a corollary, we have geometric convergence for the first 2 moments of θ k . Corollary 4.4 (Non-Asymptotic Convergence Rate). Under the setting of Theorem 4.1, for any initialization of θ 0 ∈ R d , we haveE[θ k ] -E[θ (α) ∞ ] ≤ (1 -αµ) k/2 • s ′ (θ 0 , θ * , µ, L, R), and E[θ k θ ⊤ k ] -E[θ (α) ∞ (θ (α) ∞ ) ⊤ ] ≤ (1 -αµ) k/2 • s ′′ (θ 0 , θ * , µ, L, R).Moreover, the convergence rate established in Theorem 4.1 is fast enough that we can use it to prove a Central Limit Theorem for the average iterates. Corollary 4.5 (Central Limit Theorem). Under the setting of Theorem 4.1, as k → ∞ we have1 √ k k-1 t=0 θ t -E[θ ∞ ] ⇒ N (0, Σ (a)), where Σ (α) := lim k→∞1 k E k-1 t=0 θ t -E[θ (α) ∞ ] ⊗2 .Establishing the CLT sets the stage for using the SA iterates for statistical inference tasks such as confidence interval estimation. We discuss this in greater detail in Section 4.4 after characterizing the asymptotic bias, another important ingredient for using SA for inference.In this subsection, we characterize the asymptotic bias E[θ∞ ] -θ * . Understanding the bias structure has important algorithmic implications for bias reduction, which we explore in Section 4.5, as well as for more efficient statistical inference and confidence interval estimation [34] . Theorem 4.6 (Bias Characterization). Suppose Assumptions 1-4(p = 3) hold. For each stepsize α > 0 satisfying ατ α L 2 < c 3 µ, the following holds for some vector b independent of α : EQUATIONEQUATIONEQUATIONEQUATION)with A = (ḡ ′ (θ * ) ⊗ I + I ⊗ ḡ′ (θ * )) -1 and h(θ * , x) = X (I -P * + Π) -1 (P * -Π)(x, dx ′ )g(θ * , x ′ ),with the kernel P * being a regular conditional probability on X that satisfies B π(dx)P (x, C) = C π(dy)P * (y, B), for all B, C ∈ B(X ).We defer the detailed proof to Appendix I. A few remarks are in order. First, we emphasize that (4.2) is essentially an equality, indicating a non-zero bias of order α whenever b ̸ = 0 (up to higher order terms). Notably, the Polyak-Ruppert averaging of the iterates cannot eliminate this bias. Note that the bias expansion in (4.2) applies to both weakly converged projected and non-projected SA. Our analysis shows that compared with the non-projected SA, the projection operator induces an extra bias term of the order O(α 2 τ 3 α ), which is negligible relative to the main terms in in (4.2). More importantly, Theorem 4.6 provides an explicit expression of the leading bias, which decomposes into three components: the Markovian part, the nonlinearity contribution, and a compound term, which is unique in nonlinear Markovian SA. Specifically, b m in (4.3) is associated with the Markovian multiplicative noise, where the matrix P * -Π in the h function determines the mixing time of the data sequence (x k ) k≥ . The term b n in (4.4) is linked to nonlinearity, as reflected by the Hessian term ḡ′′ (θ * ) = E[g ′′ (θ * , x ∞ )], which quantifies the nonlinearity of g and is equal to zero in the case of a linear g. Lastly, b c in (4.5) is the compound term, due to its dependence on both the Markov noise (h function) and the nonlinearity measure ḡ′′ . In particular, we note the following two special cases: (1) When g is a linear function, ḡ′′ (θ * ) = 0. Hence, b n = b c = 0, and b m recovers the result in [33] ; (2) When (x k ) k≥0 is i.i.d. sampled from the stationary distribution π, we have h(θ * , x) ≡ 0 ∀x ∈ X , for P = P * = Π. As such, b m = b c = 0, recovering the result in [20] . The presence of the compound term b c suggests that as the SA structure becomes more nonlinear and the underlying Markov chain mixes more slowly, the impact on the bias is multiplicative rather than simply additive, a surprising phenomenon not unveiled in previous studies. It is possible to improve the residual order from O(α 3/2 ) to O(α 2 ) with a more refined characterization of the asymptotic second moment E[(θ ∞ -θ * ) ⊗2 ] by following a similar strategy as our current approach. We leave this refinement out of the scope of the current paper.We examine the practical implications of our weak convergence and bias characterization results, particularly for Polyak-Ruppert (PR) tail averaging and Richardson-Romberg (RR) extrapolation. In this subsection, we focus on the dependence on the stepsize α and iteration index k, and make use of the big-O notation from Section 1. Recall that b is the bias vector defined in Theorem 4.6.PR averaging [56, 59] is a classical approach for reducing the variance and accelerating the convergence of SA. Here we consider the tail-averaging variant of PR averaging, defined as θk0,k := 1 k-k0 k-1 t=k0 θ t , for k ≥ k 0 , with a user-specified burn-in period k 0 ≥ 0 (a common choice is k 0 = k/2). The following corollary, proved in Appendix J, provides a non-asymptotic bound on the mean squared error (MSE) for the averaged iterates θk0,k . Corollary 4.7 (Tail Averaging). Under the setting of Theorem 4.6, the tail-averaged iterates satisfy the following bounds for all k > k 0 + 2τ α and kEQUATIONT1: asymptotic squared biasEQUATIONT3: optimization error .Corollary 4.7 shows that the MSE can be decomposed into three terms and elucidates how these terms depend on α, k, and other problem parameters. In particular, the term T 1 corresponds to the asymptotic squared bias ∥E[θ∞ -θ * ]∥ 2 , which is not affected by averaging. The term T 2 is associated with the variance Var( θk0,k ), which decays at rate 1/k due to averaging. Lastly, the term T 3 represents the optimization error ∥E θk0,k -θ ∞ ∥ 2 , which decays geometrically in k 0 thanks to the use of a constant stepsize α and the tail-averaging procedure.Note that averaging does not affect the bias of order α. With the precise bias characterization in Theorem 4.6, we can order-wise reduce the bias to O(α 3/2 ) by employing the RR extrapolation technique [61] . Let θ(α) k0,k and θ(2α) k0,k denote the tail-averaged iterates using two stepsizes α and 2α with the same data (x k ) k≥0 . The RR extrapolated iterates are defined as θ (α) k0,k = 2 θ(α) k0,k - θ(2α) k0,k .≥ τ α + 1 αµ log 1 ατα , E ∥ θ k0,k -θ * ∥ 2 = O (ατ α ) 3 + O τ α k -k 0 + O (1 -αµ) k0/2 α (k -k 0 ) 2 .Backed by the CLT in Corollary 4.5, the iterates of constant-stepsize SA can be used to construct confidence intervals of θ * . For i.i.d. data or linear SA, this approach has been explored in [34, 47, 66, 67] along with an appropriate variance estimator [26, 66] . In our Markovian nonlinear setting, where the iterates are biased, it is crucial to use RR extrapolation for bias reduction. Once the bias is accounted for, the power of using constant stepsizes reveals itself as it leads to rapid mixing and low correlation of the iterates. Together, they lead to efficient confidence interval estimation schemes using nonlinear Markovian SA; see the empirical results in [34] showing its efficacy. In contrast, the classical diminishing stepsize paradigm often suffers from high correlation [13] and in turn inaccurate variance estimation, resulting in unsatisfactory coverage probability with finite data [34] .Generalized linear models (GLM) extend linear regression to the modelE[Y |W ] = σ(W ⊤ θ * ),where W is the covariate, Y the response variable, and σ is called the mean function. For any monotone (and potentially nonlinear) σ, the powerful framework developed in [19, 36, 37, 64] allows one to formulate the estimation of θ * as minimizing an appropriate convex (surrogate) loss function.Applying SGD to this loss leads to a nonlinear SA update, to which our results are applicable. Below we discuss their applications in two concrete examples of GLMs.Logistic Regression Logistic regression uses a sigmoid mean function σ(x) = 1 1+exp(-x) . Suppose the covariate w k is sequentially sampled from a uniformly ergodic Markov chain with a bounded state space W ⊂ R d , and conditioned on w k the response y k is Bernoulli distributed with parameter (1 + exp(-w ⊤ k θ * )) -1 . SGD applied to the L 2 -regularized negative log-likelihood function takes the form of the SA updateθ k+1 = θ k + αg(θ k , x k ), where x k = (w k , y k ) ∈ W × {0, 1} and g(θ k , x k ) = -w k σ(-w ⊤ k θ k ) -y k -λθ k .For simplicity, we do not consider ξ-perturbation, i.e., ξ k+1 (θ k ) ≡ 0. It is easy to verify that this g is strongly monotone and sufficiently smooth with at most linear growth in |θ|, hence satisfying Assumption 1-3. Therefore, all the results in Sections 4.1-4.5 apply to logistic regression with constant stepsizes and Markovian data.The mean function σ can be interpreted as playing a similar role as the activation function in neural networks. Widely adopted is ReLU activation σ(x) = max(0, x) as well as its various smooth approximations [7, 30] . The problem of learning θ * in this setting, sometimes called ReLU Regression, has been studied in the last decade and recently regained attention [19, 36, 37, 64] . Unlike linear or logistic regression, the least squares and maximum likelihood formulation associated with such nonlinear mean functions σ is non-convex. Nevertheless, the convex surrogate loss framework in [19, 64] still applies. As an example, we focus on the SoftPlus activation σ(x) = log(1 + exp(ιx))/ι with a temperature parameter ι > 0 [30] . With L 2 -regularization the resulting SGD iteration isθ k+1 = θ k -α w k 1 ι log(1 + exp(ιw ⊤ k θ k )) -y k + λθ k ,where the covariate-response pair (θ k , x k ) is as before. This problem can again be cast as nonlinear SA with a strongly monotone and smooth g, satisfying Assumptions 1-3. All results in Sections 4.1-4.5 apply. 
numerical experiments: 
 In this section, we provide numerical experiment results to verify our theoretical results. We run SGD on L 2 -regularized logistic regression with Markovian data and constant stepsizes, where the covariate x t ∈ R is sequentially sampled from an autoregressive (AR) model of order 1; specifically, x t+1 = 0.9x t + ζ t+1 with ζ t i.i.d. following N (0, 1), and the stationary distribution is x ∞ ∼ N (0, 1/(1-0.9 2 )). The binary dependent variable y t is sampled from Bernoulli(1/(1+exp(-w * x t )) with w * = 1. The regularized parameter is set to λ = 0.0001. 1 To examine the asymptotic bias, we run the experiment for an episode length of 10 7 , with Markovian data as well as i.i.d. data sampled from N (0, 1/(1 -0.9 2 )). We plot the errors (distance to θ * ) of the PR averaged iterates and the RR extrapolated iterates, for different stepsizes α. Figure 1 (a) verifies the presence of an asymptotic bias approximately proportional to the stepsize α, and illustrates the effectiveness of RR extrapolation in reducing this bias. In Figure 1 (b), we compare the bias under Markovian data (x t+1 ∼ P (•|x t )) and i.i.d. data (x t ∼ x ∞ ). Interestingly, Figure 1 (b) reveals that Markovian data does not necessarily lead to a larger bias than i.i.d. data. This is consistent with our theory, as the three bias terms b m , b n , b c may have opposite signs leading to cancellation. This result suggests that in the presence of nonlinearity, one should not avoid Markovian data simply for the sake of reducing bias. Rather, RR extrapolation may be more effective for bias reduction. To verify the CLT in Corollary 4.5, we repeat the experiment 1000 times with an episode length of 10 6 and stepsize α = 0.8. We compute the PR averaged iterates and plot the histogram and the quantile-quantile (QQ) plot in Figure 2 in Appendix C. The close alignment between the histogram and the normal curve in Figure 2 (a) and the linearity of the points along the 45-degree reference line in the QQ plot in Figure 2 (b) confirm that the empirical distribution follows a normal distribution. 
related work: 
 General SA and SGD. SA and SGD can be traced back to the seminal work of [58] . Classical work assumes a diminishing stepsize sequence, and has shown almost sure asymptotic convergence to θ * [8, 58] . Subsequent works propose the iterate averaging technique, now known as Polyak-Ruppert (PR) averaging, to reduce variance and accelerate convergence [56, 59] , and also establish a Central Limit Theorem for the asymptotic normality of the averaged iterates [57] . The asymptotic convergence theory of SA and SGD is well developed and extensively addressed in many exemplary textbooks, see [3, 40, 65] . There are also recent works studying the non-asymptotic convergence with diminishing stepsizes [12, 14] . The recent work [15] establishes the high probability bound on the estimation error of contractive SA with diminishing stepsize.Stepsizes. There has been an increasing interest in studying SA with constant stepsize. Many works in this line provide non-asymptotic upper bounds on mean squared error [25, 42, 51] study linear SA (LSA) under i.i.d. data. Recent works extend the analysis of the MSE to LSA with Markovian data, such as [24, 52, 60] . There are also works providing upper bounds of MSE for general contractive SA with Markovian noise [14, 17] .(MSE) E[∥θ t -θ * ∥ 2 ]. Works inIn addition to non-asymptotic guarantees, some works focus on the asymptotic behavior of SA iterates. Recent works have shown that when using constant stepsize, one loses the almost sure convergence guarantee in the diminishing stepsize sequence regime, and at best can achieve distributional convergence, as demonstrated in [16, 20, 25, 33, 66, 67, 69] . The presence of asymptotic bias is also a recurring theme in recent literature, with precise characterization given in [20] for strongly-convex SGD with i.i.d. data and in [33] for LSA with Markovian data. Works in [34, 51, 66, 67, 69 ] also establish Central Limit Theorems for averaged SA iterates with constant stepsizes. 
conclusion: 
 We provide the first weak convergence and steady-state analysis for constant-stepsize SA with both nonlinear update and Markovian data. Our analysis elucidates the compound effect of nonlinearity and memory, which leads to new analytical challenges and behaviors. A limitation of our results is the use of a projection step or the noise minorization assumption. Whether they can be removed is worth investigating. Other future directions include refining the dimension dependence in our results, as well as a theoretical investigation of statistical inference. 
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: The computation resources are disclosed in Section 5.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?Answer: [No] Justification: This paper is theoretical and has no potential negative social impacts.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: This paper poses no such risks.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [NA] Justification: This paper does not use existing assets.Guidelines:• The answer NA means that the paper does not use existing assets.• 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 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: This 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: This paper does not involve crowdsourcing nor research with 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] Justification: This paper does not involve crowdsourcing nor research with human subjects.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.• We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: General Probability We write z 1 ⊥ ⊥ z 2 | z 3 if random variables z 1 and z 2 are conditionally independent given z 3 .Recall that we define the metric d (x, θ), (x ′ , θ ′ ) := 1{x ̸ = x ′ } + ∥θ -θ ′ ∥ 2 for the space X ×R d . Thus, for μ and ν in P 2 (X × R d ), the Wasserstein-2 distance w.r.t. d is computed asW2 (μ, ν) = inf E[1{x ̸ = x ′ } + ∥θ -θ ′ ∥ 2 ] 1 2 : L (x, θ) = μ, L (x ′ , θ ′ ) = ν .General State Space Markov Chains Throughout the paper, we assume that X is a Borel space. Let P denote the transition kernel. We call π the stationary distribution ofP if it satisfies X π(dx)P (x, B) = π(B), for B ∈ B(X ). Define the π-weighted inner product ⟨f, g⟩ L 2 (π) = X π(dx)f ⊤ (x)g(x) and the induced norm ∥f ∥ L 2 (π) = (⟨f, f ⟩ L 2 (π) ) 1/2 . Let L 2 (π) = {f : ∥f ∥ L 2 (π)< ∞} denote the corresponding Hilbert space of R d -valued, squareintegrable and measurable functions on X . For an operator T :L 2 (π) → L 2 (π), its operator norm is defined as ∥T ∥ L 2 (π) = sup ∥f ∥ L 2 (π) =1 ∥T f ∥ L 2 (π). The transition kernel is a bounded linear operator on L 2 (π), in particular with norm ∥P ∥ L 2 (π)=1 . Also, we define the kernel/operatorΠ = 1 ⊗ π by Π(x, •) = π.Throughout the paper, we assume that X is a Borel space. Let P denote the transition kernel. We call π the stationary distribution of P if it satisfies X π(dx)P (x, B) = π(B), for B ∈ B(X ). There exists a kernel P * as a regular conditional probability that satisfies A π(ds)P (x, B) = B π(dy)P * (y, A), for A, B ∈ B(X ) [27, Chapter 21.4, Theorem 19], and P * defines the probability law for the time-reversed chain of (x k ) k≥0 .In this section, we provide a more detailed discussion of the assumptions taken in this work.Projection and Minorization Projection steps have a longstanding presence in SA literature for tractability in convergence theory, as seen in many analyses of SGD [6, 11, 39, 41, 54] . Although not an algorithmic proposal, this additional projection step does not incur much computational cost in practice, as it only involves rescaling the iterates, and the projection radius can be estimated a priori. Before our work, no studies had proven weak convergence for non-linear SA with Markovian data and constant stepsize, with or without the projection. Thus, our result is the first to prove detailed weak convergence in this setting. Nonetheless, we provide an alternative proof of weak convergence in Theorem 4.3 using the Drift and Minorization technique, which does not require a projection.Differentiability The differentiability condition of the SA update operator g in Assumption 2 ensures controlled evolution of the iterates θ. This differentiability assumption supports a third-order Taylor expansion of g with a bounded remainder, which is crucial for both analyzing the ♠ term in the convergence proof as discussed in Section 3 and for bias characterization. Moreover, some form of differentiability assumption is standard in SA literature, such as [1, 20, 45] , particularly when one seeks a fine-grained characterization of the iterates' distributional property. Such an assumption is satisfied by many GLMs, such as logistic regression and Poisson regression. When g is not differentiable, the bias of SA iterates behaves drastically different [68] , which is beyond the scope of this paper.The strong monotonicity assumption is common in SA literature. Together with smoothness, it allows us to establish geometric distributional convergence. While some GLMs by themselves do not satisfy this condition, it is a common practice to apply L 2 -regularization (equivalently, weight decay) to ensure strong convexity and improve statistical performance. It is a standard calculation that one can appropriately choose the regularization parameter to derive tight results for non-strongly-convex functions.In this section, we present the plots from numerical experiments in Section 5 that verify the Central Limit Theorem (CLT) in Corollary 4.5. In these plots, we plot the centered and scaled PR-averaged iterates, i.e,√ T θ(k) T - 1000 l=1 θ(l) T, where θ(k) = T t=1 θ t /T with T = 10 6 is the PR-averaged iterates for the k-th repeat. In this section, we prove the pilot result, namely Proposition 4.2. We prove the desired moments for β = ∞, i.e., without any projection. It is easy to see that when the projection radiusβ ∈ [2∥θ * ∥, ∞], E[∥θ t+1 -θ * ∥ 2p ] ≤ E[∥θ t+1/2 -θ * ∥ 2p ],where θ t+1/2 denotes the iterate before projection. The term on the right hand side can be further bounded by the moment bounds for iteration without projection. Therefore, it suffices for us to prove the respective moment bounds without any projection.Given Assumption 4 hold for 2p-th moment, with p ≥ 1, we prove the moment bound in Proposition 4.2 for n with 1 ≤ n ≤ p by induction.In this section, we prove the base case of Proposition 4.2, i.e., with n = 1. The base case gives the desired mean squared error (MSE) convergence bound, which will subsequently be used in the proof of weak convergence.We start by noting the following decomposition,E[∥θ k+1 -θ * ∥ 2 ] -E[∥θ k -θ * ∥ 2 ] = 2αE[⟨θ k -θ * , g(θ k , x k )⟩] + α 2 E[∥g(θ k , x k )∥ 2 ] + α 2 E[∥ξ k+1 (θ k )∥ 2 ] = 2αE[⟨θ k -θ * , g(θ k , x k ) -ḡ(θ k )⟩] + 2αE[⟨θ k -θ * , ḡ(θ k )⟩] + α 2 E[∥g(θ k , x k )∥ 2 ] + α 2 E[∥ξ k+1 (θ k )∥ 2 ].It is easy to see that under Assumption 3, we have⟨θ k -θ * , ḡ(θ k )⟩ = ⟨θ k -θ * , ḡ(θ k ) -ḡ(θ * )⟩ ≤ -µ∥θ k -θ * ∥ 2 . (D.1)Additionally, under Assumption 2 and 4, we have the following upper boundα 2 E[∥g(θ k , x k )∥ 2 ] + E[∥ξ k+1 (θ k )∥ 2 ] ≤ α 2 L 2 1 E[(∥θ k -θ * ∥ + 1) 2 ] + L 2 2 E[(∥θ k -θ * ∥ + 1) 2 ] ≤ 2α 2 L 2 E[∥θ k -θ * ∥ 2 ] + 1 . (D.2)Therefore, the key to analyze the remaining inner product⟨θ k -θ * , g(θ k , x k ) -ḡ(θ k )⟩.Consider the following decomposition⟨θ k -θ * , g(θ k , x k ) -ḡ(θ k )⟩ = ⟨θ k -θ k-τ , g(θ k , x k ) -ḡ(θ k )⟩ (D.3) + ⟨θ k-τ -θ * , g(θ k-τ , x k ) -ḡ(θ k-τ )⟩ (D.4) + ⟨θ k-τ -θ * , g(θ k , x k ) -g(θ k-τ , x k )⟩ (D.5) + ⟨θ k-τ -θ * , ḡ(θ k ) -ḡ(θ k-τ )⟩. (D.6)Hence, we need some upper bound on ∥θ k -θ k-τ ∥.We next note the following technical Lemma, which is adapted from [17, 60] E[∥θ k -θ k-τ ∥|F k-τ ] ≤ 2ατ L∥θ k-τ -θ * ∥ + 2ατ L (D.7) E[∥θ k -θ k-τ ∥|F k-τ ] ≤ 4ατ LE[∥θ k -θ * ∥|F k-τ ] + 4ατ L (D.8) E[∥θ k -θ k-τ ∥ 2 |F k-τ ] ≤ 8α 2 τ 2 L 2 ∥θ k-τ -θ * ∥ 2 + 8α 2 τ 2 L 2 (D.9) E[∥θ k -θ k-τ ∥ 2 |F k-τ ] ≤ 32α 2 τ 2 L 2 E[∥θ k -θ * ∥ 2 |F k-τ ] + 32α 2 τ 2 L 2 . (D.10)Given (D.10), we additionally note that∥θ k-τ -θ * ∥ 2 + 1 = E[∥θ k-τ -θ * ∥ 2 |F k-τ ] + 1 ≤ 2 E[∥θ k -θ k-τ ∥ 2 |F k-τ ] + E[∥θ k -θ * ∥ 2 |F k-τ ] + 1 ≤ 2 32α 2 τ 2 L 2 (E[∥θ k -θ * ∥ 2 |F k-τ ] + 1) + E[∥θ k -θ * ∥ 2 |F k-τ ] + 1 ≤ 4(E[∥θ k -θ * ∥ 2 |F k-τ ] + 1). (D.11)We next use the above four technical inequalities to analyze the four terms in (D.3)-(D.4).To bound (D.3), we first note that∥E[⟨θ k -θ k-τ , g(θ k , x k ) -ḡ(θ k )⟩|F k-τ ]∥ ≤ E[∥θ k -θ k-τ ∥ • 2L(∥θ k -θ * ∥ + 1)|F k-τ ] (i) ≤ 2L E[∥θ k -θ k-τ ∥ 2 |F k-τ ] E[(∥θ k -θ * ∥ + 1) 2 |F k-τ ] (ii) ≤ 2L 32α 2 τ 2 L 2 (E[∥θ k -θ * ∥ 2 |F k-τ ] + 1)] 2(E[∥θ k -θ * ∥ 2 |F k-τ ] + 1) ≤ 16ατ L 2 (E[∥θ k -θ * ∥ 2 |F k-τ ] + 1),where (i) holds for the Cauchy-Schwarz inequality and (ii) holds for (D.10).To bound (D.4), we next note that∥E[⟨θ k-τ -θ * , g(θ k-τ , x k ) -ḡ(θ k-τ )⟩|F k-τ ]∥ = ∥⟨θ k-τ -θ * , E[g(θ k-τ , x k ) -ḡ(θ k-τ )|F k-τ ]⟩∥ ≤ ∥θ k-τ -θ * ∥∥E[g(θ k-τ , x k ) -ḡ(θ k-τ )|F k-τ ]∥ (iii) ≤ ∥θ k-τ -θ * ∥ • αL(∥θ k-τ -θ * ∥ + 1) ≤ 2αL(∥θ k-τ -θ * ∥ 2 + 1) (iv) ≤ 8αL(E[∥θ k -θ * ∥ 2 |F k-τ ] + 1) ≤ 8ατ L 2 (E[∥θ k -θ * ∥ 2 |F k-τ ] + 1),where (iii) holds due to the mixing property of Markov chain (x k ) k≥0 and (iv) holds for (D.11).To bound (D.5), we have∥E[⟨θ k-τ -θ * , g(θ k , x k ) -g(θ k-τ , x k )⟩|F k-τ ]∥ = ∥⟨θ k-τ -θ * , E[g(θ k , x k ) -g(θ k-τ , x k )|F k-τ ]⟩∥ ≤ L∥θ k-τ -θ * ∥ • E[∥θ k -θ k-τ ∥|F k-τ ] (v) ≤ L∥θ k-τ -θ * ∥ • 2ατ L(∥θ k-τ -θ * ∥ + 1) ≤ 4ατ L 2 (∥θ k-τ -θ * ∥ 2 + 1) (vi) ≤ 16ατ L 2 (E[∥θ k -θ * ∥ 2 |F k-τ ] + 1),where (v) holds for (D.7) and (vi) holds for (D.11).Lastly, to bound (D.6), we apply the similar technique used in bounding the third term in (D.5) and obtain a similar result∥E[⟨θ k-τ -θ * , ḡ(θ k ) -ḡ(θ k-τ )⟩|F k-τ ]∥ ≤ 16ατ L 2 (E[∥θ k -θ * ∥ 2 |F k-τ ] + 1).Combining all analyses above, we have∥2αE[⟨θ k -θ * ⟩, g(θ k , x k ) -ḡ(θ k )|F k-τ ]∥ ≤ 2α(16ατ L 2 + 8ατ L 2 + 32ατ L 2 )(E[∥θ k -θ * ∥ 2 |F k-τ ] + 1) ≤ 112α 2 τ L 2 (E[∥θ k -θ * ∥ 2 |F k-τ ] + 1). (D.12)Hence, making use of (D.1), (D.2), and (D.12), we obtain the followingE[∥θ k+1 -θ * ∥ 2 |F k-τ ] -E[∥θ k -θ * ∥ 2 |F k-τ ] ≤ -2αµE[∥θ -θ * ∥ 2 |F k-τ ] + 112α 2 τ L 2 (E[∥θ k -θ * ∥ 2 |F k-τ ] + 1) + 2α 2 L 2 (E[∥θ k -θ * ∥ 2 |F k-τ ] + 1) ≤ -2αµE[∥θ -θ * ∥ 2 |F k-τ ] + 114α 2 τ L 2 E[∥θ k -θ * ∥ 2 |F k-τ ] + 114α 2 τ L 2 = -2α(µ -57ατ L 2 )E[∥θ k -θ * ∥ 2 |F k-τ ] + 114α 2 τ L 2 .Therefore, when we have α satisfying the constraint, i.e., ατ L 2 < c 2,1 µ, we obtainE[∥θ k+1 -θ * ∥ 2 |F k-τ ] ≤ (1 -αµ)E[∥θ k -θ * ∥ 2 |F k-τ ] + 114α 2 τ L 2 .Recursively, we getE[∥θ k -θ * ∥ 2 ] ≤ (1 -αµ) k-τ E[∥θ τ -θ * ∥ 2 ] + 114ατ L 2 µ ≤ 2(1 -αµ) k-τ E[∥θ 0 -θ * ∥ 2 ] + E[∥θ τ -θ 0 ∥ 2 ] + 114ατ L 2 µ ≤ 2(1 -αµ) k-τ E[∥θ 0 -θ * ∥ 2 ] + 8α 2 τ 2 L 2 E[∥θ 0 -θ * ∥ 2 ] + 1 + 114ατ L 2 µ ≤ 4(1 -αµ) k-τ E[∥θ 0 -θ * ∥ 2 ] + 122ατ L 2 µ .Lastly, we note that1 (1 -αµ) τ (i) ≤ 1 1 -ατ µ (ii) ≤ 1 1 -ατ L (iii) ≤ 2, (D.13)where (i) holds by the Bernoulli inequality, that (1 + x) r ≥ 1 + rx for x ≥ -1 and r ≥ 1; (ii) holds for µ ≤ L; (iii) holds for ατ L < µ/(114L) < 1 2 .Hence, for k ≥ τ , we haveE[∥θ k -θ * ∥ 2 ] ≤ c 2,1 (1 -αµ) k ∥θ 0 -θ * ∥ 2 + c 2,2 ατ α L 2 µ ,for c 2,1 and c 2,2 some universal constants. As such, we have completed the proof of base case for Proposition 4.2.D.1.1 Proof of Lemma D.1In this section, we provide the proofs of the four technical inequalities in Lemma D.1.Proof of (D.7).E[∥θ k -θ k-τ ∥|F k-τ ] ≤ 2ατ L∥θ k-τ -θ * ∥ + 2ατ L.Proof. Note that∥θ k -θ k-τ ∥ ≤ k-1 t=k-τ ∥θ t+1 -θ t ∥, so we start with analyzing ∥θ t+1 -θ t ∥. ∥θ t+1 -θ * ∥ -∥θ t -θ * ∥ ≤ ∥θ t+1 -θ t ∥ = α∥g(θ t , x t ) + ξ t+1 (θ t )∥ ≤ α∥g(θ t , x t )∥ + α∥ξ t+1 (θ t )∥ ≤ αL 1 (∥θ t -θ * ∥ + 1) + α∥ξ t+1 (θ t )∥ ∥θ t+1 -θ * ∥ ≤ (1 + αL 1 )∥θ t -θ * ∥ + αL 1 + α∥ξ t+1 (θ t )∥.Recall that we assumeE 1/2 [∥ξ t+1 (θ t )∥ 2 |F t ] ≤ L 2 (∥θ t ∥ + 1), then we have for k -τ ≤ t ≤ k, E[∥ξ t+1 (θ t )∥|F k-τ ] = E[E[∥ξ t+1 (θ k )|F t ]|F k-τ ] ≤ E[L 2 (∥θ k ∥ + 1)|F k-τ ] E[∥θ t+1 -θ t ∥|F k-τ ] ≤ αL(E[∥θ t -θ * ∥|F k-τ ] + 1) E[∥θ k+1 -θ * ∥|F k-τ ] ≤ (1 + αL)E[∥θ k -θ * ∥|F k-τ ] + αL. Hence, for 0 ≤ n ≤ τ , E[∥θ k-τ +n -θ * ∥|F k-τ ] ≤ (1 + αL) n E[∥θ k-τ -θ * ∥|F k-τ ] + αL n-1 l=0 (1 + αL) l = (1 + αL) n ∥θ k-τ -θ * ∥ + ((1 + αL) n -1).We next note that(1 + x) y = e y log(1+x) ≤ e xy ≤ 1 + 2xy, xy ∈ [0, 1/2].Hence, at this stage, if we require ατ L < µ/(4L) < 1/4, we have the following upper bound(1 + αL) n ≤ (1 + αL) τ ≤ 1 + 2ατ L ≤ 2.Therefore, for 0≤ n ≤ τ , E[∥θ k-τ +n -θ * ∥|F k-τ ] ≤ (1 + 2ατ L)∥θ k-τ -θ * ∥ + 2ατ L ≤ 2∥θ k-τ -θ * ∥ + 2ατ L.As such, we haveE[∥θ k -θ k-τ ∥|F k-τ ] ≤ k-1 t=k-τ E[∥θ t+1 -θ t ∥|F k-τ ] ≤ αL k-1 t=k-τ E[∥θ t -θ * ∥|F k-τ ] + ατ L ≤ ατ L(2∥θ k-τ -θ * ∥ + 2ατ L) + ατ L ≤ 2ατ L∥θ k-τ -θ * ∥ + 2ατ L,and prove the desired inequality.Proof of (D.8).E[∥θ k -θ k-τ ∥|F k-τ ] ≤ 4ατ LE[∥θ k -θ * ∥|F k-τ ] + 4ατ L.Proof. We prove this inequality based on the claim that we have just shown,E[∥θ k -θ k-τ ∥|F k-τ ] ≤ 2ατ L∥θ k-τ -θ * ∥ + 2ατ L.We simply note that∥θ k-τ -θ * ∥ = E[∥θ k-τ -θ * ∥|F k-τ ] ≤ E[∥θ k -θ k-τ ∥|F k-τ ] + E[∥θ k -θ * ∥|F k-τ ].Hence,E[∥θ k -θ k-τ ∥|F k-τ ] ≤ 2ατ L(E[∥θ k -θ k-τ ∥|F k-τ ] + E[∥θ k -θ * ∥|F k-τ ] + 1) (1 -2ατ L)E[∥θ k -θ k-τ ∥|F k-τ ] ≤ 2ατ LE[∥θ k -θ * ∥|F k-τ ] + 2ατ L.Therefore, we obtainE[∥θ k -θ k-τ ∥|F k-τ ] ≤ 4ατ LE[∥θ k -θ * ∥|F k-τ ] + 4ατ L.Proof of (D.9).E[∥θ k -θ k-τ ∥ 2 |F k-τ ] ≤ 32α 2 τ 2 L 2 E[∥θ k -θ * ∥ 2 |F k-τ ] + 32α 2 τ 2 L 2 . Proof. To analyze E[∥θ k -θ k-τ ∥ 2 |F k-τ ],we consider the following attempt.E[∥θ k -θ k-τ ∥ 2 |F k-τ ] ≤ τ k-1 t=k-τ E[∥θ t+1 -θ t ∥ 2 |F k-τ ] = α 2 τ k-1 t=k-τ E[(∥g(θ t , x t )∥ + ∥ξ t+1 (θ t )∥) 2 |F k-τ ] ≤ 2α 2 τ k-1 t=k-τ E[∥g(θ t , x t )∥ 2 |F k-τ ] + E[∥ξ t+1 (θ t )∥ 2 |F k-τ ] ≤ 2α 2 τ L 2 k-1 t=k-τ E[(∥θ t -θ * ∥ + 1) 2 |F k-τ ] ≤ 4α 2 τ L 2 k-1 t=k-τ E[∥θ t -θ * ∥ 2 |F k-τ ] + 4α 2 τ 2 L 2 . Next, we study E[∥θ t -θ * ∥ 2 |F k-τ ]. We start with the following, for k -τ ≤ t < k, E[∥θ t+1 -θ * ∥ 2 |F k-τ ] = E[∥θ t -θ * ∥ 2 |F k-τ ] + 2αE[⟨θ t -θ * , g(θ t , x t )⟩|F k-τ ] + α 2 E[∥g(θ t , x t ) + ξ t+1 (θ t )∥ 2 |F k-τ ] ≤ E[∥θ t -θ * ∥ 2 |F k-τ ] + 2α 2 E[∥g(θ t , x t )∥ 2 |F k-τ ] + E[∥ξ t+1 (θ t )∥ 2 |F k-τ ] + 2αE[∥θ t -θ * ∥∥g(θ t , x t )∥|F k-τ ].We note that2E[∥θ t -θ * ∥∥g(θ t , x t )∥|F k-τ ] ≤ 2 E[∥θ t -θ * ∥ 2 |F k-τ ]E[∥g(θ t , x t )∥ 2 |F k-τ ] ≤ 2 E[∥θ t -θ * ∥ 2 |F k-τ ]E[(L(∥θ t -θ * ∥ + 1)) 2 |F k-τ ] ≤ 2 E[∥θ t -θ * ∥ 2 |F k-τ ]2L 2 E[∥θ t -θ * ∥ 2 + 1|F k-τ ] ≤ 4L(E[∥θ t -θ * ∥ 2 |F k-τ ] + 1).Substituting the above inequality back, we obtainE[∥θ t+1 -θ * ∥ 2 |F k-τ ] ≤ E[∥θ t -θ * ∥ 2 |F k-τ ] + 4α 2 L 2 E[∥θ t -θ * ∥ 2 |F k-τ ] + 1 + 4αLE[∥θ t -θ * ∥ 2 |F k-τ ] + 4αL ≤ (1 + 4α 2 L 2 + 4αL)E[∥θ t -θ * ∥ 2 |F k-τ ] + (4α 2 L 2 + 4αL).We further recall that 4α 2 L 2 ≤ 4αL(ατ L) ≤ αL, and hence we obtain the following upper boundE[∥θ t+1 -θ * ∥ 2 |F k-τ ] ≤ (1 + 5αL)E[∥θ t -θ * ∥ 2 |F k-τ ] + 5αL.Then, recursively, for 0 ≤ n ≤ τ , we haveE[∥θ k-τ +n -θ * ∥ 2 |F k-τ ] ≤ (1 + 5αL) n ∥θ k-τ -θ * ∥ 2 + 5αL n-1 l=0 (1 + 5αL) l .As such, under the assumption that 4ατ L < µ/(4L) < 1/4, then for k -τ ≤ t ≤ k, we haveE[∥θ t -θ * ∥ 2 |F k-τ ] ≤ (1 + 10ατ L)∥θ k-τ -θ * ∥ 2 + 10ατ L ≤ 2∥θ k-τ -θ * ∥ 2 + 10ατ L.Combining all the analyses above, we haveE[∥θ k -θ k-τ ∥ 2 |F k-τ ] ≤ 4α 2 τ L 2 k-1 t=k-τ E[∥θ t -θ * ∥ 2 |F k-τ ] + 4α 2 τ 2 L 2 ≤ 4α 2 τ 2 L 2 2∥θ k-τ -θ * ∥ 2 + 10ατ L + 4α 2 τ 2 L 2 ≤ 8α 2 τ 2 L 2 ∥θ k-τ -θ * ∥ 2 + 8α 2 τ 2 L 2 .Proof of (D.10).E[∥θ k -θ k-τ ∥ 2 |F k-τ ] ≤ 32α 2 τ 2 L 2 E[∥θ k -θ * ∥ 2 |F k-τ ] + 32α 2 τ 2 L 2 .Proof. This inequality simply extends the result from (D.9),i.e.,E[∥θ k -θ k-τ ∥ 2 |F k-τ ] ≤ 8α 2 τ 2 L 2 ∥θ k-τ -θ * ∥ 2 + 8α 2 τ 2 L 2 .We first note that∥θ k-τ -θ * ∥ 2 = E[∥θ k-τ -θ * ∥ 2 |F k-τ ] ≤ 2E[∥θ k -θ k-τ ∥ 2 |F k-τ ] + 2E[∥θ k -θ * ∥ 2 |F k-τ ].Hence,E[∥θ k -θ k-τ ∥ 2 |F k-τ ] ≤ 8α 2 τ 2 L 2 (2E[∥θ k -θ k-τ ∥ 2 |F k-τ ] + 2E[∥θ k -θ * ∥ 2 |F k-τ ] + 1) (1 -16α 2 τ 2 L 2 )E[∥θ k -θ k-τ ∥ 2 |F k-τ ] ≤ 16α 2 τ 2 L 2 E[∥θ k -θ * ∥ 2 |F k-τ ] + 8α 2 τ 2 L 2 .Again, under the assumption that 16ατ L 2 < µ/4, we can conclude thatE[∥θ k -θ k-τ ∥ 2 |F k-τ ] ≤ 32α 2 τ 2 L 2 E[∥θ k -θ * ∥ 2 |F k-τ ] + 32α 2 τ 2 L 2 .In this step, assume that the moment bound in Proposition 4.2 has been proven for k ≤ n -1, we now proceed to show that the desired moment convergence holds for n with 2 ≤ n ≤ p.We start with the following decomposition of ∥θ k+1 -θ * ∥ 2n∥θ k+1 -θ * ∥ 2n = ∥θ k -θ * ∥ 2 + 2α⟨θ k -θ * , g(θ k , x k ) + ξ k+1 (θ k )⟩ + α 2 ∥g(θ k , x k ) + ξ k+1 (θ t )∥ 2 n = i,j,l i+j+l=n n i, j, l ∥θ k -θ * ∥ 2i 2α⟨θ k -θ * , g(θ k , x k ) + ξ k+1 (θ k )⟩ j α∥g(θ k , x k ) + ξ k+1 (θ k )∥ 2lWe note the following cases.1. i = n, j = l = 0. In this case, the summand is simply ∥θ k -θ * ∥ 2i . 2. When i = n -1, j = 1 and l = 0. In this case, the summand is of order α, i.e., α2n⟨θ k - 1) . We can further compose it asθ * , g(θ k , x k ) + ξ k+1 (θ k )⟩ j ∥θ k -θ * ∥ 2(n-2nα⟨θ k -θ * , g(θ k , x k ) + ξ k+1 (θ k )⟩∥θ k -θ * ∥ 2(n-1) = 2nα⟨θ k -θ * , g(θ k , x k ) -ḡ(θ k ) + ξ k+1 (θ k )⟩∥θ k -θ * ∥ 2(n-1) T1 + 2nα⟨θ k -θ * , ḡ(θ k )⟩∥θ k -θ * ∥ 2(n-1).Note that, when (x k ) is i.i.d. or from a martingale noise sequence, we haveE[T 1 |θ k ] = 0.However, when (x k ) is Markovian, the above equality then does not hold and T 1 requires a careful analysis. Nonetheless, under the strong monotonicity assumption, we haveT 2 ≤ -2nαµ∥θ k -θ * ∥ 2n .3. For the remaining terms, we see that they are of higher orders of α. Therefore, when α is selected sufficiently small, these terms do not raise concern.Therefore, to prove the desired moment bound, we spend the remaining section analyzing T 1 .Immediately, we note thatE[T 1 |F k-τ ] = E 2nα⟨θ k -θ * , g(θ k , x k ) -ḡ(θ k ) + E[ξ k+1 (θ k )|θ k ]⟩∥θ k -θ * ∥ 2(n-1) |F k-τ = E 2nα⟨θ k -θ * , g(θ k , x k ) -ḡ(θ k )⟩∥θ k -θ * ∥ 2(n-1) T ′ 1 |F k-τ .Subsequently, we focus on analyzing T ′ 1 . We start with the following decomposition ofT ′ 1 . 2nα⟨g(θ k , x k ) -ḡ(θ k ), θ k -θ * ⟩∥θ k -θ * ∥ 2(n-1) ≤ 2nα∥g(θ k-τ , x k ) -ḡ(θ k-τ )∥∥θ k-τ -θ * ∥ 2n-1 (D.14) + 2nα∥g(θ k , x k ) -g(θ k-τ , x k )∥∥θ k-τ -θ * ∥ 2n-1 (D.15) + 2nα∥ḡ(θ k-τ ) -ḡ(θ k )∥∥θ k-τ -θ * ∥ 2n-1 (D.16) + 2nα∥g(θ k , x k ) -ḡ(θ k )∥∥θ k -θ * ∥ 2(n-1) ∥θ k -θ k-τ ∥ (D.17) + 2nα∥g(θ k , x k ) -ḡ(θ k )∥∥θ k-τ -θ * ∥ • ∥θ k -θ * ∥ 2(n-1) -∥θ k-τ -θ * ∥ 2(n-1) . (D.18)We note the following technical lemma, which will offer significant help in the analysis of T ′ 1 . We postpone the proof of the lemma to the end of this subsection.Lemma D.2. For cn ατ ≤ µ/(4L 2 ), where cn denotes some constant dependent of the higher-moment 2n, we haveE[∥θ k -θ k-τ ∥ 2n |F k-τ ] ≤ c n α 2n τ 2n L 2n (∥θ k-τ -θ * ∥ 2n + 1).Following the lemma, we observe that a natural consequence is for any m ≤ 2n, we haveE[∥θ k -θ k-τ ∥ m |F k-τ ] ≤ E[∥θ k -θ k-τ ∥ 2n |F k-τ ] m 2n ≤ c n α 2n τ 2n L 2n (∥θ k-τ -θ * ∥ 2n + 1) m 2n ≤ c m α m τ m L m ∥θ k-τ -θ * ∥ m + 1 ,where we use the inequality a p + b p > (a + b) p for a, b > 0, p ∈ (0, 1) to obtain the final inequality. Now, we are ready to analyze (D.14)-(D.18). Firstly, for (D.14), we make use of the mixing assumption of τ , and have thatE[|(D.14)||F k-τ ] ≤ 2nα∥θ k-τ -θ * ∥ 2n-1 E[∥g(θ k-τ , x k ) -ḡ(θ k-τ )∥|F k-τ ] ≤ 2nα 2 L∥θ k-τ -θ * ∥ 2n-1 (∥θ k-τ -θ * ∥ + 1) ≤ 2nα 2 L∥θ k-τ -θ * ∥ 2n + 2nα 2 L∥θ k-τ -θ * ∥ 2n-1 ≤ 3nα 2 L∥θ k-τ -θ * ∥ 2n + nα 2 L∥θ k-τ -θ * ∥ 2(n-1) ,where we make use of the inequality 2|x| 3 ≤ x 2 + x 4 to obtain the final step.Next, we proceed to analyze (D.15). It is easy to see that 1) .E[|(D.15)||F k-τ ] = 2nα∥θ k-τ -θ * ∥ 2n-1 E[∥g(θ k , x k ) -g(θ k-τ , x k )|F k-τ ] ≤ 2nα∥θ k-τ -θ * ∥ 2n-1 E[∥θ k -θ k-τ ∥|F k-τ ] ≤ 2nα∥θ k-τ -θ * ∥ 2n-1 2ατ L(∥θ k-τ -θ * ∥ + 1) ≤ 4nα 2 τ L∥θ k-τ -θ * ∥ 2n + 4nα 2 τ L∥θ k-τ -θ * ∥ 2n-1 ≤ 6nα 2 τ L∥θ k-τ -θ * ∥ 2n + 2nα 2 τ L∥θ k-τ -θ * ∥ 2(n-The term in (D. 16 ) can be analyzed in a similar fashion as the (D.15).For (D.17), we first derive the followingE[|(D.17)||F k-τ ] ≤ 2nαE 2L ∥θ k -θ * ∥ + 1 ∥θ k -θ k-τ ∥∥θ k -θ * ∥ 2(n-1) |F k-τ = 4nαLE[∥θ k -θ k-τ ∥∥θ k -θ * ∥ 2n-1 |F k-τ ] Ta + 4nαLE[∥θ k -θ k-τ ∥∥θ k -θ * ∥ 2(n-1) |F k-τ ] T b .We next analyze the two terms T a and T b respectively. Starting with T a , we have4nαLE[∥θ k -θ k-τ ∥∥θ k -θ * ∥ 2n-1 |F k-τ ] (D.19) ≤ 4nαLE ∥θ k -θ k-τ ∥ ∥θ k -θ k-τ ∥ + ∥θ k-τ -θ * ∥ 2n-1 |F k-τ (D.20) ≤ 2 2(n-1) 4nαLE[∥θ k -θ k-τ ∥(∥θ k -θ k-τ ∥ 2n-1 + ∥θ k-τ -θ * ∥ 2n-1 |F k-τ ] (D.21) = 4 n nαL E[∥θ k -θ k-τ ∥ 2n |F k-τ ] by Lemma D.2 +∥θ k-τ -θ * ∥ 2n-1 E[∥θ k -θ k-τ ∥|F k-τ ] (D.22) ≤ 4 n nαL c n α 2n τ 2n L 2n ≤2ατ L (∥θ k-τ -θ * ∥ 2n + 1) + ∥θ k-τ -θ * ∥ 2n-1 (2ατ L(∥θ k-τ -θ * ∥ + 1)) (D.23) ≤ 4 n nαL 4ατ L∥θ k-τ -θ * ∥ 2n + 2ατ L∥θ k-τ -θ * ∥ 2n-1 + c n α 2n τ 2n L 2n (D.24) ≤ 4 n nαL 5ατ L∥θ k-τ -θ * ∥ 2n + ατ L∥θ k-τ -θ * ∥ 2(n-1) + c ′ n α 2n-1 τ 2n-1 L 2n-1 . (D.25)For T b , we have4nαLE[∥θ k -θ k-τ ∥∥θ k -θ * ∥ 2(n-1) |F k-τ ] ≤ 4nαLE ∥θ k -θ k-τ ∥ ∥θ k -θ k-τ ∥ + ∥θ k-τ -θ * ∥ 2(n-1) |F k-τ ≤ 2 2n-1 nαLE[∥θ k -θ k-τ ∥(∥θ k -θ k-τ ∥ 2(n-1) + ∥θ k-τ -θ * ∥ 2(n-1) |F k-τ ] = 2 2n-1 nαL E[∥θ k -θ k-τ ∥ 2n-1 |F k-τ ] by Lemma D.2 +∥θ k-τ -θ * ∥ 2(n-1) E[∥θ k -θ k-τ ∥|F k-τ ] ≤ 2 2n-1 nαL c n-1 α 2n-1 τ 2n-1 L 2n-1 ≤2ατ L (∥θ k-τ -θ * ∥ 2n-1 + 1) + ∥θ k-τ -θ * ∥ 2(n-1) (2ατ L(∥θ k-τ -θ * ∥ + 1)) ≤ 2 2n-1 nαL 4ατ L∥θ k-τ -θ * ∥ 2n-1 + 2ατ L∥θ k-τ -θ * ∥ 2(n-1) + c n-1 α 2n-1 τ 2n-1 L 2n-1 ≤ 2 2n-1 nαL 2ατ L∥θ k-τ -θ * ∥ 2n + 4ατ L∥θ k-τ -θ * ∥ 2(n-1) + c n-1 α 2n-1 τ 2n-1 L 2n-1 .Combining the analyses of the two terms, we get the following upper bound to (D.17)E[|(D.17)||F k-τ ] ≤ 4 n nαL 5ατ L∥θ k-τ -θ * ∥ 2n + ατ L∥θ k-τ -θ * ∥ 2(n-1) + c ′ n α 2n-1 τ 2n-1 L 2n-1 + 2 2n-1 nαL 2ατ L∥θ k-τ -θ * ∥ 2n + 4ατ L∥θ k-τ -θ * ∥ 2(n-1) + c n-1 α 2n-1 τ 2n-1 L 2n-1 = 2 2n-1 nαL 12ατ L∥θ k-τ -θ * ∥ 2n + 6ατ L∥θ k-τ -θ * ∥ 2(n-1) + c ′′ n-1 α 2n-1 τ 2n-1 L 2n-1 .Lastly, we analyze (D.18). We first make use of the mean-value theorem, with a ∈ [0, 1], we have∥θ k -θ * ∥ 2(n-1) -∥θ k-τ -θ * ∥ 2(n-1) = ∥θ k -θ k-τ ∥ • 2(n -1)∥a(θ k -θ * ) + (1 -a)(θ k-τ -θ * )∥ 2n-3 = ∥θ k -θ k-τ ∥ • 2(n -1)∥a(θ k -θ k-τ ) + θ k-τ -θ * ∥ 2n-3 ≤ 2 2n-3 (n -1)∥θ k -θ k-τ ∥ ∥θ k -θ k-τ ∥ 2n-3 + ∥θ k-τ -θ * ∥ 2n-3Substituting the above upper bound back into (D.18), we obtainE[|(D.18)||F k-τ ] ≤ 2 2n-1 n(n -1)αL∥θ k-τ -θ * ∥ E (∥θ k -θ * ∥ + 1)∥θ k -θ k-τ ∥ ∥θ k -θ k-τ ∥ 2n-3 + ∥θ k-τ -θ * ∥ 2n-3 |F k-τ ≤ 2 2n-1 n(n -1)αL ∥θ k-τ -θ * ∥E[∥θ k -θ k-τ ∥ 2n-1 |F k-τ ] + ∥θ k-τ -θ * ∥ 2 E[∥θ k -θ k-τ ∥ 2n-2 |F k-τ ] + ∥θ k-τ -θ * ∥E[∥θ k -θ k-τ ∥ 2n-2 |F k-τ ] + ∥θ k-τ -θ * ∥ 2n-2 E[∥θ k -θ k-τ ∥ 2 |F k-τ ] + ∥θ k-τ -θ * ∥ 2n-1 E[∥θ k -θ k-τ ∥|F k-τ ] + ∥θ k-τ -θ * ∥ 2n-2 E[∥θ k -θ k-τ ∥|F k-τ ] ≤ 2 2n-1 n(n -1)αL c n ατ L∥θ k-τ -θ * ∥ 2n + c n-1 ατ L∥θ k-τ -θ * ∥ 2(n-1) + c n-1 α 2n-1 τ 2n-1 L 2n-1 .Combining the analyses above, we have the following bound for T 1 ,E[|T 1 ||F k-τ ] ≤ E[∥(D.14)||F k-τ ] + E[∥(D.15)||F k-τ ] + E[∥(D.16)||F k-τ ] + E[∥(D.17)||F k-τ ] + E[∥(D.18)||F k-τ ] ≤ c n,1 α 2 τ L 2 ∥θ k-τ -θ * ∥ 2n + c n,2 α 2 τ L 2 ∥θ k-τ -θ * ∥ 2(n-1) + c n,3 α 2n τ 2n-1 L 2n ,where c n,1 , c n,2 and c n,3 are some constants that depend on n.Additionally, we note that∥θ k-τ -θ * ∥ 2n = E[∥θ k-τ -θ * ∥ 2n |F k-τ ] ≤ E ∥θ k -θ k-τ ∥ + ∥θ k -θ * ∥ 2n |F k-τ ≤ 2 2n-1 E[∥θ k -θ k-τ ∥ 2n |F k-τ ] + 2 2n-1 E[∥θ k -θ * ∥ 2n |F k-τ ] ≤ c n α 2n τ 2n L 2n (∥θ k-τ -θ * ∥ 2n + 1) + 2 2n-1 E[∥θ k -θ * ∥ 2n |F k-τ ]. Therefore, for sufficiently small ατ L < µ/(c ′ n L), we have (1 -c ′ n α 2n τ 2n L 2n )∥θ k-τ -θ * ∥ 2n ≤ c ′′ n E[∥θ k -θ * ∥ 2n |F k-τ ] + c n α 2n τ 2n L 2n ⇒ ∥θ k-τ -θ * ∥ 2n ≤ 2c ′′ n E[∥θ k -θ * ∥ 2n |F k-τ ] + 2c n α 2n τ 2n L 2n .As such, for sufficiently small α, we haveE[|T 1 ||F k-τ ] ≤ c n,1 α 2 τ L 2 c n E[∥θ k -θ * ∥ 2n |F k-τ ] + c ′ n α 2n τ 2n L 2n + c n,2 α 2 τ L 2 c n-1 E[∥θ k -θ * ∥ 2(n-1) |F k-τ ] + c ′ n-1 α 2(n-1) τ 2(n-1) L 2(n-1) + c n,3 α 2n τ 2n-1 L 2n = c n,1 α 2 τ L 2 E[∥θ k -θ * ∥ 2n |F k-τ ] + c n,2 α 2 τ L 2 E[∥θ k -θ * ∥ 2(n-1) |F k-τ ] + c n,3 α 2n τ 2n-1 L 2n .Hence, up til this point, we have obtainedE[∥θ k+1 -θ * ∥ 2n |F k-τ ] ≤ (1 -2nαµ)E[∥θ k -θ * ∥ 2n |F k-τ ] + c n,1 α 2 τ L 2 E[∥θ k -θ * ∥ 2n |F k-τ ] + c n,2 α 2 τ L 2 E[∥θ k -θ * ∥ 2(n-1) |F k-τ ] + c n,3 α 2n τ 2n-1 L 2n ≤ (1 -2nα(µ -c ′ n,1 ατ L 2 ))E[∥θ k -θ * ∥ 2n |F k-τ ] + c n,2 α 2 τ L 2 E[∥θ k -θ * ∥ 2(n-1) ]|F k-τ ] + c n,3 α 2n τ 2n-1 L 2n .Following the induction hypothesis, when k is sufficiently large, we haveE[∥θ k -θ * ∥ 2(n-1) |F k-τ ] ≤ c n-1 α n-1 τ n-1 s(θ 0 , L, µ).Substituting the above upper bound back into our analysis of the 2n-th moment bound, we obtainE[∥θ k+1 -θ * ∥ 2n |F k-τ ] ≤ (1 -2nα(µ -c ′ n,1 ατ L 2 ))E[∥θ k -θ * ∥ 2n |F k-τ ] + α n+1 τ n L 2 c n,2 • c n-1 s(θ 0 , L, µ) + c n,3 α 2n τ 2n-1 L 2n .Subsequently, if we set α sufficiently small, such that ατ L 2 < c n • µ, we obtainE[∥θ k+1 -θ * ∥ 2n |F k-τ ] ≤ (1 -αµ)E[∥θ k -θ * ∥ 2n |F k-τ ] + α n+1 τ n c ′ n,2 • s(θ 0 , L, µ), where s(θ 0 , L, µ) is some constant that may depend on the initialization θ 0 and the problem primitives µ and L but is independent of α.E[∥θ k -θ * ∥ 2n ] ≤ (1 -αµ) k-τ E[∥θ τ -θ * ∥ 2n ] + α n τ n • s(θ 0 , L, µ).Lastly, we recall thatE[∥θ τ -θ * ∥ 2n ] ≤ 2 2n-1 E[∥θ τ -θ 0 ∥ 2n ] + 2 2n-1 E[∥θ 0 -θ * ∥ 2n ] ≤ c n,1 α 2n τ 2n L 2n (E[∥θ 0 -θ * ∥ 2n ] + 1) + c n,2 ∥θ 0 -θ * ∥ 2n ≤ c n,1 E[∥θ 0 -θ * ∥ 2n ] + c n,2 α 2n τ 2n L 2n. Substituting back, we obtain for sufficiently large k,E[∥θ k -θ * ∥ 2n ] ≤ c n,1 (1 -αµ) k-τ E[∥θ 0 -θ * ∥ 2n ] + α 2n τ 2n s(θ 0 , L, µ).As such, we have proven the desired n-th moment bound.We now come back to Lemma D.2 and provide the complete proof.Proof. The proof follows a similar strategy as (D.9) and (D.10) in Section D.1.1.We start with the following relaxation and obtain thatE[∥θ k -θ k-τ ∥ 2n |F k-τ ] ≤ E k-1 t=k-τ ∥θ t+1 -θ t ∥ 2n |F k-τ ≤ τ 2n-1 k-1 t=k-τ E[∥θ t+1 -θ t ∥ 2n |F k-τ ] = α 2n τ 2n-1 k-1 t=k-τ E[∥g(θ t , x t ) + ξ t+1 (θ t )∥ 2n |F k-τ ] ≤ 2 2n-1 α 2n τ 2n-1 k-1 t=k-τ E[∥g(θ t , x t )∥ 2n |F k-τ ] + E[∥ξ t+1 (θ t )∥ 2n |F k-τ ] ≤ 2 2n-1 α 2n τ 2n-1 k-1 t=k-τ L 2n 1 E[(∥θ t -θ * ∥ + 1) 2n |F k-τ ] + L 2n 2 (E[∥θ t -θ * ∥|F k-τ ] + 1) 2n ≤ 4 2n-1 α 2n τ 2n-1 L 2n k-1 t=k-τ E[∥θ t -θ * ∥ 2n |F k-τ ] + 4 2n-1 α 2n τ 2n L 2n .Next, in order to obtain a bound on ∥θ t -θ * ∥ 2n , we study the following term.E[∥θ t+1 -θ * ∥ 2n |F k-τ ] ≤ E ∥θ t+1 -θ t ∥ + ∥θ t -θ * ∥ 2n |F k-τ = 2n i=0 2n i E[∥θ t+1 -θ t ∥ i ∥θ t -θ * ∥ 2n-i |F k-τ ] = E[∥θ t -θ * ∥ 2n |F k-τ ] + 2n i=1 α i E[∥g(θ t , x t ) + ξ t+1 (θ t )∥ i ∥θ t -θ * ∥ 2n-i |F k-τ ]Note thatE[∥g(θ t , x t ) + ξ t+1 (θ t )∥ i ∥θ t -θ * ∥ 2n-i |F k-τ ] ≤ 2 i-1 E ∥g(θ t , x t )∥ i + ∥ξ t+1 (θ t )∥ i ∥θ t -θ * ∥ 2n-i |F k-τ = 2 i-1 E E[(∥g(θ t , x t )∥ i + ∥ξ t+1 (θ t )∥ i )|θ t ]∥θ t -θ * ∥ 2n-i |F k-τ ≤ 2 i-1 L i E (∥θ t -θ * ∥ + 1) 2n |F k-τ ≤ 2 2(n-1) 2 i L i E[∥θ t -θ * ∥|F k-τ ] + 1Substituting back, we obtain2n i=1 α i E[∥g(θ t , x t ) + ξ t+1 (θ t )∥ i ∥θ t -θ * ∥ 2n-i |F k-τ ] ≤ 2 2(n-1) 2n i=1 2 i α i L i E[∥θ t -θ * ∥|F k-τ ] + 1 = 2 2(n-1) E[∥θ t -θ * ∥|F k-τ ] + 1 • 2αL(1 + 2αL) 2n-1 ≤ 4 2n-1 αL E[∥θ t -θ * ∥|F k-τ ] + 1Consolidating the terms, we haveE[∥θ t+1 -θ * ∥ 2n |F k-τ ] ≤ (1 + 4 2n-1 αL) E[∥θ t -θ * ∥|F k-τ ] + 1Recursively, for 0 ≤ l ≤ τ , we haveE[∥θ k-τ +l -θ * ∥ 2n |F k-τ ] ≤ (1 + 4 2n-1 αL) l ∥θ k-τ -θ * ∥ 2n + 4 2n-1 αL l-1 i=0 (1 + 4 2n-1 αL) i = (1 + 4 2n-1 αL) l ∥θ k-τ -θ * ∥ 2n + (1 + 4 2n-1 αL) lThen, for4 2n-1 ατ L ≤ µ/4L < 1/4, we have for k -τ ≤ t ≤ k, E[∥θ t -θ * ∥ 2n |F k-τ ] ≤ (1 + 2 4n-1 ατ L)∥θ k-τ -θ * ∥ 2n + 2 4n-1 ατ L ≤ 2∥θ k-τ -θ * ∥ 2n + 2 4n-1 ατ LFinally, we haveE[∥θ k -θ k-τ ∥ 2n |F k-τ ] ≤ 4 2n-1 α 2n τ 2n-1 L 2n k-1 t=k-τ E[∥θ t -θ * ∥ 2n |F k-τ ] + 4 2n-1 α 2n τ 2n L 2n ≤ 4 2n-1 α 2n τ 2n L 2n 2∥θ k-τ -θ * ∥ 2n + 2 4n-1 ατ L + 1 ≤ 2 4n-1 α 2n τ 2n L 2n ∥θ k-τ -θ * ∥ 2n + 1 .As such, we have completed the proof.In this section, we prove the weak convergence result in Theorem 4.1. In fact, the proof of the projected SA weak convergence result can be seen as a special case of unprojected SA with the asymptotic linearity condition, which we have briefly discussed in Section 4. Therefore, the proof proceeds in the following two subsections. First, we formally define the asymptotic linearity condition and present our weak convergence result for unprojected SA under this additional assumption. Next, we relate this result for unprojected SA to projected SA and specialize the proof to obtain Theorem 4.1.In this subsection, we formally introduce the asymptotic linearity condition, which is crucial for establishing weak convergence in the context of unprojected SA (β = ∞). Additionally, we explore the implications of this condition. Assumption 6 (Asymptotic Linearity). The noise sequence (ξ k ) k≥1 is a collection of i.i.d. random fields satisfying the following conditions: (1)E[ξ k+1 (θ)|F k ] = 0, (2) there exists a constant L 3 > 0 such that ξ 1 is L 3 -Lipschitz, i.e., ∥ξ 1 (θ) -ξ 1 (θ ′ )∥ ≤ L 3 ∥θ -θ ′ ∥, for all θ, θ ′ ∈ R d , and (3) ∥ξ 1 (0)∥ ≤ L 3 .Moreover, there exists a function G(•) : X → R d×d such that given ϵ > 0, define δ(ϵ) := min δ : ∥g ′ (θ, x) -G(x)∥ ≤ ϵ, ∀x ∈ X and ∀θ ∈ {θ : ∥θ∥ ≥ δ} , and we have lim ϵ→0 ϵδ(ϵ) = 0.The first part of Assumption 6 states that the random field grows at most linearly in θ. The second part of Assumption 6 implies that g ′ (θ, x) converges to a limit G(x) when ∥θ∥ → ∞ for all x ∈ X , which shows the asymptotic linearity of g(θ, x). Furthermore, Assumption 6 also requires how fast g ′ (θ, x) converges to G(x). A sufficient condition under which the second part of Assumption 6 holds is that there exists ω > 0 such that ∥θ∥ 1+ω ∥g ′ (θ, x) -G(x)∥ < ∞ for ∀θ ∈ R d and x ∈ X . We can verify that to ensure ∥g ′ (θ, x) -G(x)∥ < ϵ, we can set ∥θ∥ ∈ Θ(ϵ -1 1+ω ), which can ensure ϵδ(ϵ) ∈ O(ϵ ω 1+ω ) → 0 as ϵ → 0. This sufficient condition implies that g ′ (θ, x) uniformly converge to G(x) with convergence rate of O(∥θ∥ -(1+w) ). By definition, we conclude that the structure of linear SA is also asymptotic linear. Besides that, the 1-dimensional logistic regression also satisfies Assumption 6. For 1-dimensional logistic regression, we have g(θ, x, y) = x 1 1+e -θx -y + λθ, where (x, y) presents the data. Therefore, we have g ′ (θ, x, y) = x 2 e -θx (1+e -θx ) 2 + λ and g ′ (θ, x, y) uniformly converges to λ with geometric convergence rate, thereby satisfies the Assumption 6.With the asymptotic linearity condition now formally defined, we proceed to prove the weak convergence for unprojected SA. For convenient reference, we state the theorem below. Theorem E.1 (Ergodicity of SA-Asymptotic Linearity). Suppose that Assumption 1-Assumption 4 hold. Additionally, assume 6. For stepsize α > 0 that satisfies the constraint ατ α L 2 < min(c 2 µ, κ µ ), with c 2 formalized in Proposition 4.2 and κ µ defined in (E.1), the Markov chain (x k , θ k ) k≥0 converges to a unique stationary distribution να ∈ P 2 (X × R d ).Moreover, there exist κ µ > 0 and some universal constant c ′ such thatϵδ(ϵ) ≤ c ′ µ, ∀ϵ ≤ κ µ . (E.1)We let ν α := L(θ ∞ ) be the second marginal of να . For k ≥ 2τ α , it holds thatW 2 (L(θ k ), ν α ) ≤ W2 (L(x k , θ k ), να ) ≤ (1 -αµ) k/2 • s(θ 0 , L, µ). (E.2)The proof of Theorem E.1 consists of two major steps. Firstly, we assume that x 0 ∼ π, and show that (x k , θ k ) k≥0 converges to a unique limiting invariant distribution. Next, we relax the assumption of x 0 ∼ π, and prove that for arbitrary initialization (x 0 , θ 0 ) ∈ X × R d , the Markov chain will converge to the same limit.Step 1: Initialization with x 0 ∼ π. To prove the convergence of the Markov chain, we consider the following coupling construction. We have a pair of Markov chains (x k , θk ) k≥0 and (x k , θk ) k≥0 sharing the same underlying process and noise, i.e., (x k , ξ k+1 ) k≥0 , i.e., θ[1] k+1 = θ [1] k + α(g(θ [1] k , x k ) + ξ k+1 (θ [1] k )), θ [2] k+1 = θ [2] k + α(g(θ [2] k , x k ) + ξ k+1 (θ [2] k )). (E.3)We assume that the initial iterates θ [1] 0 and θ [2] 0 may depend on each other and on x 0 , but are independent of subsequent (x k ) k≥1 given x 0 . For the iterates difference θ[1] k -θ [2]k , we have the following Proposition E.2, whose proof is given at the end of this subsection. Proposition E.2. ∀k ≥ τ and ατ ≤ min( µ 908L 2 , κµ L 2 ), we haveE[∥θ [1] k -θ [2] k ∥ 2 ] ≤ 4(1 -µα) k-τ E[∥θ [1] 0 -θ [2] 0 ∥ 2 ], where κ µ > 0 and ϵδ(ϵ) ≤ µ 768 , ∀ϵ ≤ κ µ .By Proposition E.2 and the definition of W 2 and W2 , we haveW 2 2 L(θ [1] k ), L(θ [2] k ) (i) ≤ W 2 2 L(x k , θ [1] k ), L(x k , θ [2] k ) (ii) ≤ E[∥θ [1] k -θ [2] k ∥ 2 ] (iii) ≤ 4(1 -µα) k-τ E[∥θ [1] 0 -θ [2] 0 ∥ 2 ], (E.4)where (i) and (ii) hold by the definition of W 2 and W2 and (iii) holds by applying Proposition E.2.Note that equation (E.4) always holds for any joint distribution of initial iterates (x 0 , θ[1] 0 , θ [2]0 ). Recall that P * represents the transition kernel for the time-reversed Markov chain of {x k } k≥0 , and the initial distribution of x 0 is assumed to be mixed already. Given a specific x 0 , we sample x -1 from P * (• | x 0 ). Additionally, we use θ [2] -1 to denote the random varible that satisfies θ[2] -1 d = θ [1]0 and is independent of {x k } k≥0 . Finally, we set θ [2] 0 as θ[2] 0 = θ [2] -1 + α(g(x -1 , θ [2] -1 ) + ξ 0 (θ [2] -1 )).By the property of time-reversed Markov chain, we have{x k } k≥-1 d = {x k } k≥0 . Given that θ [2] -1 d = θ [1] 0 and θ [2] -1 is independent with {x k } k≥-1 , we can prove (x k , θ [2] k ) d = (x k+1 , θ [1]k+1 ) for k ≥ 0. We thus have for all k ≥ τ :W 2 2 L x k , θ [1] k , L x k+1 , θ [1] k+1 = W 2 2 L x k , θ [1] k , L x k , θ [2] k (i) ≤ 4(1 -µα) k-τ E[∥θ [1] 0 -θ [2] 0 ∥ 2 ], where (i) holds by inequality (E.4). Then, we haveEQUATIONk )} k≥0 forms a Cauchy sequence w.r.t. the metric W2 . Since the space P 2 (X × R d ) endowed with W2 is a Polish space, every Cauchy sequence converges [63, Theorem 6.18] . Furthermore, convergence in Wasserstein 2-distance also implies weak convergence [63, Theorem 6.9]. Therefore, we conclude that the sequence (L(x k , θ [1] k )) k≥0 converges weakly to a limit distribution μ ∈ P 2 (X × R d ).Now that we have established the existence of a limiting distribution, we next proceed to show the uniqueness. We prove this by contradiction. Note that we currently assume that x 0 ∼ π, hence to show that the limit (x ∞ , θ ∞ ) is unique, we only need to show that the limit is independent of the initial distribution of θ 0 , which can be correlated to x 0 .Consider two Markov chains (x k , θ k ) k≥0 and (x k , θ ′ k ) k≥0 , sharing (x k , ξ k+1 ) k≥0 but with arbitrary initialization of θ 0 and θ ′ 0 . For the sake of contradiction, we assume that (x 0 , θ 0 ) ⇒ (x ∞ , θ ∞ ) and (x 0 , θ ′ 0 ) ⇒ (x ∞ , θ ′ ∞ ) respectively. Then, by the triangle inequality, we have thatW2 (x ∞ , θ ∞ ), (x ∞ , θ ′ ∞ ) ≤ W2 (x ∞ , θ ∞ ), (x k , θ k ) + W2 (x k , θ k ), (x k , θ ′ k ) + W2 (x k , θ ′ k ), (x ∞ , θ ′ ∞ ) → 0.As such, we have shown that the limit ν is unique.Lastly, we prove that ν is invariant. Suppose that we initialize the joint process at its limit, i.e., (x 0 , θ 0 ) ∼ ν. We first apply the triangle inequality, and we obtain W2 (x 1 , θ 1 ), (x 0 , θ 0 ) ≤ W2 (x 1 , θ 1 ), (x k+1 , θ k+1 ) + W2 (x k+1 , θ k+1 ), (x 0 , θ 0 ) .Clearly, as k → ∞, W2 (x k+1 , θ k+1 ), (x 0 , θ 0 ) → 0. To bound W2 (x 1 , θ 1 ), (x k+1 , θ k+1 ) , we need the following lemma. Lemma E.3. Consider two copies of the SA trajectory, where L(x 0 , θ 0 ) = ν and L(x ′ 0 , θ ′ 0 ) is allowed to be arbitrary.W 2 2 L(x 1 , θ 1 ), L(x ′ 1 , θ ′ 1 ) ≤ ρ 1 • W 2 2 L(x 0 , θ 0 ), L(x ′ 0 , θ ′ 0 ) + ρ 2 • W 2 2 L(x 0 , θ 0 ), L(x ′ 0 , θ ′ 0 ) ,Lemma E.4. Under x 0 ∼ π, Assumption 1-4 and 6 and the same setting as Proposition E.2,W 2 2 (L(θ k ), ν α ) ≤ W 2 2 (L(x k , θ k ), να ) ≤ 16(1 -µα) k • E ∥θ [1] 0 -θ * ∥ 2 + c ′ 2,2 .Proof. Let us consider the coupled processes defined as equation (E.3). Suppose that the initial iterate (x 0 , θ [2] 0 ) follows the stationary distribution ν, thus L(x k , θk ) = ν and L(θ[2]k ) = ν for all k ≥ 0. By equation (E.4), we have for all k ≥ τ :W 2 2 L(θ [1] k ), µ = W 2 2 L(θ [1] k ), L(θ [2] k ) ≤ W 2 2 L(x k , θ [1] k ), L(x k , θ [2] k ) ≤ 4(1 -µα) k-τ E[∥θ [1] 0 -θ [2] 0 ∥ 2 ] ≤ 8(1 -µα) k-τ • E ∥θ [1] 0 -θ * ∥ 2 + E ∥θ ∞ -θ * ∥ 2 ≤ 16(1 -µα) k • E ∥θ [1] 0 -θ * ∥ 2 + E ∥θ ∞ -θ * ∥ 2 , (E.5)where we make use of the derivation in (D.13) to obtain the last inequality.We note thatE[∥θ ∞ -θ * ∥ 2 ] ≤ lim inf k→∞ E[∥θ k -θ * ∥ 2 ∞ ] ≤ c 2,2 • ατ L 2 /µ ≤ c ′ 2,2 ,the last inequality holds for ατ ≤ min( µ 908L 2 , κµ L 2 ). Therefore, we prove the desired inequalityW 2 2 L(θ [1] k ), µ ≤ 16(1 -µα) k • E ∥θ [1] 0 -θ * ∥ 2 + c ′ 2,2 .Step 2: Arbitrary Initialization for (x 0 , θ 0 ). In this step, we remove the assumption of x 0 ∼ π needed in the previous step. We need the following lemma to prove our result. Lemma E.5. Consider two trajectories (x k , θ k ) k≥0 and (x ′ k , θ ′ k ) k≥0 . Suppose that θ 0 = θ ′ 0 , x ′ 0 ∼ π and x 0 is initialized from some arbitrary distribution that satisfies ∥L(x 0 ) -π∥ TV = ϵ. Then for k ≥ τ, we haveW2 (L(x k , θ k ), L(x ′ k , θ ′ k )) ≤ ϵ 4c 2,1 (1 -αµ) k E[∥θ 0 -θ * ∥ 2 ] + 4c s,2 ατ • L 2 µ + 1 1/2 .Proof. We consider the following coupling between two joint processes (x k , θ k ) k≥0 and (x ′ k , θ ′ k ) k≥0 . We first apply the maximal coupling on x 0 and x ′ 0 such thatP(x 0 ̸ = x ′ 0 ) = ∥L(x 0 ) -L(x ′ 0 )∥ TV = ϵ.For the case x 0 = x ′ 0 , we can couple the two Markov chains {x k } k≥0 and {x ′ k } k≥0 such thatx k ≡ x ′ k , ∀k ≥ 0. Under this coupling, we have θ k ≡ θ ′ k , ∀k ≥ 0.For the case x 0 ̸ = x ′ 0 , we let the two processes (x k , θ k ) k≥1 and (x ′ k , θ ′ k ) k≥1 evolve independently. Given the above coupling, we first observe thatW2 (L(x k , θ k ), L(x ′ k , θ ′ k )) = E[ W2 (L(x k , θ k ), L(x ′ k , θ ′ k ))|x 0 = x ′ 0 ]P(x 0 = x ′ 0 ) + E[ W2 (L(x k , θ k ), L(x ′ k , θ ′ k ))|x 0 ̸ = x ′ 0 ]P(x 0 ̸ = x ′ 0 ) = ϵE[ W2 (L(x k , θ k ), L(x ′ k , θ ′ k ))|x 0 ̸ = x ′ 0 ]. The second equality holds since E[ W2 (L(x k , θ k ), L(x k , θ ′ k ))|x 0 = x ′ 0 ] = 0 and P(x 0 ̸ = x ′ 0 ) = ϵ.Next, we note the following upper bound of the Wasserstein distance,W 2 2 (L(x k , θ k ), L(x ′ k , θ ′ k )) = inf E 1{x k ̸ = x ′ k } + ∥θ k -θ ′ k ∥ 2 E.2.1 Proof of Proposition E.2First, we present the following lemma that is similar to [17, Lemma 2.3] . Lemma E.6. For any k 1 < k 2 satisfying α(k 2 -k 1 ) ≤ 1 8L , the following six inequalities hold: ∥θ[1] k2 -θ [2] k2 -θ [1] k1 + θ [2] k1 ∥ ≤ 8αL(k 2 -k 1 )∥θ [1] k2 -θ [2] k2 ∥ ∥θ [1] k2 -θ [2] k2 -θ [1] k1 + θ [2] k1 ∥ ≤ 8αL(k 2 -k 1 )∥θ [1] k1 -θ [2] k1 ∥ ∥θ [1] k2 -θ [1] k1 ∥ ≤ 8αL(k 2 -k 1 ) ∥θ [1] k2 ∥ + 1 ∥θ [1] k2 -θ [1] k1 ∥ ≤ 8αL(k 2 -k 1 ) ∥θ [1] k1 ∥ + 1 ∥θ [2] k2 -θ [2] k1 ∥ ≤ 8αL(k 2 -k 1 ) ∥θ [2] k2 ∥ + 1 ∥θ [2] k2 -θ [2] k1 ∥ ≤ 8αL(k 2 -k 1 ) ∥θ [2] k1 ∥ + 1 .Proof. Consider the coupling given by equation (E.3), by Assumption 2 and 6, we have∥θ [1] k+1 -θ [2] k+1 ∥ -∥θ [1] k -θ [2] k ∥ ≤ ∥θ [1] k+1 -θ [2] k+1 -θ [1] k + θ [2] k ∥ = α∥g(θ [1] k , x k ) + ξ k+1 (θ [1] k ) -g(θ [1] k , x k ) -ξ k+1 (θ [1] k )∥ ≤ 2αL∥θ [1] k -θ [2] k ∥. Given k 1 < k 2 , for ∀t ∈ [k 1 , k 2 ], since 1 + x ≤ e x for ∀x ∈ R, we have ∥θ [1] t -θ [2] t ∥ ≤ t-1 j=k1 (1 + 2αL)∥θ [1] k1 -θ [2] k1 ∥ ≤ exp(2α(k 2 -k 1 )L)∥θ [1] k1 -θ [2] k1 ∥ (i) ≤ (1 + 4α(k 2 -k 1 )L)∥θ [1] k1 -θ [2] k1 ∥ ≤ 2∥θ [1] k1 -θ [2]k1 ∥, where (i) holds for e x ≤ 1 + 2x ∀x ∈ [0, 1 2 ] and α(k 2 -k 1 ) ≤ 1 8L . Then, we have∥θ [1] k2 -θ [2] k2 -θ [1] k1 + θ [2] k1 ∥ ≤ k2-1 t=k1 ∥θ [1] t+1 -θ [2] t+1 -θ [1] t + θ [2] t ∥ ≤ 4α(k 2 -k 1 )L∥θ [1] k1 -θ [2] k1 ∥. Therefore, following α(k 2 -k 1 ) ≤ 1 8L , we have ∥θ [1] k2 -θ [2] k2 -θ [1] k1 + θ [2] k1 ∥ ≤ 4α(k 2 -k 1 )L∥θ [1] k1 -θ [2] k1 ∥ ≤ 4α(k 2 -k 1 )L(∥θ [1] k2 -θ [2] k2 -θ [1] k1 + θ [2] k1 ∥ + ∥θ [1] k2 -θ [2] k2 ∥) ≤ 1 2 ∥θ [1] k2 -θ [2] k2 -θ [1] k1 + θ [2] k1 ∥ + 4α(k 2 -k 1 )L∥θ [1] k2 -θ [2]k2 ∥. Then, by rearranging the terms, we have∥θ [1] k2 -θ [2] k2 -θ [1] k1 + θ [2] k1 ∥ ≤ 8α(k 2 -k 1 )L∥θ [1] k2 -θ [2]k2 ∥, thereby we have proved the first two inequalities of Lemma E.6.By Assumption 2 and 6 and [17, Lemma 2.3], we can prove the last four inequalities and we omit the details here. Now, we are ready to prove Proposition E.2. We start with the following decomposition. By equation (E.3), we first haveE[∥θ [1] k+1 -θ [2] k+1 ∥ 2 ] =E ∥θ [1] k -θ [2] k + α g(θ [1] k , x k ) -g(θ [2] k , x k ) + ξ k+1 (θ [1] k ) -ξ k+1 (θ [2] k ) ∥ 2 =E ∥θ [1] k -θ [2] k ∥ 2 + 2α E ⟨θ [1] k -θ [2] k , g(θ [1] k , x k ) -g(θ [2] k , x k )⟩ T1 + α 2 E ∥g(θ [1] k , x k ) -g(θ [2] k , x k ) + ξ k+1 (θ [1] k ) -ξ k+1 (θ [2] k )∥ 2 T2 .For T 2 , by Assumption 2 and 6, we haveT 2 ≤ 4L 2 E ∥θ [1] k -θ [2] k ∥ 2 .We denote ε(•, x k ) := g(•, x k ) -ḡ(•) to be the noise function. Then, by Assumption 2, we conclude that ε(•, x k ) is 2L-Lipschitz continuous. Therefore, T 1 can be rewritten as:T 1 = E ⟨θ [1] k -θ [2] k , ε(θ [1] k , x k ) -ε(θ [2] k , x k )⟩ + E ⟨θ [1] k -θ [2] k , ḡ(θ [1] k ) -ḡ(θ [2] k )⟩ . For E ⟨θ [1] k -θ [2] k , ḡ(θ [1] k ) -ḡ(θ [2]k )⟩ , by Assumption 3, we haveE ⟨θ [1] k -θ [2] k , ḡ(θ [1] k ) -ḡ(θ [2] k )⟩ ≤ -µE ∥θ [1] k -θ [2] k ∥ 2 . Let F k := σ (θ [1] t , θ [2] t , x t ) | t ≤ k . For E ⟨θ [1] k -θ [2] k , ε(θ [1] k , x k ) -ε(θ [2] k , x k )⟩ , we have E ⟨θ [1] k -θ [2] k , ε(θ [1] k , x k ) -ε(θ [2] k , x k )⟩ =E E ⟨θ [1] k -θ [2] k , ε(θ [1] k , x k ) -ε(θ [2] k , x k )⟩ | F k-τ =E E ⟨θ [1] k-τ -θ [2] k-τ , ε(θ [1] k-τ , x k ) -ε(θ [2] k-τ , x k )⟩ | F k-τ (T 3 ) + E ⟨θ [1] k -θ [2] k -θ [1] k-τ -θ [2] k-τ , ε(θ [1] k-τ , x k ) -ε(θ [2] k-τ , x k )⟩ (T 4 ) + E ⟨θ [1] k -θ [2]k , ε(θ[1] k , x k ) -ε(θ [2] k , x k ) -ε(θ [1] k-τ , x k ) + ε(θ [2] k-τ , x k )⟩ ♠ .We assume ατ ≤ 1 8L . For T 3 , by definition of mixing time τ , we obtainT 3 = E E ⟨θ [1] k-τ -θ [2] k-τ , ε(θ [1] k-τ , x k ) -ε(θ [2] k-τ , x k )⟩ | θ [1] k-τ , θ [2] k-τ , x k-τ ≤ 2αLE ∥θ [1] k-τ -θ [2] k-τ ∥ 2 ≤ 4αLE ∥θ [1] k -θ [2] k ∥ 2 + 4αLE ∥θ [1] k-τ -θ [2] k-τ -θ [1] k + θ [2] k ∥ 2 ≤ (4αL + 256α 3 τ 2 L 3 )E ∥θ [1] k -θ [2] k ∥ 2 ,where the last inequality holds by Lemma E.6.For T 4 , we obtainT 4 ≤ E ∥θ [1] k -θ [2] k -θ [1] k-τ + θ [2] k-τ ∥∥ε(θ [1] k-τ , x k ) -ε(θ [2] k-τ , x k )∥ ≤ 16ατ L 2 E ∥θ [1] k -θ [2] k ∥∥θ [1] k-τ -θ [2] k-τ ∥ ≤ (16ατ L 2 + 128α 2 τ 2 L 3 )E ∥θ [1] k -θ [2] k ∥ 2 .Below, we bound term ♠ by two different Taylor expansions. One the one hand, there existλ 1 , λ 2 ∈ [0, 1] such that h k = λ 1 θ [1] k + (1 -λ 1 )θ [2] k , h k-τ = λ 2 θ [1] k-τ + (1 -λ 2 )θ [2]k-τ and|♠| =|⟨θ [1] k -θ [2] k , ε ′ (h k , x k )(θ [1] k -θ [2] k ) -ε ′ (h k-τ , x k )(θ [1] k-τ -θ [2] k-τ )⟩| =|⟨θ [1] k -θ [2] k , ε ′ (h k-τ , x k )(θ [1] k -θ [2] k -θ [1] k-τ + θ [2] k-τ )⟩ + ⟨θ [1] k -θ [2] k , (ε ′ (h k , x k ) -ε ′ (h k-τ , x k ))(θ [1] k -θ [2] k )⟩| ≤16ατ L 2 ∥θ [1] k -θ [2] k ∥ 2 + |⟨θ [1] k -θ [2] k , (ε ′ (h k , x k ) -ε ′ (h k-τ , x k ))(θ [1] k -θ [2] k )⟩| (E.6) ≤16ατ L 2 ∥θ [1] k -θ [2] k ∥ 2 + 2L∥θ [1] k -θ [2] k ∥ 2 ∥h k -h k-τ ∥ ≤16ατ L 2 ∥θ [1] k -θ [2] k ∥ 2 + 2L∥θ [1] k -θ [2] k ∥ 2 (∥h k -θ [2] k ∥ + ∥θ [2] k -θ [2] k-τ ∥ + ∥θ [2] k-τ -h k-τ ∥) (E.7) ≤16ατ L 2 ∥θ [1] k -θ [2] k ∥ 2 + 2L∥θ [1] k -θ [2] k ∥ 2 (∥θ [1] k -θ [2] k ∥ + ∥θ [2] k -θ [2] k-τ ∥ + ∥θ [1] k-τ -θ [2] k-τ ∥) ≤16ατ L 2 ∥θ [1] k -θ [2] k ∥ 2 + 2L∥θ [1] k -θ [2] k ∥ 2 (2∥θ [1] k -θ [2] k ∥ + 8ατ L(∥θ [2] k ∥ + 1) + 8ατ L∥θ [1] k -θ [2] k ∥) (E.8) ≤16ατ L 2 ∥θ [1] k -θ [2] k ∥ 2 + 2L∥θ [1] k -θ [2] k ∥ 2 3∥θ [1] k -θ [2] k ∥ + 8ατ L(∥θ [2] k ∥ + 1) A := (Bound1) .where we note that choosing the second iterates for triangle inequality in equation (E.7) and choosing θ [2] k for bounding θ [2] k -θ [2] k-τ in equation (E.8) are both symmetric, which implies that we can replace the θ [2] k in term A with arbitrary one in {θ [1] k-τ , θ [1] k , θ [2] k-τ , θ [2] k }.On the other hand, there exist λ1 , λ2 ∈ [0, 1] such that p k = λ1 θ[1] k + (1 -λ1 )θ [1] k-τ , q k = λ2 θ [2] k + (1 -λ2 )θ [2] k-τ and |♠| =|⟨θ [1] k -θ [2] k , ε ′ (p k , x k )(θ [1] k -θ [1] k-τ ) -ε ′ (q k , x k )(θ [2] k -θ [2] k-τ )⟩| =|⟨θ [1] k -θ [2] k , ε ′ (p k , x k )(θ [1] k -θ [2] k -θ [1] k-τ + θ [2] k-τ )⟩ (E.9) + ⟨θ [1] k -θ [2] k , (ε ′ (p k , x k ) -ε ′ (q k , x k ))(θ [2] k -θ [2] k-τ )⟩| ≤16ατ L 2 ∥θ [1] k -θ [2] k ∥ 2 + 2L∥θ [1] k -θ [2] k ∥∥p k -q k ∥∥θ [2] k -θ [2] k-τ ∥ ≤16ατ L 2 ∥θ [1] k -θ [2] k ∥ 2 + 16ατ L 2 ∥θ [1] k -θ [2] k ∥∥p k -q k ∥∥(∥θ [2] k ∥ + 1), (E.10)where adding and subtracting the second iterates in equation (E.7) and choosing θ [2] k for bounding θ [2] k -θ [2] k-τ in equation (E.8) are both symmetric. We have∥p k -q k ∥ = ∥ λ1 (θ [1] k -θ [2] k ) + (1 -λ1 )(θ [1] k-τ -θ [2] k-τ ) + ( λ1 -λ2 )(θ [2] k -θ [2] k-τ )∥ ≤ (2 + 8ατ L)∥θ [1] k -θ [2] k ∥ + 8ατ L(∥θ [2] k ∥ + 1) ≤ 3∥θ [1] k -θ [2]k ∥ + 8ατ L(∥θ [2] k ∥ + 1).|♠| ≤16ατ L 2 ∥θ [1] k -θ [2] k ∥ 2 + 2L∥θ [1] k -θ [2] k ∥ 3∥θ [1] k -θ [2] k ∥ + 8ατ L(∥θ [2] k ∥ + 1) B 8ατ L(∥θ [2] k ∥ + 1) C := (Bound2),and we can replace the θ [2] k in term B and C with arbitrary one in {θ[1] k-τ , θ [1] k , θ [2] k-τ , θ [2]k }. By (Bound1) and (Bound2), we have|♠| ≤ min(Bound1, Bound2) ≤ 16ατ L 2 ∥θ [1] k -θ [2] k ∥ 2 + 48ατ L 2 ∥θ [1] k -θ [2]k ∥ 2 (∥θ [2] k ∥ + 1). Below, we discuss the upper bound for |♠| by three cases.Case 1: If one of {θ [1] k-τ , θ [1] k , θ [2] k-τ , θ [2] k } has norm that is less than 4δ ατ L 2 , where we define δ ατ L 2 := δ(ατ L 2 ), without loss of generality, we assume ∥θ [2] k ∥ ≤ 4δ ατ L 2 . Multiply 1(∥θ [2] k ∥ ≤ 4δ ατ L 2 ) to both sides of the above inequality, and we have|♠|1(∥θ [2] k ∥ ≤ 4δ ατ L 2 ) ≤ 16ατ L 2 ∥θ [1] k -θ [2] k ∥ 2 + 48ατ L 2 ∥θ [1] k -θ [2]k ∥ 2 (∥θ [2] k ∥ + 1) 1(∥θ [2] k ∥ ≤ 4δ ατ L 2 )≤(64 + 192δ ατ L 2 )ατ L 2 ∥θ [1] k -θ [2] k ∥ 2 ,where we can replace the θ [2] k in the left hand of the above inequality with an arbitrary one in {θ [1] k-τ , θ [1] k , θ [2] k-τ , θ [2] k } by the similar argument under (Bound1) and (Bound2).Case 2: If ∥θ [1] k ∥ < 2∥θ [1] k -θ [2] k ∥, we obtain ∥θ [2] k ∥ ≤ ∥θ [1] k -θ [2] k ∥ + ∥θ [1] k ∥ ≤ 3∥θ [1] k -θ [2] k ∥.Therefore, by definition of ♠ and Lemma E.6, we obtain|♠|1(∥θ [1] k ∥ ≤ 2∥θ [1] k -θ [2] k ∥) ≤ 2L∥θ [1] k -θ [2] k ∥(∥θ [1] k -θ [1] k-τ ∥ + ∥θ [2] k -θ [2] k-τ ∥)1(∥θ [1] k ∥ ≤ 2∥θ [1] k -θ [2] k ∥) ≤ 16ατ L 2 ∥θ [1] k -θ [2] k ∥(∥θ [1] k ∥ + ∥θ [2] k ∥)1(∥θ [1] k ∥ ≤ 2∥θ [1] k -θ [2] k ∥) ≤ 80ατ L 2 ∥θ [1] k -θ [2] k ∥. Case 3: If all four variables in {θ [1] k-τ , θ [1] k , θ [2] k-τ , θ [2]k } have a norm larger than 4δ ατ L 2 , and ∥θ[1] k ∥ ≥ 2∥θ [1] k -θ [2]k ∥, we obtain∥h k ∥ = ∥λ 1 θ [1] k + (1 -λ 1 )θ [2] k ∥ ≥ ∥θ [1] k ∥ -∥θ [1] k -θ [2] k ∥ ≥ ∥θ [1] k ∥ 2 ≥ 2δ ατ L 2 ∥h k-τ ∥ = ∥λ 2 θ [1] k-τ + (1 -λ 2 )θ [2] k-τ ∥ ≥ ∥θ [1] k-τ ∥ -∥θ [1] k-τ -θ [2] k-τ ∥ (i) ≥ δ ατ L 2 -(1 + 8ατ L)∥θ [1] k -θ [2] k ∥ ≥ δ ατ L 2 ,where (i) holds by choosing ατ ≤ 1 16L . Therefore, we have ∥h k ∥ ≥ δ ατ L 2 and ∥h k-τ ∥ ≥ δ ατ L 2 . Therefore, by equation (E.6) and Assumption 6, we obtain|♠|1(Case 3) ≤16ατ L 2 ∥θ [1] k -θ [2] k ∥ 2 + ∥θ [1] k -θ [2] k ∥ 2 (∥ε ′ (h k , x k ) -G(x k )∥ + ∥ε ′ (h k-τ , x k ) -G(x k )∥) ≤18ατ L 2 ∥θ [1] k -θ [2] k ∥ 2 .Therefore, we obtain|E[♠]| ≤ E[|♠|] ≤E[|♠|1(∥θ [1] k-τ ∥ ≤ δ ατ L 2 )] + E[|♠|1(∥θ [1] k ∥ ≤ δ ατ L 2 )] + E[|♠|1(∥θ [2] k-τ ∥ ≤ δ ατ L 2 )] + E[|♠|1(∥θ [2] k ∥ ≤ δ ατ L 2 )] + E[|♠|1(∥θ [1] k ∥ ≤ 2∥θ [1] k -θ [2] k ∥)] + E[|♠|1(Case 3)] ≤(354 + 768δ ατ L 2 )ατ L 2 E[∥θ [1] k -θ [2] k ∥ 2 ].By Assumption 6, there exists κ µ > 0 such that ϵδ(ϵ) ≤ µ 3072 for ∀ϵ ≤ κ µ . By the above bounds, when ατ ≤ min( cµ L 2 , κµ L 2 ), where we specify the constant c < 1 later, we obtainE[∥θ [1] k+1 -θ [2] k+1 ∥ 2 ] ≤ 1 + 2α(-µ + 4αL + 256α 3 τ 2 L 3 + 128α 2 τ 2 L 3 + 784ατ L 2 + 768δ ατ L 2 ατ L 2 ) + 4α 2 L 2 E[∥θ [1] k -θ [2] k ∥ 2 ] = 1 + α(-2µ + 1568ατ L 2 + 256α 2 τ 2 L 3 + 1536δ ατ L 2 ατ L 2 ) + α 2 (8L + 512α 2 τ 2 L 3 + 4L 2 ) E[∥θ [1] k -θ [2] k ∥ 2 ] ≤ 1 + αµ(-2 + 1568c + 256c + 1 2 + 524c) E[∥θ [1] k -θ [2] k ∥ 2 ] ≤(1 -µα)E[∥θ [1] k -θ [2] k ∥ 2 ],where the last inequality holds by choosing c = 1 4696 . Therefore, for k ≥ τ and ατ ≤ min( µ 906L 2 , κµ L 2 ) we haveE[∥θ [1] k -θ [2] k ∥ 2 ] ≤ (1 -µα) k-τ E[∥θ [1] τ -θ [2] τ ∥ 2 ] ≤ 2(1 -µα) k-τ E[∥θ [1] 0 -θ [2]0 ∥ 2 + ∥θ [1] τ -θ [2] τ -θ[1] 0 + θ [2] 0 ∥ 2 ] ≤ 2(1 + 8ατ L)(1 -µα) k-τ E[∥θ [1] 0 -θ [2] 0 ∥ 2 ] ≤ 4(1 -µα) k-τ E[∥θ [1] 0 -θ [2] 0 ∥ 2 ].Now, we specialize the proof of Theorem E.1 to the projected SA iterates.We consider the same coupling:θ [1] k+ 1 2 = θ [1] k + α g(θ [1] k , x k ) + ξ k+1 (θ [1] k ) , θ [1] k+1 = Π B(β) θ [1] k+ 1 2 , θ [2] k+ 1 2 = θ [2] k + α g(θ [2] k , x k ) + ξ k+1 (θ [2] k ) , θ [2] k+1 = Π B(β) θ [2] k+ 1 2 .We first need to verify that Proposition E.2 holds for projected SA. By the non-expansion property of Π B(β) with respect to ∥ • ∥, we obtain the following inequality:∥θ [1] k+1 -θ [2] k+1 ∥ -∥θ [1] k -θ [2] k ∥ ≤ ∥θ [1] k+ 1 2 -θ [2] k+ 1 2 ∥ -∥θ [1] k -θ [2]k ∥, which implies Lemma E.6 still holds for projected SA. For Proposition E.2, we can notice that the iterates of projected SA will always satisfy Case 1 with a finite bound β and we do not need to discuss Cases 2 and 3.Therefore, when ατ ≤ min( cµ L 2 , 1 8L ), where we specify the constant c < 1 later, we obtainE[∥θ [1] k+1 -θ [2] k+1 ∥ 2 ] ≤ 1 + α(-2µ + (160 + 96β))ατ L 2 + 256α 2 τ 2 L 3 ) + α 2 (8L + 512α 2 τ 2 L 3 + 4L 2 ) E[∥θ [1] k -θ [2] k ∥ 2 ] ≤ (1 + αµ(-2 + (160 + 96β)c + 256c + 524c)) E[∥θ [1] k -θ [2] k ∥ 2 ] =(1 -µα)E[∥θ [1] k -θ [2] k ∥ 2 ],where we set c =1 940+96β . Therefore, ∀k ≥ τ and ατ ≤ µ (940+96β)L 2 , we have E[∥θ [1] k -θ [2] k ∥ 2 ] ≤ 4(1 -µα) k-τ E[∥θ [1] 0 -θ [2] 0 ∥ 2 ].Then, still by the non-expansion of Π B(β) with respect to ∥ • ∥, the rest of the proof simply follows the same proof for non-projected SA. As such, we have proven Theorem 4.1.In this section, we present the proof of Corollary 4.4.Recall that by Theorem 4.3, we obtain for k ≥ 2τ ,W 2 2 L(θ k ), ν α ≤ (1 -αµ) k • s(θ 0 , θ * , µ, L, R).By [63, Theorem 4.1] , there exists a coupling between θ k and θ ∞ such thatW 2 2 (L(θ k ), ν α ) = E[∥θ k -θ ∞ ∥ 2 ]. Applying Jensen's inequality twice, we obtain that∥E[θ k -θ ∞ ]∥ 2 ≤ (E[∥θ k -θ ∞ ∥]) 2 ≤ E[∥θ k -θ ∞ ∥ 2 ] ≤ (1 -αµ) k • s(θ 0 , θ * , µ, L, R).We thus have for all k ≥ 2τ ,∥E[θ k ] -E[θ ∞ ]∥ ≤ E[∥θ k -θ ∞ ∥] ≤ (1 -αµ) k/2 • s ′ (θ 0 , θ * , µ, L, R).For the second moment, we first note that∥E[θ k θ ⊤ k ] -E[θ ∞ θ ⊤ ∞ ]∥ = ∥E[(θ k -θ ∞ )(θ k -θ ∞ ) ⊤ ] + E[θ ∞ (θ k -θ ∞ ) ⊤ ] + E[(θ k -θ ∞ )θ ⊤ ∞ ]∥ ≤ ∥E[(θ k -θ ∞ )(θ k -θ ∞ ) ⊤ ]∥ + ∥E[θ ∞ (θ k -θ ∞ ) ⊤ ]∥ + ∥E[(θ k -θ ∞ )θ ⊤ ∞ ]∥ ≤ E[∥θ k -θ ∞ ∥ 2 ]∥ + 2E[∥θ ∞ ∥∥θ k -θ ∞ ∥]∥ ≤ E[∥θ k -θ ∞ ∥ 2 ]∥ + 2 E[∥θ ∞ ∥ 2 ∥θ k -θ ∞ ∥ 2 ]∥, (F.1)where we apply Cauchy-Schwarz to obtain the last inequality.Meanwhile, we haveE ∥θ k -θ ∞ ∥ 2 ≤ (1 -αµ) k • s(θ 0 , θ * , µ, L, R) and E ∥θ ∞ ∥ 2 = O(1).Substituting the above bounds into the right-hand side of inequality (F.1) yieldsE θ k θ ⊤ k -E θ ∞ θ ⊤ ∞ ≤ (1 -αµ) k/2 • s ′′ (θ 0 , θ * , µ, L, R).G Proof of Corollary 4.5In this section, we prove the CLT result.Proof. Consider the following centered test function h :X × R d → R d defined as h(x, θ) = θ -E[θ ∞ ].To prove that the CLT for function h, we need to verify the Maxwell-Woodroofe condition [49] , i.e.,∞ n=1 n -3/2 n-1 t=0 Q t h L 2 (ν) < ∞,where Q denotes the transition kernel of the joint Markov chain. If we can show the followingn-1 t=0 Q t h L 2 (ν) = O(n r ) (G.1)with r ∈ [0, 1/2), then the Maxwell-Woodroofe condition is verified, as∞ n=1 n -3/2 n-1 t=0 Q t h L 2 (ν) = ∞ n=1 n -3/2 O(n r ) < ∞.We now proceed to prove the desired order in (G.1). For sufficiently large n ≥ 2τ α , we observen-1 t=0 Q t h L 2 (ν) = E ν n-1 t=0 Q t h 2 ≤ n-1 t=0 E ν ∥Q t h∥ 2 = 2τα-1 t=0 E ν ∥Q t h∥ 2 T1 + n-1 t=2τα E ν ∥Q t h∥ 2.We now show that both terms T 1 and T 2 are of order O(1) with respect to the parameter n.For T 1 , since Q is a transition kernel, so its ∥Q∥ L 2 (µ) operator norm equals to 1. Hence, T 1 can be upper bounded asT 1 ≤ τ α E μ[∥h(θ, x)∥ 2 2 ] = τ α Tr(Var(θ ∞ )) < C 1 ,where the last inequality follows from Tr(Var(θ ∞ )) ≤ ατ established in (H.4) and ατ 2 α → 0 by Definition 2.1.Before proceeding to analyze the summation in T 2 , we first recall (E.2), that for t ≥ 2τ α , W2 (L(x t , θ t ), ν) = O((1 -αµ) t/2 ), which holds for any (x, θ) ∈ X × R d . Hence, by the property of Wasserstein distance [63] , there always exists a coupling that attains the optimality, i.e.,E Γ((xt,θt),ν) ∥θ t -θ ′ ∥ 2 2 + δ 0 (x t ̸ = x ′ ) = O((1 -αµ) t). Making use of this relationship, we can therefore bound T 2 ,T 2 = n-1 t=τα E ν ∥Q t h∥ 2 ≤ ∞ t=τα E ν ∥Q t h∥ 2 = O 1 1 -(1 -αµ) 1/2 = O(1),where the last O(•) is asymptotic in n.Combining the analysis of T 1 and T 2 , we have shown the desired order in (G.1). Therefore, the Maxwell-Woodroofe condition has been verified and we establish the CLT for averaged nonlinear iterates with constant stepsize and Markovian data.When assuming the perturbed continuous noise condition in Assumption 5, one takes the alternative route to prove weak convergence. This is achieved by establishing the satisfaction of both a minorization condition and a drift condition. In this section, we prove the weak convergence result in Theorem 4.3 by following this alternative approach. The subsequent corollaries of weak convergence, namely the non-asymptotic convergence rate in Corollary 4.4 and the Central Limit Theorem (CLT) in Corollary 4.5, also hold, and we will provide the proofs for these results as well.In this section, we prove the weak convergence under Assumption 5(a). The proof consists of two major steps. Firstly, built upon the MSE convergence established in Proposition 4.2, we derive a multi-state drift condition. Subsequently, we show that under the minorization condition, the Markov chain is (x k , θ k ) k≥0 is φ-irreducible. Then, follow [50, Theorem 19.1.3] , we can conclude that the Markov chain (x k , θ k ) k≥0 is geometrically ergodic.For completeness, we include the Theorem 19.1.3 from [50] below. Theorem H.1. Suppose that Φ is a φ-irreducible chain on X , and let n(x) be a measurable function from X → Z + . The chain is geometrically ergodic if it is aperiodic and there exists some petite set C, a nonnegative function V ≥ 1 and bounded on C, and positive constants λ < 1 and b satisfyingP n(x) (x, dy)V (y) ≤ λ n(x) [V (x) + b1 C (x)]. (H.1)We note that the function n(x) can be interpreted as the number of steps we must wait, starting from any x, for the drift to become negative.Step 1: Deriving the Drift Condition Given the iteration stepθ k+1 = θ k + α(g(θ k , x k ) + ξ t+1 (θ k )),we have already shown the following convergence rate on the MSE in Proposition 4.2, thatE[∥θ k -θ * ∥ 2 ] ≤ c 2,1 (1 -αµ) k ∥θ 0 -θ * ∥ 2 + c 2,2 ατ α L 2 µ .Inspired by the MSE convergence bound, we the Lyapunov function V :X × R d → [1, ∞] V (x, θ) = ∥θ -θ * ∥ 2 + 1. (H.2)Therefore, the major goal in this step is to obtain the desired drift condition as shown in (H.1).Therefore, from the above MSE convergence rate, we first obtain thatQ k V (x 0 , θ 0 ) ≤ c 2,1 (1 -αµ) k V (x 0 , θ 0 ) + c 2,2 L 2 µ ατ + 1. (H.3) Consider k = min{t ≥ 0 : c 2,1 (1 -αµ) t < 1}. Then, set η = 8(1 -αµ) k , β = 1-η 2 , and m = c 2,2 L 2µ ατ + 1, and we consider the following bounded sublevel set,C Θ = {θ : (∥θ -θ * ∥ 2 + 1) ≤ m/β}.From (H.3), we derive thatQ k V (x 0 , θ 0 ) -V (x 0 , θ 0 ) ≤ -βV (x 0 , θ 0 ) + m1 C (x 0 , θ 0 ),where C = X × C Θ . Rewriting the above Lyapunov drift condition, with b = m/(1 -β), we obtainQ k V (x 0 , θ 0 ) ≤ (1 -β) V (x 0 , θ 0 ) + b1 C (x 0 , θ 0 ) . Setting λ such that λ k = 1 -β, we have Q k V (x 0 , θ 0 ) ≤ λ k V (x 0 , θ 0 ) + b1 C (x 0 , θ 0 ) ,which gives the desired multi-step drift condition.Step 2: Proving the Minorization Condition Now that we have established the desired multi-step drift condition, it remains for us to show that C is accessible, small, and aperiodic.Under this setup of ξ t (θ) in Assumption 5(a), it is straightforward to verify the accessibility of C. For any (x, θ) ∈ X × R d , we haveQ((x, θ), C) = P((x ′ , θ ′ ) ∈ X × C Θ |(x, θ)) = P(θ ′ ∈ C Θ |(x, θ)) ≥ θ ′ ∈CΘ 1 α d p θ θ ′ -θ α -g(x, θ) dθ ′ > 0.As such, we have shown that C is accessible.Assuming C is small, we can directly conclude aperiodicity following the definition of period of an accessible small set, d(C) = g.c.d. n ∈ N * : inf x∈C P n (x, C) > 0 .Therefore, what remains to show is that C is small. For (x, θ) ∈ C and Ā ∈ B(X ) × B(R d ), we define the following projection sets,A x = {θ ∈ R d |(x, θ) ∈ Ā} and A θ = {x ∈ X |(x, θ) ∈ Ā}.Therefore, 1) , θ (m-1) ))Q m ((x, θ), Ā) = {(x (k) ,θ (k) )} m-1 k=1 ∈(X ×R d ) m-1 P((x m , θ m ) ∈ Ā|(x m-1 , θ m-1 ) = (x (m-• • • P((x 1 , θ 1 ) = d(x (1) , θ (1) )|(x 0 , θ 0 ) = (x, θ)) 1) , θ (m-1) ))≥ {(x (k) ,θ (k) )} m-1 k=1 ∈(X ×CΘ) m-1 P((x m , θ m ) ∈ Ā|(x m-1 , θ m-1 ) = (x (m-• • • P((x 1 , θ 1 ) = d(x (1) , θ (1) )|(x 0 , θ 0 ) = (x, θ)) 1) , θ (m-1) )) 2) , θ (m-2) ))= {(x (k) ,θ (k) )} m-1 k=1 ∈(X ×CΘ) m-1 x ′ ∈X P(x m = dx ′ |x m-1 = x (m-1) )P(θ m ∈ A x ′ |(x m-1 , θ m-1 ) = (x (m-P(x m-1 = dx (m-1) |x m-2 x (m-2) )P(θ m-1 = dθ (m-1) |(x m-2 , θ m-2 ) = (x (m-• • • P(x 1 = dx (1) |x 0 = x)P(θ 1 = dθ (1) |(x 0 , θ 0 ) = (x, θ)) 1) , θ (m-1) ))= x ′ ∈X {x (k) } m-1 k=1 ∈X m-1 P(x m = dx ′ |x m-1 = x (m-1) )P(x m-1 = dx (m-1) |x m-2 = x (m-2) ) • • • P(x 1 = dx (1) |x 0 = x) {θ (k) } m-1 k=1 ∈C m-1 Θ P(θ m ∈ A x ′ |(x m-1 , θ m-1 ) = (x (m-P(θ m-1 = dθ (m-1) |(x m-2 , θ m-2 ) = (x (m-2) , θ (m-2) )) • • • P(θ 1 = dθ (1) |(x 0 , θ 0 ) = (x, θ)) .Next, for (x, θ) ∈ C, we observe thatP(θ k = dθ ′ |(x k-1 , θ k-1 ) = (x, θ)) = P ξ k (θ) = d θ ′ -θ α -g(x, θ) |(x k , θ k ) = (x, θ) ≥ 1 α d p θ θ ′ -θ α -g(x, θ) dθ ′ ≥ 1 α d inf θ∈CΘ p θ θ ′ -θ α -g(x, θ) dθ ′ .We next recall the linear growth assumption in Assumption 2 that ∥g(x, θ)∥ ≤ L(∥θ∥ + 1). Hence, given θ ∈ C Θ , ∀x ∈ X , we have∥θ + αg(x, θ)∥ ≤ ∥θ∥ + α∥g(x, θ)∥ ≤ (1 + αL)(∥θ∥ + 1) ≤ (1 + αL)( M/β + ∥θ * ∥ + 1) ≤ B,for some bounded value B.We now define the measure ς † on R d as ς† (A) = 1 α d inf ∥z∥≤B θ ′ ∈A∩CΘ inf θ∈CΘ p θ θ ′ -z α dθ ′ .Hence, it is easy to see that ς † (C c Θ ) = 0. Moreover, we note the following property of measure ς † . Claim 1. For A ⊆ C Θ and λ(A) > 0, ς † (A) > 0, where λ denotes the Lebesgue measure.We delay the proof to the end. Taking the claim as true, we derive that (1) , θ (1) ))P(θ 1 = dθ (1) |(x 0 , θ 0 ) = (x, θ)){θ (k) } m-1 k=1 ∈C m-1 Θ P(θ m ∈ A x ′ |(x m-1 , θ m-1 ) = (x (m-1) , θ (m-1) )) P(θ m-1 = dθ (m-1) |(x m-2 , θ m-2 ) = (x (m-2) , θ (m-2) )) • • • P(θ 1 = dθ (1) |(x 0 , θ 0 ) = (x, θ)) = θ (m) ∈A x ′ θ (m-1) ∈CΘ P(θ m dθ (m) |(x m-1 , θ m-1 ) = (x (m-1) , θ (m-1) )) θ (m-2) ∈CΘ P(θ m-1 = dθ (m-1) |(x m-2 , θ m-2 ) = (x (m-2) , θ (m-2) )) • • • θ (1) ∈CΘ P(θ 2 = dθ (2) |(x 1 , θ 1 ) = (x≥ θ (m) ∈A x ′ θ (m-1) ∈CΘ 1 α d inf θ∈CΘ p θ θ (m) -θ (m-1) α -g(x (m-1) , θ (m-1) ) dθ (m) θ (m-2) ∈CΘ P(θ m-1 = dθ (m-1) |(x m-2 , θ m-2 ) = (x (m-2) , θ (m-2) ))Therefore, by combining all the analyses, we obtain 1) , θ (m-1) ))Q m ((x, θ), Ā) ≥ x ′ ∈X {x (k) } m-1 k=1 ∈X m-1 P(x m = dx ′ |x m-1 = x (m-1) )P(x m-1 = dx (m-1) |x m-2 = x (m-2) ) • • • P(x 1 = dx (1) |x 0 = x) {θ (k) } m-1 k=1 ∈C m-1 Θ P(θ m ∈ A x ′ |(x m-1 , θ m-1 ) = (x (m-P(θ m-1 = dθ (m-1) |(x m-2 , θ m-2 ) = (x (m-2) , θ (m-2) )) • • • P(θ 1 = dθ (1) |(x 0 , θ 0 ) = (x, θ)) ≥ ς † (C Θ ) m-1 x ′ ∈X ς † (A x ′ ) {x (k) } m-1 k=1 ∈X m-1 P(x m = dx ′ |x m-1 = x (m-1) ) • • • P(x 1 = dx (1) |x 0 = x) ≥ ς † (C Θ ) m-1 x ′ ∈X ς † (A x ′ )ϕ(dx ′ ) = ζ • (ϕ × ς † )( Ā),where ζ = ς † (C Θ ) m-1 and ϕ × ς † being the unique induced product measure on X × R d .As such, we have proven that C is (m, ϕ × ς † )-small, and hence (δ m , ϕ × ς † )-petite, and subsequently shown that (x k , θ k ) k≥0 is geometrically ergodic.By [50, Theorem 16 .0.1 (iv)], we can further conclude that the ergodic (θ t , x t ) t≥0 is also V -uniformly ergodic with the same V as defined in (H.2). Therefore, we have the following convergence rate in the V -norm L(x k , θ k ), -ν α V ≤ κρ k , where κ and ρ implicitly depend on the stepsize α.Lastly, we provide the proof of Claim 1.Proof. We first recall that for any set B ⊆ R d such that λ(B) > 0, where λ refers to the Lebesgue measure, then for θ ∈ C Θ , we know that t∈B p θ (t)dt ≥ t∈B inf θ∈CΘ p θ (t)dt > 0. Moreover, for a given (translation) z ∈ R d and θ ∈ C Θ , we havet∈B p θ (t -z)dt = t ′ ∈B-z p θ (t ′ )dt ′ ≥ t ′ ∈B-z inf θ∈C p θ (t ′ )dt ′ > 0.Following the properties stated above, we define h(z) = t∈B inf θ∈C p θ (t -z)dt, and it is easy to verify that h(z) is a continuous function. Hence, for a bounded set D, we have inf z∈D h(z) > 0.Subsequently, we define the measure ς † induced by h and we verify thatς † (A) = inf ∥z∥≤B t∈A∩CΘ inf θ∈C p θ (t -z)dt > 0.H.2 Proof of Corollary 4.4As shown in the previous section, the joint process (x k , θ k ) k≥0 is V -uniformly ergodic and hence also exhibits a geometric non-asymptotic convergence rate under the V -weighted norm. Subsequently, this corresponds to a version of Corollary 4.4 resulting in a different set of convergence rate coefficients. Thus, we restate Corollary 4.4 in the context of the minorization setting and provide the proof below. Corollary H.2 (Non-Asymptotic Convergence Rate). For any initialization of θ 0 ∈ R d , under the setting of Theorem 4.3, we haveE[θ k ]-E[θ (α) ∞ ] ≤ κ•ρ k •s ′ (θ 0 , L, µ), and E[θ k θ ⊤ k ]-E[θ (α) ∞ (θ (α) ∞ ) ⊤ ] ≤ κ•ρ k •s ′′ (θ 0 , L, µ), where κ and ρ are defined in (4.1) and implicitly depend on α.Proof. For V -uniformly ergodic Markov chain (x k , θ k ) k≥0 , when functions f : X × R d → R d is dominated by the Lyapunov function, i.e., ∥f ∥ ≤ V , it enjoys the following convergence property,∥Q n f (x 0 , θ 0 ) -πf ∥ ≤ κρ n V (θ 0 ).Consider test function f (θ, x) = θ -θ * . It is easy to see that ∥f ∥ ≤ V . Hence, we obtain∥E[θ n ] -E[θ ∞ ]∥ = ∥Q n f (x, θ) -πf ∥ ≤ κρ n V (θ 0 ). Next, consider f ′ (x) = (θ -θ * )(θ -θ * ) ⊤ . Clearly, ∥f ′ ∥ ≤ V . Therefore, ∥E[(θ t -θ * )(θ t -θ * ) ⊤ ] -E[(θ ∞ -θ * )(θ ∞ -θ * ) ⊤ ]∥ ≤ κρ n V (θ 0 ).For the LHS, we have∥E[(θ t -θ * )(θ t -θ * ) ⊤ ] -E[(θ ∞ -θ * )(θ ∞ -θ * ) ⊤ ]∥ = ∥E[θ t θ ⊤ t ] -E[θ ∞ θ ⊤ ∞ ] -E[θ t -θ ∞ ](θ * ) ⊤ -θ * E[(θ t -θ ∞ ) ⊤ ]∥ ≥ ∥E[θ t θ ⊤ t ] -E[θ ∞ θ ⊤ ∞ ]∥ -2∥E[θ t -θ ∞ ]∥∥θ * ∥. Subsequently, ∥E[θ t θ ⊤ t ] -E[θ ∞ θ ⊤ ∞ ]∥ ≤ (2∥θ * ∥ + 1)κρ n V (θ 0 ).The above results imply the convergence of the first two moments. Moreover, we conclude thatVar(θ ∞ ) = Var(θ ∞ -θ * ) ≤ E[∥θ ∞ -θ * ∥ 2 ] = lim t→∞ E[∥θ t -θ * ∥ 2 ] ≲ ατ.Additionally,E[∥θ ∞ -θ * ∥] 2 ≤ E[∥θ ∞ -θ * ∥ 2 ] ≲ ατ. (H.4)H.3 Proof of Corollary 4.5After establishing the V -uniform ergodicity the joint process (x k , θ k ) k≥0 , the central limit theorem for averaged iterates follows as a straightforward consequence.Proof. For any test function h : X × R d → R d that satisfies ∥h∥ 2 ≤ V , by Theorem 17.0.1 in [50] , it has the following CLT results,1 √ k k-1 t=0 h t -E[h ∞ ] ⇒ N (0, Σ (a) ).Therefore, consider h(x, θ) = θ -θ * , it is easy to see that ∥h∥ 2 ≤ V , and hence we naturally obtain the desired CLT result, 1√ k k-1 t=0 θ t -E[θ ∞ ] ⇒ N (0, Σ (a) ) as k → ∞.I Proof of Theorem 4.6We now provide the proof of Theorem 4.6 on characterizing the asymptotic bias of nonlinear SA.The proof utilizes the basic adjoint relationship (BAR) approach to study the stationary distribution E ν [(P -I)h(θ, x)] = 0, via carefully designed test functions h. We refer readers to [33] for the derivation of the following properties of Markovian SA at stationarity,E[1{θ ∞ ∈ S} | x ∞ ](x) = E[1{θ ∞+1 ∈ S} | x ∞+1 ](x), (I.1) E[θ ∞ | x ∞ ](x) = E[θ ∞+1 | x ∞+1 ](x), (I.2) E[(θ ∞ ) ⊗2 | x ∞ ](x) = E[(θ ∞+1 ) ⊗2 | x ∞+1 ](x). (I.3) Following the Borel state space assumption in 1, E[1{θ ∞ ∈ •}|x ∞ = x] induce a regular conditional probability measure, which we denote as ν(•, x ∞ = x), and hence (I.1) can be reformulated as ν(θ ∞ ∈ S, x ∞ = x) = ν(θ ∞+1 ∈ S, x ∞+1 = x).Before proceeding to the proof, we introduce the following shorthands and notations. For x ∈ X ,z i (x) := E[(θ ∞ -θ * ) ⊗i |x ∞ = x],δ i (x) := z i (x) -πz i , Following the differentiability assumption of g in Assumption 3, and we can apply Taylor expansion to g and we note the following notation on residuals.g(θ, x) = g(θ * , x) + g ′ (θ * , x)(θ -θ * ) + 1 2 g ′′ (θ * , x)(θ -θ * ) ⊗2 + R 3 (θ, x) (I.4) = g(θ * , x) + g ′ (θ * , x)(θ -θ * ) + R 2 (θ, x).(I.5) By Assumption 3 and results from Proposition 4.2 and 4.2, we note that the residual R n (θ, x) satisfiessup x∈X ,θ∈R d ∥R n (θ, x)∥/∥θ -θ * ∥ n < +∞.Hence, we have∥R n (θ, x ∞ )∥ L 2 (π) ≲ E[∥θ ∞ -θ * ∥ n ] = O((ατ ) n/2 ), n = 2, 3, 4. Lastly, we denote ḡ(θ) := E x∼π [g(θ, x)], ḡ2 (θ) := E x∼π [(g(θ, x)) ⊗2 ] ḡ(1) (θ) := E x∼π [g ′ (θ, x)], ḡ(1) 2 (θ) := E x∼π [(g ′ (θ, x)) ⊗2 ], ḡ(2) (θ) := E x∼π [g ′′ (θ, x)].We are now ready to present our proof. The proof consists of two major steps. For the complexity of this problem, we first set aside the projection constraint and focus on the BAR analysis, and we shall present the asymptotic bias characterization without the projection step. The analysis in this step thus shall work for the minorization proof technique as well. Then, in the second step, we elaborate on the impact brought along by the projection analysis and conclude our proof.Substituting the above analyses of the terms and consolidating the terms, we obtain-(ḡ (θ * ) ⊗ I + I ⊗ ḡ(1) (θ * ))E[(θ ∞ -θ * ) ⊗2 ] = αḡ 2 (θ * ) + αE[(ξ ∞+1 (θ * )) ⊗2 ] + αḡ (1) 2 (θ * )E[(θ ∞ -θ * ) ⊗2 ] + E[g(θ * , x ∞ ) ⊗ δ 1 (x ∞ )] + E[δ 1 (x ∞ ) ⊗ g(θ * , x ∞ )] + E[g ′ (θ * , x ∞ )δ 2 (x ∞ )] + E[δ 2 (x ∞ )g ′ (θ * , x ∞ )] + αE[g(θ * , x ∞ ) ⊗ (g ′ (θ * , x ∞ )δ 1 (x ∞ ))] + αE[g(θ * , x ∞ ) ⊗ (g ′ (θ * , x ∞ )E[θ ∞ -x ∞ ])] + αE[(g ′ (θ * , x ∞ )δ 1 (x ∞ )) ⊗ g(θ * , x ∞ )] + αE[(g ′ (θ * , x ∞ )E[θ ∞ -x ∞ ]) ⊗ g(θ * , x ∞ )] + αE[g ′ (θ * , x ∞ )δ 2 (x ∞ )g ′ (θ * , x ∞ )] + E[R 2 (θ ∞ , x ∞ ) ⊗ (θ ∞ -θ * )] + E[(θ ∞ -θ * ) ⊗ R 2 (θ ∞ , x ∞ )] + αC ′ (θ * )E[θ ∞ -θ * ] + O(α 2 τ ).(I.17) By far, we observe that the remaining terms all contain δ 1 and δ 2 . Therefore, we conclude our analysis of E[(θ ∞ -θ * ) ⊗2 ] at this step and leave the analysis of δ 1 and δ 2 to the next section.In this section, we analyze δ 1 and δ 2 .Analysis of δ 1 . Starting with δ 1 , we first consider the following recursive relationship induced by (I.2).E[θ ∞+1 -θ * |x ∞+1 = s] = R P * (s, ds ′ )E[θ ∞+1 -θ * |x ∞+1 = s, x ∞ = s ′ ] (i) = R P * (s, ds ′ )E[θ ∞ -θ * + αg(θ ∞ , x ∞ )|x ∞ = s ′ ] (ii) = R P * (s, ds ′ )E θ ∞ -θ * + α g(θ * , x ∞ ) + g ′ (θ * , x ∞ )(θ ∞ -θ * ) + R 2 (θ ∞ , x ∞ ) |x ∞ = s ′ = R P * (s, ds ′ )E θ ∞ -θ * |x ∞ = s ′ + α R P * (s, ds ′ ) g(θ * , s ′ ) + g ′ (θ * , s ′ )E[θ ∞ -θ * |x ∞ = s ′ ] + E[R 2 (θ ∞ , x ∞ )|x ∞ = s ′ ] ,where in (i) we make use of the update rule in (2.1), E[ξ ∞+1 (θ ∞ )] = 0 and conditional independence x ∞+1 ⊥ ⊥ θ ∞ |x ∞ . Next, we substitute the Taylor expansion (I.5) to obtain (ii).Writing with notation shorthands z and δ, we havez 1 (s) = R P * (s, ds ′ )z 1 (s ′ ) + α R P * (s, ds ′ ) g(θ * , s ′ ) + g ′ (θ * , s ′ )z 1 (s ′ ) + E[R 2 (θ ∞ , s ′ )|x ∞ = s ′ ] . (I.18)If we apply π to both sides of (I.18), we obtainπ(ds)P * (s, ds ′ ) g(θ * , s ′ ) + g ′ (θ * , s ′ )z 1 (s ′ ) + E[R 2 (θ ∞ , s ′ )|x ∞ = s ′ ] = 0.Analyzing the three terms closely, we observe that π(ds) P * (s, ds ′ )g(θ * , s ′ ) = g(θ * , s ′ ) π(ds)P * (s, ds ′ )=π(ds ′ ) = ḡ(θ * ) = 0 π(ds) P * (s, ds ′ )g ′ (θ * , s ′ )z 1 (s ′ ) = E[g ′ (θ * , x ∞ )z 1 (x ∞ )] π(ds) P * (s, ds ′ )E[R 2 (θ ∞ , s ′ )|x ∞ = s ′ ] = E[R 2 ∞ , x ∞ )],and hence we first obtainE[g ′ (θ * , x ∞ )z 1 (x ∞ )] + E[R 2 (θ ∞ , x ∞ )] = 0. (I.19)We now subtract πz 1 on both sides of (I.18), and obtainδ 1 (s) = R P * (s, ds ′ )δ 1 (s ′ ) + α R P * (s, ds ′ ) g(θ * , s ′ ) + g ′ (θ * , s ′ )z 1 (s ′ ) + E[R 2 (θ ∞ , s ′ )|x ∞ = s ′ ] = R P * (s, ds ′ ) -π(ds ′ ) δ 1 (s ′ ) + α R P * (s, ds ′ ) -π(ds ′ ) g(θ * , s ′ ) + g ′ (θ * , s ′ )z 1 (s ′ ) + E[R 2 (θ ∞ , s ′ )|x ∞ = s ′ ] .Consolidating the terms, we have(I -P * + Π)δ 1 (s) = α R P * (s, ds ′ ) -π(ds ′ ) g(θ * , s ′ ) + g ′ (θ * , s ′ )z 1 (s ′ ) + E[R 2 (θ ∞ , s ′ )|x ∞ = s ′ ] . (I.20)We next note the following properties,∥E[R 2 (θ ∞ , x ∞ )|x ∞ = s]∥ 2 L 2 (π) ≤ E[∥R 2 (θ ∞ , x ∞ ) 2 ∥]] = O((ατ ) 2 ), ∥z 1 (s)∥ 2 L 2 (π) ≤ E[∥θ ∞ -θ * ∥ 2 ] = O(ατ ) = O(1).Therefore, we can first conclude that ∥δ∥ L 2 (π) = O(α).Subsequently, from (I.19), we can derive that0 = E[g ′ (θ * , x ∞ )z 1 (x ∞ )] + E[R 2 (θ ∞ , x ∞ )] = E[g ′ (θ * , x ∞ )δ 1 (x ∞ )] + ḡ(1) (θ * )E[θ ∞ -θ * ] + E[R 2 (θ ∞ , x ∞ )] E[θ ∞ -θ * ] = -(ḡ (1) (θ * )) -1 E[g ′ (θ * , x ∞ )δ 1 (x ∞ )] + E[R 2 (θ ∞ , x ∞ )] . (I.21)Hence, together with the relationship between δ 1 and z 1 , we havez 1 = δ 1 -(ḡ (1) (θ * )) -1 E[g ′ (θ * , x ∞ )δ 1 (x ∞ )] + E[R 2 (θ ∞ , x ∞ )] . (I.22)Lastly, we substitute (I.22) back into (I.20) and obtain(I -P * + Π)δ 1 (s) = α R P * (s, ds ′ ) -π(ds ′ ) g(θ * , s ′ ) + α R P * (s, ds ′ ) -π(ds ′ ) g ′ (θ * , s ′ ) δ 1 (s ′ ) -(ḡ (1) (θ * )) -1 E[g ′ (θ * , x ∞ )δ 1 (x ∞ )] -α R P * (s, ds ′ ) -π(ds ′ ) g ′ (θ * , s ′ ) (ḡ (1) (θ * )) -1 E[R 2 (θ ∞ , x ∞ )] + α R P * (s, ds ′ ) -π(ds ′ ) E[R 2 (θ ∞ , x ∞ )|x ∞ = s ′ ] = αv(θ * , s) + O(α 2 τ ), where v(θ * , s) = (I -P * + Π) -1 (P * -Π)g θ * (s) = R (I -P * + Π) -1 (P * -Π)(s, ds ′ )g(θ * , s ′ ).To analyze the above system of δ 2 , we make use of the Taylor expansion of g(θ ∞ , s) and ∞+1 (θ ∞ ).Starting with the first term, we haveE[(θ ∞ -θ * ) ⊗ g(θ ∞ , s)|x ∞ = s] = E θ ∞ -θ * ⊗ g(θ * , s) + g ′ (θ * , s)(θ ∞ -θ * ) + R 2 (θ, s) |x ∞ = s = z 1 (s) ⊗ g(θ * , s) + z 2 (s)g ′ (θ * , s) + E[(θ ∞ -θ * ) ⊗ R 2 (θ, s)|x ∞ = s] = δ 1 (s) -(ḡ (1) (θ * )) -1 E[g ′ (θ * , x ∞ )δ 1 (x ∞ )] + E[R 2 (θ ∞ , x ∞ )] ⊗ g(θ * , s) + z 2 (s)g ′ (θ * , s) + E[(θ ∞ -θ * ) ⊗ R 2 (θ, s)|x ∞ = s].Similarly, we have the following relationship for the second term,E[g(θ ∞ , s) ⊗ (θ ∞ -θ * )|x ∞ = s] = g(θ * , s) ⊗ δ 1 (s) -(ḡ (1) (θ * )) -1 E[g ′ (θ * , x ∞ )δ 1 (x ∞ )] + E[R 2 (θ ∞ , x ∞ )] + g ′ (θ * , s)z 2 (s) + E[R 2 (θ, s) ⊗ (θ ∞ -θ * )|x ∞ = s].Next, we proceed to analyze the third term.E[(g(θ ∞ , s)) ⊗2 |x ∞ = s] = E g(θ * , s) + g ′ (θ * , s)(θ ∞ -θ * ) + R 2 (θ ∞ , s) ⊗2 |x ∞ = s = E[(g(θ * , s)) ⊗2 |x ∞ = s] + E[(g ′ (θ * , s)(θ ∞ -θ * )) ⊗2 |x ∞ = s] + E[g(θ * , s) ⊗ (g ′ (θ * , s)(θ ∞ -θ * ))|x ∞ = s] + E[(g ′ (θ * , s)(θ ∞ -θ * )) ⊗ g(θ * , s)|x ∞ = s] + E[g(θ * , s) ⊗ R 2 (θ ∞ , s)|x ∞ = s] + E[R 2 (θ ∞ , s) ⊗ g(θ * , s)|x ∞ = s] + E[(g ′ (θ * , s)(θ ∞ -θ * )) ⊗ R 2 (θ ∞ , s)|x ∞ = s] + E[R 2 (θ ∞ , s) ⊗ (g ′ (θ * , s)(θ ∞ -θ * ))|x ∞ = s] + E[(R 2 (θ ∞ , s)) ⊗2 |x ∞ = s].Lastly, for the noise term, we derive thatE[(ξ ∞+1 (θ ∞ )) ⊗2 |x ∞ = s] = E[(ξ ∞+1 (θ * )) ⊗2 ]+C ′ (θ * )E[θ ∞ -θ * |x ∞ = s]+E[R ′ 2 (θ ∞ )|x ∞ = s].Leveraging on the respective orders, we can conclude that∥δ 2 ∥ L 2 (π) = O(α 2 τ ).Therefore, for second-moment cross-terms, we have the following ordersE[g ′ (θ * , x ∞ )δ 2 (x ∞ )] = O(α 2 τ ) and E[δ 2 (x ∞ )g ′ (θ * , x ∞ )] = O(α 2 τ ), E[g ′′ (θ * , x ∞ )δ 2 (x ∞ )] = O(α 2 τ ). (I.24)Finally, we are ready to consolidate the above analyses and conclude the characterization of the asymptotic bias.We recall that we have already shown the following expansion of the asymptotic biasE[θ ∞ -θ * ] = -(ḡ (1) (θ * )) -1 E[g ′ (θ * , x ∞ )δ 1 (x ∞ )] + 1 2 E[g ′′ (θ * , x ∞ )δ 2 (x ∞ )] + ḡ(2) (θ * )E[(θ ∞ -θ * ) ⊗2 ] + O((ατ ) 3/2 ).By our analyses above, we have shown thatδ 1 = α(I -P * + Π) -1 (P * -Π)g * θ + O(α 2 τ ). Hence, we derive that E[g ′ (θ * , x ∞ )δ 1 (x ∞ )] = αE[g ′ (θ * , x ∞ )(I -P * + Π) -1 (P * -Π)g * θ (x ∞ )] + O(α 2 τ ).From the analyses above, we can also conclude that∥E[θ ∞+1 -θ ∞+1/2 |x ∞ L 2 (π) = O((ατ ) 3/2 ) ∥E[(θ ∞+1 -θ * ) ⊗2 -(θ ∞+1/2 -θ * ) ⊗2 |x ∞ ]∥ L 2 (π) = O((ατ ) 3/2 ).Therefore, combining the analyses above, we see that the projection only introduces error terms of order O(α 2 τ 3 ). Hence, combining the analysis from Section I.2.1, we can conclude the same desired order thatE[θ (α) ∞ -θ * ] = αb + O((ατ ) 3/2 ).In this section, we present more detailed results that characterize the first and second moment of Polayk-Ruppert (PR) tail-averaged iterates and Richardson-Romberg (RR) extrapolated iterates.The following corollary provides non-asymptotic characterization for the first two moments of PR tail-averaged iterates θk0,k .Corollary J.1 (Tail Averaging). Under the setting of Theorem 4.6, the tail-averaged iterates satisfy the following bounds for all k > k 0 + 2τ and k 0 ≥ τ + 1 αµ log 1 ατα :E[ θ(α) k0,k -θ * ] = αb + O (ατ α ) 3/2 + O (1 -αµ) k0/2 α(k -k 0 ) and (J.1) E ( θ(α) k0,k -θ * )( θ(α) k0,k -θ * ) ⊤ = α 2 bb T + O(α • (ατ α ) 3/2 ) + O τ α k -k 0 + (1 -αµ) k0/2 α (k -k 0 ) 2 . (J.2)With this result, taking the trace on both sides of (J.2) recovers Corollary 4.7.Proof. First, we haveE[ θ(α) k0,k -θ * ] = (E [θ ∞ ] -θ * ) + 1 k -k 0 k-1 t=k0 E [θ t -θ ∞ ] .By Corollary 4.4, we obtain∥E[θ t ] -E[θ ∞ ]∥ ≤ (1 -αµ) t 2 • s ′ (θ 0 , L, µ). Hence, it follows that k-1 t=k0 E [θ t -θ ∞ ] ≤ k-1 t=k0 ∥E [θ t ] -E[θ ∞ ]∥ ≤ s ′ (θ 0 , L, µ) • (1 -αµ) k 0 2 1 1 -(1 -αµ) ≤ s ′ (θ 0 , L, µ) • (1 -αµ) k 0 2 2 αµ .Together with Theorem 4.6, we haveE θ(α) k0,k -θ * = αb + O (ατ α ) 3/2 + O (1 -αµ) k0/2 α(k -k 0 ) ,thereby finishing the proof of the first moment.To bound the second moment of the tail-averaged iterate, we follow the proof technique in [33, Section A.6.2] . We notice that θ(α) k0,k -θ * θ(α) k0,k -θ * ⊤ =E θ(α) k0,k -E [θ ∞ ] + E [θ ∞ ] -θ * θ(α) k0,k -E [θ ∞ ] + E [θ ∞ ] -θ * ⊤ = E θ(α) k0,k -E [θ ∞ ] θ(α) k0,k -E [θ ∞ ] ⊤ T1 + E θ(α) k0,k -E [θ ∞ ] (E [θ ∞ ] -θ * ) ⊤ T2 + E (E [θ ∞ ] -θ * ) θ(α) k0,k -E [θ ∞ ] ⊤ T3 + E (E [θ ∞ ] -θ * ) (E [θ ∞ ] -θ * ) ⊤ T4 .For T 2 , we haveT 2 = 1 k -k 0 k-1 t=k0 E [θ t -θ ∞ ] (E[θ ∞ ] -θ * ) ⊤ = O (1 -αµ) k0/2 α(k -k 0 ) • (αb + O (ατ α ) 3/2 ) = O (1 -αµ) k0/2 (k -k 0 ) .The term T 3 is similar to T 2 and obeys the same bound.For T 4 , we haveT 4 = (αb + O (ατ α ) 3/2 )(αb + O (ατ α ) 3/2 ) T = α 2 bb T + O(α • (ατ α ) 3/2 ).For T 1 , we haveT 1 = 1 (k -k 0 ) 2 E k-1 t=k0 θ t -E[θ ∞ ] k-1 t=k0 θ t -E[θ ∞ ] ⊤ = 1 (k -k 0 ) 2 k-1 t=k0 E θ t -E[θ ∞ ] θ t -E[θ ∞ ] ⊤ (J.3) + 1 (k -k 0 ) 2 k-1 t=k0 k-1 l=t+1 E θ t -E[θ ∞ ] θ l -E[θ ∞ ] ⊤ (J.4) + 1 (k -k 0 ) 2 k-1 t=k0 k-1 l=t+1 E θ l -E[θ ∞ ] θ t -E[θ ∞ ] ⊤ . (J.5)By Corollary 4.4 and Proposition 4.2 we haveE θ t -E[θ ∞ ] θ t -E[θ ∞ ] ⊤ = E θ t θ ⊤ t -E θ ∞ θ ∞ ⊤ + E θ ∞ θ ∞ ⊤ -E[θ ∞ ]E θ ∞ ⊤ -E [θ t ] E θ ∞ ⊤ + E[θ ∞ ]E θ ⊤ t -2E[θ ∞ ]E θ ∞ ⊤ = E θ t θ ⊤ t -E θ ∞ θ ∞ ⊤ + Var θ ∞ -E θ t -θ ∞ E θ ∞ ⊤ -E[θ ∞ ]E (θ t -θ ∞ ) ⊤ =O (1 -αµ) t 2 + ατ α ,where we bound Var θ ∞ with O(ατ α ) by Proposition 4.2 and Fatou's lemma.Then, for (J.3), we have(J.3) = 1 (k -k 0 ) 2 k-1 t=k0 (1 -αµ) t 2 + ατ α = O 1 (k -k 0 ) 2 ∞ t=k0 (1 -αµ) t 2 + O ατ α k -k 0 = O (1 -αµ) k0/2 α (k -k 0 ) 2 + ατ α k -k 0 .We restate the following claim, whose proof closely resembles Claim 4 in [33] .Claim 2. For t ≥ τ + 1 αµ log 1 ατα and l ≥ t + 2τ α , we haveE θ t -E[θ (α) ] θ l -E[θ (α) ] ⊤ = O (ατ α ) • (1 -αµ) (l-t) 2.Then, by [33, Claim 4] , we have term (J.4) = O τα k-k0 . Similarly, we have term (J.5)= O τα k-k0 . Therefore, we haveT 1 = O (1 -αµ) k0/2 α (k -k 0 ) 2 + τ α k -k 0 . (J.6)By adding T 1 -T 4 together, we obtainE θ(α) k0,k -θ * θ(α) k0,k -θ * ⊤ =α 2 bb T + O(α • (ατ α ) 3/2 ) + O (1 -αµ) k0/2 (k -k 0 ) + O (1 -αµ) k0/2 α (k -k 0 ) 2 + τ α k -k 0 =α 2 bb T + O(α • (ατ α ) 3/2 ) + O τ α k -k 0 + (1 -αµ) k0/2 α (k -k 0 ) 2 .Next, we present the following corollary formalizes the non-asymptotic characterization for the first two moments of the RR-extrapolated iterate θ (α) k0,k . Corollary J.2 (Richardson-Romberg Extrapolation). Under the setting of Theorem 4.6, the RR extrapolated iterates with stepsizes α and 2α satisfy the following bounds for all k > k 0 + 2τ α and k0 ≥ τ α + 1 αµ log 1 ατα : E[ θ(α) k0,k -θ * ] = O (ατ α ) 3/2 + O (1 -αµ) k0/2 α(k -k 0 ), and (J.7) With these notations, θk0,k -θ * = 2u 1 -u 2 + v. We then have the following boundE ( θ(α) k0,k -θ * )( θk0,k -θ * ) ⊤ = O (ατ α ) 3 + O τ α k -k 0 + (1 -αµ) k0/2 α (k -k 0 )E θ(α) k0,k -θ * θ(α) k0,k -θ * ⊤ ≤E ∥2u 1 -u 2 + v∥ 2 ≤3E ∥2u 1 ∥ 2 + 3E ∥u 2 ∥ 2 + 3∥v∥ 2 .By equation (J.6), we haveE ∥u 1 ∥ 2 = Tr E u 1 u ⊤ 1 = O (1 -αµ) k0/2 α (k -k 0 ) 2 + τ α k -k 0 .Similarly, we haveE ∥u 2 ∥ 2 2 = O (1 -2αµ) k0/2 α (k -k 0 ) 2 + τ 2α k -k 0 .By Theorem 4.6, we have ∥v∥ 2 2 = O (ατ α ) 3 . Combining these bounds, we haveE θk0,k -θ * θk0,k -θ * ⊤ = O (ατ α ) 3 + O τ α k -k 0 + (1 -αµ) k0/2 α (k -k 0 ) 2 .• 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 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.Justification: The research conducted in the paper conforms, in every respect, with the NeurIPS Code of Ethics.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).The experiments were conducted on two sockets of Intel(R) Xeon(R) Gold 6154 CPU @ 3.00GHz with 566Gb of RAM. Implementation details and code are available at https://github.com/lucyhuodongyan/ nonlinear-sa-bias. 


Decision: The reviewers were in unanimous agreement that the paper solves a challenging problem caused by Markovian data and nonlinear updates. Also the interplay between memory and nonlinearity in stochastic approximations is particularly striking and good for the community to know. A clear accept.

Reviews: review: paper_summary: paper_summary: This paper investigates stochastic approximation (SA) with Markovian data and nonlinear updates under constant stepsize, and establishes the weak convergence of the joint process $(x_t, \theta_t)$. It also presents a precise characterization of the asymptotic bias of the SA iterates.
strength_weakness: strength_weakness: I find the presentation of this paper a bit technical for people that are not very familiar with this area. In addition, perhaps some numerical experiments should be performed to better illustrate the theory.
questions: questions: 1. About Assump. 3: for many GLMs, we actually do not have strong convexity. I feel this assumption is a bit strong. 
2. Page 4, line 147: is the notation superscript "\cross 2" defined anywhere? I assume this denotes the outer product of a vector. 
3. Similar to Assump. 3, Assump. 4 should be justified further as well.
limitations: limitations: NA
paper_summary: paper_summary: This paper considers a nonlinear stochastic approximation (SA) problem with Markov noise (MC). It is assumed that the MC is uniformly geometrically ergodic. Instead of the standard iterative procedure, a projection onto a bounded set is additionally introduced (the latter can be relaxed under the additional assumption of the existence of a positive density). Under these assumptions the authors manage to write a decomposition that characterizes the bias. At the same time, they manage to identify three factors influencing the bias: the factor of MC, the factor of nonlinearity of the procedure, and the factor of interaction between MC and nonlinearity. In addition, bounds for the Polyak-Ruppert averaging and the Richardson Romberg procedure are given.
strength_weakness: strength_weakness: It would be good to obtain:
- high probability bounds instead of the MSE
- remove additional projection step and assumption on the density (it could be useful to consider convergence in the weighted W distance instead of V norm)
-  explicit dependence on the asymptotic variance of MC in the $O(\tau/(k - k_0))$.
questions: questions: Could you please comment on the fact that there is no dependence between stepsize and number of iterates in the Polyak Ruppert avaraging?
limitations: limitations: -
paper_summary: paper_summary: The present paper obtains a representation for the asymptotic bias of constant step-size nonlinear stochastic approximation with Markovian noise. In particular, this characterization makes the hindering effect of memory and nonlinearity explicit in terms of the algorithm's performance. Moreover, the authors establish ergodicity of the parameter-noise joint process (with and without projection of estimates), obtain finite-time bounds on L_p moments of the estimation error and establish a central limit theorem for the constant gain algorithm. Finally, a bias attenuation technique based upon the Richardson-Romberg extrapolation is proposed for the nonlinear algorithm.
strength_weakness: strength_weakness: One of the weaknesses of this paper is the fact that no numerical experiments are provided to illustrate any of the main results. Although a discussion on how the theory fits within Generalized Linear Models is given before the conclusions, I encourage the authors to include a simple toy example to illustrate some of their main results such as the CLT or the bias attenuation technique.
questions: questions: -	Could the authors clarify if assuming strong monotonicity and uniform boundedness of g and its Lipschitz constant over x are indeed needed for the result in Thm. 4.6? I can see their importance for finite-time bounds, but if there are needed for the asymptotic bias bound, I believe that the authors should mention the employment of stronger assumptions when comparing their work with previous research on asymptotic results.

-	I am not sure I understand what the authors mean by ``fine-grained’’  when talking about the result in Thm. 4.6. Upon further inspection, it seems that equation (4.2) is an extension of the result in [40] for nonlinear recursions: it consists a representation for the dominant bias term plus an upper bound as in this previous work. Is the fine-grained part related to the upper bound for \alpha?

-	I might have missed this in the text, but could the authors clarify if the results in Sections 4.3 and 4.5 require projection?
limitations: limitations: The authors addressed the limitations in their work through a clear list of assumptions and discussions after presenting the main results.
paper_summary: paper_summary: This paper studies constant step-size stochastic approximation algorithms with Markovian noise. It is known that in this case, the error of a stochastic does not vanish asymptotically, even when consider averaging. This is \theta_n has a bias.  Previous work have show how to study the bias of such algorithm for linear stochastic approximation. In this paper, the authors tackle the more challenging case of non-linear stochastic approximation. Most of the paper is devoted to proving that the bias is of order O(\alpha). Some discussion on the applicability of results is presented.
strength_weakness: strength_weakness: Most of the assumptions needed to obtain the results are relatively mild but two conditions are quite strong:
1. strong monotonicity (A3)
2. Smoothness (A2)
I think that the paper should discuss in more details these assumption and highlight that they are really limiting the applicability of the result.

The assumption developed in part 4.2 to avoid assuming that the iterates are bounded seems quite strong, not verified in many practical cases (for instance any instance of stochastic gradient descent where the noise has a finite support would not satisfy this assumption), and very artificial (is the only purpose of this assumption to

The paper is extremely long (54 pages including proofs and references). To me, this says that there is either too much content or that the results are too diluted. As a result, the paper is very technical and hard to read. There should be more effort to make the paper more readable. Some suggestions:
- consider a slightly less general setting
- avoid considering sub-cases (like section 4.2) that are a bit orthogoal to the paper

The practical applications of the results are unclear. Some potential applications are presented in Section 4.5 / 4.6 but there are no experiments or simulations to confirm that this actually work (these sections are interesting, though).

There is a lot of papers on the subject. This can be seen as a good sign (this is an active area of research), but at the same time, it is hard for me to really assess the novelty of the results. Note that there are some (recent) papers that seem related to these work and that might be cited. This last point is not a weakness since some of them are extremely recent (available online after the submission deadline):
- Computing the Bias of Constant-step Stochastic Approximation with Markovian Noise. S Allmeier, N Gast (2024)
- Bias in Stochastic Approximation Cannot Be Eliminated With Averaging, Caio Kalil Lauand, Sean P. Meyn (2022)
- Revisiting Step-Size Assumptions in Stochastic Approximation. Caio Kalil Lauand, Sean Meyn (2024)
A comparison with these papers might be useful (even of this is not mandatory for the two 2024 papers).
questions: questions: There are some questions / comments in the limit that would deserve some comments from the authors.

The Markovian noise (x_k) is assumed exogeneous (it does not depend on \theta). For some application (like Q-learning with a navigating policy derived from \theta), this would not be satisfied. Could this assumption be lifted?

The next order term is a O(\alpha^{3/2}): is it the sharpest bound or could O(\alpha^2) be obtained?
limitations: limitations: The limitation of the assumptions (see "weaknesses") should be better discussed.

Comments: comment: I thank the authors for their response. I retain my current score.
Thank you for the clarification and response. We will include those discussions in our revised paper.
I thank the authors for their clarifying and detailed responses. I also thank them for including an experiment that validates their theorems. Also, it is interesting to see the experiment with i.i.d. data performing worst than the experiment with Markovian data.

Here a few more observations:

**On a bias equation** I apologize If I was not clear in my first response. When I mentioned [40] I was referring to their equation (10) where there are no limsups and not (8).

However, it is very clear to me that this does not affect the value of our work since they do not incorporate the effects of nonlinearity. Still, I think I would be beneficial to mention that a similar fine-grained expression was obtained for the more restrictive case of linear $g$ in the final version.

**Strong Monotonicity** Maybe it would be beneficial to provide a bit more discussion regarding Assumption 3 in the final version like in the responses provided (e.g. when it is satisfied/ or if it could be lifted).
Thank you for this detailed answer. I have no further questions.
The authors have done a good job in responding to my previous comments. I find their responses clear and detailed. I have no further questions.