1 Introduction↩︎

Root-finding problems, in particular the task of solving fixed-point equations, provide a unifying framework for formulating and analyzing a broad spectrum of computational and analytical tasks that arise across the mathematical sciences. For example, minimizing a convex objective function is equivalent to finding the zeros of its gradient operator [1]. In reinforcement learning (RL) [2], [3], which has seen a surge of interest in recent years [4], [5], computing an optimal policy reduces to solving a fixed-point equation known as the Bellman equation [6]. Additional applications of fixed-point equations arise in game theory and differential equations [7], [8].

Formally, let \(\bar{F}:\mathbb{R}^d \to \mathbb{R}^d\) be an operator. The fixed-point equation of \(\bar{F}(\cdot)\) takes the form \(\bar{F}(x) = x\), and any solution is called a fixed point of \(\bar{F}(\cdot)\). A natural question is: when does \(\bar{F}(x) = x\) admit a (unique) solution? Classical results such as the Banach fixed-point theorem [9] provide sufficient conditions. In particular, it states that if \(\bar{F}(\cdot)\) is a contraction mapping with respect to some norm, denoted by \(\|\cdot\|_c\), then \(\bar{F}(x) = x\) has a unique solution \(x^*\). Moreover, the fixed-point iteration \(x_{k+1}=\bar{F}(x_k)\) converges to \(x^*\) at a geometric rate in \(\|\cdot\|_c\).

While such fixed-point iterations are appealing in theory, in practice one often lacks the information or computational resources required to evaluate \(\bar{F}(\cdot)\) exactly. For instance, in large-scale optimization, computing the full gradient may be prohibitively expensive; in RL, the underlying stochastic model is typically unknown, preventing direct evaluation of the Bellman operator. In these situations, one must design algorithms that rely on noisy estimates of \(\bar{F}(\cdot)\), which leads naturally to the stochastic approximation (SA) method [10]. SA forms a foundational framework underlying many modern large-scale methods, including stochastic gradient descent [11], temporal-difference (TD) learning [12], Q-learning [13], and fictitious play [14].

The generic SA scheme for solving the fixed-point equation \(\bar{F}(x) = x\) takes the form \[\begin{align} \label{algo:sa1} x_{k+1} = x_k + \alpha_k \big( F(x_k, Y_k) - x_k \big),\;\forall\, k \ge 0, \end{align}\tag{1}\] where \(\{Y_k\}\) is a stochastic process taking values in a set \(\mathcal{Y}\), \(F : \mathbb{R}^d \times \mathcal{Y} \to \mathbb{R}^d\) is an operator, and \(\alpha_k > 0\) is the stepsize. For simplicity, we assume that \(\mathcal{Y}\) is finite, that \(\{Y_k\}\) is an i.i.d.sequence with distribution \(\mu \in \Delta(\mathcal{Y})\),¹ and that the operator \(F(\cdot,\cdot)\) satisfies \(\mathbb{E}_{Y \sim \mu}[F(x,Y)] = \bar{F}(x)\), so that \(F(x_k, Y_k)\) provides a conditionally unbiased estimate of \(\bar{F}(x_k)\). In later sections, we will discuss several extensions, including settings where the noise process \(\{Y_k\}\) is not i.i.d.but instead forms a Markov chain, and settings where \(\bar{F}(\cdot)\) is not a contraction but admits other useful structural properties.

Early work on SA focused on asymptotic convergence. A powerful analytical tool in this setting is the ordinary differential equation (ODE) method [15]–[18], which relates the discrete-time stochastic process defined by the SA iteration 1 to the continuous-time ODE \[\begin{align} \label{eq:ODE} \dot{x}(t) = \bar{F}(x(t)) - x(t). \end{align}\tag{2}\] Under mild conditions on the operators \(F(\cdot,\cdot)\) and \(\bar{F}(\cdot)\), as well as on the decay rate of the stepsize sequence \(\{\alpha_k\}\), if the equilibrium point \(x^*\) of ODE 2 is globally asymptotically stable [19], [20], then the SA iterates satisfy \(\lim_{k\rightarrow\infty}x_k=x^*\) almost surely [16].

Although such asymptotic theory is powerful, it does not quantify the rate of convergence. For many practical applications, particularly in large-scale learning systems, finite-time or non-asymptotic guarantees are more relevant. These guarantees characterize how quickly an algorithm approaches the desired limit and, in data-driven settings, provide concrete guidance on the number of samples or iterations required to achieve a prescribed accuracy [21]. In finite-time analysis, the objective is to understand how the error \(\|x_k - x^*\|_c\) decays with \(k\). Two common performance metrics are the mean squared error and high-probability bounds. This survey emphasizes mean squared error while also reviewing results that provide high-probability error guarantees.

A central tool for the finite-time analysis we review is the Lyapunov-drift approach. Originally developed to study stability in ODEs [19], [20], [22], this method constructs a scalar function \(M:\mathbb{R}^d\to\mathbb{R}\) (a Lyapunov or potential function) whose evolution tracks the decay of the ODE towards equilibrium. For example, consider the ODE \(\dot{x}(t) = -x(t)\), which is a special case of 2 with \(\bar{F}(x) \equiv 0\). For this ODE, the choice \(M(x) = \|x\|_2^2\) yields \(\dot{M}(x(t)) = -2 M(x(t))\), implying the geometric decay \(M(x(t)) = M(x(0)) e^{-2t}\) for any \(t \ge 0\). By adapting such Lyapunov arguments to the discrete and stochastic setting, one expects to obtain convergence rates for the SA iteration 1 . However, constructing a valid Lyapunov function is well known to be challenging, and this difficulty is amplified in the context of SA algorithms. The Lyapunov function must, on the one hand, induce a negative drift, and on the other hand, enable control of errors arising from stochasticity and discretization. In this survey, we present a detailed roadmap for the systematic construction of valid Lyapunov functions tailored to SA.

The remainder of this survey is organized as follows. Section 2 presents the full mean-square analysis of the SA recursion 1 , highlighting the construction of a valid Lyapunov function. Section 3 discusses applications in RL. Section 4 reviews several generalizations, such as SA with Markovian noise, SA with seminorm-contractive operators, and SA with dissipative operators, as well as their applications in RL. Section 5 surveys high-probability bounds. We conclude in Section 6 with open research directions.

2 Stochastic Approximation: Mean-Square Bounds↩︎

Consider the SA recursion (1 ). We begin by stating the assumptions regarding the operators \(\bar{F}(\cdot)\), \(F(\cdot,\cdot)\), and the noise process \(\{Y_k\}\).

Under Assumption 1, the fixed-point equation \(\bar{F}(x)=x\) admits a unique solution, denoted by \(x^*\) [9]. As a remark, Assumption 1 can be relaxed to the setting where \(\bar{F}(\cdot)\) is a pseudo-contraction mapping, meaning that \(\|\bar{F}(x)-x^*\|_c \leq \gamma_c \|x-x^*\|_c\) for all \(x \in \mathbb{R}^d\) [3]. In that case, however, the existence of \(x^*\) does not follow automatically and must be imposed separately. Since our goal is to review the Lyapunov approach for analyzing SA algorithms, we adopt Assumption 2 for ease of presentation. Extensions to more general noise models will be discussed in Section 4.1. Assumption 3 is standard even in the classical literature on asymptotic convergence of SA [16], and it is satisfied in a wide range of applications, such as TD-learning and Q-learning in RL.

In the remainder of this section, we present our detailed approach for establishing the mean-square bounds of \(\{x_k\}\) generated by the SA recursion 1 . Adopting the Lyapunov-drift approach, we begin by illustrating the main challenges in constructing a valid Lyapunov function in Section 2.1. Building on these insights, we then describe a systematic method for constructing a smooth Lyapunov function using the generalized Moreau envelope in Section 2.2. With a valid Lyapunov function in hand, we proceed to derive a one-step Lyapunov-drift inequality for the SA recursion in Section 2.3, and finally obtain the finite-time bound by solving the resulting recursion in Section 2.4.

2.1 Fundamental Challenges↩︎

To provide intuition, assume for now that the contraction norm \(\|\cdot\|_c\) is the \(\ell_p\)-norm for some \(p \in [2,\infty)\). Consider the ODE 2 associated with the SA algorithm. It was shown in [16] that \(W(x) = \|x - x^*\|_p\) satisfies \(\dot{W}(x(t)) \le -\kappa W(x(t))\) for some \(\kappa > 0\). This inequality implies that the solution \(x(t)\) of ODE 2 converges to its equilibrium point \(x^*\) at a geometric rate, which in turn ensures the asymptotic convergence of the SA algorithm via the ODE approach [16], [18].

The ODE approach yields asymptotic convergence but does not provide finite-time guarantees. To obtain finite-time bounds, we analyze the SA algorithm directly rather than its ODE counterpart. Observe that the SA algorithm in recursion 1 can be equivalently written as \[\begin{align} \label{eq:sa3} x_{k+1} - x_k &= \underbrace{\alpha_k \left( \bar{F}(x_k) - x_k \right)}_{\text{Noiseless Update}} + \underbrace{\alpha_k \left( F(x_k, Y_k) - \bar{F}(x_k) \right)}_{\text{Effective Noise}}. \end{align}\tag{3}\] Due to discretization and stochastic errors, the ODE Lyapunov function \(W(x)\) generally cannot be directly used to analyze the SA algorithm. Suppose there exists a function \(M(x)\) that induces a negative drift for the ODE in the sense that \(\dot{M}(x(t)) \leq -\kappa' M(x(t))\) and, in addition, is \(L\)-smooth with respect to some norm, denoted by \(\|\cdot\|_s\). Then we can handle both discretization and stochastic errors to obtain \[\begin{align} \label{eq:contraction} \mathbb{E}[M(x_{k+1} - x^*)] \le \left( 1 - \mathcal{O}(\alpha_k) + o(\alpha_k) \right) \mathbb{E}[M(x_k - x^*)] + o(\alpha_k), \end{align}\tag{4}\] which implies an almost contraction in \(\mathbb{E}[M(x_{k+1} - x^*)]\). A finite-time bound then follows by recursively applying 4 . The key point is that the smoothness of \(M(\cdot)\), together with its negative drift in the ODE dynamics, yields a contraction factor of the form \(\left( 1 - \mathcal{O}(\alpha_k) + o(\alpha_k) \right)\) for \(\{x_k\}\). In the case of \(\ell_p\)-norm contraction, a natural choice is the norm-square function \(M(x) = \|x - x^*\|_p^2\), which is known to be smooth [1].

In the case where the contraction norm \(\|\cdot\|_c\) is arbitrary, the norm-square function \(f(x) = \|x - x^*\|_c^2\) is not necessarily smooth, and the key difficulty lies in constructing a smooth Lyapunov function that also induces a negative drift. An important special case is when \(\|\cdot\|_c = \|\cdot\|_\infty\), which is particularly relevant in RL, since the Bellman operator is a contraction mapping with respect to \(\|\cdot\|_\infty\).

We address this challenge by constructing a smoothed convex envelope \(M(x)\), referred to as the generalized Moreau envelope, that is smooth and closely approximates the original norm-square function \(f(x)\) in the sense that \(a M(x) \le f(x) \le b M(x)\) for some constants \(a, b\) close to one. The approximation property ensures that \(M(\cdot)\) serves as a valid Lyapunov function with a negative drift, while the smoothness property enables us to control both the discretization error and the stochastic error in the SA algorithm.

2.2 Construction of the Lyapunov Function↩︎

To present the construction of a valid Lyapunov function, we need the following definitions.

Intuitively, the \(L\)-smoothness ensures that the function does not curve too sharply: all higher-order terms in the Taylor expansion of \(g(y)\) around \(x\) can be uniformly bounded by a constant multiple of the quadratic term \(\|x-y\|^2\).

Next, we introduce the generalized Moreau envelope.

In the existing literature, \(M_{h_1}^{\theta,h_2}(\cdot)\) in Definition 2 is also referred to as the infimal convolution [1] between \(h_1(\cdot)\) and \(h_2(\cdot)/\theta\), and is sometimes denoted by \(h_1 \square (h_2/\theta)\) [1], [23]. In the special case where \(h_2(x) = \|x\|_2^2\) and \(\theta = 1\), \(M_{h_1}^{\theta,h_2}(\cdot)\) reduces to the classical Moreau envelope [1].

Let \(f(x) = \frac{1}{2}\|x\|_c^2\), where \(\|\cdot\|_c\) is the contraction norm of \(\bar{F}(\cdot)\). Let \(\|\cdot\|_s\) be an arbitrary norm on \(\mathbb{R}^d\) such that \(g(x) := \frac{1}{2}\|x\|_s^2\) is \(L\)-smooth with respect to \(\|\cdot\|_s\) for some \(L > 0\). For example, \(\|\cdot\|_s\) may be chosen as the \(\ell_p\)-norm for any \(p \in [2,\infty)\) [1]. We use the generalized Moreau envelope of \(f(\cdot)\) with respect to \(g(\cdot)\), i.e., \[\begin{align} M_f^{\theta,g}(x) = \min_{u \in \mathbb{R}^d} \left\{ f(u) + \frac{1}{\theta} g(x-u) \right\}, \end{align}\] as our Lyapunov function for analyzing the SA recursion 1 , where \(\theta>0\) is a tunable parameter to be specified. The following proposition shows that \(M_f^{\theta,g}(\cdot)\) serves as a smooth approximation of the norm-square function \(f(\cdot)\). The proof of Proposition [prop:Moreau] is provided in Appendix 8.1.

Proposition [prop:Moreau] (1) states that this generalized Moreau envelope is itself a norm-square function. Propositions [prop:Moreau] (2) and (3) together show that \(M_f^{\theta,g}(\cdot)\) provides a smooth approximation of \(f(\cdot)\).

2.3 Establishing a One-Step Contractive Inequality↩︎

For simplicity of notation, we henceforth write \(M(x)\) for \(M_f^{\theta,g}(x)\). For any \(k \ge 0\), by the smoothness of \(M(\cdot)\) and the SA recursion 1 , we have \[\begin{align} \label{eq:composition195contractive} M(x_{k+1} - x^*) \leq \;& M(x_k - x^*) + \langle \nabla M(x_k - x^*), x_{k+1} - x_k \rangle+ \frac{L}{2\theta} \|x_{k+1} - x_k\|_s^2 \nonumber\\ =\; & M(x_k - x^*) + \alpha_k \langle \nabla M(x_k - x^*), F(x_k, Y_k) - x_k \rangle \nonumber\\ &+ \frac{L \alpha_k^2}{2\theta} \|F(x_k, Y_k) - x_k\|_s^2 \nonumber\\ =\; & M(x_k - x^*) + \underbrace{\alpha_k \langle \nabla M(x_k - x^*), \bar{F}(x_k) - x_k \rangle}_{E_1:\;\text{The expected update}} \nonumber\\ &+ \underbrace{\alpha_k \langle \nabla M(x_k - x^*), F(x_k, Y_k) - \bar{F}(x_k) \rangle}_{E_2:\;\text{The error due to noise } Y_k} \nonumber\\ &+ \underbrace{\frac{L \alpha_k^2}{2\theta} \|F(x_k, Y_k) - x_k\|_s^2}_{E_3:\;\text{The error due to discretization and noise}}. \end{align}\tag{5}\]

The term \(E_1\) on the right-hand side of Inequality 5 captures the deterministic part of the SA iteration 1 , and is bounded in the following lemma. The proof of Lemma 1 is provided in Appendix 8.2.

By choosing the tunable parameter \(\theta\) such that \(1 - \gamma_c\, u_{cm} / \ell_{cm} \in (0,1)\), which is always feasible since \(\lim_{\theta \to 0} u_{cm} / \ell_{cm} = 1\) (cf.Proposition [prop:Moreau]) and \(\gamma_c \in (0,1)\), Lemma 1 ensures a negative drift.

We next consider the term \(E_2\) in Inequality 5 . Since \(\{Y_k\}\) is an i.i.d.sequence, we have by the tower property of conditional expectations that \[\begin{align} \mathbb{E}[E_2] = \alpha_k \mathbb{E} \left[ \left\langle \nabla M(x_k - x^*),\, \mathbb{E}[F(x_k, Y_k) \mid \mathcal{F}_k] - \bar{F}(x_k) \right\rangle \right] = 0, \label{eq:T952} \end{align}\tag{6}\] where \(\mathcal{F}_k\) denotes the \(\sigma\)-algebra generated by \(\{Y_n\}_{0 \le n \le k-1}\). Equation 6 follows from the fact that \(x_k\) is measurable with respect to \(\mathcal{F}_k\) and \(Y_k\) is independent of \(\mathcal{F}_k\). In the more general setting where \(\{Y_k\}\) forms a Markov chain, more sophisticated arguments are required to control \(\mathbb{E}[E_2]\), which will be discussed in Section 4.1.

Lastly, we bound the term \(E_3\) from Inequality 5 in the following lemma. The proof of Lemma 2 relies primarily on the Lipschitz continuity of \(F(\cdot,\cdot)\) (cf.Assumption 3) and is deferred to Appendix 8.3.

Using Lemma 1, Equation 6 , and Lemma 2 together in Inequality 5 , we obtain the following one-step contractive inequality, whose proof is presented in Appendix 8.4. For simplicity of notation, we denote \(M_k = \mathbb{E}[M(x_k - x^*)]\).

Lemma 3 has the desired form 4 , showing that the SA recursion 1 is overall contractive with respect to the Lyapunov function \(M(\cdot)\) up to higher-order error terms.

2.4 Solving the Recursion↩︎

The last step is to solve the recursion ?? and evaluate the resulting bounds under stepsizes with different decay rates. Let \(\varphi_2 = 1 - \gamma_c u_{cm} / \ell_{cm}\). By repeatedly applying Inequality ?? , we have for all \(k \ge 0\) that \[\begin{align} M_k \leq \prod_{j=0}^{k-1}\left(1-\varphi_2\alpha_j\right)M_0 + \frac{L}{\theta \ell_{cs}^2} \big( (L_1+1)\|x^*\|_c + L_2 \big)^2 \sum_{i=0}^{k-1}\alpha_i^2 \prod_{j=i+1}^{k-1}\left(1-\varphi_2\alpha_j\right). \end{align}\] After using Proposition [prop:Moreau] (2) to translate \(M_k\) back to the mean squared error \(\mathbb{E}[\|x_k-x^*\|_c^2]\), we obtain for all \(k \ge 0\) that \[\begin{align} \label{eq:before95solving} \mathbb{E}[\|x_k - x^*\|_c^2] \leq\;& \varphi_1 \prod_{j=0}^{k-1}\left(1-\varphi_2\alpha_j\right) + \varphi_3 \sum_{i=0}^{k-1}\alpha_i^2 \prod_{j=i+1}^{k-1}\left(1-\varphi_2\alpha_j\right), \end{align}\tag{7}\] where \[\begin{align} \varphi_1 = \frac{\|x_0 - x^*\|_c^2 u_{cm}^2}{\ell_{cm}^2}, \qquad \varphi_3 = \frac{2 u_{cm}^2 L}{\theta \ell_{cs}^2} \big( (L_1+1)\|x^*\|_c + L_2 \big)^2. \end{align}\] The final step is to evaluate the right-hand side of Inequality 7 for different choices of stepsizes. Specifically, we consider stepsizes of the form \(\alpha_k = \alpha/(k+h)^z\), where \(\alpha, h > 0\) and \(z \in [0,1]\). Note that \(z = 0\) corresponds to a constant stepsize, while \(z \in (0,1]\) corresponds to diminishing stepsizes. Evaluating such terms under this type of stepsize schedule is well studied in the convergence-rate analysis of iterative algorithms [1], [24] and is omitted here; see [25] for additional details.

2.5 Finite-Time Mean-Square Bounds↩︎

After evaluating Inequality 7 under the specified stepsize rule, we obtain the finite-time bounds for the SA recursion 1 stated in the following theorem.

In view of Theorem 1, it is clear that the choice of stepsizes plays an important role in determining the convergence rate of the SA algorithm. Before illustrating each case of Theorem 1 in detail, we provide a high-level intuition. For ease of illustration, consider the special case where the contraction norm \(\|\cdot\|_c\) is the Euclidean norm \(\|\cdot\|_2\). In this case, we have the following standard decomposition: \[\begin{align} \underbrace{\mathbb{E}[\|x_k - x^*\|_2^2]}_{\text{Mean Squared Error}} = \underbrace{\|\mathbb{E}[x_k] - x^*\|_2^2}_{\text{Squared bias}} + \underbrace{\mathbb{E}[\|x_k - \mathbb{E}[x_k]\|_2^2]}_{\text{Variance}}. \end{align}\]

Suppose there is no noise, i.e., \(F(\cdot, y) \equiv \bar{F}(\cdot)\) for any \(y \in \mathcal{Y}\), so the variance term vanishes and the mean-square error reduces to the squared bias. In this case, we should choose a large stepsize to maximize the per-step improvement. This is indeed the case: choosing \(\alpha_k \equiv 1\) recovers the fixed-point iteration and yields geometric convergence. When there is noise, however, a trade-off arises in choosing the stepsizes. The bias term favors large stepsizes, as discussed above, while the variance term favors small stepsizes. This is because at each iteration the effective noise injected into \(x_{k+1}\) is \(\alpha_k (F(x_k, Y_k) - \bar{F}(x_k))\), whose variance is proportional to \(\alpha_k^2\). With this high-level intuition regarding the bias–variance trade-off in choosing the stepsizes in mind, we next discuss the three cases presented in Theorem 1.

In all cases of Theorem 1, the convergence bound consists of a bias term (the first term on the right-hand side of the inequality) and a variance term (the second term on the right-hand side). A constant stepsize yields geometric decay of the bias with a constant variance that is proportional to the stepsize. When using diminishing stepsizes of the form \(\alpha_k = \alpha/(k+h)\), and when the parameter \(\alpha\) is above a certain threshold, namely \(1/\varphi_2\), both terms decay to zero at a rate of \(\mathcal{O}(1/k)\). However, when \(\alpha\) is not chosen properly, the convergence rate can be arbitrarily bad. We will return to this issue in Section 6 when discussing open problems in SA. Finally, when using diminishing stepsizes of the form \(\alpha_k = \alpha/(k+h)^z\) with \(z \in (0,1)\), the bias term decays to zero at an almost geometric rate, while the variance term decays at a suboptimal rate of \(\mathcal{O}(1/k^z)\), so the overall convergence rate is \(\mathcal{O}(1/k^z)\). Importantly, the convergence rate in this case, while being sub-optimal, is independent of \(\alpha\) (which only appears as a multiplicative constant), and is therefore more robust.

2.6 Connection to SGD and linear SA↩︎

In this section, we illustrate that while we present the finite-time mean squared bound for a generic contractive SA framework, the result also captures several other popular SA settings, including SGD and linear SA.

2.6.0.1 Stochastic Gradient Descent.

Consider the problem \(\min_{x\in\mathbb{R}^d} J(x)\). With access to a noisy oracle that returns a stochastic estimate of the gradient \(\nabla J(x)\), the SGD updates \(x_k\) according to \[\begin{align} x_{k+1} = x_k + \beta_k \big( -\nabla J(x_k) + w_k \big), \end{align}\] where \(w_k\) denotes the mean-zero noise and \(\beta_k\) is the stepsize. For any \(\eta>0\), since the SGD update rule can be equivalently written as \[\begin{align} x_{k+1} = x_k + \frac{\beta_k}{\eta} \left( -\eta \nabla J(x_k) + x_k + \eta w_k - x_k \right), \end{align}\] by defining \(\alpha_k = \beta_k / \eta\), \(Y_k = \eta w_k\), and \(F(x,y) = -\eta \nabla J(x) + x + y\), the SGD update can be formulated as an instance of the SA recursion 1 . Moreover, under the assumption that \(J(\cdot)\) is both smooth and strongly convex, with parameters \(L_J\) and \(\sigma_J\), respectively, the operator \(\bar{F}(x) = -\eta \nabla J(x) + x\) is Lipschitz with respect to the Euclidean norm \(\|\cdot\|_2\), with Lipschitz constant \(L_{\mathrm{SGD}} = \max(|1 - \eta \sigma_J|,\, |1 - \eta L_J|)\) [26]. Therefore, when \(\eta \in (0, 2/L_J)\), we have \(L_{\mathrm{SGD}} < 1\), and hence the operator \(\bar{F}(\cdot)\) is a contraction with respect to \(\|\cdot\|_2\). As a result, Theorem 1 applies to SGD in this setting.

2.6.0.2 Linear Stochastic Approximation.

A linear SA is a recursion of the form \[\begin{align} \label{eq:linear95SA} x_{k+1} = x_k + \beta_k (A_k x_k - b_k), \end{align}\tag{8}\] where \(\{A_k\}\) and \(\{b_k\}\) are sequences of i.i.d.random matrices and vectors, respectively, and \(\beta_k\) is the stepsize. Letting \(\bar{A} = \mathbb{E}[A_k]\) and \(\bar{b} = \mathbb{E}[b_k]\), the linear SA can be viewed as iteratively solving the linear system of equations \(\bar{A}x = \bar{b}\). Popular algorithms such as linear regression and a large class of TD-learning-based algorithms in RL, e.g., TD\((0)\), \(n\)-step TD, TD\((\lambda)\), off-policy TD-learning, and TD-learning with linear function approximation, can all be formulated as linear SA of the form presented above.

Assuming that \(\bar{A}\) is Hurwitz, that is, all of its eigenvalues have negative real parts, such a linear SA algorithm can be reformulated as a contractive SA in the form of 1 . Specifically, for any \(\eta > 0\), observe that the update rule 8 is equivalent to \[\begin{align} x_{k+1} = x_k + \frac{\beta_k}{\eta} \big( (I + \eta A_k)x_k - \eta b_k - x_k \big). \end{align}\] Therefore, by defining \(\alpha_k = \beta_k / \eta\), \(Y_k = (A_k, b_k)\), and \(F(x,y)=F(x,A,b) = (I + \eta A)x - \eta b\), the linear SA presented above is a special instance of the generic SA recursion 1 . Moreover, since \(\bar{A}\) is Hurwitz, there exists a unique positive definite matrix \(P\) that solves the Lyapunov equation \(\bar{A}^\top P + P \bar{A} + I = 0\) [20]. As a result, with an appropriately chosen \(\eta\), the operator \(\bar{F}(x) = (I + \eta \bar{A})x - \eta \bar{b}\) is a contraction mapping with respect to the weighted norm \(\|x\|_P = \sqrt{x^\top P x}\). To see this, note that for any \(x_1,x_2\in\mathbb{R}^d\), we have \[\begin{align} \|\bar{F}(x_1) - \bar{F}(x_2)\|_P^2 =\,& (x_1 - x_2)^\top (I + \eta \bar{A})^\top P (I + \eta \bar{A})(x_1 - x_2)\\ =\,& (x_1 - x_2)^\top \big( P + \eta \bar{A}^\top P + \eta P \bar{A} + \eta^2 \bar{A}^\top P \bar{A} \big)(x_1 - x_2)\\ =\,& (x_1 - x_2)^\top \big( P - \eta I + \eta^2 \bar{A}^\top P \bar{A} \big)(x_1 - x_2)\\ \le\,&\left( 1 - \frac{\eta}{\lambda_{\max}(P)} + \frac{\eta^2 \lambda_{\max}(\bar{A}^\top P \bar{A})}{\lambda_{\min}(P)} \right) \|x_1 - x_2\|_P^2\\ =\,&\left( 1 - \frac{\lambda_{\min}(P)}{4 \lambda_{\max}(\bar{A}^\top P \bar{A}) \lambda_{\max}(P)^2} \right) \|x_1 - x_2\|_P^2, \end{align}\] where the last inequality follows by choosing \[\begin{align} \eta = \frac{\lambda_{\min}(P)}{2 \lambda_{\max}(\bar{A}^\top P \bar{A}) \lambda_{\max}(P)}. \end{align}\] Since the operator \(\bar{F}(\cdot)\) is a contraction with respect to \(\|\cdot\|_P\), Theorem 2 applies to linear SA.

2.7 Related Work↩︎

The SA method was originally proposed in [10] for solving root-finding problems under a noisy oracle.

2.7.0.1 The Asymptotic Results.

The asymptotic convergence of SA has been studied extensively in the literature. A central line of work is based on the ODE method, which compares the interpolated trajectory of the discrete-time recursion with the flow of the mean ODE. Under standard stability and noise assumptions, if the stepsizes are non-summable and square summable, and if the relevant equilibrium point of the ODE is globally asymptotically stable, then the SA iterates converge almost surely [16]–[18], [27]. Such almost-sure convergence results are quite general: the ODE can be replaced by a differential inclusion, which is useful for set-valued or nonsmooth mean fields [28], [29], and the classical square-summability condition on the stepsizes can be relaxed in several directions [30], [31].

Beyond almost-sure convergence, a parallel line of work studies asymptotic rates and distributional limits. For smooth mean fields, after linearizing the mean drift \(\bar F(x)-x\) around a locally stable equilibrium \(x^*\), the suitably normalized error process satisfies central-limit-theorem-type results; for instance, under stepsizes of order \(1/k\), the last iterate typically has \(\sqrt{k}\)-scale fluctuations whose covariance is determined by the Jacobian of the mean field and the covariance of the martingale noise [27], [32]–[34]. These results provide asymptotic convergence rates and form the basis for statistical inference based on SA iterates. A major refinement is Polyak–Ruppert averaging, which shows that averaging the iterates can achieve the optimal asymptotic covariance under substantially less delicate stepsize tuning than is required for the last iterate [35], [36]. More recent CLT and asymptotic efficiency results have further extended this theory to controlled Markovian dynamics, parameter-dependent Markovian noise, nonsmooth problems, and decision-dependent data distributions [37]–[40].

2.7.0.2 The Non-Asymptotic Results.

Finite-time analysis has received a lot of attention in the last decade [24], [41], [42]. A central approach for establishing mean-square bounds is the Lyapunov method: one constructs a scalar function that simultaneously captures the stability of the limiting dynamics and permits control of the stochastic and discretization errors. In the strongly stable case, this typically leads to a one-step recursion consisting of a contracting bias term and a variance term. For example, with a constant stepsize, such recursions often yield a geometrically decaying transient plus an \(O(\alpha)\) steady-state error, whereas with a properly tuned diminishing stepsize they yield an \(O(1/k)\) mean-square rate.

For linear SA [24], [43], [44], the Lyapunov function is usually constructed from the solution of a Lyapunov equation associated with the Hurwitz mean matrix [19], [20]. This quadratic structure has led to sharp finite-time mean-square bounds for both last iterates and averaged iterates, including results that clarify the role of constant stepsizes and Polyak–Ruppert averaging. For nonlinear SA, finite-time analysis is more delicate, and existing results typically exploit additional structure, such as stochastic gradient methods and their variants [41], [42], [45]–[47], contractive SA [25], [48]–[50], or SA under dissipativity assumptions [51]. These structural conditions play the role of producing a negative drift. Beyond standard mean-square bounds, finite-time analysis can be carried out in a more refined way by separating the optimization bias, the asymptotic bias, and the variance [52], [53].

Another important direction concerns weakly stable operators, especially nonexpansive operators. The study of SA under nonexpansive maps has a long history in reinforcement learning [54], and recent work has developed finite-time analyses for stochastic approximation with non-expansive operators. In this setting, because strict contraction and uniqueness of the fixed point may fail, the guarantees are often stated in terms of fixed-point residuals rather than mean-square distance to a distinguished fixed point [55]. More recently, finite-sample analyses have also been developed for nonexpansive SA with Markovian noise [56].

3 Applications in Reinforcement Learning↩︎

In this section, we show how to use the SA results presented in the previous section to establish finite-time bounds of popular RL algorithms. RL provides a principled framework for sequential decision-making under uncertainty [2], with broad applications in game playing [57], robotics [58], recommendation systems [59], and large language models (LLMs) [5], inspiring a surge of theoretical research aimed at deepening the mathematical foundations of RL and guiding its practical deployment.

More formally, an RL problem is usually modeled as a Markov decision process (MDP) [60]. In this survey, we consider a discounted infinite-horizon MDP consisting of a finite state space \(\mathcal{S}\), a finite action space \(\mathcal{A}\), a transition kernel \(\{p(\cdot \mid s,a)\}_{(s,a)\in\mathcal{S}\times \mathcal{A}}\), a reward function \(\mathcal{R}:\mathcal{S}\times\mathcal{A}\to[0,1]\), and a discount factor \(\gamma\in(0,1)\). Since we work with finite MDPs, assuming bounded rewards is indeed without loss of generality. Notably, the parameters of the stochastic model, such as the reward function and the transition kernel, are unknown to the agent.

Given a stationary policy \(\pi:\mathcal{S}\to\Delta(\mathcal{A})\), the associated value function \(V^\pi:\mathcal{S}\to\mathbb{R}\) and Q-function \(Q^\pi:\mathcal{S}\times\mathcal{A}\to\mathbb{R}\) are defined as \[\begin{align} V^\pi(s)=\,&\mathbb{E}_\pi\!\left[\sum_{k=0}^\infty\gamma^k\mathcal{R}(S_k,A_k)\,\middle|\, S_0=s\right],\;\forall\,s\in\mathcal{S},\\ Q^\pi(s,a)=\,&\mathbb{E}_\pi\!\left[\sum_{k=0}^\infty\gamma^k\mathcal{R}(S_k,A_k)\,\middle|\, S_0=s, A_0=a\right],\;\forall\,(s,a)\in\mathcal{S}\times \mathcal{A}, \end{align}\] respectively, where we use \(\mathbb{E}_\pi[\cdot]\) to indicate that the actions are selected according to the policy \(\pi\). Alternatively, the value function and the Q-function can be viewed as vectors in \(\mathbb{R}^{|\mathcal{S}|}\) and \(\mathbb{R}^{|\mathcal{S}||\mathcal{A}|}\).

The goal is to find an optimal policy \(\pi^*\) so that its value function \(V^*\), or equivalently its Q-function \(Q^*\), is uniformly maximized [2], [3]; that is, \(V^*(s)\ge V^\pi(s)\) for all \(\pi\) and \(s\), respectively, \(Q^*(s,a)\ge Q^\pi(s,a)\) for all \(\pi\) and \((s,a)\). It is known that such an optimal policy always exists [3].

In RL, the problem of finding an optimal policy is called the control problem. One of the most popular algorithms for solving the control problem is Q-learning [13]. While the ultimate goal is to find an optimal policy, there is usually a smaller goal of finding the value function of a given policy, which is called the prediction problem and is typically solved with TD-learning and its variants [12]. Both Q-learning and TD-learning are, in nature, SA algorithms for solving their associated Bellman equations. Therefore, the results of SA unify the finite-time analysis of value-based RL algorithms. We next present a detailed case study using Q-learning.

3.1 Q-Learning↩︎

The Q-learning algorithm proposed in [13] finds an optimal policy \(\pi^*\) by computing the optimal Q-function \(Q^* = Q^{\pi^*}\). The motivation for Q-learning is based on the following fact: \(\pi\) is an optimal policy if and only if \(\{a \mid \pi(a\mid s)>0\} \subseteq \arg\max_{a\in\mathcal{A}} Q^*(s,a)\) for all \(s\in\mathcal{S}\) [2], [3]. Therefore, finding the optimal Q-function is sufficient for finding an optimal policy.

We next present the Bellman equation for the optimal Q-function \(Q^*\), which serves as the foundation for the design of Q-learning. Let \(\mathcal{H}:\mathbb{R}^{|\mathcal{S}||\mathcal{A}|}\to\mathbb{R}^{|\mathcal{S}||\mathcal{A}|}\) be the Bellman operator defined as \[\begin{align} [\mathcal{H}(Q)](s,a) = \mathcal{R}(s,a) + \gamma\,\mathbb{E}\!\left[\max_{a'\in\mathcal{A}} Q(S_{k+1},a') \;\middle|\; S_k=s, A_k=a\right] \end{align}\] for all \(Q\in\mathbb{R}^{|\mathcal{S}||\mathcal{A}|}\) and \((s,a)\). Then, the Bellman equation is simply the fixed-point equation of the Bellman operator \(\mathcal{H}(\cdot)\): \[\begin{align} Q = \mathcal{H}(Q). \end{align}\] It is well known that \(\mathcal{H}(\cdot)\) is a contraction mapping with respect to \(\|\cdot\|_\infty\), with \(Q^*\) as its unique fixed point [60]. Therefore, a natural approach to finding \(Q^*\) is to perform the fixed-point iteration \(Q_{k+1} = \mathcal{H}(Q_k)\). However, carrying out such a fixed-point iteration requires complete knowledge of the transition probabilities of the underlying MDP, and hence is not feasible in RL when the stochastic model is unknown. The Q-learning algorithm is then proposed as a data-driven stochastic variant of this fixed-point iteration.

Let \(\pi_b\) be the behavior policy, i.e., the policy used to collect samples. Assuming that the Markov chain \(\{S_k\}\) induced by \(\pi_b\) is irreducible, it admits a unique stationary distribution \(\kappa_b \in \Delta(\mathcal{S})\) satisfying \(\kappa_b(s) > 0\) for all \(s \in \mathcal{S}\) [61]. Then, with a sequence of i.i.d.samples \(\{(S_k, A_k, S_k')\}\) generated as \(S_k \sim \kappa_b(\cdot)\), \(A_k \sim \pi_b(\cdot \mid S_k)\), and \(S_k' \sim p(\cdot \mid S_k, A_k)\) for all \(k \ge 0\), the Q-learning algorithm is given in Algorithm 1. In this survey paper, we assume i.i.d.sampling for ease of exposition. Finite-time bounds for Q-learning under Markovian sampling can be established using results on Markovian SA, to be presented in later sections (see Theorem 3).

To establish the finite-time bounds for Q-learning, we first reformulate Q-learning in the form of the SA recursion (1 ). Let \(Y_k=(S_k,A_k,S_k')\in\mathcal{Y}:=\mathcal{S}\times\mathcal{A}\times\mathcal{S}\) for all \(k\ge 0\). It is clear that \(\{Y_k\}\) is an i.i.d.sequence with distribution \(\mu\in\Delta(\mathcal{Y})\) satisfying \(\mu(s,a,s')=\kappa_b(s)\pi_b(a\mid s)p(s'\mid s,a)\). Let \(F:\mathbb{R}^{|\mathcal{S}||\mathcal{A}|}\times\mathcal{Y}\to\mathbb{R}^{|\mathcal{S}||\mathcal{A}|}\) be an operator defined such that, given inputs \(Q\in\mathbb{R}^{|\mathcal{S}||\mathcal{A}|}\) and \(y_0=(s_0,a_0,s_1)\), the \((s,a)\)-th component of the output is \[\begin{align} [F(Q,y_0)](s,a) = \begin{dcases} \mathcal{R}(s_0,a_0)+\gamma\max_{a'\in\mathcal{A}}Q(s_1,a'),&(s,a)=(s_0,a_0),\\ Q(s,a),&(s,a)\neq (s_0,a_0).\\ \end{dcases} \end{align}\] Then, the update equation of the Q-learning algorithm presented in Line \(4\) of Algorithm 1 can be equivalently written as \[\begin{align} Q_{k+1} = Q_k + \alpha_k \left(F(Q_k,Y_k)-Q_k\right),\;\forall\,k\geq 0, \end{align}\] which is in the same form as the SA recursion (1 ).

The next step is to verify that all assumptions required to apply Theorem 1 are satisfied in the context of Q-learning. Assumption 2 is clearly satisfied since \(\{Y_k\}\) is a sequence of i.i.d.samples. We next show that Assumptions 1 and 3 are satisfied. To present the result, let \(D_b\) be an \(|\mathcal{S}||\mathcal{A}| \times |\mathcal{S}||\mathcal{A}|\) diagonal matrix with diagonal components \(\{\kappa_b(s)\pi_b(a \mid s)\}_{(s,a)\in\mathcal{S}\times \mathcal{A}}\). Denote the minimum diagonal entry of \(D_b\) by \(D_{b,\min}\). The proof of the following proposition can be found in [25].

We next present the mean-square bounds for Q-learning. For simplicity, we focus on the case of constant stepsize; the extension to diminishing stepsizes is straightforward in view of Theorem 1.

Beyond Q-learning, most value-based RL algorithms can be modeled as SA algorithms under contractive operators. Consequently, the SA results presented in this survey paper apply to a broad class of RL methods. We briefly discuss a few examples below.

3.2 The Efficiency of Bootstrapping in TD-Learning↩︎

Given a target policy \(\pi\), consider the problem of predicting its value function \(V^\pi\) (or its Q-function \(Q^\pi\)). Such problems are typically solved with TD-learning [12] and its variants. At a high level, TD-learning estimates \(V^\pi\) iteratively using a sample trajectory collected by applying the target policy \(\pi\) to interact with the environment. An important idea for potentially improving the performance of TD-learning is to use bootstrapping, with typical examples including \(n\)-step TD and TD\((\lambda)\) [2]. Intuitively, the degree of bootstrapping is determined by how much the current estimate depends on previous estimates. For instance, in \(n\)-step TD, \(n=1\) corresponds to extreme bootstrapping (the next iterate depends solely on the previous one), whereas \(n=\infty\) corresponds to no bootstrapping at all, in which case the method reduces to Monte Carlo.

Despite strong empirical evidence, theoretically analyzing the effect of bootstrapping remains challenging. It is listed as one of the four major open problems in [65], which states that “While it remains unclear exactly what should or could be proved here, it is clear that this [referring to the efficiency of bootstrapping] is a key open question at the heart of current and future RL.” The results presented here for general contractive SA algorithms allow us to establish finite-time bounds for \(n\)-step TD (and TD\((\lambda)\)) as explicit functions of \(n\) (and \(\lambda\)), thereby providing theoretical insight into the efficiency of bootstrapping. We emphasize that the open problem is not resolved here: to provide a definitive answer, upper bounds alone are insufficient; a matching lower bound is needed.

As an aside, since TD-learning is often used within an actor–critic framework [66] to ultimately compute an optimal policy, the finite-time analyses of TD-learning also enable a deeper understanding of the performance of a variety of actor–critic algorithms.

3.3 The Bias-Variance Trade-Off in Off-Policy TD-Learning↩︎

Given a target policy \(\pi\), the TD-learning algorithm for predicting \(V^\pi\) (or \(Q^\pi\)) can be divided into two categories: on-policy TD and off-policy TD. In on-policy TD-learning, the samples are collected under the target policy \(\pi\), which was discussed in the previous section. In off-policy TD, samples are generated by a behavior policy \(\pi_b \neq \pi\). Off-policy sampling is used for three important reasons: (1) it is typically necessary to incorporate exploration in the behavior policy \(\pi_b\), which makes it differ from the target policy \(\pi\); (2) it is widely used in multi-agent training, where different agents collect rewards using behavior policies that lag behind the target policies in an actor–critic framework [67]; and (3) it enables learning from historical data, thereby improving sample efficiency.

In practice, off-policy TD-learning is implemented via importance sampling to obtain an unbiased estimate of the desired value function. However, the product of importance sampling ratios induces high variance in the estimate [68], which is a fundamental issue in off-policy learning. To mitigate this problem, many variants of off-policy TD-learning have been proposed, including the \(Q^\pi(\lambda)\) algorithm [69], the \(\mathrm{TB}(\lambda)\) algorithm [70], the \(\mathrm{Retrace}(\lambda)\) algorithm [71], and the Q-trace algorithm [25]. These methods modify the importance sampling ratios in various ways to reduce variance, albeit at the cost of introducing bias in the limit. While numerical experiments can compare their empirical performance, the lack of analytical tools has made the bias–variance trade-off in tuning the importance sampling ratios unclear.

The results for generic contractive SA presented in this survey paper enable a unified theoretical analysis of all these algorithms and provide quantitative characterizations of the associated bias–variance trade-offs; see [72] for further details.

4 Beyond Contractive Stochastic Approximation with i.i.d.Noise↩︎

So far, we have focused on SA under a contractive operator with i.i.d. noise sequences. In this section, we discuss several extensions of the SA framework.

4.1 Stochastic Approximation under Markovian Noise↩︎

In RL, the samples, namely trajectories of state-action pairs, are typically collected by an agent implementing a behavior policy to interact with the environment. In this case, the resulting sample trajectory forms a Markov chain rather than an i.i.d.sequence [3]. In this section, we discuss how to extend the finite-time analysis to the Markovian noise setting.

To begin, we state the assumption that formalizes the notion of Markovian noise.

Assumption 4 implies that \(\{Y_k\}\) has a unique stationary distribution \(\mu\in\Delta(\mathcal{Y})\). In this case, the expected operator \(\bar{F}(\cdot):=\mathbb{E}_{Y\sim \mu}[F(\cdot,Y)]\) is defined with respect to the stationary distribution \(\mu\). In addition, \(\{Y_k\}\) mixes at a geometric rate; that is, there exist \(C>0\) and \(\rho\in(0,1)\) such that \(\max_{y\in\mathcal{Y}}\|P^k(y,\cdot)-\mu(\cdot)\|_{\text{TV}} \le C\rho^k\) for all \(k\ge 0\), where \(\|\cdot\|_{\text{TV}}\) denotes the total variation distance and \(P\in\mathbb{R}^{|\mathcal{Y}|\times|\mathcal{Y}|}\) is the transition probability matrix of \(\{Y_k\}\) [61]. Let \(t_\alpha := \min\left\{k \ge 0 : \max_{y \in \mathcal{Y}} \|P^{k}(y,\cdot) - \mu(\cdot)\|_{\text{TV}} \le \alpha\right\}\), which captures how fast the distribution of \(Y_k\) converges to its stationary distribution \(\mu\) and is called the mixing time of the Markov chain \(\{Y_k\}\) with precision \(\alpha\). Note that under Assumption 4, we have \(t_\alpha = \mathcal{O}(\log(1/\alpha))\).

We next state the finite-time bounds for Markovian SA under contractive operators.

For ease of exposition, we consider only constant stepsizes in Theorem 3. See [25] for results under diminishing stepsizes.

Compared with the finite-time bound of SA under i.i.d.noise (cf.Theorem 1 (1)), there are two differences in Theorem 3. The first is that the bound holds only when \(k \ge t_\alpha\) instead of for all \(t \ge 0\). Intuitively, this is because the distribution of \(Y_k\) is not “stable” until \(\{Y_k\}\) is sufficiently mixed. Moreover, the second term on the right-hand side of Inequality (?? ) includes an additional factor of \(t_\alpha = \mathcal{O}(\log(1/\alpha))\), which also arises from the Markovian noise. Other than these two points, Theorem 3 is qualitatively similar to Theorem 1, as the mean squared error consists of a geometrically decaying bias and a constant variance of order \(\mathcal{O}(\alpha \log(1/\alpha))\).

To establish Theorem 3, the high-level idea is identical to that of Theorem 1 in Section 2, namely to use the generalized Moreau envelope as a Lyapunov function and to establish a one-step contractive inequality. After the error decomposition in Inequality 5 , the main challenge under Markovian noise is to control the expected value of the term \(E_2\) from Inequality 5 , which we rewrite below as \(\hat{E}_2\) to avoid confusion: \[\begin{align} \hat{E}_2 = \alpha_k \mathbb{E}[\langle \nabla M(x_k - x^*),\, F(x_k, Y_k) - \bar{F}(x_k) \rangle]. \end{align}\] In view of Lemma 1, since the negative drift is \(\mathcal{O}(\alpha_k)\), we need to show that \(\hat{E}_2 = o(\alpha_k)\) so that it is dominated by the negative drift. When \(\{Y_k\}\) is an i.i.d.sequence, as considered in Section 2, we immediately have \(\hat{E}_2 = 0\) by the tower property of conditional expectations. When \(\{Y_k\}\) is a Markov chain, controlling \(\hat{E}_2\) requires additional effort.

We next present the high-level idea of an approach for bounding \(\hat{E}_2\) using a conditioning argument, which was first developed in [24] for linear SA and later extended in [51] to nonlinear SA. The first step is to rewrite the term \(\hat{E}_2\) as \[\begin{align} \label{eq:markov95decomposition} \hat{E}_2 = \alpha_k \mathbb{E}[\langle \nabla M(x_{k-t} - x^*),\, F(x_{k-t}, Y_k) - \bar{F}(x_{k-t}) \rangle] + \textit{residual terms}, \end{align}\tag{9}\] where \(t\) is a positive integer yet to be chosen, and the residual terms arise from replacing \(x_k\) with \(x_{k-t}\) from the original expression of \(\hat{E}_2\).

For the first term on the right-hand side of (9 ), by applying the tower property of conditional expectations, we have \[\begin{align} &\alpha_k \mathbb{E}[\langle \nabla M(x_{k-t} - x^*),\, F(x_{k-t}, Y_k) - \bar{F}(x_{k-t}) \rangle] \\ =\;& \alpha_k \mathbb{E}[\langle \nabla M(x_{k-t} - x^*),\, \mathbb{E}[F(x_{k-t}, Y_k) \mid \mathcal{F}_{k-t}] - \bar{F}(x_{k-t}) \rangle] \\ \leq \;& \alpha_k \mathcal{O}(C \rho^t), \end{align}\] where \(\mathcal{F}_{k-t}\) is the \(\sigma\)-algebra generated by \(\{Y_n\}_{0 \le n \le k-t-1}\), and the last inequality can be established using the geometric mixing of \(\{Y_k\}\). As for the residual terms in (9 ), since there are \(t\) updates required to obtain \(x_k\) from \(x_{k-t}\), and each update contributes an amount proportional to the stepsize, we can show that \[\begin{align} \textit{residual terms} = \mathcal{O}\!\left(\alpha_k \sum_{\ell = k-t}^k \alpha_\ell \right). \end{align}\] Finally, by choosing \(t = \mathcal{O}(\log(1/\alpha_k))\), we obtain \(\hat{E}_2 = o(\alpha_k)\) as desired. See [24], [51] for more details regarding the conditioning argument for dealing with Markovian noise.

4.1.0.1 Relaxing the Aperiodicity Assumption.

To relax the aperiodicity requirement, the most straightforward approach is to apply an aperiodic transformation to the data samples. Specifically, for a finite Markov chain that is irreducible, the corresponding lazy chain is both irreducible and aperiodic. Moreover, with an aperiodicity transformation, also known as the Schweitzer transformation [73], one can construct a trajectory from the original sample path as if it were generated by the lazy transition matrix.

Specifically, let \(\{Y_k\}\) be a finite and irreducible Markov chain with transition matrix \(P\). Fix some \(\eta\in(0,1)\) and define the lazy transition matrix \(\tilde{P}=(1-\eta)P+\eta I\). Given a realized trajectory \(\{Y_k\}\), we construct a new trajectory \(\{Y_k'\}\) as follows. Let \(Y_0'=Y_0\). For each \(k\ge 0\), insert a geometrically distributed number of copies of \(Y_k\) (with success probability \(1-\eta\)) before moving to \(Y_{k+1}\). In other words, between \(Y_k\) and \(Y_{k+1}\) in the original trajectory, we repeat the state \(Y_k\) a random number of times, and then append \(Y_{k+1}\). The resulting process \(\{Y_k'\}\) is a Markov chain with transition matrix \(\tilde{P}\). Since \(P\) is irreducible and \(\eta>0\), the matrix \(\tilde{P}\) is irreducible and aperiodic. Therefore, by working with the modified trajectory \(\{Y_k'\}\) constructed from the original sample path, one may assume without loss of generality that the underlying Markov chain is aperiodic.

While the aforementioned approach is theoretically sound, it is practically inconvenient, since one needs to repeatedly insert identical samples to justify the analysis. A more direct but technically more challenging approach is to use the Poisson equation to decompose the Markovian noise. For each fixed \(x\), define the forcing function \(h_x(y)=F(x,y)-\bar{F}(x)\). The Poisson equation is then given by \(\hat{h}_x-P\hat{h}_x=h_x\), where \(\hat{h}_x\) is the solution (unique up to an additive constant) [74]. Using the Poisson equation, we write \(F(x_k,Y_k)-\bar{F}(x_k)=h_{x_k}(Y_k)\) from \(\hat{E}_2\) in 9 as \[\begin{align} F(x_k,Y_k)-\bar{F}(x_k)=h_{x_k}(Y_k)=\hat{h}_{x_k}(Y_k)-P\hat{h}_{x_k}(Y_k). \end{align}\] Adding and subtracting \(\hat{h}_{x_k}(Y_{k+1})\), we obtain \[\begin{align} F(x_k,Y_k)-\bar{F}(x_k) =\big(\hat{h}_{x_k}(Y_k)-\hat{h}_{x_k}(Y_{k+1})\big) +\big(\hat{h}_{x_k}(Y_{k+1})-P\hat{h}_{x_k}(Y_k)\big). \end{align}\] The second term satisfies \(\mathbb{E}[\hat{h}_{x_k}(Y_{k+1})-P\hat{h}_{x_k}(Y_k)\mid \mathcal{F}_k]=0\), so it forms a martingale difference sequence, while the first term telescopes across iterations. This representation separates the Markovian dependence from the stochastic fluctuations, allowing the finite-time analysis to proceed by combining martingale arguments with bounds on the Poisson solution \(\hat{h}_x\). A detailed illustration of the Poisson equation-based approach for dealing with Markovian noise can be found in [75]–[77].

4.2 Stochastic Approximation under Seminorm Contractive Operators↩︎

In this section, instead of imposing the contraction assumption, we consider SA under seminorm contractive operators. This setting is motivated by average-reward RL, where the Bellman operator is not contractive in any norm but can be seminorm contractive under certain structural assumptions on the underlying MDP [60], and by control problems in which only partial stability can be established [19]. We begin by defining what we mean by a seminorm.

Recall that a norm \(\|\cdot\|\) must also satisfy \(\|x\| = 0\) if and only if \(x = 0\), a property that is not required for a seminorm. Before proceeding, we present two representative examples of seminorms.

• The span seminorm \(\textit{sp}(\cdot)\) is defined as \(\textit{sp}(x) = \max_i x_i - \min_j x_j\) for all \(x \in \mathbb{R}^d\). One can easily verify that \(\textit{sp}(\cdot)\) is indeed a seminorm. The span seminorm is widely used in studying the convergence behavior of learning algorithms for solving average-reward MDPs [60].

• Let \(P \in \mathbb{R}^{d \times d}\) be a positive semidefinite matrix. Then the function \(p(x) = \sqrt{x^\top P x}\) is a seminorm. Such seminorms are widely used in the analysis of least-squares methods and in control, where quadratic Lyapunov functions naturally induce such seminorm structures [19].

Returning to the SA recursion (1 ), we impose the following assumptions on the operators \(\bar{F}(\cdot)\) and \(F(\cdot,\cdot)\).

Similar to the case of norm-contractive operators, seminorm-contractive operators also enjoy several useful properties. In particular, the set of solutions to the seminorm fixed-point equation \(p_c(\bar{F}(x) - x) = 0\) must be of the form \(x^* + \ker(p_c)\), where \(x^*\) is a particular solution and \(\ker(p_c) = \{x \mid p_c(x) = 0\}\). Moreover, the fixed-point iteration \(x_{k+1} = \bar{F}(x_k)\) converges geometrically in the seminorm: \(p_c(x_k - x^*) \le \gamma_c^k p_c(x_0 - x^*)\). See [78] for more details. We next state the finite-time bounds of SA under seminorm-contractive operators.

Theorem 4 is qualitatively similar to Theorem 1, except that the norm is replaced by the seminorm. The proof of Theorem 4 also follows a Lyapunov approach, where we choose the following variant of the generalized Moreau envelope, \[\begin{align} \label{average95Lyapunov} \hat{M}(x) := \min_{u \in \mathbb{R}^d} \left\{ \frac{1}{2} p_c^2(u) + \frac{1}{2\theta} \|x - u\|_s^2 \right\}, \end{align}\tag{10}\] as the Lyapunov function. The resulting analysis is similar to that of Theorem 1; see [78] for more details.

4.3 Stochastic Approximation under a Dissipativity Assumption↩︎

The following assumption introduces a dissipative operator.

Assumption 7 can be viewed as an exponential dissipativity property of the ODE 2 with the quadratic storage function \(x^\top \bar{P} x\) [19], which motivates referring to Assumption 7 as the dissipativity assumption. A special case in which Assumption 7 holds is linear SA with a Hurwitz matrix, i.e., \[\begin{align} x_{k+1} =\;& x_k + \alpha_k \bigl(A(Y_k)x_k - b(Y_k)\bigr) \\ =\;& x_k + \alpha_k\bigl(\underbrace{(A(Y_k) + I_d)x_k - b(Y_k)}_{F(x_k,Y_k)} - x_k\bigr), \end{align}\] where \(A:\mathcal{Y} \to \mathbb{R}^{d \times d}\) and \(b:\mathcal{Y} \to \mathbb{R}^d\) are deterministic functions. Denote \(\bar{A} = \mathbb{E}_{Y \sim \mu}[A(Y)]\), assumed to be Hurwitz, and \(\bar{b} = \mathbb{E}_{Y \sim \mu}[b(Y)]\). In this case, \(\bar{F}(x) = (\bar{A} + I_d)x - \bar{b}\) and the unique fixed point of \(\bar{F}(\cdot)\) is \(x^* = \bar{A}^{-1}\bar{b}\). Since \(\bar{A}\) is Hurwitz, there exists a unique positive definite matrix \(\bar{P}\) such that \[\begin{align} \bar{A}^\top \bar{P} + \bar{P}\bar{A} + I_d = 0 \end{align}\] [20]. Consequently, for any \(x \in \mathbb{R}^d\), \[\begin{align} (x - x^*)^\top \bar{P}\bigl(\bar{F}(x) - x\bigr) =\;& (x - x^*)^\top \bar{P}(\bar{A}x - \bar{b}) \\ =\;& (x - x^*)^\top \bar{P}\bar{A}(x - x^*) \\ =\;& \tfrac{1}{2}(x - x^*)^\top\bigl(\bar{A}^\top \bar{P} + \bar{P}\bar{A}\bigr)(x - x^*) \\ =\;& -\tfrac{1}{2}(x - x^*)^\top(x - x^*) \\ \le\;& -\tfrac{1}{2\lambda_{\max}(\bar{P})} (x - x^*)^\top \bar{P}(x - x^*), \end{align}\] where \(\lambda_{\max}(\bar{P})\) is the largest eigenvalue of \(\bar{P}\). Hence, Assumption 7 is satisfied.

Under Assumptions 2, 3, and 7, using \(\tilde{M}(x) = x^\top \bar{P} x\) as the Lyapunov function and following the approach in Section 2, we can establish finite-time guarantees for the SA algorithm with update 1 , achieving a qualitatively similar rate of convergence as in Theorem 1. See [51] for more details.

4.4 Other Applications in Reinforcement Learning↩︎

The extensions (SA with Markovian noise, semi-norm contractive SA, and SA under a dissipativity assumption) presented in earlier sections enables studying many other RL algorithms. We will briefly talk about a few in the following.

4.4.1 Q-Learning with Linear Function Approximation↩︎

In practice, RL algorithms are typically combined with function approximation to overcome the curse of dimensionality, which refers to the computational bottleneck as the size of the state-action space increases. One of the most successful practical algorithms is the deep Q-network, which is essentially Q-learning with neural network approximation [79]. On the other hand, the behavior of Q-learning with function approximation is not yet fully understood theoretically and was identified in [65] as an open problem. In fact, the infamous deadly triad [2] arises in Q-learning with function approximation, and hence even in the basic setting with linear function approximation, the algorithm can be unstable in general [80].

The SA results in the previous sections not only allow us to analyze existing RL algorithms but also enable the development of new RL algorithms with theoretical performance guarantees. Specifically, they lead to a variant of Q-learning with linear function approximation that incorporates a target network and truncation, is provably stable, and achieves a \(\tilde{\mathcal{O}}(\epsilon^{-2})\) sample complexity up to a function approximation error [81].

4.4.2 Robust Reinforcement Learning↩︎

Robust RL studies sequential decision-making when the transition model used for training may differ from the transition model encountered at deployment. This is often formulated through robust MDPs, where the goal is to optimize the worst-case value over an ambiguity set of transition kernels. In the discounted setting, under rectangular uncertainty sets, the robust Bellman operator retains a contraction structure. This contraction property is the key reason why value-based robust RL algorithms can be analyzed through the same contractive-SA perspective.

This connection is made explicit in [82]. Although the robust Bellman update involves an inner worst-case optimization problem, the induced robust Bellman operator remains contractive in \(\|\cdot\|_\infty\). Hence, the analysis uses the generalized Moreau-envelope Lyapunov function to establish the one-step drift bound. The same SA viewpoint is used in robust policy optimization. In [83], stochastic robust policy mirror descent relies on repeated robust policy evaluation from samples. The resulting robust TD subroutine is a Markovian contractive SA recursion, so Theorem 3 provides the finite-time critic bound, which is then propagated through the mirror-descent analysis to obtain the overall sample-complexity guarantee. More recently, [84] applied this Lyapunov approach to robust average-reward policy evaluation. Since discounting is absent, the key structural step is to show that the robust average-reward Bellman operator is contractive under a suitable seminorm. The robust TD recursion is then analyzed with a Moreau-envelope-type Lyapunov function, analogous to \(\hat{M}(\cdot)\) in Equation 10 .

4.4.3 Federated Reinforcement Learning↩︎

Federated RL studies settings where multiple agents collaboratively learn a shared value function or \(Q\)-function while keeping their local trajectories, data, and possibly behavior policies private. This setting is natural in RL because sampling is often the computational bottleneck, and parallel data collection across \(N\) agents can in principle reduce the sample complexity.

In [85], the authors extend the Markovian-SA Lyapunov framework to this federated setting. They introduce a general federated stochastic approximation scheme with Markovian noise, which includes federated on-policy TD-learning, off-policy TD-learning, and Q-learning as special cases. The analysis again makes use of the generalized Moreau-envelope Lyapunov function to obtain a one-step drift inequality for the global averaged iterate.

5 Stochastic Approximation: High Probability Bounds↩︎

In this section, we survey results on high-probability bounds for the SA recursion 1 , i.e., bounds on the quantity \(\mathbb{P}(\|x_k - x^*\|_c \le \epsilon)\) for some \(\epsilon > 0\). High-probability bounds are often preferable to mean-square bounds because they provide both a convergence rate and a confidence level. A naive way to obtain such bounds from a mean-square guarantee is to apply Markov’s inequality [86]; however, this yields only polynomially decaying tails. Establishing high-probability bounds with super-polynomial tails (i.e., tails lighter than polynomial tails) is substantially more challenging. For example, even in the classical setting of the law of large numbers,² deriving exponential tail bounds, such as Hoeffding’s inequality, the Chernoff bound, or Bernstein’s inequality, requires significantly more work.

To establish high-probability bounds with super-polynomial tails for the SA recursion 1 , the tail behavior of the effective noise term \(F(x_k, Y_k) - \bar{F}(x_k)\) plays a central role. Although Assumption 2 ensures that this noise is conditionally mean zero, its tail behavior depends on the interplay between the random process \(\{Y_k\}\) and the iterates \(\{x_k\}\), in particular, whether \(Y_k\) is additive or multiplicative with respect to \(x_k\). In the remainder of this section, we focus on the additive (yet unbounded) noise setting. The Lyapunov technique introduced here can be extended to the multiplicative noise case, but the latter requires substantially more technical work.

5.1 Stochastic Approximation under Sub-Gaussian Additive Noise↩︎

Consider the SA recursion 1 . In addition to Assumptions 1 and 2, we impose the following condition on the tail behavior of the effective noise \(F(x_k, Y_k) - \bar{F}(x_k)\). Let \(\mathcal{F}_k\) denote the \(\sigma\)-algebra generated by \(\{Y_n\}_{0 \le n \le k-1}\).

Assumption 8 generalizes the standard notion of a random vector being norm sub-Gaussian [87] to the setting where an arbitrary norm \(\|\cdot\|_c\) is used in place of \(\|\cdot\|_2\). Indeed, when \(\|\cdot\|_c = \|\cdot\|_2\) and \(c_d = d\), Inequalities ?? and ?? specialize exactly to the equivalent definitions of a sub-Gaussian random vector [87], [88]. Since we allow an arbitrary norm, we correspondingly allow a (possibly different) dimension-dependent constant \(c_d\). A representative setting in which Assumptions 2, 1, and 8 all hold is when the noise \(Y_k\) is purely additive and forms an i.i.d.mean-zero sequence with sub-Gaussian tails.

We next state the high-probability bound of the SA recursion (1 ).

The first term captures the deterministic contraction of the initial error, while the remaining terms capture the stochastic fluctuations. The logarithmic factors arise from controlling the martingale terms uniformly over time. Next, we discuss the implications of Theorem 5 in terms of its dependence on \(\delta\), \(K\), and \(k\). Since the tolerance level \(\delta\) appears as \(\log(1/\delta)\) in the norm-square bound, the norm error \(\|x_k - x^*\|_c\) has a sub-Gaussian tail. As for the dependence on \(K\) and \(k\), Theorem 5 implies that, with probability at least \(1-\delta\), all the iterates lie in a cone with radius \(\tilde{\Theta}((\log(1/\delta)^{1/2}+\log(k/K^{1/2}))k^{-1/2})\) for all \(k \ge K\). As byproducts of Theorem 5, one can also establish high-probability bounds for a fixed time \(K\) (rather than for all \(k \ge K\)), as well as the full tail bound on the quantity \(\mathbb{P}(\|x_k - x^*\|_c > \epsilon)\). Consequently, using the formula \[\begin{align} \mathbb{E}[\|x_k - x^*\|_c^r] = \int_0^\infty \mathbb{P}(\|x_k - x^*\|_c^r > x)\, dx \end{align}\] for any positive integer \(r\), we can integrate the tail bound to obtain moment bounds for the error at any fixed time. These results are omitted here.

5.2 High Probability Bounds of Q-Learning↩︎

Theorem 5 enables us to establish the concentration bounds of Q-learning (cf. Algorithm 1) for solving the RL problem, which is presented in the following.

Theorem 6 is qualitatively similar to Theorem 5 in that Q-learning achieves a \(1/k\) convergence rate with a sub-Gaussian tail.

While Theorem 5 is useful for studying Q-learning, it requires the noise to appear in an additive manner. In the context of RL, this limits the applicability of the results to the on-policy and/or tabular setting (i.e., without function approximation). To study RL in the off-policy and/or function approximation setting, we need to analyze SA with multiplicative noise, which is presented in full detail in Appendix 9 for interested readers.

5.3 Proof Sketch of Theorem 5↩︎

Recall that in the establishment of classical concentration inequalities such as Hoeffding’s inequality, the main idea is to bound the moment-generating function (MGF) of the random variable and then use Markov’s inequality to derive the concentration results. Inspired by this approach, we develop a Lyapunov method to prove Theorem 5, where we first bound the MGF of the generalized Moreau envelope and then use Ville’s maximal inequality to establish the maximal concentration bounds.

5.3.1 Bounding the Log-MGF of the Generalized Moreau Envelope↩︎

Let \(\lambda_k = \theta / \alpha_k\), where \(\alpha_k\) is the stepsize and \(\theta > 0\) is a tunable parameter. For any \(k \ge 0\), let \[\begin{align} Z_k = \log\left(\mathbb{E}\left[\exp\left(\lambda_k M(x_k - x^*)\right)\right]\right) \end{align}\] be the log-MGF, where \(M(\cdot)\) is the generalized Moreau envelope introduced in Section 2 as a smooth approximation of the norm-square function \(f(x) = \frac{1}{2}\|x\|_c^2\). The properties of \(M(\cdot)\) were summarized in Proposition [prop:Moreau]. We view \(Z_k\) as our Lyapunov function, and the key step is to derive a bound for \(Z_k\).

Working with \(\mathbb{E}[\exp(\lambda_k M(x_k - x^*))]\) presents new challenges, as the nonlinear nature of the log-MGF prevents us from exploiting the linearity of expectation, which was used extensively in deriving mean-square bounds. Instead, after expressing \(Z_{k+1}\) in terms of \(Z_k\), we obtain the expectation of a product of random variables. To overcome this difficulty, we use a conditioning argument along with Assumption 8 to obtain the following inequality: \[\begin{align} Z_k \le\;& W_1\left(\frac{h}{k+h}\right)^{\alpha D_0/2 - 1} + W_2, \quad \forall\, k \ge 0, \label{sketch:overall} \end{align}\tag{11}\] where \(D_0\), \(W_1\), and \(W_2\) are problem-dependent constants.

5.3.2 An Exponential Supermartingale and Ville’s Maximal Inequality↩︎

Let \(\overline{M}_k = \exp(\lambda_k M(x_k - x^*) - W_3 \sum_{i=0}^{k-1} \alpha_i)\) for all \(k \ge 0\), where \(W_3\) is a constant. Using Inequality (11 ), we show that \(\{\overline{M}_k\}\) is a supermartingale with respect to the filtration \(\{\mathcal{F}_k\}\) [90], which enables us to use Ville’s maximal inequality [90] together with Inequality 11 to establish a maximal concentration inequality. Specifically, for any \(K \ge 0\), we have \[\begin{align} &\mathbb{P}\left(\sup_{k \ge K} \left\{\lambda_k M(x_k - x^*) - W_3 \sum_{i=0}^{k-1} \alpha_i\right\} > \epsilon \right)\\ = \;&\mathbb{P}\left(\sup_{k \ge K} \left\{\exp\left(\lambda_k M(x_k - x^*) - W_3 \sum_{i=0}^{k-1} \alpha_i\right)\right\} > e^\epsilon \right)\\ \leq \;&\exp\left(W_1\left(\frac{h}{K+h}\right)^{\alpha D_0/2 - 1} + W_2 - W_3 \sum_{i=0}^{K-1} \alpha_i - \epsilon\right), \end{align}\] where the last line follows from Ville’s maximal inequality and Inequality 11 . The result follows by rearranging terms in the previous inequality and then connecting \(M(\cdot)\) with \(\|\cdot\|_c^2\) using Proposition [prop:Moreau] (2).

5.4 Related Work↩︎

There is now a substantial literature on high-probability bounds for SA-type algorithms. For stochastic optimization, such results have been developed for SGD and mirror-descent-type methods under various assumptions on convexity, boundedness, and tail behavior of the noise [45], [91]–[95]. For general nonlinear SA, one line of work uses the ODE method and Alekseev’s formula to obtain lock-in or concentration bounds near stable equilibria [96], [97]. Another line focuses on contractive SA and derives maximal concentration bounds for martingale- difference or Markovian noise, with applications to asynchronous Q-learning and TD-learning [75], [89], [98]. More recent work has also studied time-uniform concentration for general iterative algorithms and high-probability bounds for SA with Markovian and heavy-tailed noise [99], [100].

For linear SA, high-probability bounds have been obtained for both last iterates and Polyak–Ruppert averaged iterates, including settings with i.i.d. and Markovian data [44], [101]–[103]. In reinforcement learning, high-probability analyses have been developed for Q-learning and its asynchronous variants [49], [104]–[107], as well as for TD-learning and policy evaluation with linear function approximation [96], [101], [108]–[111]. These results differ along several important dimensions: whether the iterates are projected or otherwise known to remain in a compact set, whether the noise is additive or multiplicative, whether the noise is bounded, sub-Gaussian, Markovian, or heavy-tailed, whether the stepsize depends on the target confidence level or time horizon, and whether the guarantee is fixed-time, finite-horizon uniform, or anytime. See [89] for a more detailed comparison of these assumptions and guarantee formats.

6 Future Research Directions↩︎

There are many potential future directions for this line of work. We highlight three of them below.

6.1 Multiple Time-Scale Stochastic Approximation↩︎

The SA algorithm considered in earlier sections maintains a single sequence of iterates \(\{x_k\}\), updated with stepsizes \(\{\alpha_k\}\). In many applications, however, the algorithm must update several coupled sequences on different time scales. A canonical example is actor–critic learning in RL, where the critic estimates the value function on a faster time scale while the actor updates the policy on a slower time scale. Similar multi-time-scale structures also arise in gradient TD methods with auxiliary variables, constrained RL, bilevel optimization, and multi-agent learning.

The asymptotic theory of multiple time-scale SA has a long history and is typically based on the ODE method; see [16] and the references therein. More recently, general stability and convergence conditions have also been developed for \(N\)-time-scale stochastic recursions [112]. In comparison, finite-time theory is more fragmented. Existing results cover several important special cases, including linear two-time-scale SA, two-timescale TD-learning, and actor–critic-type algorithms [101], [113]–[116]. Recent work has further extended the Lyapunov/Moreau-envelope approach to nonlinear two-time-scale SA with arbitrary-norm contractions and Markovian noise [117], while another line develops accelerated multi-time-scale schemes that achieve \(\widetilde{\mathcal{O}}(1/k)\) rates under strong monotonicity-type assumptions [118].

Despite the progress, it remains unclear whether the obtained convergence bounds are tight in terms of their dependence on rates and other important problem-dependent parameters. For example, despite the extensive existing analyses of actor-critic methods, there are no results showing that two-timescale actor-critic can achieve the statistical lower bound [63]. As another example, consider the simple two-time-scale recursion \[\begin{align} x_{k+1}=\;&x_k+\alpha_k(F(x_k,Y_k)-x_k),\tag{12}\\ \bar{x}_{k+1}=\;&\bar{x}_k+\frac{1}{k+1}(-\bar{x}_k+x_k),\tag{13} \end{align}\] where 13 computes the Polyak average \(\bar{x}_k=k^{-1}\sum_{i=0}^{k-1}x_i\). Polyak averaging is known to have strong asymptotic efficiency and robustness properties [36]; for example, it can achieve the optimal asymptotic \(\mathcal{O}(1/k)\) mean-square rate under stepsizes such as \(\alpha_k=\alpha/(k+h)^z\) with \(z\in(0,1)\). However, finite-time explanations of this robustness phenomenon were obtained only relatively recently in linear SA [43], [44], [103].

In general, a sharp and unified finite-time theory for general multiple time-scale SA remains incomplete.

6.2 Stochastic Approximation with Rapidly Time-Inhomogeneous Markovian Noise↩︎

The Markovian-noise theory discussed in Section 4.1 assumes that the driving process \(\{Y_k\}\) is a time-homogeneous Markov chain with a fixed transition matrix \(P\). In many adaptive algorithms, however, the sampling process is itself affected by the current iterate or by another evolving decision variable. A generic model takes the form \[\begin{align} x_{k+1}=x_k+\alpha_k\big(F(x_k,Y_k)-x_k\big),\qquad \mathbb{P}(Y_{k+1}=y'\mid \mathcal{F}_k,Y_k=y)=P_{\theta_k}(y,y'), \end{align}\] where \(\theta_k\) may be a function of \(x_k\), e.g., a policy in RL, a strategy profile, or a belief vector in learning in games. Thus, the noise process is Markovian only conditionally on the current value of \(\theta_k\), and the transition kernel changes over time.

When the map \(\theta\mapsto P_\theta\) is Lipschitz and \(\theta_k\) evolves slowly, existing controlled-Markov-noise and Poisson-equation techniques can often be adapted. Indeed, the Lipschitz dependence of \(P_\theta\), the invariant distribution, and the corresponding Poisson-equation solution allows one to control the error caused by replacing \(P_{\theta_k}\) with an earlier kernel \(P_{\theta_{k-t}}\). This is the slowly time-inhomogeneous regime, and it is closely related to classical SA with parameter-dependent Markov chains [16], [27], [119].

A substantially more difficult regime arises when the transition kernel changes rapidly with the iterate. In this case, \(\|P_{\theta_{k+1}}-P_{\theta_k}\|\) need not be of order \(\|x_{k+1}-x_k\|\), and hence need not be of order \(\alpha_k\). The usual mixing-time or Poisson-equation residual terms can then contain order-one jumps, rather than small perturbations accumulated over a mixing window. This phenomenon is not specific to reinforcement learning. It is a general SA issue that appears whenever the sampling kernel is generated by a nonsmooth, discontinuous, or set-valued decision rule.

Two motivating examples illustrate the difficulty. In RL with \(\epsilon\)-greedy exploration, the behavior policy is a discontinuous function of the current \(Q\)-estimate: a small perturbation of \(Q\) can change the identity of the greedy action and hence produce an order-one change in the state-action transition matrix. This discontinuity is closely related to the chattering and policy-oscillation phenomena observed in value-based RL with \(\epsilon\)-greedy exploration [120]. A similar issue arises in learning in games. In fictitious play and best-response-type dynamics, strategies may switch abruptly when a belief or payoff estimate crosses an indifference surface. When the play of the agents determines the sampling distribution or the state transition law, such strategy switches induce rapidly time-inhomogeneous Markovian noise in the underlying SA recursion. This issue is especially natural in stochastic games, where fictitious-play-type algorithms combine belief updates, value-function estimates, and greedy or best-response decisions [121]–[123].

Recent work has started to address special cases of this problem. For example, [77] develops a finite-time analysis of Q-learning with time-varying learning policies using a Poisson-equation decomposition and sensitivity bounds, while [124] studies SA with fast-changing transition functions in the analysis of linear \(Q\)-learning with \(\epsilon\)-softmax exploration. High-probability bounds for time-inhomogeneous Markovian SA have also begun to appear in the analysis of online Q-learning [125]. Nevertheless, these results do not yet constitute a general Lyapunov theory for SA with rapidly time-inhomogeneous Markovian noise. Developing such a theory, especially for discontinuous or set-valued kernel maps such as \(\epsilon\)-greedy policies and best-response dynamics, remains an important open direction.

6.3 Stochastic Approximation under Nonexpansive Operators↩︎

Another natural direction is to go beyond contractive or seminorm-contractive operators and consider SA driven by nonexpansive operators. An operator \(\bar F(\cdot)\) is nonexpansive with respect to a norm \(\|\cdot\|_c\) if \[\begin{align} \|\bar F(x_1)-\bar F(x_2)\|_c\leq \|x_1-x_2\|_c, \qquad \forall x_1,x_2\in\mathbb{R}^d . \end{align}\] This setting is substantially more delicate than the contractive case. Fixed points need not exist; even when they exist, they need not be unique; and the plain fixed-point iteration \(x_{k+1}=\bar F(x_k)\) need not converge to a fixed point. For example, \(\bar F(x)=x+1\) has no fixed point, \(\bar F(x)=x\) has every point as a fixed point, and a reflection map may cycle rather than converge.

Suppose that the fixed-point set \(\mathcal{X}=\{x:\bar F(x)=x\}\) is nonempty. In the Euclidean setting, a standard way to restore convergence is to use the averaged operator \[\begin{align} \bar F_\eta(x)=(1-\eta)x+\eta \bar F(x),\qquad \eta\in(0,1), \end{align}\] which has the same fixed-point set as \(\bar F\). The resulting Krasnosel’skii–Mann iteration converges to a fixed point under suitable assumptions, and deterministic rates are typically stated in terms of the fixed-point residual, e.g., \[\begin{align} \min_{0\leq i\leq k}\|\bar F_\eta(x_i)-x_i\|_2^2=\mathcal{O}(1/k). \end{align}\] Finite-sample analyses of stochastic fixed-point iterations for nonexpansive maps have been developed for martingale-difference noise and arbitrary finite dimensional normed spaces [55]. More recently, [56] studied nonexpansive SA with Markovian noise, using Poisson-equation decompositions to obtain both asymptotic and finite-sample results, with applications to tabular average-reward TD-learning.

Despite this progress, the convergence rates obtained in [55], [56] imply a sample complexity of \(\tilde{\mathcal{O}}(\epsilon^{-10})\) in terms of the expected residual. An important future direction is to determine whether this sample complexity can be improved, identify the corresponding fundamental lower bounds, and develop sharp, possibly optimal finite-time rates under general norms, Markovian or controlled Markovian noise, multiplicative or unbounded noise, and high-probability criteria.

7 Conclusion↩︎

In this tutotial paper, we surveyed Lyapunov-based techniques for the finite-time analysis of stochastic approximation algorithms under contractive operators. We explained how generalized Moreau envelopes yield smooth Lyapunov functions for arbitrary contraction norms, and how the resulting drift arguments lead to mean-square and high-probability guarantees. We also discussed extensions to Markovian noise, seminorm-contractive operators, and dissipative operators, with applications to reinforcement learning algorithms such as Q-learning and TD-learning. We conclude by highlighting several open directions.

8 Proof of All Technical Results in Section 2↩︎

8.1 Proof of Proposition [prop:Moreau]↩︎

(1) It is clear from the definition of \(M_f^{\theta,g}(\cdot)\) that it is non-negative and is equal to zero if and only if \(x=0\). Now for any \(c\neq 0\) (the case \(c=0\) is immediate), we have \[\begin{align} M_f^{\theta,g}(c x)&=\min_{u}\left\{\frac{1}{2}\|u\|_c^2+\frac{1}{2\theta}\|c x-u\|_s^2\right\}\\ &= \min_{v}\left\{\frac{1}{2}\|c v\|_c^2+\frac{1}{2\theta}\|c x-c v\|_s^2\right\}\\ &=|c|^2M_f^{\theta,g}(x). \end{align}\] Therefore, we have \(\sqrt{M_f^{\theta,g}(c x)}=|c|\sqrt{M_f^{\theta,g}(x)}\). We next show the triangle inequality. For any \(x_1,x_2\in\mathbb{R}^d\), let \[\begin{align} u_1\in\,&\arg\min_{u\in\mathbb{R}^d}\left\{\frac{1}{2}\|u\|_c^2+\frac{1}{2\theta}\|x_1-u\|_s^2\right\},\\ u_2\in\,&\arg\min_{u\in\mathbb{R}^d}\left\{\frac{1}{2}\|u\|_c^2+\frac{1}{2\theta}\|x_2-u\|_s^2\right\}. \end{align}\] Then we have \[\begin{align} &M_f^{\theta,g}(x_1+x_2)\\ =\;&\min_{u}\left\{\frac{1}{2}\|u\|_c^2+\frac{1}{2\theta}\|x_1+x_2-u\|_s^2\right\}\\ \leq \;& \frac{1}{2}\|u_1+u_2\|_c^2+\frac{1}{2\theta}\|x_1+x_2-u_1-u_2\|_s^2\\ \leq \;&\frac{1}{2}(\|u_1\|_c+\|u_2\|_c)^2+\frac{1}{2\theta}(\|x_1-u_1\|_s+\|x_2-u_2\|_s)^2\\ =\;&M_f^{\theta,g}(x_1)+M_f^{\theta,g}(x_2)+ \|u_1\|_c\|u_2\|_c+\frac{1}{\theta}\|x_1-u_1\|_s\|x_2-u_2\|_s\\ \leq\;& M_f^{\theta,g}(x_1)+M_f^{\theta,g}(x_2)\\ &+2\sqrt{\frac{1}{2}\|u_1\|_c^2+\frac{1}{2\theta}\|x_1-u_1\|_s^2}\sqrt{\frac{1}{2}\|u_2\|_c^2+\frac{1}{2\theta}\|x_2-u_2\|_s^2}\\ =\;&M_f^{\theta,g}(x_1)+M_f^{\theta,g}(x_2)+2\sqrt{M_f^{\theta,g}(x_1)M_f^{\theta,g}(x_2)}. \end{align}\] It follows that \(\sqrt{M_f^{\theta,g}(x_1+x_2)}\leq \sqrt{M_f^{\theta,g}(x_1)}+\sqrt{M_f^{\theta,g}(x_2)}\) for any \(x_1,x_2\in\mathbb{R}^d\). Therefore, \(M_f^{\theta,g}(\cdot)\) is a norm-square function and we can write \(M_f^{\theta,g}(x)=\frac{1}{2}\|x\|_m^2\) for some norm \(\|\cdot\|_m\).

(2) We first derive the upper bound. By definition of \(M_f^{\theta,g}(\cdot)\), we have \[\begin{align} M_f^{\theta,g}(x)&=\min_{u\in\mathbb{R}^d}\left\{\frac{1}{2}\|u\|_c^2+\frac{1}{2\theta}\|x-u\|_s^2\right\}\\ &\geq \min_{u\in\mathbb{R}^d}\left\{\frac{1}{2}\|u\|_c^2+\frac{1}{2\theta u_{cs}^2}\|x-u\|_c^2\right\}\|\cdot\|_s}\\ &\geq \min_{u\in\mathbb{R}^d}\left\{\frac{1}{2}\|u\|_c^2+\frac{1}{2\theta u_{cs}^2}(\|x\|_c-\|u\|_c)^2\right\}\\ &=\min_{y\in\mathbb{R}}\left\{\frac{1}{2}y^2+\frac{1}{2\theta u_{cs}^2}(\|x\|_c-y)^2\right\}\\ &=\min_{y\in\mathbb{R}}\left\{\left(\frac{1}{2}+\frac{1}{2\theta u_{cs}^2}\right)y^2-\frac{1}{\theta u_{cs}^2}\|x\|_cy+\frac{1}{2\theta u_{cs}^2}\|x\|_c^2\right\}\\ &=\frac{1}{2}\|x\|_c^2\frac{1}{\theta u_{cs}^2+1}\\ &=\frac{1}{\theta u_{cs}^2+1}f(x). \end{align}\] It follows that \(f(x)\leq \left(1+\theta u_{cs}^2\right) M_f^{\theta,g}(x)\) for all \(x\), which implies \(\|\cdot\|_c\leq (1+\theta u_{cs}^2)^{1/2}\|\cdot\|_m\).

Next we show the lower bound. Similarly, by definition we have for
any $x\in\mathbb{R}^d$ that
$$\begin{align} M_f^{\theta,g}(x)&=\min_{u\in\mathbb{R}^d}\left\{\frac{1}{2}\|u\|_c^2+\frac{1}{2\theta}\|x-u\|_s^2\right\}\\ &\leq \min_{\alpha\in(0,1)}\left\{\frac{1}{2}\|\alpha x\|_c^2+\frac{1}{2\theta}\|x-\alpha x\|_s^2\right\}\\ &\leq \frac{1}{2}\|x\|_c^2\min_{\alpha\in(0,1)}\left\{\alpha^2+\frac{(1-\alpha)^2}{\theta \ell_{cs}^2}\right\}\|\cdot\|_s\leq \|\cdot\|_c}\\ &=\frac{1}{1+\theta\ell_{cs}^2}\frac{1}{2}\|x\|_c^2\\ &= \frac{1}{1+\theta\ell_{cs}^2}f(x).
\end{align}$$ It follows that
$f(x)\geq \left(1+\theta \ell_{cs}^2\right)M_f^{\theta,g}(x)$ for
all $x$, which implies
$\|\cdot\|_c\geq (1+\theta \ell_{cs}^2)^{1/2}\|\cdot\|_m$.

(3) The convexity of \(M_f^{\theta,g}(\cdot)\) follows from [1]. Since \(f(\cdot)\) is proper, closed, and convex, and \(g(\cdot)\) is \(L\) – smooth with respect to \(\|\cdot\|_s\), we have by [1] that \(M_f^{\theta,g}(\cdot)\) is \(\frac{L}{\theta}\) – smooth with respect to \(\|\cdot\|_s\).

8.2 Proof of Lemma 1↩︎

Using the fact that \(\bar{F}(x^*)=x^*\), we have \[\begin{align} \label{eq:split-E-terms} &\langle \nabla M_f^{\theta,g}(x_k-x^*),\bar{F}(x_k)-x_k\rangle\nonumber\\ =\;&\underbrace{\langle \nabla M_f^{\theta,g}(x_k-x^*),\bar{F}(x_k)-\bar{F}(x^*)\rangle}_{E_{1,1}}-\underbrace{\langle \nabla M_f^{\theta,g}(x_k-x^*),x_k-x^*\rangle}_{E_{1,2}}. \end{align}\tag{14}\] For the gradient of \(M_f^{\theta,g}(x)\), since \(M_f^{\theta,g}(x)=\frac{1}{2}\|x\|_m^2\) is a smooth function, we have by the chain rule of calculus that \(\nabla M_f^{\theta,g}(x)=\|x\|_m \nabla \|x\|_m\).

Now consider the term \(E_{1,1}\). Using Hölder’s inequality, we have \[\begin{align} E_{1,1}&=\|x_k-x^*\|_m\langle \nabla \|x_k-x^*\|_m,\bar{F}(x_k)-\bar{F}(x^*)\rangle\nonumber\\ &\leq \|x_k-x^*\|_m\|\nabla \|x_k-x^*\|_m\|_m^*\|\bar{F}(x_k)-\bar{F}(x^*)\|_m,\label{eq:dual-norm} \end{align}\tag{15}\] where \(\|\cdot\|_m^*\) is the dual norm of \(\|\cdot\|_m\). To further control \(E_{1,1}\), the following result from [126] is needed.

Since \(\|x\|_m\) as a function of \(x\) is \(1\) – Lipschitz with respect to \(\|\cdot\|_m\), we have by Lemma 4 that \(\|\nabla \|x_k-x^*\|_m\|_m^*\leq 1\). For the term \(\|\bar{F}(x_k)-\bar{F}(x^*)\|_m\) in Inequality (15 ), using Proposition [prop:Moreau] (2) and the contraction property of \(\bar{F}(\cdot)\) with respect to \(\|\cdot\|_c\), we have \[\begin{align} \|\bar{F}(x_k)-\bar{F}(x^*)\|_m\leq \,&\frac{1}{\ell_{cm}}\|\bar{F}(x_k)-\bar{F}(x^*)\|_c\\ \leq \,&\frac{\gamma_c}{\ell_{cm}}\|x_k-x^*\|_c\\ \leq \,&\frac{\gamma_c u_{cm}}{\ell_{cm}}\|x_k-x^*\|_m. \end{align}\] Substituting the upper bounds we obtained for \(\|\nabla \|x_k-x^*\|_m\|_m^*\) and \(\|\bar{F}(x_k)-\bar{F}(x^*)\|_m\) into Inequality (15 ), we have \[\begin{align} E_{1,1} \leq \frac{\gamma_c u_{cm}}{\ell_{cm}}\|x_k-x^*\|_m^2=\frac{2\gamma_c u_{cm}}{\ell_{cm}} M_f^{\theta,g}(x_k-x^*). \end{align}\]

Now consider the term \(E_{1,2}\) in Inequality (14 ). Since the norm \(\|\cdot\|_m\) is a convex function of \(x\), we have by definition of convexity that \[\begin{align} \|0\|_m-\|x_k-x^*\|_m\geq \langle \nabla \|x_k-x^*\|_m,-(x_k-x^*)\rangle, \end{align}\] which implies \[\begin{align} E_{1,2}=\|x_k-x^*\|_m\langle \nabla \|x_k-x^*\|_m,x_k-x^*\rangle\geq \|x_k-x^*\|_m^2=2M_f^{\theta,g}(x_k-x^*). \end{align}\] Combining the bounds on \(E_{1,1}\) and \(E_{1,2}\), we obtain \[\begin{align} E_1=\alpha_k (E_{1,1}-E_{1,2})\leq -2\left(1-\gamma_c\frac{u_{cm}}{\ell_{cm}}\right)\alpha_k M_f^{\theta,g}(x_k-x^*). \end{align}\]

8.3 Proof of Lemma 2↩︎

We first show that \(\|F(x,y)\|_c\) can grow at most affinely with respect to \(\|x\|_c\). Indeed, by the triangle inequality, we have \[\begin{align} \label{eq:affine95growth} \|F(x,y)\|_c \leq \|F(x,y) - F(0,y)\|_c + \|F(0,y)\|_c \leq L_1 \|x\|_c + L_2. \end{align}\tag{16}\] Further, by Jensen’s inequality, we have \[\begin{align} \|\bar{F}(x)\|_c\leq L_1\|x\|_c+L_2. \end{align}\] To bound the term \(E_3\), for any \(k\geq 0\), we have \[\begin{align} E_3&=\frac{L\alpha_k^2}{2\theta}\|F(x_k,Y_k)-x_k\|_s^2\\ &\leq \frac{L\alpha_k^2}{2\theta\ell_{cs}^2}\|F(x_k,Y_k)-x_k\|_c^2 (2)}\\ &\leq \frac{L\alpha_k^2}{2\theta\ell_{cs}^2}\left(\|F(x_k,Y_k)\|_c+\|x_k\|_c\right)^2\\ &\leq \frac{L\alpha_k^2}{2\theta\ell_{cs}^2}\left((L_1+1)\|x_k\|_c+L_2\right)^2)}\\ &\leq \frac{L\alpha_k^2}{2\theta\ell_{cs}^2}\left((L_1+1)\|x_k-x^*\|_c+(L_1+1)\|x^*\|_c+L_2\right)^2\\ &\leq \frac{L\alpha_k^2}{\theta\ell_{cs}^2}(L_1+1)^2\|x_k-x^*\|_c^2+\frac{L\alpha_k^2}{\theta\ell_{cs}^2}\left((L_1+1)\|x^*\|_c+L_2\right)^2\\ &\leq \frac{2L(L_1+1)^2u_{cm}^2\alpha_k^2}{\theta\ell_{cs}^2}M_f^{\theta,g}(x_k-x^*)+\frac{L\alpha_k^2}{\theta\ell_{cs}^2}\left((L_1+1)\|x^*\|_c+L_2\right)^2, \end{align}\] where the last line follows from Proposition [prop:Moreau] (2).

8.4 Proof of Lemma 3↩︎

Using Lemma 1, Equation (6 ), and Lemma 2 in Inequality (5 ), we have for all \(k\geq 0\) that \[\begin{align} M_{k+1}\leq\,& \left(1- 2\left(1 - \frac{\gamma_c u_{cm}}{\ell_{cm}}\right)\alpha_k+\frac{2L(L_1+1)^2u_{cm}^2\alpha_k^2}{\theta\ell_{cs}^2}\right)M_k\\ &+\frac{L\alpha_k^2}{\theta\ell_{cs}^2}\left((L_1+1)\|x^*\|_c+L_2\right)^2\\ \leq \,&\left(1-\left(1 - \frac{\gamma_c u_{cm}}{\ell_{cm}}\right)\alpha_k \right)M_k+\frac{L\alpha_k^2}{\theta\ell_{cs}^2}\left((L_1+1)\|x^*\|_c+L_2\right)^2, \end{align}\] where the last inequality follows from choosing \[\begin{align} \alpha_0\leq \varphi_0=\frac{\theta \ell_{cs}^2\left(1 - \gamma_c u_{cm}/\ell_{cm}\right)}{2 L (L_1+1)^2 u_{cm}^2 }. \end{align}\]

9 Stochastic Approximation under Bounded Multiplicative Noise↩︎

The following assumption explains what we mean by multiplicative noise.

Assumption 9 is in fact a relaxed version of Assumption 3. To see this, suppose that Assumption 3 is satisfied. Then we have the affine growth property (cf. Equation (16 )), which in turn implies that \[\begin{align} \| F(x_k,Y_k) - \Bar{F}(x_k) \|_c\leq \;&\| F(x_k,Y_k)\|_c +\| \Bar{F}(x_k) \|_c\\ = \;&\| F(x_k,Y_k)\|_c +\| \mathbb{E}_{Y\sim \mu}[F(x_k,Y)] \|_c\\ \leq \;&\| F(x_k,Y_k)\|_c +\mathbb{E}_{Y\sim \mu}[\| F(x_k,Y) \|_c]\\ \leq\;& 2L_1\|x_k\|_c+2L_2. \end{align}\] Therefore, Assumption 9 is satisfied with \(\sigma=2\max(L_1,L_2)\).

The almost sure boundedness in Assumption 9 is necessary for us to establish high probability bounds with super-polynomial tail. To see this, consider scalar-valued linear SA \(x_{k+1}=x_k+\alpha_k\tilde{a}_kx_k\), and suppose that \(\tilde{a}_k\sim \mathcal{N}(-1/2,1)\). Note that \(\tilde{a}_k\) is negative in expectation but is not bounded. It is easy to see that \(x_k=x_0\prod_{j=0}^{k-1}(1+\alpha_j \tilde{a}_j)\). Since \(1+\alpha_k \tilde{a}_k\sim \mathcal{N}(1-\alpha_k/2,\alpha_k^2)\), the tail of \(x_k\) gets heavier and heavier as \(k\) increases. For example, \(x_1\) has sub-Gaussian tail, \(x_2\) has a \(\mathcal{X}^2\)-distribution (which only has sub-exponential tail), \(x_3\) has no MGF, etc. One special case where Assumption 9 is satisfied is when the operator \(F(x,y)\) is Lipschitz continuous in \(x\) (cf. Assumption 3). The implication easily follows from using triangle inequality.

Observe that under Assumptions 1 and 9, we have by triangle inequality that \[\begin{align} \|x_{k+1}\|_c\leq\;& (1-\alpha_k)\|x_k\|_c+\alpha_k\|F(x_k,Y_k)\|_c\\ \leq \;&(1-\alpha_k)\|x_k\|_c+\alpha_k\|F(x_k,Y_k)-\bar{F}(x_k)\|_c\\ &+\alpha_k\|\bar{F}(x_k)-\bar{F}(0)\|_c+\alpha_k\|\bar{F}(0)\|_c\\ \leq \;&(1+(\sigma+\gamma_c-1)\alpha_k)\|x_k\|_c+\alpha_k (\sigma+\|\bar{F}(0)\|_c). \end{align}\] Note that the parameter \(\gamma_c\) captures the contraction effect of the expected operator and the parameter \(\sigma\) captures the expansive effect of the noise. The combined effect is captured by the parameter \(D:=\sigma+\gamma_c-1\). When using \(\alpha_k=\alpha/(k+h)\), the previous inequality implies that \(\|x_k\|_c\) is either uniformly bounded by a deterministic constant or grow at most logarithmically when \(D\leq 0\). However, when \(D>0\), \(\|x_k\|_c\) can grow at a polynomial rate \(\mathcal{O}(k^{\alpha D})\). See Proposition [prop:worst95case95bound] for more details.

We next state the result. For ease of exposition, we only state the result for \(D>0\) and \(D<0\). The case \(D=0\) is a straightforward extension. We use \(\alpha_k=\alpha/(k+h)\) as our stepsize, and assume that \(2\alpha D\) is an integer, which is in fact without loss of generality because \(D=\sigma+\gamma_c-1\) and if Assumption 3 holds with some \(\sigma>0\) it also holds for any \(\sigma'>\sigma\). Let \(m=2\alpha D+1\). The parameters \(\{c_i\}_{1\leq i\leq 4}\), \(c_1'\), and \(D_0\in (0,1)\) used to present the following theorem are constants.

Several remarks are in order. We first discuss the convergence rate in \(k\). Since the dominant polynomial term in \(k\) is \(\mathcal{O}(1/k)\) in both cases of Theorem 7 (i.e., \(D<0\) and \(D>0\)), the norm-square error achieves \(\mathcal{O}(1/k)\) rate of convergence. This agrees with the mean-square bound presented in Section 2, and existing results on SGD [41].

We next discuss the tail. In the case \(D<0\), Theorem 7 (1) states that Algorithm (1 ) achieves \(\mathcal{O}(1/k)\) rate of convergence with sub-Gaussian tail, as \(\delta\) appears as \(\log(1/\delta)\) in the norm-square bound. This is not surprising as \(D<0\) implies that \(\|x_k\|_c\) is uniformly bounded by a deterministic constant (cf. Proposition [prop:worst95case95bound]). Therefore, we in fact only have bounded (hence sub-Gaussian) additive noise.

The case where \(D>0\) is more complicated. In this case, the tail depends on the parameter \(m\). Since \(m=2\alpha D+1\) and \(D>0\), we in general only have super-polynomial tail since the probability tolerance level \(\delta\) appears as \([\log(m/\delta)]^m\) in the bound. The fact that \(m\) is affine in \(D\) makes intuitive sense as larger \(D\) implies noisier updates, which in turn implies heavier tail. In terms of the dependency on \(k\) and \(K\), Theorem 7 (2) states that, with probability at least \(1-\delta\), all the iterates lie in a cone that starts with radius \(\tilde{\Theta}((1+\log^{m/2}(1/\delta))K^{-1/2})\) when \(k=K \geq 0\), and then (for all \(k>K\)) its radius is of order \(\tilde{\Theta}((\log^{m/2}(1/\delta)+\log^{1/2}(k/K))k^{-1/2})\). Moreover, for small values of \(k\), this bound can be tightened by an a.s. bound that is polynomial in \(k\) (cf. Proposition [prop:worst95case95bound]). Note that, as a function of \(k\), the initial radius is always of order at most \(K^{-1}\), matching the rate obtained for the mean-square error in [25]. On the other hand, the radius decays at only a slightly slower rate than the initial radius as a function of \(k\).

Anytime concentration bounds immediately imply concentration bounds for a fixed iteration number (cf. Corollary 1), which in turn gives the full tail bound (cf. Corollary 2) and the sample complexity result (cf. Corollary 3). For ease of exposition, we only state the results when \(D>0\).

Corollary 1 follows by setting \(K=k\) in Theorem 7 (2).

Corollary 2 is a direct implication of Corollary 1, and provides an upper bound for the whole complementary cumulative distribution function (CDF) of the error \(\|x_k-x^*\|_c\) for any iterate \(k\geq 0\), which can be integrated to obtain bounds for any moment of the error at any point in time.

As we see from Corollary 3, the sample complexity dependency on \(\epsilon\) is \(\tilde{\mathcal{O}}(\epsilon^{-2})\), which is known to be optimal (up to a logarithmic factor). In addition, we have super-polynomial tail as \(\delta\) appears as \(\log^m(1/\delta)\) in the bound.

Theorem 7 shows that SA with multiplicative noise in general is able to achieve an \(\tilde{\mathcal{O}}(1/k)\) rate of convergence with a super-polynomial tail. One may ask if sub-Gaussian (or sub-exponential) tail is achievable. In fact, for SA with multiplicative noise, it is in general not possible to obtain sub-exponential tail bound. See [89] for more details.

9.1 Proof Sketch of Theorem 7↩︎

The main challenge of obtaining super-polynomial high probability bounds is due to the combination of unbounded iterates and multiplicative noise. While having unbounded iterates and multiplicative noise are not too problematic in isolation, the combination of both creates a setting where the variance of the noise is unbounded. In this case, since we allow the multiplicative noise to be large enough so that the “noisy” operator can be expansive with positive probability, the error can grow extremely fast with a significant probability. This creates a challenge that no approach in the literature can deal with in general. To overcome this challenge, we develop a novel bootstrapping argument. The high level ideas are presented in the following.

9.1.1 Initialization: Time-Varying Worst-Case Bounds↩︎

While the iterates of SA with multiplicative noise are not uniformly bounded by a constant, we show in the following proposition that they do admit a time-varying a.s. bound. The behavior of such time-varying bound depends on the contraction effect in the expected operator and the expansive effect in the multiplicative noise. In general, the bound can be polynomially increasing with time.

9.1.2 Bootstrapping: An Iterative Framework to Improve the Bound↩︎

The key in the bootstrapping argument is to start with a non-decreasing sequence \(\{T_k(\delta)\}_{k\geq 0}\) such that \[\begin{align} \mathbb{P}(\|x_k-x^*\|_c^2\leq T_k(\delta),\forall\;k\geq 0)\geq 1-\delta, \end{align}\] and obtain a sequence \(\{T_k(\delta,\delta')\}_{k\geq 0}\), with \(T_k(\delta,\delta')=\tilde{\mathcal{O}}(T_k(\delta)/k)\), such that \[\begin{align} \mathbb{P}(\|x_i-x^*\|_c^2\leq T_k(\delta,\delta'),\forall\;k\geq 0)\geq 1-\delta-\delta'. \end{align}\] This blueprint enables us to start with the time-varying worst-case bound for the error (which can be polynomially increasing) and iteratively improve it to obtain our super-polynomial concentration bound with the desired convergence rate. To establish this blueprint, we develop a refined variant of the two-step Lyapunov approach presented in Section 5.3.

References↩︎

---

1. Given a finite set \(\mathcal{X}\), we use \(\Delta(\mathcal{X})\) to denote the set of probability distributions supported on \(\mathcal{X}\).↩︎

2. The sample average \(\frac{1}{k}\sum_{i=0}^{k-1} Y_i\) can be computed iteratively as \(x_{k+1} = x_k + \frac{1}{k+1}(-x_k + Y_k)\) with \(x_0 = 0\), which is a special case of Algorithm 1 .↩︎